NAMD
ComputeHomeTuples.h
Go to the documentation of this file.
1 
7 #ifndef COMPUTEHOMETUPLES_H
8 #define COMPUTEHOMETUPLES_H
9 
10 #ifdef USE_HOMETUPLES
11 #include <vector>
12 #endif
13 #include "NamdTypes.h"
14 #include "common.h"
15 #include "structures.h"
16 #include "Compute.h"
17 #include "HomePatch.h"
18 
19 #include "Box.h"
20 #include "OwnerBox.h"
21 #include "UniqueSet.h"
22 
23 #include "Node.h"
24 #include "SimParameters.h"
25 #include "PatchMap.inl"
26 #include "AtomMap.h"
27 #include "PatchMgr.h"
28 #include "ProxyMgr.h"
29 #include "HomePatchList.h"
30 #include "Molecule.h"
31 #include "Parameters.h"
32 #include "ReductionMgr.h"
33 #include "UniqueSet.h"
34 #include "UniqueSetIter.h"
35 #include "Priorities.h"
36 #include "LdbCoordinator.h"
37 
39  public:
41  Patch *p;
49  Force *f;
51 
52  int hash() const { return patchID; }
53 
54  TuplePatchElem(PatchID pid = -1) {
55  patchID = pid;
56  p = NULL;
57  positionBox = NULL;
58  avgPositionBox = NULL;
59  forceBox = NULL;
60  x = NULL;
61  xExt = NULL;
62  x_avg = NULL;
63  r = NULL;
64  f = NULL;
65  af = NULL;
66  }
67 
68  TuplePatchElem(Patch *p_param, Compute *cid) {
69  patchID = p_param->getPatchID();
70  p = p_param;
71  positionBox = p_param->registerPositionPickup(cid);
73  forceBox = p_param->registerForceDeposit(cid);
74  x = NULL;
75  xExt = NULL;
76  x_avg = NULL;
77  r = NULL;
78  f = NULL;
79  af = NULL;
80  }
81 
83 
84  int operator==(const TuplePatchElem &elem) const {
85  return (elem.patchID == patchID);
86  }
87 
88  int operator<(const TuplePatchElem &elem) const {
89  return (patchID < elem.patchID);
90  }
91 };
92 
95 
96 class AtomMap;
97 class ReductionMgr;
98 
99 #ifdef MEM_OPT_VERSION
100 template <class T> struct ElemTraits {
101  typedef AtomSignature signature;
102  static signature* get_sig_pointer(Molecule *mol) { return mol->atomSigPool; }
103  static int get_sig_id(const CompAtomExt &a) { return a.sigId; }
104 };
105 
106 template <> struct ElemTraits <ExclElem> {
107  typedef ExclusionSignature signature;
108  static signature* get_sig_pointer(Molecule *mol) { return mol->exclSigPool; }
109  static int get_sig_id(const CompAtomExt &a) { return a.exclId; }
110 };
111 #endif
112 
113 #ifdef USE_HOMETUPLES
114 //
115 // Simple base class for HomeTuples and SelfTuples that stores the type of the tuple
116 //
117 class Tuples {
118 private:
119  int type;
120 protected:
121  Tuples(int type) : type(type) {}
122 public:
123  // Tuple types
124  enum {BOND=0, ANGLE, DIHEDRAL, IMPROPER, EXCLUSION, CROSSTERM, NUM_TUPLE_TYPES};
125 
126  int getType() {return type;}
127  virtual void submitTupleCount(SubmitReduction *reduction, int tupleCount)=0;
128  // virtual void copyTupleData(void* tupleData)=0;
129  virtual int getNumTuples()=0;
130  virtual void* getTupleList()=0;
131  virtual void loadTuples(TuplePatchList& tuplePatchList, const char* isBasePatch, AtomMap *atomMap,
132  const std::vector<int>& pids = std::vector<int>())=0;
133 };
134 
135 //
136 // HomeTuples class. These are created and stored in ComputeBondedCUDA::registerCompute()
137 // e.g.: new HomeTuples<BondElem, Bond, BondValue>(BOND)
138 //
139 template <class T, class S, class P> class HomeTuples : public Tuples {
140  protected:
141  std::vector<T> tupleList;
142 
143  public:
144 
145  HomeTuples(int type=-1) : Tuples(type) {}
146 
147 #if __cplusplus < 201103L
148 #define final
149 #endif
150 
151  virtual void* getTupleList() final {
152  return (void*)tupleList.data();
153  }
154 
155  virtual void submitTupleCount(SubmitReduction *reduction, int tupleCount) final {
156  reduction->item(T::reductionChecksumLabel) += (BigReal)tupleCount;
157  }
158 
159  // virtual void copyTupleData(void* tupleData) final {
160  // for (int i=0;i < tupleList.size();i++) {
161  // tupleData[i] =
162  // }
163  // T::loadTupleData(tupleData);
164  // }
165 
166  virtual int getNumTuples() final {
167  return tupleList.size();
168  }
169 
170  virtual void loadTuples(TuplePatchList& tuplePatchList, const char* isBasePatch, AtomMap *atomMap,
171  const std::vector<int>& pids = std::vector<int>()) {
172 
173  if (isBasePatch == NULL) {
174  NAMD_bug("NULL isBasePatch detected in HomeTuples::loadTuples()");
175  }
176 
177  int numTuples;
178 
179 #ifdef MEM_OPT_VERSION
180  typename ElemTraits<T>::signature *allSigs;
181 #else
182  int32 **tuplesByAtom;
183  /* const (need to propagate const) */ S *tupleStructs;
184 #endif
185 
186  const P *tupleValues;
187  Node *node = Node::Object();
188  PatchMap *patchMap = PatchMap::Object();
189  // AtomMap *atomMap = AtomMap::Object();
190 
191 #ifdef MEM_OPT_VERSION
192  allSigs = ElemTraits<T>::get_sig_pointer(node->molecule);
193 #else
194  T::getMoleculePointers(node->molecule,
195  &numTuples, &tuplesByAtom, &tupleStructs);
196 #endif
197 
198  T::getParameterPointers(node->parameters, &tupleValues);
199 
200  tupleList.clear();
201 
202  LocalID aid[T::size];
203  int partition[T::size];
204 
205  const int lesOn = node->simParameters->lesOn;
206  const int soluteScalingOn = node->simParameters->soluteScalingOn;
207  const int fepOn = node->simParameters->singleTopology;
208  const int sdScaling = node->simParameters->sdScaling;
209  Real invLesFactor = lesOn ?
210  1.0/node->simParameters->lesFactor :
211  1.0;
212  const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
213  const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
214  BigReal OneMinusLambda = 1.0 - node->simParameters->alchLambda;
215  BigReal Lambda = node->simParameters->alchLambda;
216  const int num_unpert_bonds = node->molecule->num_alch_unpert_Bonds;
217  const int num_unpert_angles = node->molecule->num_alch_unpert_Angles;
218  const int num_unpert_dihedrals = node->molecule->num_alch_unpert_Dihedrals;
219  Bond *unpert_bonds = node->molecule->alch_unpert_bonds;
220  Angle *unpert_angles = node->molecule->alch_unpert_angles;
221  Dihedral *unpert_dihedrals = node->molecule->alch_unpert_dihedrals;
222 
223  // cycle through each patch and gather all tuples
224  TuplePatchListIter ai(tuplePatchList);
225  if (pids.size() == 0) ai = ai.begin();
226 
227  int numPid = (pids.size() == 0) ? tuplePatchList.size() : pids.size();
228 
229  for (int ipid=0;ipid < numPid;ipid++) {
230  // Patch *patch;
231  int numAtoms;
232  CompAtomExt *atomExt;
233  // Take next patch
234  if (pids.size() == 0) {
235  Patch* patch = (*ai).p;
236  numAtoms = patch->getNumAtoms();
237  atomExt = (*ai).xExt;
238  ai++;
239  } else {
240  TuplePatchElem *tpe = tuplePatchList.find(TuplePatchElem(pids[ipid]));
241  Patch* patch = tpe->p;
242  numAtoms = patch->getNumAtoms();
243  atomExt = tpe->xExt;
244  }
245 
246  // cycle through each atom in the patch and load up tuples
247  for (int j=0; j < numAtoms; j++)
248  {
249  /* cycle through each tuple */
250 #ifdef MEM_OPT_VERSION
251  typename ElemTraits<T>::signature *thisAtomSig =
252  &allSigs[ElemTraits<T>::get_sig_id(atomExt[j])];
253  TupleSignature *allTuples;
254  T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
255  for(int k=0; k<numTuples; k++) {
256  T t(atomExt[j].id, &allTuples[k], tupleValues);
257 #else
258  /* get list of all tuples for the atom */
259  int32 *curTuple = tuplesByAtom[atomExt[j].id];
260  for( ; *curTuple != -1; ++curTuple) {
261  T t(&tupleStructs[*curTuple],tupleValues);
262 #endif
263  register int i;
264  aid[0] = atomMap->localID(t.atomID[0]);
265  int homepatch = aid[0].pid;
266  int samepatch = 1;
267  partition[0] = fepOn ? node->molecule->get_fep_type(t.atomID[0]) : 0; //using atom partition to determine if a bonded term to be scaled by lambda or 1-lambda in single topology relative FEP.
268  int has_les = lesOn ? node->molecule->get_fep_type(t.atomID[0]) : 0;
269  int has_ss = soluteScalingOn ? node->molecule->get_ss_type(t.atomID[0]) : 0;
270  int is_fep_ss = partition[0] > 2;
271  int is_fep_sd = 0;
272  int fep_tuple_type = 0;
273  for (i=1; i < T::size; i++) {
274  aid[i] = atomMap->localID(t.atomID[i]);
275  samepatch = samepatch && ( homepatch == aid[i].pid );
276  partition[i] = fepOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
277  has_les |= lesOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
278  has_ss |= soluteScalingOn ? node->molecule->get_ss_type(t.atomID[i]) : 0;
279  if (fepOn) {
280  is_fep_ss &= partition[i] > 2;
281  is_fep_sd |= (abs(partition[i] - partition[0]) == 2);
282  fep_tuple_type = partition[i]; }
283  }
284  if (sdScaling && is_fep_sd) {
285  // check if this bonded term is one of Shobana term.
286  // This segment looks ugly and not GPU friendly,
287  // and might not appear in GPU code.
288  //
289  // XXX Could optimize in a number of ways:
290  // - could switch on T::size, then loop for that sized tuple
291  // - could use hash table to look up unpert_*[] elements
292  // - could add flag field to BondElem, et al., classes
293  for (i=0; i < num_unpert_bonds; i++) {
294  if (T::size == 2
295  && t.atomID[0]==unpert_bonds[i].atom1
296  && t.atomID[1]==unpert_bonds[i].atom2) is_fep_sd = 0;
297  }
298  for (i=0; i < num_unpert_angles; i++) {
299  if (T::size == 3
300  && t.atomID[0]==unpert_angles[i].atom1
301  && t.atomID[1]==unpert_angles[i].atom2
302  && t.atomID[2]==unpert_angles[i].atom3) is_fep_sd = 0;
303  }
304  for (i=0; i < num_unpert_dihedrals; i++) {
305  if (T::size == 4
306  && t.atomID[0]==unpert_dihedrals[i].atom1
307  && t.atomID[1]==unpert_dihedrals[i].atom2
308  && t.atomID[2]==unpert_dihedrals[i].atom3
309  && t.atomID[3]==unpert_dihedrals[i].atom4) is_fep_sd = 0;
310  }
311  }
312  if (T::size < 4 && !soluteScalingAll) has_ss = false;
313  if ( samepatch ) continue;
314  t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
315  if (is_fep_ss) t.scale = (fep_tuple_type == 4) ? OneMinusLambda : Lambda;
316  if (is_fep_sd && sdScaling) t.scale = (fep_tuple_type == 4 || fep_tuple_type == 2) ? OneMinusLambda : Lambda;
317 
318  for (i=1; i < T::size; i++) {
319  homepatch = patchMap->downstream(homepatch,aid[i].pid);
320  }
321  if ( homepatch != notUsed && isBasePatch[homepatch] ) {
322  TuplePatchElem *p;
323  for (i=0; i < T::size; i++) {
324  t.p[i] = p = tuplePatchList.find(TuplePatchElem(aid[i].pid));
325  if ( ! p ) {
326 #ifdef MEM_OPT_VERSION
327  iout << iWARN << "Tuple with atoms ";
328 #else
329  iout << iWARN << "Tuple " << *curTuple << " with atoms ";
330 #endif
331  int erri;
332  for( erri = 0; erri < T::size; erri++ ) {
333  iout << t.atomID[erri] << "(" << aid[erri].pid << ") ";
334  }
335  iout << "missing patch " << aid[i].pid << "\n" << endi;
336  break;
337  }
338  t.localIndex[i] = aid[i].index;
339  }
340  if ( ! p ) continue;
341 #ifdef MEM_OPT_VERSION
342  //avoid adding Tuples whose atoms are all fixed
344  int allfixed = 1;
345  for(i=0; i<T::size; i++){
346  CompAtomExt *one = &(t.p[i]->xExt[aid[i].index]);
347  allfixed = allfixed & one->atomFixed;
348  }
349  if(!allfixed) tupleList.push_back(t);
350  }else{
351  tupleList.push_back(t);
352  }
353 #else
354  tupleList.push_back(t);
355 #endif
356  }
357  }
358  }
359  }
360  }
361 
362 };
363 #endif
364 
365 template <class T, class S, class P> class ComputeHomeTuples : public Compute {
366 
367  protected:
368 
369 #ifndef USE_HOMETUPLES
370  virtual void loadTuples(void) {
371  int numTuples;
372 
373  #ifdef MEM_OPT_VERSION
374  typename ElemTraits<T>::signature *allSigs;
375  #else
376  int32 **tuplesByAtom;
377  /* const (need to propagate const) */ S *tupleStructs;
378  #endif
379 
380  const P *tupleValues;
381  Node *node = Node::Object();
382 
383  #ifdef MEM_OPT_VERSION
384  allSigs = ElemTraits<T>::get_sig_pointer(node->molecule);
385  #else
386  T::getMoleculePointers(node->molecule,
387  &numTuples, &tuplesByAtom, &tupleStructs);
388  #endif
389 
390  T::getParameterPointers(node->parameters, &tupleValues);
391 
392  tupleList.resize(0);
393 
394  LocalID aid[T::size];
395  int partition[T::size];
396 
397  const int lesOn = node->simParameters->lesOn;
398  const int soluteScalingOn = node->simParameters->soluteScalingOn;
399  const int fepOn = node->simParameters->singleTopology;
400  const int sdScaling = node->simParameters->sdScaling;
401  Real invLesFactor = lesOn ?
402  1.0/node->simParameters->lesFactor :
403  1.0;
404  const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
405  const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
406  BigReal OneMinusLambda = 1.0 - node->simParameters->alchLambda;
407  BigReal Lambda = node->simParameters->alchLambda;
408  const int num_unpert_bonds = node->molecule->num_alch_unpert_Bonds;
409  const int num_unpert_angles = node->molecule->num_alch_unpert_Angles;
410  const int num_unpert_dihedrals = node->molecule->num_alch_unpert_Dihedrals;
411  Bond *unpert_bonds = node->molecule->alch_unpert_bonds;
412  Angle *unpert_angles = node->molecule->alch_unpert_angles;
413  Dihedral *unpert_dihedrals = node->molecule->alch_unpert_dihedrals;
414 
415  // cycle through each patch and gather all tuples
416  TuplePatchListIter ai(tuplePatchList);
417 
418  for ( ai = ai.begin(); ai != ai.end(); ai++ )
419  {
420  // CompAtom *atom = (*ai).x;
421  Patch *patch = (*ai).p;
422  int numAtoms = patch->getNumAtoms();
423  CompAtomExt *atomExt = (*ai).xExt; //patch->getCompAtomExtInfo();
424 
425  // cycle through each atom in the patch and load up tuples
426  for (int j=0; j < numAtoms; j++)
427  {
428  /* cycle through each tuple */
429  #ifdef MEM_OPT_VERSION
430  typename ElemTraits<T>::signature *thisAtomSig =
431  &allSigs[ElemTraits<T>::get_sig_id(atomExt[j])];
432  TupleSignature *allTuples;
433  T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
434  for(int k=0; k<numTuples; k++) {
435  T t(atomExt[j].id, &allTuples[k], tupleValues);
436  #else
437  /* get list of all tuples for the atom */
438  int32 *curTuple = tuplesByAtom[atomExt[j].id];
439  for( ; *curTuple != -1; ++curTuple) {
440  T t(&tupleStructs[*curTuple],tupleValues);
441  #endif
442  register int i;
443  aid[0] = atomMap->localID(t.atomID[0]);
444  int homepatch = aid[0].pid;
445  int samepatch = 1;
446  partition[0] = fepOn ? node->molecule->get_fep_type(t.atomID[0]) : 0;
447  int has_les = lesOn ? node->molecule->get_fep_type(t.atomID[0]) : 0;
448  int has_ss = soluteScalingOn ? node->molecule->get_ss_type(t.atomID[0]) : 0;
449  int is_fep_ss = partition[0] > 2;
450  int is_fep_sd = 0;
451  int fep_tuple_type = 0;
452  for (i=1; i < T::size; i++) {
453  aid[i] = atomMap->localID(t.atomID[i]);
454  samepatch = samepatch && ( homepatch == aid[i].pid );
455  partition[i] = fepOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
456  has_les |= lesOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
457  has_ss |= soluteScalingOn ? node->molecule->get_ss_type(t.atomID[i]) : 0;
458  if (fepOn) {
459  is_fep_ss &= partition[i] > 2;
460  is_fep_sd |= (abs(partition[i] - partition[0]) == 2);
461  fep_tuple_type = partition[i];
462  }
463  }
464  if (sdScaling && is_fep_sd) {
465  for (i=0; i < num_unpert_bonds; i++) {
466  if (T::size == 2
467  && t.atomID[0]==unpert_bonds[i].atom1
468  && t.atomID[1]==unpert_bonds[i].atom2) is_fep_sd = 0;
469  }
470  for (i=0; i < num_unpert_angles; i++) {
471  if (T::size == 3
472  && t.atomID[0]==unpert_angles[i].atom1
473  && t.atomID[1]==unpert_angles[i].atom2
474  && t.atomID[2]==unpert_angles[i].atom3) is_fep_sd = 0;
475  }
476  for (i=0; i < num_unpert_dihedrals; i++) {
477  if (T::size == 4
478  && t.atomID[0]==unpert_dihedrals[i].atom1
479  && t.atomID[1]==unpert_dihedrals[i].atom2
480  && t.atomID[2]==unpert_dihedrals[i].atom3
481  && t.atomID[3]==unpert_dihedrals[i].atom4) is_fep_sd = 0;
482  }
483  }
484  if (T::size < 4 && !soluteScalingAll) has_ss = false;
485  if ( samepatch ) continue;
486  t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
487  if (is_fep_ss) t.scale = (fep_tuple_type == 4) ? OneMinusLambda : Lambda;
488  if (is_fep_sd && sdScaling) t.scale = (fep_tuple_type == 4 || fep_tuple_type == 2) ? OneMinusLambda : Lambda;
489 
490  for (i=1; i < T::size; i++) {
491  homepatch = patchMap->downstream(homepatch,aid[i].pid);
492  }
493  if ( homepatch != notUsed && isBasePatch[homepatch] ) {
494  TuplePatchElem *p;
495  for (i=0; i < T::size; i++) {
496  t.p[i] = p = tuplePatchList.find(TuplePatchElem(aid[i].pid));
497  if ( ! p ) {
498  #ifdef MEM_OPT_VERSION
499  iout << iWARN << "Tuple with atoms ";
500  #else
501  iout << iWARN << "Tuple " << *curTuple << " with atoms ";
502  #endif
503  int erri;
504  for( erri = 0; erri < T::size; erri++ ) {
505  iout << t.atomID[erri] << "(" << aid[erri].pid << ") ";
506  }
507  iout << "missing patch " << aid[i].pid << "\n" << endi;
508  break;
509  }
510  t.localIndex[i] = aid[i].index;
511  }
512  if ( ! p ) continue;
513  #ifdef MEM_OPT_VERSION
514  //avoid adding Tuples whose atoms are all fixed
515  if(node->simParameters->fixedAtomsOn &&
517  int allfixed = 1;
518  for(i=0; i<T::size; i++){
519  CompAtomExt *one = &(t.p[i]->xExt[aid[i].index]);
520  allfixed = allfixed & one->atomFixed;
521  }
522  if(!allfixed) tupleList.add(t);
523  }else{
524  tupleList.add(t);
525  }
526  #else
527  tupleList.add(t);
528  #endif
529  }
530  }
531  }
532  }
533  }
534 #endif
535 
537 
538  protected:
539 
540 #ifdef USE_HOMETUPLES
541  HomeTuples<T, S, P>* tuples;
542  TuplePatchList tuplePatchList;
543 #else
546 #endif
547 
555  char *isBasePatch;
556 
558  patchMap = PatchMap::Object();
559  atomMap = AtomMap::Object();
561 
563  accelMDdoDihe=false;
564  if (params->accelMDOn) {
565  if (params->accelMDdihe || params->accelMDdual) accelMDdoDihe=true;
566  }
567  if (params->pressureProfileOn) {
568  pressureProfileSlabs = T::pressureProfileSlabs =
569  params->pressureProfileSlabs;
570  int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
571  pressureProfileReduction = ReductionMgr::Object()->willSubmit(
572  REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
573  int numAtomTypePairs = n*n;
574  pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
575  } else {
576  pressureProfileReduction = NULL;
577  pressureProfileData = NULL;
578  }
579  doLoadTuples = false;
580  isBasePatch = 0;
581 #ifdef USE_HOMETUPLES
582  tuples = NULL;
583 #endif
584  }
585 
587  patchMap = PatchMap::Object();
588  atomMap = AtomMap::Object();
591  accelMDdoDihe=false;
592  if (params->accelMDOn) {
593  if (params->accelMDdihe || params->accelMDdual) accelMDdoDihe=true;
594  }
595  if (params->pressureProfileOn) {
596  pressureProfileSlabs = T::pressureProfileSlabs =
597  params->pressureProfileSlabs;
598  int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
599  pressureProfileReduction = ReductionMgr::Object()->willSubmit(
600  REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
601  int numAtomTypePairs = n*n;
602  pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
603  } else {
604  pressureProfileReduction = NULL;
605  pressureProfileData = NULL;
606  }
607  doLoadTuples = false;
608  int nPatches = patchMap->numPatches();
609  isBasePatch = new char[nPatches];
610  int i;
611  for (i=0; i<nPatches; ++i) { isBasePatch[i] = 0; }
612  for (i=0; i<pids.size(); ++i) { isBasePatch[pids[i]] = 1; }
613 #ifdef USE_HOMETUPLES
614  tuples = NULL;
615 #endif
616  }
617 
618  public:
619 
620  virtual ~ComputeHomeTuples() {
621  delete reduction;
622  delete [] isBasePatch;
623  delete pressureProfileReduction;
624  delete pressureProfileData;
625 #ifdef USE_HOMETUPLES
626  if (tuples != NULL) delete tuples;
627 #endif
628  }
629 
630  //======================================================================
631  // initialize() - Method is invoked only the first time
632  // atom maps, patchmaps etc are ready and we are about to start computations
633  //======================================================================
634  virtual void initialize(void) {
635 
636 #ifdef NAMD_CUDA
637  ProxyMgr *proxyMgr = ProxyMgr::Object();
638 #endif
639 
640 #ifdef USE_HOMETUPLES
641  tuples = new HomeTuples<T, S, P>();
642 #endif
643 
644  // Start with empty list
645  tuplePatchList.clear();
646 
647  int nPatches = patchMap->numPatches();
648  int pid;
649  for (pid=0; pid<nPatches; ++pid) {
650  if ( isBasePatch[pid] ) {
651 #ifdef NAMD_CUDA
652  proxyMgr->createProxy(pid);
653 #endif
654  Patch *patch = patchMap->patch(pid);
655  tuplePatchList.add(TuplePatchElem(patch, this));
656  }
657  }
658 
659  // Gather all proxy patches (neighbors, that is)
661 
662  for (pid=0; pid<nPatches; ++pid) if ( isBasePatch[pid] ) {
663  int numNeighbors = patchMap->upstreamNeighbors(pid,neighbors);
664  for ( int i = 0; i < numNeighbors; ++i ) {
665  if ( ! tuplePatchList.find(TuplePatchElem(neighbors[i])) ) {
666 #ifdef NAMD_CUDA
667  proxyMgr->createProxy(neighbors[i]);
668 #endif
669  Patch *patch = patchMap->patch(neighbors[i]);
670  tuplePatchList.add(TuplePatchElem(patch, this));
671  }
672  }
673  }
674  setNumPatches(tuplePatchList.size());
675  doLoadTuples = true;
676 
677  basePriority = COMPUTE_PROXY_PRIORITY; // no patch dependence
678  }
679 
680  //======================================================================
681  // atomUpdate() - Method is invoked after anytime that atoms have been
682  // changed in patches used by this Compute object.
683  //======================================================================
684  void atomUpdate(void) {
685  doLoadTuples = true;
686  }
687 
688 //-------------------------------------------------------------------
689 // Routine which is called by enqueued work msg. It wraps
690 // actualy Force computation with the apparatus needed
691 // to get access to atom positions, return forces etc.
692 //-------------------------------------------------------------------
693  virtual void doWork(void) {
694 
695  LdbCoordinator::Object()->startWork(ldObjHandle);
696 
697  // Open Boxes - register that we are using Positions
698  // and will be depositing Forces.
699  UniqueSetIter<TuplePatchElem> ap(tuplePatchList);
700  for (ap = ap.begin(); ap != ap.end(); ap++) {
701  ap->x = ap->positionBox->open();
702  ap->xExt = ap->p->getCompAtomExtInfo();
703  if ( ap->p->flags.doMolly ) ap->x_avg = ap->avgPositionBox->open();
704  ap->r = ap->forceBox->open();
705  ap->f = ap->r->f[Results::normal];
706  if (accelMDdoDihe) ap->af = ap->r->f[Results::amdf]; // for dihedral-only or dual-boost accelMD
707  }
708 
709  BigReal reductionData[T::reductionDataSize];
710  int tupleCount = 0;
711  int numAtomTypes = T::pressureProfileAtomTypes;
712  int numAtomTypePairs = numAtomTypes*numAtomTypes;
713 
714  for ( int i = 0; i < T::reductionDataSize; ++i ) reductionData[i] = 0;
715  if (pressureProfileData) {
716  memset(pressureProfileData, 0, 3*pressureProfileSlabs*numAtomTypePairs*sizeof(BigReal));
717  // Silly variable hiding of the previous iterator
718  UniqueSetIter<TuplePatchElem> newap(tuplePatchList);
719  newap = newap.begin();
720  const Lattice &lattice = newap->p->lattice;
721  T::pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
722  T::pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
723  }
724 
725  if ( ! Node::Object()->simParameters->commOnly ) {
726  if ( doLoadTuples ) {
727 #ifdef USE_HOMETUPLES
728  tuples->loadTuples(tuplePatchList, isBasePatch, AtomMap::Object());
729 #else
730  loadTuples();
731 #endif
732  doLoadTuples = false;
733  }
734  // take triplet and pass with tuple info to force eval
735 #ifdef USE_HOMETUPLES
736  T *al = (T *)tuples->getTupleList();
737  const int ntuple = tuples->getNumTuples();
738 #else
739  T *al = tupleList.begin();
740  const int ntuple = tupleList.size();
741 #endif
742  if ( ntuple ) T::computeForce(al, ntuple, reductionData, pressureProfileData);
743  tupleCount += ntuple;
744  }
745 
746  LdbCoordinator::Object()->endWork(ldObjHandle);
747 
748  T::submitReductionData(reductionData,reduction);
749  reduction->item(T::reductionChecksumLabel) += (BigReal)tupleCount;
750  reduction->submit();
751 
752  if (pressureProfileReduction) {
753  // For ease of calculation we stored interactions between types
754  // i and j in (ni+j). For efficiency now we coalesce the
755  // cross interactions so that just i<=j are stored.
756  const int arraysize = 3*pressureProfileSlabs;
757  const BigReal *data = pressureProfileData;
758  for (int i=0; i<numAtomTypes; i++) {
759  for (int j=0; j<numAtomTypes; j++) {
760  int ii=i;
761  int jj=j;
762  if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
763  const int reductionOffset =
764  (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
765  for (int k=0; k<arraysize; k++) {
766  pressureProfileReduction->item(reductionOffset+k) += data[k];
767  }
768  data += arraysize;
769  }
770  }
771  pressureProfileReduction->submit();
772  }
773 
774  // Close boxes - i.e. signal we are done with Positions and
775  // AtomProperties and that we are depositing Forces
776  for (ap = ap.begin(); ap != ap.end(); ap++) {
777  ap->positionBox->close(&(ap->x));
778  if ( ap->p->flags.doMolly ) ap->avgPositionBox->close(&(ap->x_avg));
779  ap->forceBox->close(&(ap->r));
780  }
781  }
782 };
783 
784 
785 #endif
786 
static Node * Object()
Definition: Node.h:86
Elem * find(const Elem &elem)
Definition: UniqueSet.h:60
UniqueSet< TuplePatchElem > TuplePatchList
#define COMPUTE_PROXY_PRIORITY
Definition: Priorities.h:71
int num_alch_unpert_Dihedrals
Definition: Molecule.h:574
Box< Patch, CompAtom > * registerAvgPositionPickup(Compute *cid)
Definition: Patch.C:134
static ProxyMgr * Object()
Definition: ProxyMgr.h:394
Definition: Node.h:78
short int32
Definition: dumpdcd.c:24
int ComputeID
Definition: NamdTypes.h:183
unsigned int atomFixed
Definition: NamdTypes.h:96
void clear(void)
Definition: UniqueSet.h:62
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:41
int32 atom3
Definition: structures.h:65
static PatchMap * Object()
Definition: PatchMap.h:27
Angle * alch_unpert_angles
Definition: Molecule.h:576
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
int num_alch_unpert_Bonds
Definition: Molecule.h:572
float Real
Definition: common.h:107
BigReal & item(int i)
Definition: ReductionMgr.h:312
BigReal * pressureProfileData
void startWork(const LDObjHandle &handle)
BigReal z
Definition: Vector.h:66
virtual ~ComputeHomeTuples()
int downstream(int pid1, int pid2)
Definition: PatchMap.inl:51
int upstreamNeighbors(int pid, PatchID *neighbor_ids)
Definition: PatchMap.C:669
TuplePatchList tuplePatchList
ComputeHomeTuples(ComputeID c)
SubmitReduction * pressureProfileReduction
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
int num_alch_unpert_Angles
Definition: Molecule.h:573
int add(const Elem &elem)
Definition: UniqueSet.h:52
int size(void) const
Definition: UniqueSet.h:58
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define iout
Definition: InfoStream.h:87
int pressureProfileSlabs
Patch * patch(PatchID pid)
Definition: PatchMap.h:235
virtual void initialize(void)
Vector origin() const
Definition: Lattice.h:262
BigReal alchLambda
Dihedral * alch_unpert_dihedrals
Definition: Molecule.h:577
Definition: Patch.h:35
UniqueSetIter< T > begin(void) const
Definition: UniqueSetIter.h:55
int32 atom4
Definition: structures.h:66
int32 atom1
Definition: structures.h:48
TuplePatchElem(Patch *p_param, Compute *cid)
void NAMD_bug(const char *err_msg)
Definition: common.C:123
int index
Definition: NamdTypes.h:195
int32 atom2
Definition: structures.h:56
TuplePatchElem(PatchID pid=-1)
BigReal soluteScalingFactor
int32 atom1
Definition: structures.h:55
int Bool
Definition: common.h:131
unsigned char get_ss_type(int anum) const
Definition: Molecule.h:1351
CompAtomList p
Definition: Patch.h:146
LocalID localID(AtomID id)
Definition: AtomMap.h:74
int32 atom3
Definition: structures.h:57
Bond * alch_unpert_bonds
Definition: Molecule.h:575
int PatchID
Definition: NamdTypes.h:182
Box< Patch, CompAtom > * positionBox
void createProxy(PatchID pid)
Definition: ProxyMgr.C:493
virtual void doWork(void)
static LdbCoordinator * Object()
PatchID pid
Definition: NamdTypes.h:194
static AtomMap * Object()
Definition: AtomMap.h:36
int32 atom2
Definition: structures.h:64
unsigned char get_fep_type(int anum) const
Definition: Molecule.h:1345
Parameters * parameters
Definition: Node.h:177
void endWork(const LDObjHandle &handle)
ResizeArray< T > tupleList
PatchID getPatchID()
Definition: Patch.h:114
UniqueSetIter< T > end(void) const
Definition: UniqueSetIter.h:64
int pressureProfileAtomTypes
int numPatches(void) const
Definition: PatchMap.h:59
SubmitReduction * reduction
int getNumAtoms()
Definition: Patch.h:105
ComputeHomeTuples(ComputeID c, PatchIDList &pids)
Bool pressureProfileOn
virtual void loadTuples(void)
Box< Patch, CompAtom > * avgPositionBox
CompAtomExt * xExt
void submit(void)
Definition: ReductionMgr.h:323
int size(void) const
Definition: ResizeArray.h:127
infostream & endi(infostream &s)
Definition: InfoStream.C:38
int hash() const
Molecule * molecule
Definition: Node.h:176
UniqueSetIter< TuplePatchElem > TuplePatchListIter
Box< Patch, Results > * forceBox
int operator<(const TuplePatchElem &elem) const
int32 atom1
Definition: structures.h:63
int32 atom2
Definition: structures.h:49
Vector c() const
Definition: Lattice.h:254
Box< Patch, CompAtom > * registerPositionPickup(Compute *cid)
Definition: Patch.C:107
double BigReal
Definition: common.h:112
int operator==(const TuplePatchElem &elem) const
iterator begin(void)
Definition: ResizeArray.h:36
Box< Patch, Results > * registerForceDeposit(Compute *cid)
Definition: Patch.C:228