NAMD
ComputeSelfTuples.h
Go to the documentation of this file.
1 
7 #ifndef COMPUTESELFTUPLES_H
8 #define COMPUTESELFTUPLES_H
9 
10 #include "ComputeHomeTuples.h"
11 #include "LdbCoordinator.h"
12 
13 #ifdef USE_HOMETUPLES
14 template <class T, class S, class P> class SelfTuples : public HomeTuples<T, S, P> {
15 
16  public:
17  SelfTuples(int type=-1) : HomeTuples<T,S,P>(type) {}
18 
19  private:
20 
21  virtual void loadTuples(TuplePatchList& tuplePatchList, const char* isBasePatch, AtomMap *atomMap,
22  const std::vector<int>& pids = std::vector<int>()) {
23 
24  if (isBasePatch != NULL) {
25  iout << iWARN << "Non-NULL isBasePatch detected in SelfTuples::loadTuples()" << endi;
26  }
27 
28  int numTuples;
29 
30  #ifdef MEM_OPT_VERSION
31  typename ElemTraits<T>::signature *allSigs;
32  #else
33  int32 **tuplesByAtom;
34  /* const (need to propagate const) */ S *tupleStructs;
35  #endif
36 
37  const P *tupleValues;
38  Node *node = Node::Object();
39  PatchMap *patchMap = PatchMap::Object();
40  // AtomMap *atomMap = AtomMap::Object();
41 
42  #ifdef MEM_OPT_VERSION
43  allSigs = ElemTraits<T>::get_sig_pointer(node->molecule);
44  #else
45  T::getMoleculePointers(node->molecule,
46  &numTuples, &tuplesByAtom, &tupleStructs);
47  #endif
48 
49  T::getParameterPointers(node->parameters, &tupleValues);
50 
51  this->tupleList.clear();
52 
53  LocalID aid[T::size];
54  int partition[T::size];
55 
56  const int lesOn = node->simParameters->lesOn;
57  const int soluteScalingOn = node->simParameters->soluteScalingOn;
58  const int fepOn = node->simParameters->singleTopology;
59  const int sdScaling = node->simParameters->sdScaling;
60  Real invLesFactor = lesOn ?
61  1.0/node->simParameters->lesFactor :
62  1.0;
63  const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
64  const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
65  BigReal OneMinusLambda = 1.0 - node->simParameters->alchLambda;
66  BigReal Lambda = node->simParameters->alchLambda;
67  const int num_unpert_bonds = node->molecule->num_alch_unpert_Bonds;
68  const int num_unpert_angles = node->molecule->num_alch_unpert_Angles;
69  const int num_unpert_dihedrals = node->molecule->num_alch_unpert_Dihedrals;
70  Bond *unpert_bonds = node->molecule->alch_unpert_bonds;
71  Angle *unpert_angles = node->molecule->alch_unpert_angles;
72  Dihedral *unpert_dihedrals = node->molecule->alch_unpert_dihedrals;
73 
74  // cycle through each patch and gather all tuples
75  // There should be only one!
76  TuplePatchListIter ai(tuplePatchList);
77  if (pids.size() == 0) ai = ai.begin();
78 
79  int numPid = (pids.size() == 0) ? tuplePatchList.size() : pids.size();
80 
81  for (int ipid=0;ipid < numPid;ipid++) {
82  Patch *patch;
83  int numAtoms;
84  CompAtomExt *atomExt;
85  // Take next patch
86  if (pids.size() == 0) {
87  patch = (*ai).p;
88  numAtoms = patch->getNumAtoms();
89  atomExt = (*ai).xExt;
90  ai++;
91  } else {
92  TuplePatchElem *tpe = tuplePatchList.find(TuplePatchElem(pids[ipid]));
93  patch = tpe->p;
94  numAtoms = patch->getNumAtoms();
95  atomExt = tpe->xExt;
96  }
97 
98  // cycle through each atom in the patch and load up tuples
99  for (int j=0; j < numAtoms; j++)
100  {
101 #ifdef MEM_OPT_VERSION
102  typename ElemTraits<T>::signature *thisAtomSig =
103  &allSigs[ElemTraits<T>::get_sig_id(atomExt[j])];
104  TupleSignature *allTuples;
105  T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
106  for(int k=0; k<numTuples; k++) {
107  T t(atomExt[j].id, &allTuples[k], tupleValues);
108 #else
109  /* get list of all tuples for the atom */
110  int32 *curTuple = tuplesByAtom[atomExt[j].id];
111  /* cycle through each tuple */
112  for( ; *curTuple != -1; ++curTuple) {
113  T t(&tupleStructs[*curTuple],tupleValues);
114 #endif
115  register int i;
116  aid[0] = atomMap->localID(t.atomID[0]);
117  int homepatch = aid[0].pid;
118  int samepatch = 1;
119  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.
120  int has_les = lesOn ? node->molecule->get_fep_type(t.atomID[0]) : 0;
121  int has_ss = soluteScalingOn ? node->molecule->get_ss_type(t.atomID[0]) : 0;
122  int is_fep_ss = partition[0] > 2;
123  int is_fep_sd = 0;
124  int fep_tuple_type = 0;
125  for (i=1; i < T::size; i++) {
126  aid[i] = atomMap->localID(t.atomID[i]);
127  samepatch = samepatch && ( homepatch == aid[i].pid );
128  partition[i] = fepOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
129  has_les |= lesOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
130  has_ss |= soluteScalingOn ? node->molecule->get_ss_type(t.atomID[i]) : 0;
131  if (fepOn) {
132  is_fep_ss &= partition[i] > 2;
133  is_fep_sd |= (abs(partition[i] - partition[0]) == 2);
134  fep_tuple_type = partition[i]; }
135  }
136  if (sdScaling && is_fep_sd) { //check if this bonded term is one of Shobana term. This segment looks ugly and not GPU friendly, and might not appear in GPU code.
137  for (i=0; i < num_unpert_bonds; i++) {
138  if (T::size == 2
139  && t.atomID[0]==unpert_bonds[i].atom1
140  && t.atomID[1]==unpert_bonds[i].atom2) is_fep_sd = 0;
141  }
142  for (i=0; i < num_unpert_angles; i++) {
143  if (T::size == 3
144  && t.atomID[0]==unpert_angles[i].atom1
145  && t.atomID[1]==unpert_angles[i].atom2
146  && t.atomID[2]==unpert_angles[i].atom3) is_fep_sd = 0;
147  }
148  for (i=0; i < num_unpert_dihedrals; i++) {
149  if (T::size == 4
150  && t.atomID[0]==unpert_dihedrals[i].atom1
151  && t.atomID[1]==unpert_dihedrals[i].atom2
152  && t.atomID[2]==unpert_dihedrals[i].atom3
153  && t.atomID[3]==unpert_dihedrals[i].atom4) is_fep_sd = 0;
154  }
155  }
156  if (T::size < 4 && !soluteScalingAll) has_ss = false;
157  if ( samepatch ) {
158  t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
159  if (is_fep_ss) t.scale = (fep_tuple_type == 4) ? OneMinusLambda : Lambda;
160  if (is_fep_sd && sdScaling) t.scale = (fep_tuple_type == 4 || fep_tuple_type == 2) ? OneMinusLambda : Lambda;
161  TuplePatchElem *p;
162  p = tuplePatchList.find(TuplePatchElem(homepatch));
163  for(i=0; i < T::size; i++) {
164  t.p[i] = p;
165  t.localIndex[i] = aid[i].index;
166  }
167 #ifdef MEM_OPT_VERSION
168  //avoid adding Tuples whose atoms are all fixed
169  if(node->simParameters->fixedAtomsOn &&
170  !node->simParameters->fixedAtomsForces) {
171  int allfixed = 1;
172  for(i=0; i<T::size; i++){
173  CompAtomExt *one = &(p->xExt[aid[i].index]);
174  allfixed = allfixed & one->atomFixed;
175  }
176  if(!allfixed) this->tupleList.push_back(t);
177  } else {
178  this->tupleList.push_back(t);
179  }
180 #else
181  this->tupleList.push_back(t);
182 #endif
183  }
184  }
185  }
186  }
187  }
188 
189 };
190 #endif
191 
192 template <class T, class S, class P> class ComputeSelfTuples :
193  public ComputeHomeTuples<T,S,P> {
194 
195 #ifndef USE_HOMETUPLES
196  private:
197 
198  virtual void loadTuples(void) {
199  int numTuples;
200 
201  #ifdef MEM_OPT_VERSION
202  typename ElemTraits<T>::signature *allSigs;
203  #else
204  int32 **tuplesByAtom;
205  /* const (need to propagate const) */ S *tupleStructs;
206  #endif
207 
208  const P *tupleValues;
209  Node *node = Node::Object();
210 
211  #ifdef MEM_OPT_VERSION
212  allSigs = ElemTraits<T>::get_sig_pointer(node->molecule);
213  #else
214  T::getMoleculePointers(node->molecule,
215  &numTuples, &tuplesByAtom, &tupleStructs);
216  #endif
217 
218  T::getParameterPointers(node->parameters, &tupleValues);
219 
220  this->tupleList.resize(0);
221 
222  LocalID aid[T::size];
223  int partition[T::size];
224 
225  const int lesOn = node->simParameters->lesOn;
226  const int soluteScalingOn = node->simParameters->soluteScalingOn;
227  const int fepOn = node->simParameters->singleTopology;
228  const int sdScaling = node->simParameters->sdScaling;
229  Real invLesFactor = lesOn ?
230  1.0/node->simParameters->lesFactor :
231  1.0;
232  const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
233  const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
234  BigReal OneMinusLambda = 1.0 - node->simParameters->alchLambda;
235  BigReal Lambda = node->simParameters->alchLambda;
236  const int num_unpert_bonds = node->molecule->num_alch_unpert_Bonds;
237  const int num_unpert_angles = node->molecule->num_alch_unpert_Angles;
238  const int num_unpert_dihedrals = node->molecule->num_alch_unpert_Dihedrals;
239  Bond *unpert_bonds = node->molecule->alch_unpert_bonds;
240  Angle *unpert_angles = node->molecule->alch_unpert_angles;
241  Dihedral *unpert_dihedrals = node->molecule->alch_unpert_dihedrals;
242 
243  // cycle through each patch and gather all tuples
244  // There should be only one!
245  TuplePatchListIter ai(this->tuplePatchList);
246 
247  for ( ai = ai.begin(); ai != ai.end(); ai++ )
248  {
249 
250  // CompAtomExt *atom = (*ai).x;
251  Patch *patch = (*ai).p;
252  int numAtoms = patch->getNumAtoms();
253  CompAtomExt *atomExt = (*ai).xExt; //patch->getCompAtomExtInfo();
254 
255  // cycle through each atom in the patch and load up tuples
256  for (int j=0; j < numAtoms; j++)
257  {
258  #ifdef MEM_OPT_VERSION
259  typename ElemTraits<T>::signature *thisAtomSig =
260  &allSigs[ElemTraits<T>::get_sig_id(atomExt[j])];
261  TupleSignature *allTuples;
262  T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
263  for(int k=0; k<numTuples; k++) {
264  T t(atomExt[j].id, &allTuples[k], tupleValues);
265  #else
266  /* get list of all tuples for the atom */
267  int32 *curTuple = tuplesByAtom[atomExt[j].id];
268  /* cycle through each tuple */
269  for( ; *curTuple != -1; ++curTuple) {
270  T t(&tupleStructs[*curTuple],tupleValues);
271  #endif
272  register int i;
273  aid[0] = this->atomMap->localID(t.atomID[0]);
274  int homepatch = aid[0].pid;
275  int samepatch = 1;
276  partition[0] = fepOn ? node->molecule->get_fep_type(t.atomID[0]) : 0;
277  int has_les = lesOn ? node->molecule->get_fep_type(t.atomID[0]) : 0;
278  int has_ss = soluteScalingOn ? node->molecule->get_ss_type(t.atomID[0]) : 0;
279  int is_fep_ss = partition[0] > 2;
280  int is_fep_sd = 0;
281  int fep_tuple_type = 0;
282  for (i=1; i < T::size; i++) {
283  aid[i] = this->atomMap->localID(t.atomID[i]);
284  samepatch = samepatch && ( homepatch == aid[i].pid );
285  partition[i] = fepOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
286  has_les |= lesOn ? node->molecule->get_fep_type(t.atomID[i]) : 0;
287  has_ss |= soluteScalingOn ? node->molecule->get_ss_type(t.atomID[i]) : 0;
288  if (fepOn) {
289  is_fep_ss &= partition[i] > 2;
290  is_fep_sd |= (abs(partition[i] - partition[0]) == 2);
291  fep_tuple_type = partition[i]; }
292  }
293  if (T::size < 4 && !soluteScalingAll) has_ss = false;
294  if (sdScaling && is_fep_sd) {
295  for (i=0; i < num_unpert_bonds; i++) {
296  if (T::size == 2
297  && t.atomID[0]==unpert_bonds[i].atom1
298  && t.atomID[1]==unpert_bonds[i].atom2) is_fep_sd = 0;
299  }
300  for (i=0; i < num_unpert_angles; i++) {
301  if (T::size == 3
302  && t.atomID[0]==unpert_angles[i].atom1
303  && t.atomID[1]==unpert_angles[i].atom2
304  && t.atomID[2]==unpert_angles[i].atom3) is_fep_sd = 0;
305  }
306  for (i=0; i < num_unpert_dihedrals; i++) {
307  if (T::size == 4
308  && t.atomID[0]==unpert_dihedrals[i].atom1
309  && t.atomID[1]==unpert_dihedrals[i].atom2
310  && t.atomID[2]==unpert_dihedrals[i].atom3
311  && t.atomID[3]==unpert_dihedrals[i].atom4) is_fep_sd = 0;
312  }
313  }
314  if ( samepatch ) {
315  t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
316  if (is_fep_ss) t.scale = (fep_tuple_type == 4) ? OneMinusLambda : Lambda;
317  if (is_fep_sd && sdScaling) t.scale = (fep_tuple_type == 4 || fep_tuple_type == 2) ? OneMinusLambda : Lambda;
318  TuplePatchElem *p;
319  p = this->tuplePatchList.find(TuplePatchElem(homepatch));
320  for(i=0; i < T::size; i++) {
321  t.p[i] = p;
322  t.localIndex[i] = aid[i].index;
323  }
324  #ifdef MEM_OPT_VERSION
325  //avoid adding Tuples whose atoms are all fixed
326  if(node->simParameters->fixedAtomsOn &&
327  !node->simParameters->fixedAtomsForces) {
328  int allfixed = 1;
329  for(i=0; i<T::size; i++){
330  CompAtomExt *one = &(p->xExt[aid[i].index]);
331  allfixed = allfixed & one->atomFixed;
332  }
333  if(!allfixed) this->tupleList.add(t);
334  }else{
335  this->tupleList.add(t);
336  }
337  #else
338  this->tupleList.add(t);
339  #endif
340  }
341  }
342  }
343  }
344  }
345 #endif
346 
347  PatchID patchID;
348 
349  public:
350 
351  ComputeSelfTuples(ComputeID c, PatchID p) : ComputeHomeTuples<T,S,P>(c) {
352  patchID = p;
353  }
354 
355  virtual ~ComputeSelfTuples() {
356  UniqueSetIter<TuplePatchElem> ap(this->tuplePatchList);
357  for (ap = ap.begin(); ap != ap.end(); ap++) {
358  ap->p->unregisterPositionPickup(this,&(ap->positionBox));
359  ap->p->unregisterAvgPositionPickup(this,&(ap->avgPositionBox));
360  ap->p->unregisterForceDeposit(this,&(ap->forceBox));
361  }
362  }
363 
364 
365  //======================================================================
366  // initialize() - Method is invoked only the first time
367  // atom maps, patchmaps etc are ready and we are about to start computations
368  //======================================================================
369  virtual void initialize(void) {
370 #ifdef USE_HOMETUPLES
371  this->tuples = new SelfTuples<T, S, P>();
372 #endif
373  // Start with empty list
374  this->tuplePatchList.clear();
375 
376  this->tuplePatchList.add(TuplePatchElem(ComputeHomeTuples<T,S,P>::patchMap->patch(patchID), this));
377 
378  this->setNumPatches(this->tuplePatchList.size());
379 
380  this->doLoadTuples = true;
381 
382  int myNode = CkMyPe();
383  if ( PatchMap::Object()->node(patchID) != myNode )
384  {
385  this->basePriority = COMPUTE_PROXY_PRIORITY + PATCH_PRIORITY(patchID);
386  }
387  else
388  {
389  this->basePriority = COMPUTE_HOME_PRIORITY + PATCH_PRIORITY(patchID);
390  }
391  }
392 
393  void doWork(void) {
394 // LdbCoordinator::Object()->startWork(this->ldObjHandle);
395 
396 #ifdef TRACE_COMPUTE_OBJECTS
397  double traceObjStartTime = CmiWallTimer();
398 #endif
399 
400  ComputeHomeTuples<T,S,P>::doWork();
401 
402 #ifdef TRACE_COMPUTE_OBJECTS
403  traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+this->cid, traceObjStartTime, CmiWallTimer());
404 #endif
405 
406 // LdbCoordinator::Object()->endWork(this->ldObjHandle);
407  }
408 
409 };
410 
411 
412 #endif
413 
Node::Object
static Node * Object()
Definition: Node.h:86
UniqueSet::find
Elem * find(const Elem &elem)
Definition: UniqueSet.h:60
COMPUTE_PROXY_PRIORITY
#define COMPUTE_PROXY_PRIORITY
Definition: Priorities.h:71
Molecule::num_alch_unpert_Dihedrals
int num_alch_unpert_Dihedrals
Definition: Molecule.h:575
Node
Definition: Node.h:78
int32
short int32
Definition: dumpdcd.c:24
ComputeID
int ComputeID
Definition: NamdTypes.h:183
TRACE_COMPOBJ_IDOFFSET
#define TRACE_COMPOBJ_IDOFFSET
Definition: Compute.h:77
CompAtomExt::atomFixed
unsigned int atomFixed
Definition: NamdTypes.h:96
UniqueSet::clear
void clear(void)
Definition: UniqueSet.h:62
partition
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
dihedral::atom3
int32 atom3
Definition: structures.h:65
PatchMap::Object
static PatchMap * Object()
Definition: PatchMap.h:27
Molecule::alch_unpert_angles
Angle * alch_unpert_angles
Definition: Molecule.h:577
Node::simParameters
SimParameters * simParameters
Definition: Node.h:178
Molecule::num_alch_unpert_Bonds
int num_alch_unpert_Bonds
Definition: Molecule.h:573
Real
float Real
Definition: common.h:109
endi
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
iWARN
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:82
Molecule::num_alch_unpert_Angles
int num_alch_unpert_Angles
Definition: Molecule.h:574
UniqueSet::add
int add(const Elem &elem)
Definition: UniqueSet.h:52
UniqueSet::size
int size(void) const
Definition: UniqueSet.h:58
iout
#define iout
Definition: InfoStream.h:51
SimParameters::alchLambda
BigReal alchLambda
Definition: SimParameters.h:396
Molecule::alch_unpert_dihedrals
Dihedral * alch_unpert_dihedrals
Definition: Molecule.h:578
Patch
Definition: Patch.h:35
ComputeSelfTuples::~ComputeSelfTuples
virtual ~ComputeSelfTuples()
Definition: ComputeSelfTuples.h:355
UniqueSetIter::begin
UniqueSetIter< T > begin(void) const
Definition: UniqueSetIter.h:55
dihedral::atom4
int32 atom4
Definition: structures.h:66
bond::atom1
int32 atom1
Definition: structures.h:48
COMPUTE_HOME_PRIORITY
#define COMPUTE_HOME_PRIORITY
Definition: Priorities.h:76
LocalID::index
int index
Definition: NamdTypes.h:195
angle::atom2
int32 atom2
Definition: structures.h:56
SimParameters::soluteScalingFactor
BigReal soluteScalingFactor
Definition: SimParameters.h:448
ComputeSelfTuples::ComputeSelfTuples
ComputeSelfTuples(ComputeID c, PatchID p)
Definition: ComputeSelfTuples.h:351
angle::atom1
int32 atom1
Definition: structures.h:55
Bool
int Bool
Definition: common.h:133
Molecule::get_ss_type
unsigned char get_ss_type(int anum) const
Definition: Molecule.h:1355
Patch::p
CompAtomList p
Definition: Patch.h:146
AtomMap::localID
LocalID localID(AtomID id)
Definition: AtomMap.h:74
angle::atom3
int32 atom3
Definition: structures.h:57
Molecule::alch_unpert_bonds
Bond * alch_unpert_bonds
Definition: Molecule.h:576
PatchID
int PatchID
Definition: NamdTypes.h:182
ComputeHomeTuples::doWork
virtual void doWork(void)
Definition: ComputeHomeTuples.h:693
LocalID::pid
PatchID pid
Definition: NamdTypes.h:194
dihedral::atom2
int32 atom2
Definition: structures.h:64
Molecule::get_fep_type
unsigned char get_fep_type(int anum) const
Definition: Molecule.h:1349
Node::parameters
Parameters * parameters
Definition: Node.h:177
UniqueSetIter::end
UniqueSetIter< T > end(void) const
Definition: UniqueSetIter.h:64
Patch::getNumAtoms
int getNumAtoms()
Definition: Patch.h:105
ComputeSelfTuples::initialize
virtual void initialize(void)
Definition: ComputeSelfTuples.h:369
TuplePatchElem::xExt
CompAtomExt * xExt
Definition: ComputeHomeTuples.h:46
Node::molecule
Molecule * molecule
Definition: Node.h:176
dihedral::atom1
int32 atom1
Definition: structures.h:63
bond::atom2
int32 atom2
Definition: structures.h:49
BigReal
double BigReal
Definition: common.h:114
PATCH_PRIORITY
#define PATCH_PRIORITY(PID)
Definition: Priorities.h:25
Generated by   doxygen 1.8.5