NAMD
HomePatch.C
Go to the documentation of this file.
1 
8 /*
9  HomePatch owns the actual atoms of a Patch of space
10  Proxy(s) get messages via ProxyMgr from HomePatch(es)
11  to update lists of atoms and their coordinates
12  HomePatch(es) also have a Sequencer bound to them
13 
14  superclass: Patch
15 */
16 
17 #include "time.h"
18 #include <math.h>
19 #include "charm++.h"
20 #include "qd.h"
21 
22 #include "SimParameters.h"
23 #include "HomePatch.h"
24 #include "AtomMap.h"
25 #include "Node.h"
26 #include "PatchMap.inl"
27 #include "main.h"
28 #include "ProxyMgr.decl.h"
29 #include "ProxyMgr.h"
30 #include "Migration.h"
31 #include "Molecule.h"
32 #include "PatchMgr.h"
33 #include "Sequencer.h"
34 #include "Broadcasts.h"
35 #include "LdbCoordinator.h"
36 #include "ReductionMgr.h"
37 #include "Sync.h"
38 #include "Random.h"
39 #include "Priorities.h"
40 #include "ComputeNonbondedUtil.h"
41 #include "ComputeGBIS.inl"
42 #include "Priorities.h"
43 #include "SortAtoms.h"
44 
45 #include "ComputeQM.h"
46 #include "ComputeQMMgr.decl.h"
47 
48 #include "NamdEventsProfiling.h"
49 
50 //#define PRINT_COMP
51 #define TINY 1.0e-20;
52 #define MAXHGS 10
53 #define MIN_DEBUG_LEVEL 2
54 //#define DEBUGM
55 //#define NL_DEBUG
56 #include "Debug.h"
57 
58 #include <vector>
59 #include <algorithm>
60 using namespace std;
61 
62 typedef int HGArrayInt[MAXHGS];
67 
69 
70 int average(CompAtom *qtilde,const HGArrayVector &q,BigReal *lambda,const int n,const int m, const HGArrayBigReal &imass, const HGArrayBigReal &length2, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab, const BigReal tolf, const int ntrial);
71 
72 void mollify(CompAtom *qtilde,const HGArrayVector &q0,const BigReal *lambda, HGArrayVector &force,const int n, const int m, const HGArrayBigReal &imass,const HGArrayInt &ial,const HGArrayInt &ibl,const HGArrayVector &refab);
73 
74 #define MASS_EPSILON (1.0e-35) //a very small floating point number
75 
76 
77 // DMK - Atom Separation (water vs. non-water)
78 #if NAMD_SeparateWaters != 0
79 
80 // Macro to test if a hydrogen group represents a water molecule.
81 // NOTE: This test is the same test in Molecule.C for setting the
82 // OxygenAtom flag in status.
83 // hgtype should be the number of atoms in a water hydrogen group
84 // It must now be set based on simulation parameters because we might
85 // be using tip4p
86 
87 // DJH: This will give false positive for full Drude model,
88 // e.g. O D H is not water but has hgs==3
89 #define IS_HYDROGEN_GROUP_WATER(hgs, mass) \
90  ((hgs >= 3) && ((mass >= 14.0) && (mass <= 18.0)))
91 
92 #endif
93 
94 #ifdef TIMER_COLLECTION
95 const char *TimerSet::tlabel[TimerSet::NUMTIMERS] = {
96  "kick",
97  "maxmove",
98  "drift",
99  "piston",
100  "submithalf",
101  "velbbk1",
102  "velbbk2",
103  "rattle1",
104  "submitfull",
105  "submitcollect",
106 };
107 #endif
108 
109 HomePatch::HomePatch(PatchID pd, FullAtomList &al) : Patch(pd)
110 // DMK - Atom Separation (water vs. non-water)
112  ,tempAtom()
113 #endif
114 {
115  atom.swap(al);
116  settle_initialized = 0;
117 
118  doAtomUpdate = true;
119  rattleListValid = false;
120 
121  exchange_msg = 0;
122  exchange_req = -1;
123 
124  numGBISP1Arrived = 0;
125  numGBISP2Arrived = 0;
126  numGBISP3Arrived = 0;
127  phase1BoxClosedCalled = false;
128  phase2BoxClosedCalled = false;
129  phase3BoxClosedCalled = false;
130 
131  min.x = PatchMap::Object()->min_a(patchID);
132  min.y = PatchMap::Object()->min_b(patchID);
133  min.z = PatchMap::Object()->min_c(patchID);
134  max.x = PatchMap::Object()->max_a(patchID);
135  max.y = PatchMap::Object()->max_b(patchID);
136  max.z = PatchMap::Object()->max_c(patchID);
137  center = 0.5*(min+max);
138 
139  int aAway = PatchMap::Object()->numaway_a();
140  if ( PatchMap::Object()->periodic_a() ||
141  PatchMap::Object()->gridsize_a() > aAway + 1 ) {
142  aAwayDist = (max.x - min.x) * aAway;
143  } else {
144  aAwayDist = Node::Object()->simParameters->patchDimension;
145  }
146  int bAway = PatchMap::Object()->numaway_b();
147  if ( PatchMap::Object()->periodic_b() ||
148  PatchMap::Object()->gridsize_b() > bAway + 1 ) {
149  bAwayDist = (max.y - min.y) * bAway;
150  } else {
151  bAwayDist = Node::Object()->simParameters->patchDimension;
152  }
153  int cAway = PatchMap::Object()->numaway_c();
154  if ( PatchMap::Object()->periodic_c() ||
155  PatchMap::Object()->gridsize_c() > cAway + 1 ) {
156  cAwayDist = (max.z - min.z) * cAway;
157  } else {
158  cAwayDist = Node::Object()->simParameters->patchDimension;
159  }
160 
161  migrationSuspended = false;
162  allMigrationIn = false;
163  marginViolations = 0;
164  patchMapRead = 0; // We delay read of PatchMap data
165  // to make sure it is really valid
166  inMigration = false;
167  numMlBuf = 0;
168  flags.sequence = -1;
169  flags.maxForceUsed = -1;
170 
171  numAtoms = atom.size();
172  replacementForces = 0;
173 
175  doPairlistCheck_newTolerance =
176  0.5 * ( simParams->pairlistDist - simParams->cutoff );
177 
178 
179  numFixedAtoms = 0;
180  if ( simParams->fixedAtomsOn ) {
181  for ( int i = 0; i < numAtoms; ++i ) {
182  numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
183  }
184  }
185 
186 #ifdef NODEAWARE_PROXY_SPANNINGTREE
187  ptnTree.resize(0);
188  /*children = NULL;
189  numChild = 0;*/
190 #else
191  child = new int[proxySpanDim];
192  nChild = 0; // number of proxy spanning tree children
193 #endif
194 
195 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
196  nphs = 0;
197  localphs = NULL;
198  isProxyChanged = 0;
199 #endif
200 
201 
202  // DMK - Atom Separation (water vs. non-water)
203  #if NAMD_SeparateWaters != 0
204 
205  // Create the scratch memory for separating atoms
206  tempAtom.resize(numAtoms);
207  numWaterAtoms = -1;
208 
209  // Separate the current list of atoms
210  separateAtoms();
211 
212  #endif
213 
214  // Handle unusual water models here
215  if (simParams->watmodel == WAT_TIP4) init_tip4();
216  else if (simParams->watmodel == WAT_SWM4) init_swm4();
217 
218  isNewProxyAdded = 0;
219 }
220 
221 void HomePatch::write_tip4_props() {
222  printf("Writing r_om and r_ohc: %f | %f\n", r_om, r_ohc);
223 }
224 
225 void HomePatch::init_tip4() {
226  // initialize the distances needed for the tip4p water model
227  Molecule *mol = Node::Object()->molecule;
228  r_om = mol->r_om;
229  r_ohc = mol->r_ohc;
230 }
231 
232 
233 void ::HomePatch::init_swm4() {
234  // initialize the distances needed for the SWM4 water model
235  Molecule *mol = Node::Object()->molecule;
236  r_om = mol->r_om;
237  r_ohc = mol->r_ohc;
238 }
239 
240 
241 void HomePatch::reinitAtoms(FullAtomList &al) {
242  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
243 
244  atom.swap(al);
245  numAtoms = atom.size();
246 
247  doAtomUpdate = true;
248  rattleListValid = false;
249 
250  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
251 
252  // DMK - Atom Separation (water vs. non-water)
253  #if NAMD_SeparateWaters != 0
254 
255  // Reset the numWaterAtoms value
256  numWaterAtoms = -1;
257 
258  // Separate the atoms
259  separateAtoms();
260 
261  #endif
262 }
263 
264 // Bind a Sequencer to this HomePatch
266 { sequencer=sequencerPtr; }
267 
268 // start simulation over this Patch of atoms
270 { sequencer->run(); }
271 
272 void HomePatch::readPatchMap() {
273  // iout << "Patch " << patchID << " has " << proxy.size() << " proxies.\n" << endi;
275  PatchID nnPatchID[PatchMap::MaxOneAway];
276 
277  patchMigrationCounter = numNeighbors
278  = PatchMap::Object()->oneAwayNeighbors(patchID, nnPatchID);
279  DebugM( 1, "NumNeighbors for pid " <<patchID<<" is "<< numNeighbors << "\n");
280  int n;
281  for (n=0; n<numNeighbors; n++) {
282  realInfo[n].destNodeID = p->node(realInfo[n].destPatchID = nnPatchID[n]);
283  DebugM( 1, " nnPatchID=" <<nnPatchID[n]<<" nnNodeID="<< realInfo[n].destNodeID<< "\n");
284  realInfo[n].mList.resize(0);
285  }
286 
287  // Make mapping from the 3x3x3 cube of pointers to real migration info
288  for (int i=0; i<3; i++)
289  for (int j=0; j<3; j++)
290  for (int k=0; k<3; k++)
291  {
292  int pid = p->pid(p->index_a(patchID)+i-1,
293  p->index_b(patchID)+j-1, p->index_c(patchID)+k-1);
294  if (pid < 0) {
295  DebugM(5, "ERROR, for patchID " << patchID <<" I got neigh pid = " << pid << "\n");
296  }
297  if (pid == patchID && ! (
298  ( (i-1) && p->periodic_a() ) ||
299  ( (j-1) && p->periodic_b() ) ||
300  ( (k-1) && p->periodic_c() ) )) {
301  mInfo[i][j][k] = NULL;
302  }
303  else {
304  // Does not work as expected for periodic with only two patches.
305  // Also need to check which image we want, but OK for now. -JCP
306  for (n = 0; n<numNeighbors; n++) {
307  if (pid == realInfo[n].destPatchID) {
308  mInfo[i][j][k] = &realInfo[n];
309  break;
310  }
311  }
312  if (n == numNeighbors) { // disaster!
313  DebugM(4,"BAD News, I could not find PID " << pid << "\n");
314  }
315  }
316  }
317 
318  DebugM(1,"Patch("<<patchID<<") # of neighbors = " << numNeighbors << "\n");
319 }
320 
322 {
323  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
324 #ifdef NODEAWARE_PROXY_SPANNINGTREE
325  ptnTree.resize(0);
326  #ifdef USE_NODEPATCHMGR
327  delete [] nodeChildren;
328  #endif
329 #endif
330  delete [] child;
331 }
332 
333 
334 void HomePatch::boxClosed(int box) {
335  // begin gbis
336  if (box == 5) {// end of phase 1
337  phase1BoxClosedCalled = true;
338  if (!psiSumBox.isOpen() && numGBISP1Arrived == proxy.size()) {
339  if (flags.doGBIS && flags.doNonbonded) {
340  sequencer->awaken();
341  }
342  } else {
343  //need to wait until proxies arrive before awakening
344  }
345  } else if (box == 6) {// intRad
346  //do nothing
347  } else if (box == 7) {// bornRad
348  //do nothing
349  } else if (box == 8) {// end of phase 2
350  phase2BoxClosedCalled = true;
351  //if no proxies, AfterP1 can't be called from receive
352  //so it will be called from here
353  if (!dEdaSumBox.isOpen() && numGBISP2Arrived == proxy.size()) {
354  if (flags.doGBIS && flags.doNonbonded) {
355  sequencer->awaken();
356  }
357  } else {
358  //need to wait until proxies arrive before awakening
359  }
360  } else if (box == 9) {
361  //do nothing
362  } else if (box == 10) {
363  //lcpoType Box closed: do nothing
364  } else {
365  //do nothing
366  }
367  // end gbis
368 
369  if ( ! --boxesOpen ) {
370  if ( replacementForces ) {
371  for ( int i = 0; i < numAtoms; ++i ) {
372  if ( replacementForces[i].replace ) {
373  for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
374  f[Results::normal][i] = replacementForces[i].force;
375  }
376  }
377  replacementForces = 0;
378  }
379  DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
380  << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
381  // only awaken suspended threads. Then say it is suspended.
382 
383  phase3BoxClosedCalled = true;
384  if (flags.doGBIS) {
385  if (flags.doNonbonded) {
386  sequencer->awaken();
387  } else {
388  if (numGBISP1Arrived == proxy.size() &&
389  numGBISP2Arrived == proxy.size() &&
390  numGBISP3Arrived == proxy.size()) {
391  sequencer->awaken();//all boxes closed and all proxies arrived
392  }
393  }
394  } else {//non-gbis awaken
395  sequencer->awaken();
396  }
397  } else {
398  DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
399  }
400 }
401 
403  DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
404  proxy.add(msg->node);
406 
407  isNewProxyAdded = 1;
408 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
409  isProxyChanged = 1;
410 #endif
411 
412  Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
413  delete msg;
414 }
415 
417 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
418  isProxyChanged = 1;
419 #endif
420  int n = msg->node;
421  NodeID *pe = proxy.begin();
422  for ( ; *pe != n; ++pe );
424  proxy.del(pe - proxy.begin());
425  delete msg;
426 }
427 
428 #if USE_TOPOMAP && USE_SPANNING_TREE
429 
430 int HomePatch::findSubroots(int dim, int* subroots, int psize, int* pidscopy){
431  int nChild = 0;
432  int cones[6][proxySpanDim*proxySpanDim+proxySpanDim];
433  int conesizes[6] = {0,0,0,0,0,0};
434  int conecounters[6] = {0,0,0,0,0,0};
435  int childcounter = 0;
436  nChild = (psize>proxySpanDim)?proxySpanDim:psize;
437  TopoManager tmgr;
438  for(int i=0;i<psize;i++){
439  int cone = tmgr.getConeNumberForRank(pidscopy[i]);
440  cones[cone][conesizes[cone]++] = pidscopy[i];
441  }
442 
443  while(childcounter<nChild){
444  for(int i=0;i<6;i++){
445  if(conecounters[i]<conesizes[i]){
446  subroots[childcounter++] = cones[i][conecounters[i]++];
447  }
448  }
449  }
450  for(int i=nChild;i<proxySpanDim;i++)
451  subroots[i] = -1;
452  return nChild;
453 }
454 #endif // USE_TOPOMAP
455 
456 static int compDistance(const void *a, const void *b)
457 {
458  int d1 = abs(*(int *)a - CkMyPe());
459  int d2 = abs(*(int *)b - CkMyPe());
460  if (d1 < d2)
461  return -1;
462  else if (d1 == d2)
463  return 0;
464  else
465  return 1;
466 }
467 
469 {
470 #if USE_NODEPATCHMGR
471  CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
472  NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
473  npm->sendProxyList(patchID, proxy.begin(), proxy.size());
474 #else
475  ProxyMgr::Object()->sendProxies(patchID, proxy.begin(), proxy.size());
476 #endif
477 }
478 
479 #ifdef NODEAWARE_PROXY_SPANNINGTREE
480 void HomePatch::buildNodeAwareSpanningTree(void){
481 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
482  DebugFileTrace *dft = DebugFileTrace::Object();
483  dft->openTrace();
484  dft->writeTrace("HomePatch[%d] has %d proxy on proc[%d] node[%d]\n", patchID, proxy.size(), CkMyPe(), CkMyNode());
485  dft->writeTrace("Proxies are: ");
486  for(int i=0; i<proxy.size(); i++) dft->writeTrace("%d(%d), ", proxy[i], CkNodeOf(proxy[i]));
487  dft->writeTrace("\n");
488  dft->closeTrace();
489 #endif
490 
491  //build the naive spanning tree for this home patch
492  if(! proxy.size()) {
493  //this case will not happen in practice.
494  //In debugging state where spanning tree is enforced, then this could happen
495  //Chao Mei
496  return;
497  }
498  ProxyMgr::buildSinglePatchNodeAwareSpanningTree(patchID, proxy, ptnTree);
499  //optimize on the naive spanning tree
500 
501  //setup the children
502  setupChildrenFromProxySpanningTree();
503  //send down to children
505 }
506 
507 void HomePatch::setupChildrenFromProxySpanningTree(){
508 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
509  isProxyChanged = 1;
510 #endif
511  if(ptnTree.size()==0) {
512  nChild = 0;
513  delete [] child;
514  child = NULL;
515  #ifdef USE_NODEPATCHMGR
516  numNodeChild = 0;
517  delete [] nodeChildren;
518  nodeChildren = NULL;
519  #endif
520  return;
521  }
522  proxyTreeNode *rootnode = &ptnTree.item(0);
523  CmiAssert(rootnode->peIDs[0] == CkMyPe());
524  //set up children
525  //1. add external children (the first proc inside the proxy tree node)
526  //2. add internal children (with threshold that would enable spanning
527  int internalChild = rootnode->numPes-1;
528  int externalChild = ptnTree.size()-1;
529  externalChild = (proxySpanDim>externalChild)?externalChild:proxySpanDim;
530  int internalSlots = proxySpanDim-externalChild;
531  if(internalChild>0){
532  if(internalSlots==0) {
533  //at least having one internal child
534  internalChild = 1;
535  }else{
536  internalChild = (internalSlots>internalChild)?internalChild:internalSlots;
537  }
538  }
539 
540  nChild = externalChild+internalChild;
541  CmiAssert(nChild>0);
542 
543  //exclude the root node
544  delete [] child;
545  child = new int[nChild];
546 
547  for(int i=0; i<externalChild; i++) {
548  child[i] = ptnTree.item(i+1).peIDs[0];
549  }
550  for(int i=externalChild, j=1; i<nChild; i++, j++) {
551  child[i] = rootnode->peIDs[j];
552  }
553 
554 #ifdef USE_NODEPATCHMGR
555  //only register the cores that have proxy patches. The HomePach's core
556  //doesn't need to be registered.
557  CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
558  NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
559  if(rootnode->numPes==1){
560  npm->registerPatch(patchID, 0, NULL);
561  }
562  else{
563  npm->registerPatch(patchID, rootnode->numPes-1, &rootnode->peIDs[1]);
564  }
565 
566  //set up childrens in terms of node ids
567  numNodeChild = externalChild;
568  if(internalChild) numNodeChild++;
569  delete [] nodeChildren;
570  nodeChildren = new int[numNodeChild];
571  for(int i=0; i<externalChild; i++) {
572  nodeChildren[i] = ptnTree.item(i+1).nodeID;
573  }
574  //the last entry always stores this node id if there are proxies
575  //on other cores of the same node for this patch.
576  if(internalChild)
577  nodeChildren[numNodeChild-1] = rootnode->nodeID;
578 #endif
579 
580 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
581  DebugFileTrace *dft = DebugFileTrace::Object();
582  dft->openTrace();
583  dft->writeTrace("HomePatch[%d] has %d children: ", patchID, nChild);
584  for(int i=0; i<nChild; i++)
585  dft->writeTrace("%d ", child[i]);
586  dft->writeTrace("\n");
587  dft->closeTrace();
588 #endif
589 }
590 #endif
591 
592 #ifdef NODEAWARE_PROXY_SPANNINGTREE
593 //This is not an entry method, but takes an argument of message type
595  //set up the whole tree ptnTree
596  int treesize = msg->numNodesWithProxies;
597  ptnTree.resize(treesize);
598  int *pAllPes = msg->allPes;
599  for(int i=0; i<treesize; i++) {
600  proxyTreeNode *oneNode = &ptnTree.item(i);
601  delete [] oneNode->peIDs;
602  oneNode->numPes = msg->numPesOfNode[i];
603  oneNode->nodeID = CkNodeOf(*pAllPes);
604  oneNode->peIDs = new int[oneNode->numPes];
605  for(int j=0; j<oneNode->numPes; j++) {
606  oneNode->peIDs[j] = *pAllPes;
607  pAllPes++;
608  }
609  }
610  //setup children
611  setupChildrenFromProxySpanningTree();
612  //send down to children
614 }
615 
617  if(ptnTree.size()==0) return;
619  ProxyNodeAwareSpanningTreeMsg::getANewMsg(patchID, CkMyPe(), ptnTree.begin(), ptnTree.size());
620 
621  #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
622  msg->printOut("HP::sendST");
623  #endif
624 
625  CmiAssert(CkMyPe() == msg->allPes[0]);
627 
628 }
629 #else
632 #endif
633 
634 #ifndef NODEAWARE_PROXY_SPANNINGTREE
635 // recv a spanning tree from load balancer
636 void HomePatch::recvSpanningTree(int *t, int n)
637 {
638  int i;
639  nChild=0;
640  tree.resize(n);
641  for (i=0; i<n; i++) {
642  tree[i] = t[i];
643  }
644 
645  for (i=1; i<=proxySpanDim; i++) {
646  if (tree.size() <= i) break;
647  child[i-1] = tree[i];
648  nChild++;
649  }
650 
651 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
652  isProxyChanged = 1;
653 #endif
654 
655  // send down to children
657 }
658 
660 {
661  if (tree.size() == 0) return;
663  msg->patch = patchID;
664  msg->node = CkMyPe();
665  msg->tree.copy(tree); // copy data for thread safety
667 }
668 #else
669 void HomePatch::recvSpanningTree(int *t, int n){}
671 #endif
672 
673 #ifndef NODEAWARE_PROXY_SPANNINGTREE
675 {
676  nChild = 0;
677  int psize = proxy.size();
678  if (psize == 0) return;
679  NodeIDList oldtree; oldtree.swap(tree);
680  int oldsize = oldtree.size();
681  tree.resize(psize + 1);
682  tree.setall(-1);
683  tree[0] = CkMyPe();
684  int s=1, e=psize+1;
686  int patchNodesLast =
687  ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
688  int nNonPatch = 0;
689 #if 1
690  // try to put it to the same old tree
691  for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
692  {
693  int oldindex = oldtree.find(*pli);
694  if (oldindex != -1 && oldindex < psize) {
695  tree[oldindex] = *pli;
696  }
697  }
698  s=1; e=psize;
699  for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
700  {
701  if (tree.find(*pli) != -1) continue; // already assigned
702  if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
703  while (tree[e] != -1) { e--; if (e==-1) e = psize; }
704  tree[e] = *pli;
705  } else {
706  while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
707  tree[s] = *pli;
708  nNonPatch++;
709  }
710  }
711 #if 1
712  if (oldsize==0 && nNonPatch) {
713  // first time, sort by distance
714  qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
715  }
716 #endif
717 
718  //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
719 
720 #if USE_TOPOMAP && USE_SPANNING_TREE
721 
722  //Right now only works for spanning trees with two levels
723  int *treecopy = new int [psize];
724  int subroots[proxySpanDim];
725  int subsizes[proxySpanDim];
726  int subtrees[proxySpanDim][proxySpanDim];
727  int idxes[proxySpanDim];
728  int i = 0;
729 
730  for(i=0;i<proxySpanDim;i++){
731  subsizes[i] = 0;
732  idxes[i] = i;
733  }
734 
735  for(i=0;i<psize;i++){
736  treecopy[i] = tree[i+1];
737  }
738 
739  TopoManager tmgr;
740  tmgr.sortRanksByHops(treecopy,nNonPatch);
741  tmgr.sortRanksByHops(treecopy+nNonPatch,
742  psize-nNonPatch);
743 
744  /* build tree and subtrees */
745  nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
746  delete [] treecopy;
747 
748  for(int i=1;i<psize+1;i++){
749  int isSubroot=0;
750  for(int j=0;j<nChild;j++)
751  if(tree[i]==subroots[j]){
752  isSubroot=1;
753  break;
754  }
755  if(isSubroot) continue;
756 
757  int bAdded = 0;
758  tmgr.sortIndexByHops(tree[i], subroots,
759  idxes, proxySpanDim);
760  for(int j=0;j<proxySpanDim;j++){
761  if(subsizes[idxes[j]]<proxySpanDim){
762  subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
763  bAdded = 1;
764  break;
765  }
766  }
767  if( psize > proxySpanDim && ! bAdded ) {
768  NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
769  }
770  }
771 
772 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
773 
774  for (int i=1; i<=proxySpanDim; i++) {
775  if (tree.size() <= i) break;
776  child[i-1] = tree[i];
777  nChild++;
778  }
779 #endif
780 #endif
781 
782 #if 0
783  // for debugging
784  CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
785  for (int i=0; i<psize+1; i++) {
786  CkPrintf("%d ", tree[i]);
787  }
788  CkPrintf("\n");
789 #endif
790  // send to children nodes
792 }
793 #endif
794 
795 
797 
798  numGBISP3Arrived++;
799  DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
800  int n = msg->node;
801  Results *r = forceBox.clientOpen();
802 
803  char *iszeroPtr = msg->isZero;
804  Force *msgFPtr = msg->forceArr;
805 
806  for ( int k = 0; k < Results::maxNumForces; ++k )
807  {
808  Force *rfPtr = r->f[k];
809  for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
810  if((*iszeroPtr)!=1) {
811  *rfPtr += *msgFPtr;
812  msgFPtr++;
813  }
814  }
815  }
817  delete msg;
818 }
819 
821  numGBISP3Arrived++;
822  DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
823  int n = msg->node;
824  Results *r = forceBox.clientOpen();
825  for ( int k = 0; k < Results::maxNumForces; ++k )
826  {
827  Force *f = r->f[k];
828  register ForceList::iterator f_i, f_e;
829  f_i = msg->forceList[k]->begin();
830  f_e = msg->forceList[k]->end();
831  for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
832  }
834  delete msg;
835 }
836 
838 {
839  numGBISP3Arrived++;
840  DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodeSize<<")\n");
841  Results *r = forceBox.clientOpen(msg->nodeSize);
842  register char* isNonZero = msg->isForceNonZero;
843  register Force* f_i = msg->forceArr;
844  for ( int k = 0; k < Results::maxNumForces; ++k )
845  {
846  Force *f = r->f[k];
847  int nf = msg->flLen[k];
848 #ifdef ARCH_POWERPC
849 #pragma disjoint (*f_i, *f)
850 #endif
851  for (int count = 0; count < nf; count++) {
852  if(*isNonZero){
853  f[count].x += f_i->x;
854  f[count].y += f_i->y;
855  f[count].z += f_i->z;
856  f_i++;
857  }
858  isNonZero++;
859  }
860  }
862 
863  delete msg;
864 }
865 
867 {
868  // This is used for LSS in QM/MM simulations.
869  // Changes atom labels so that we effectively exchange solvent
870  // molecules between classical and quantum modes.
871 
872  SimParameters *simParams = Node::Object()->simParameters;
873  int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
874  const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
875  const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
876  Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
877 
878  ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
879  CkpvAccess(BOCclass_group).computeQMMgr) ;
880 
881  FullAtom *a_i = atom.begin();
882 
883  for (int i=0; i<numAtoms; ++i ) {
884 
885  LSSSubsDat *subP = lssSubs(mgrP).find( LSSSubsDat(a_i[i].id) ) ;
886 
887  if ( subP != NULL ) {
888  a_i[i].id = subP->newID ;
889  a_i[i].vdwType = subP->newVdWType ;
890 
891  // If we are swappign a classical atom with a QM one, the charge
892  // will need extra handling.
893  if (qmAtomGroup[subP->newID] > 0 && simParams->PMEOn) {
894  // We make sure that the last atom charge calculated for the
895  // QM atom being transfered here is available for PME
896  // in the next step.
897 
898  // Loops over all QM atoms (in all QM groups) comparing their
899  // global indices (the sequential atom ID from NAMD).
900  for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
901 
902  if (qmAtmIndx[qmIter] == subP->newID) {
903  qmAtmChrg[qmIter] = subP->newCharge;
904  break;
905  }
906 
907  }
908 
909  // For QM atoms, the charge in the full atom structure is zero.
910  // Electrostatic interactions between QM atoms and their
911  // environment is handled in ComputeQM.
912  a_i[i].charge = 0;
913  }
914  else {
915  // If we are swappign a QM atom with a Classical one, only the charge
916  // in the full atomstructure needs updating, since it used to be zero.
917  a_i[i].charge = subP->newCharge ;
918  }
919  }
920  }
921 
922  return;
923 }
924 
925 void HomePatch::positionsReady(int doMigration)
926 {
927  SimParameters *simParams = Node::Object()->simParameters;
928 
929  flags.sequence++;
930 
931  if (!patchMapRead) { readPatchMap(); }
932 
933  if (numNeighbors && ! simParams->staticAtomAssignment) {
934  if (doMigration) {
935  rattleListValid = false;
936  doAtomMigration();
937  } else {
938  doMarginCheck();
939  }
940  }
941 
942 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
943  char prbuf[32];
944  sprintf(prbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::POSITIONS_READY], this->getPatchID());
945  NAMD_EVENT_START_EX(1, NamdProfileEvent::POSITIONS_READY, prbuf);
946 #endif
947 
948  if (doMigration && simParams->qmLSSOn)
949  qmSwapAtoms();
950 
951 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
952  #ifdef NAMD_AVXTILES
953  if ( simParams->useAVXTiles ) {
954  #endif
955  if ( doMigration ) {
956  int n = numAtoms;
957  FullAtom *a_i = atom.begin();
958  #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_AVXTILES) || \
959  (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
960  int *ao = new int[n];
961  int nfree;
962  if ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces ) {
963  int k = 0;
964  int k2 = n;
965  for ( int j=0; j<n; ++j ) {
966  // put fixed atoms at end
967  if ( a_i[j].atomFixed ) ao[--k2] = j;
968  else ao[k++] = j;
969  }
970  nfree = k;
971  } else {
972  nfree = n;
973  for ( int j=0; j<n; ++j ) {
974  ao[j] = j;
975  }
976  }
977 
978  sortAtomsForCUDA(ao,a_i,nfree,n);
979 
980  for ( int i=0; i<n; ++i ) {
981  a_i[i].sortOrder = ao[i];
982  }
983  delete [] ao;
984  #else
985  for (int i = 0; i < n; ++i) {
986  a_i[i].sortOrder = i;
987  }
988  #endif
989  }
990 
991  {
992  const double charge_scaling = sqrt(COULOMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1);
993  const Vector ucenter = lattice.unscale(center);
994  const BigReal ucenter_x = ucenter.x;
995  const BigReal ucenter_y = ucenter.y;
996  const BigReal ucenter_z = ucenter.z;
997  const int n = numAtoms;
998  #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
999  int n_16 = n;
1000  n_16 = (n + 15) & (~15);
1001  cudaAtomList.resize(n_16);
1002  CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
1003  mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
1004  mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
1005  mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
1006  mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
1007  #elif defined(NAMD_AVXTILES)
1008  int n_avx = (n + 15) & (~15);
1009  cudaAtomList.resize(n_avx);
1010  CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
1011  tiles.realloc(n, ac);
1012  #else
1013  cudaAtomList.resize(n);
1014  CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
1015  #endif
1016  const FullAtom *a = atom.begin();
1017  for ( int k=0; k<n; ++k ) {
1018  #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
1019  int j = a[k].sortOrder;
1020  atom_x[k] = a[j].position.x - ucenter_x;
1021  atom_y[k] = a[j].position.y - ucenter_y;
1022  atom_z[k] = a[j].position.z - ucenter_z;
1023  atom_q[k] = charge_scaling * a[j].charge;
1024  #else
1025  int j = a[k].sortOrder;
1026  ac[k].x = a[j].position.x - ucenter_x;
1027  ac[k].y = a[j].position.y - ucenter_y;
1028  ac[k].z = a[j].position.z - ucenter_z;
1029  ac[k].q = charge_scaling * a[j].charge;
1030  #endif
1031  }
1032  #ifdef NAMD_AVXTILES
1033  {
1034  if (n > 0) {
1035  int j = a[n-1].sortOrder;
1036  for ( int k=n; k<n_avx; ++k ) {
1037  ac[k].x = a[j].position.x - ucenter_x;
1038  ac[k].y = a[j].position.y - ucenter_y;
1039  ac[k].z = a[j].position.z - ucenter_z;
1040  }
1041  }
1042  }
1043  #endif
1044  }
1045  #ifdef NAMD_AVXTILES
1046  // If "Tiles" mode disabled, no use of the CUDA data structures
1047  } else doMigration = doMigration && numNeighbors;
1048  #endif
1049 #else
1050  doMigration = doMigration && numNeighbors;
1051 #endif
1052  doMigration = doMigration || ! patchMapRead;
1053 
1054  doMigration = doMigration || doAtomUpdate;
1055  doAtomUpdate = false;
1056 
1057  // Workaround for oversize groups:
1058  // reset nonbondedGroupSize (ngs) before force calculation,
1059  // making sure that subset of hydrogen group starting with
1060  // parent atom are all within 0.5 * hgroupCutoff.
1061  // XXX hydrogentGroupSize remains constant but is checked for nonzero
1062  // XXX should be skipped for CUDA, ngs not used by CUDA kernels
1063  // XXX should this also be skipped for KNL kernels?
1064  // ngs used by ComputeNonbondedBase.h - CPU nonbonded kernels
1065  // ngs used by ComputeGBIS.C - CPU GB nonbonded kernels
1066 
1067 #if ! (defined(NAMD_CUDA) || defined(NAMD_HIP))
1068 #if defined(NAMD_AVXTILES)
1069  if (!simParams->useAVXTiles)
1070 #endif
1071  doGroupSizeCheck();
1072 #endif
1073 
1074  // Copy information needed by computes and proxys to Patch::p.
1075  // Resize only if atoms were migrated
1076  if (doMigration) {
1077  p.resize(numAtoms);
1078  pExt.resize(numAtoms);
1079  }
1080  CompAtom *p_i = p.begin();
1081  CompAtomExt *pExt_i = pExt.begin();
1082  FullAtom *a_i = atom.begin();
1083  int i; int n = numAtoms;
1084  for ( i=0; i<n; ++i ) {
1085  p_i[i] = a_i[i];
1086  pExt_i[i] = a_i[i];
1087  }
1088 
1089  // Measure atom movement to test pairlist validity
1090  doPairlistCheck();
1091 
1092  if (flags.doMolly) mollyAverage();
1093  // BEGIN LA
1095  // END LA
1096 
1097  if (flags.doGBIS) {
1098  //reset for next time step
1099  numGBISP1Arrived = 0;
1100  phase1BoxClosedCalled = false;
1101  numGBISP2Arrived = 0;
1102  phase2BoxClosedCalled = false;
1103  numGBISP3Arrived = 0;
1104  phase3BoxClosedCalled = false;
1105  if (doMigration || isNewProxyAdded)
1107  }
1108 
1109  if (flags.doLCPO) {
1110  if (doMigration || isNewProxyAdded) {
1111  setLcpoType();
1112  }
1113  }
1114 
1115  // Must Add Proxy Changes when migration completed!
1117  int *pids = NULL;
1118  int pidsPreAllocated = 1;
1119  int npid;
1120  if (proxySendSpanning == 0) {
1121  npid = proxy.size();
1122  pids = new int[npid];
1123  pidsPreAllocated = 0;
1124  int *pidi = pids;
1125  int *pide = pids + proxy.size();
1126  int patchNodesLast =
1127  ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
1128  for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
1129  {
1130  if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
1131  *(--pide) = *pli;
1132  } else {
1133  *(pidi++) = *pli;
1134  }
1135  }
1136  }
1137  else {
1138 #ifdef NODEAWARE_PROXY_SPANNINGTREE
1139  #ifdef USE_NODEPATCHMGR
1140  npid = numNodeChild;
1141  pids = nodeChildren;
1142  #else
1143  npid = nChild;
1144  pids = child;
1145  #endif
1146 #else
1147  npid = nChild;
1148  pidsPreAllocated = 0;
1149  pids = new int[proxySpanDim];
1150  for (int i=0; i<nChild; i++) pids[i] = child[i];
1151 #endif
1152  }
1153  if (npid) { //have proxies
1154  int seq = flags.sequence;
1155  int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
1156  //begin to prepare proxy msg and send it
1157  int pdMsgPLLen = p.size();
1158  int pdMsgAvgPLLen = 0;
1159  if(flags.doMolly) {
1160  pdMsgAvgPLLen = p_avg.size();
1161  }
1162  // BEGIN LA
1163  int pdMsgVLLen = 0;
1164  if (flags.doLoweAndersen) {
1165  pdMsgVLLen = v.size();
1166  }
1167  // END LA
1168 
1169  int intRadLen = 0;
1170  if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
1171  intRadLen = numAtoms * 2;
1172  }
1173 
1174  //LCPO
1175  int lcpoTypeLen = 0;
1176  if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
1177  lcpoTypeLen = numAtoms;
1178  }
1179 
1180  int pdMsgPLExtLen = 0;
1181  if(doMigration || isNewProxyAdded) {
1182  pdMsgPLExtLen = pExt.size();
1183  }
1184 
1185  int cudaAtomLen = 0;
1186 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
1187  cudaAtomLen = numAtoms;
1188 #elif defined(NAMD_AVXTILES)
1189  if (simParams->useAVXTiles)
1190  cudaAtomLen = (numAtoms + 15) & (~15);
1191 #elif defined(NAMD_MIC)
1192  #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
1193  cudaAtomLen = numAtoms;
1194  #else
1195  cudaAtomLen = (numAtoms + 15) & (~15);
1196  #endif
1197 #endif
1198 
1199  ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
1200  lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
1201 
1202  SET_PRIORITY(nmsg,seq,priority);
1203  nmsg->patch = patchID;
1204  nmsg->flags = flags;
1205  nmsg->plLen = pdMsgPLLen;
1206  //copying data to the newly created msg
1207  memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
1208  nmsg->avgPlLen = pdMsgAvgPLLen;
1209  if(flags.doMolly) {
1210  memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
1211  }
1212  // BEGIN LA
1213  nmsg->vlLen = pdMsgVLLen;
1214  if (flags.doLoweAndersen) {
1215  memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
1216  }
1217  // END LA
1218 
1219  if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
1220  for (int i = 0; i < numAtoms * 2; i++) {
1221  nmsg->intRadList[i] = intRad[i];
1222  }
1223  }
1224 
1225  if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
1226  for (int i = 0; i < numAtoms; i++) {
1227  nmsg->lcpoTypeList[i] = lcpoType[i];
1228  }
1229  }
1230 
1231  nmsg->plExtLen = pdMsgPLExtLen;
1232  if(doMigration || isNewProxyAdded){
1233  memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
1234  }
1235 
1236 // DMK
1237 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
1238  #ifdef NAMD_AVXTILES
1239  if (simParams->useAVXTiles)
1240  #endif
1241  memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
1242 #endif
1243 
1244 #if NAMD_SeparateWaters != 0
1245  //DMK - Atom Separation (water vs. non-water)
1246  nmsg->numWaterAtoms = numWaterAtoms;
1247 #endif
1248 
1249 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
1250  nmsg->isFromImmMsgCall = 0;
1251 #endif
1252 
1253  #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
1254  DebugFileTrace *dft = DebugFileTrace::Object();
1255  dft->openTrace();
1256  dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
1257  for(int i=0; i<npid; i++) {
1258  dft->writeTrace("%d ", pids[i]);
1259  }
1260  dft->writeTrace("\n");
1261  dft->closeTrace();
1262  #endif
1263 
1264 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1265  if (isProxyChanged || localphs == NULL)
1266  {
1267 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
1268  //CmiAssert(isProxyChanged);
1269  if (nphs) {
1270  for (int i=0; i<nphs; i++) {
1271  CmiDestoryPersistent(localphs[i]);
1272  }
1273  delete [] localphs;
1274  }
1275  localphs = new PersistentHandle[npid];
1276  int persist_size = sizeof(envelope) + sizeof(ProxyDataMsg) + sizeof(CompAtom)*(pdMsgPLLen+pdMsgAvgPLLen+pdMsgVLLen) + intRadLen*sizeof(Real) + lcpoTypeLen*sizeof(int) + sizeof(CompAtomExt)*pdMsgPLExtLen + sizeof(CudaAtom)*cudaAtomLen + PRIORITY_SIZE/8 + 2048;
1277  for (int i=0; i<npid; i++) {
1278 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
1279  if (proxySendSpanning)
1280  localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
1281  else
1282 #endif
1283  localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
1284  }
1285  nphs = npid;
1286  }
1287  CmiAssert(nphs == npid && localphs != NULL);
1288  CmiUsePersistentHandle(localphs, nphs);
1289 #endif
1290  if(doMigration || isNewProxyAdded) {
1291  ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
1292  }else{
1293  ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
1294  }
1295 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1296  CmiUsePersistentHandle(NULL, 0);
1297 #endif
1298  isNewProxyAdded = 0;
1299  }
1300  isProxyChanged = 0;
1301  if(!pidsPreAllocated) delete [] pids;
1302  DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
1303 
1304 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
1305  positionPtrBegin = p.begin();
1306  positionPtrEnd = p.end();
1307 #endif
1308 
1309  if(flags.doMolly) {
1312  }
1313  // BEGIN LA
1314  if (flags.doLoweAndersen) {
1315  velocityPtrBegin = v.begin();
1316  velocityPtrEnd = v.end();
1317  }
1318  // END LA
1319 
1320  Patch::positionsReady(doMigration);
1321 
1322  patchMapRead = 1;
1323 
1324  // gzheng
1325  Sync::Object()->PatchReady();
1326 
1327  NAMD_EVENT_STOP(1, NamdProfileEvent::POSITIONS_READY);
1328 
1329 }
1330 
1332 {
1333  replacementForces = f;
1334 }
1335 
1336 void HomePatch::saveForce(const int ftag)
1337 {
1338  f_saved[ftag].resize(numAtoms);
1339  for ( int i = 0; i < numAtoms; ++i )
1340  {
1341  f_saved[ftag][i] = f[ftag][i];
1342  }
1343 }
1344 
1345 
1346 #undef DEBUG_REDISTRIB_FORCE
1347 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
1348 //#define DEBUG_REDISTRIB_FORCE
1349 
1364 void HomePatch::redistrib_colinear_lp_force(
1365  Vector& fi, Vector& fj, Vector& fk,
1366  const Vector& ri, const Vector& rj, const Vector& rk,
1367  Real distance_f, Real scale_f) {
1368  BigReal distance = distance_f;
1369  BigReal scale = scale_f;
1370  Vector r_jk = rj - rk;
1371  // TODO: Add a check for small distances?
1372  BigReal r_jk_rlength = r_jk.rlength();
1373  distance *= r_jk_rlength;
1374  BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
1375  fj += (1. + scale + distance)*fi - r_jk*fdot;
1376  fk -= (scale + distance)*fi + r_jk*fdot;
1377 }
1378 
1404 #define FIX_FOR_WATER
1405 void HomePatch::redistrib_relative_lp_force(
1406  Vector& fi, Vector& fj, Vector& fk, Vector& fl,
1407  const Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
1408  Tensor *virial, int midpt) {
1409 #ifdef DEBUG_REDISTRIB_FORCE
1410  Vector foldnet, toldnet; // old net force, old net torque
1411  foldnet = fi + fj + fk + fl;
1412  toldnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
1413 #endif
1414  Vector fja(0), fka(0), fla(0);
1415 
1416  Vector r = ri - rj;
1417  BigReal invr2 = r.rlength();
1418  invr2 *= invr2;
1419  BigReal fdot = (fi*r) * invr2;
1420  Vector fr = r * fdot;
1421 
1422  fja += fr;
1423 
1424  Vector s, t;
1425  if (midpt) {
1426  s = rj - 0.5*(rk + rl);
1427  t = 0.5*(rk - rl);
1428  }
1429  else {
1430  s = rj - rk;
1431  t = rk - rl;
1432  }
1433  BigReal invs2 = s.rlength();
1434  invs2 *= invs2;
1435 
1436  Vector p = cross(r,s);
1437 #if !defined(FIX_FOR_WATER)
1438  BigReal invp = p.rlength();
1439 #else
1440  BigReal p2 = p.length2(); // fix division by zero above
1441 #endif
1442 
1443  Vector q = cross(s,t);
1444  BigReal invq = q.rlength();
1445 
1446 #if !defined(FIX_FOR_WATER)
1447  BigReal fpdot = (fi*p) * invp;
1448  Vector fp = p * fpdot;
1449  Vector ft = fi - fr - fp;
1450 #else
1451  BigReal fpdot;
1452  Vector fp, ft;
1453  if (p2 < 1e-6) { // vector is near zero, assume no fp contribution to force
1454  fpdot = 0;
1455  fp = 0;
1456  ft = fi - fr;
1457  }
1458  else {
1459  fpdot = (fi*p) / sqrt(p2);
1460  fp = p * fpdot;
1461  ft = fi - fr - fp;
1462  }
1463 #endif
1464 
1465  fja += ft;
1466  Vector v = cross(r,ft); // torque
1467  ft = cross(s,v) * invs2;
1468  fja -= ft;
1469 
1470  if (midpt) {
1471  fka += 0.5 * ft;
1472  fla += 0.5 * ft;
1473  }
1474  else {
1475  fka += ft;
1476  }
1477 
1478  BigReal srdot = (s*r) * invs2;
1479  Vector rr = r - s*srdot;
1480  BigReal rrdot = rr.length();
1481  BigReal stdot = (s*t) * invs2;
1482  Vector tt = t - s*stdot;
1483  BigReal invtt = tt.rlength();
1484  BigReal fact = rrdot*fpdot*invtt*invq;
1485  Vector fq = q * fact;
1486 
1487  fla += fq;
1488  fja += fp*(1+srdot) + fq*stdot;
1489 
1490  ft = fq*(1+stdot) + fp*srdot;
1491 
1492  if (midpt) {
1493  fka += -0.5*ft;
1494  fla += -0.5*ft;
1495  }
1496  else {
1497  fka -= ft;
1498  }
1499 
1500  if (virial) {
1501  Tensor va = outer(fja,rj);
1502  va += outer(fka,rk);
1503  va += outer(fla,rl);
1504  va -= outer(fi,ri);
1505  *virial += va;
1506  }
1507 
1508  fi = 0; // lone pair has zero force
1509  fj += fja;
1510  fk += fka;
1511  fl += fla;
1512 
1513 #ifdef DEBUG_REDISTRIB_FORCE
1514 #define TOL_REDISTRIB 1e-4
1515  Vector fnewnet, tnewnet; // new net force, new net torque
1516  fnewnet = fi + fj + fk + fl;
1517  tnewnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
1518  Vector fdiff = fnewnet - foldnet;
1519  Vector tdiff = tnewnet - toldnet;
1520  if (fdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1521  printf("Error: force redistribution for water exceeded tolerance: "
1522  "fdiff=(%f, %f, %f)\n", fdiff.x, fdiff.y, fdiff.z);
1523  }
1524  if (tdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1525  printf("Error: torque redistribution for water exceeded tolerance: "
1526  "tdiff=(%f, %f, %f)\n", tdiff.x, tdiff.y, tdiff.z);
1527  }
1528 #endif
1529 }
1530 
1531 void HomePatch::redistrib_ap_force(Vector& fi, Vector& fj) {
1532  // final state 'atom' force must transfer to initial state atom
1533  // in single topology FEP.
1534  Vector fi_old = fi;
1535  Vector fj_old = fj;
1536  fi = fi_old + fj_old;
1537  fj = fi_old + fj_old;
1538 }
1539 
1540 /* Redistribute forces from the massless lonepair charge particle onto
1541  * the other atoms of the water.
1542  *
1543  * This is done using the same algorithm as charmm uses for TIP4P lonepairs.
1544  *
1545  * Pass by reference the forces (O H1 H2 LP) to be modified,
1546  * pass by constant reference the corresponding positions,
1547  * and a pointer to virial.
1548  */
1549 void HomePatch::redistrib_lp_water_force(
1550  Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
1551  const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
1552  const Vector& p_lp, Tensor *virial) {
1553 
1554 #ifdef DEBUG_REDISTRIB_FORCE
1555  // Debug information to check against results at end
1556 
1557  // total force and torque relative to origin
1558  Vector totforce, tottorque;
1559  totforce = f_ox + f_h1 + f_h2 + f_lp;
1560  tottorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
1561  //printf("Torque without LP is %f/%f/%f\n",
1562  // tottorque.x, tottorque.y, tottorque.z);
1563  tottorque += cross(f_lp, p_lp);
1564  //printf("Torque with LP is %f/%f/%f\n",
1565  // tottorque.x, tottorque.y, tottorque.z);
1566 #endif
1567 
1568  // accumulate force adjustments
1569  Vector fad_ox(0), fad_h(0);
1570 
1571  // Calculate the radial component of the force and add it to the oxygen
1572  Vector r_ox_lp = p_lp - p_ox;
1573  BigReal invlen2_r_ox_lp = r_ox_lp.rlength();
1574  invlen2_r_ox_lp *= invlen2_r_ox_lp;
1575  BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
1576  Vector f_rad = r_ox_lp * rad_factor;
1577 
1578  fad_ox += f_rad;
1579 
1580  // Calculate the angular component
1581  Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
1582  Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5; // half of r_h2_h1
1583 
1584  // deviation from collinearity of charge site
1585  //Vector r_oop = cross(r_ox_lp, r_hcom_ox);
1586  //
1587  // vector out of o-h-h plane
1588  //Vector r_perp = cross(r_hcom_ox, r_h2_h1_2);
1589 
1590  // Here we assume that Ox/Lp/Hcom are linear
1591  // If you want to correct for deviations, this is the place
1592 
1593 // printf("Deviation from linearity for ox %i: %f/%f/%f\n", oxind, r_oop.x, r_oop.y, r_oop.z);
1594 
1595  Vector f_ang = f_lp - f_rad; // leave the angular component
1596 
1597  // now split this component onto the other atoms
1598  BigReal len_r_ox_lp = r_ox_lp.length();
1599  BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
1600  BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
1601  BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
1602 
1603  fad_ox += (f_ang * oxcomp);
1604  fad_h += (f_ang * hydcomp); // adjustment for both hydrogens
1605 
1606  // Add virial contributions
1607  if (virial) {
1608  Tensor vir = outer(fad_ox, p_ox);
1609  vir += outer(fad_h, p_h1);
1610  vir += outer(fad_h, p_h2);
1611  vir -= outer(f_lp, p_lp);
1612  *virial += vir;
1613  }
1614 
1615  //Vector zerovec(0.0, 0.0, 0.0);
1616  f_lp = 0;
1617  f_ox += fad_ox;
1618  f_h1 += fad_h;
1619  f_h2 += fad_h;
1620 
1621 #ifdef DEBUG_REDISTRIB_FORCE
1622  // Check that the total force and torque come out right
1623  Vector newforce, newtorque;
1624  newforce = f_ox + f_h1 + f_h2;
1625  newtorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
1626  Vector fdiff = newforce - totforce;
1627  Vector tdiff = newtorque - tottorque;
1628  BigReal error = fdiff.length();
1629  if (error > 0.0001) {
1630  printf("Error: Force redistribution for water "
1631  "exceeded force tolerance: error=%f\n", error);
1632  }
1633 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1634  printf("Error in net force: %f\n", error);
1635 #endif
1636 
1637  error = tdiff.length();
1638  if (error > 0.0001) {
1639  printf("Error: Force redistribution for water "
1640  "exceeded torque tolerance: error=%f\n", error);
1641  }
1642 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1643  printf("Error in net torque: %f\n", error);
1644 #endif
1645 #endif /* DEBUG */
1646 }
1647 
1666 void HomePatch::reposition_colinear_lonepair(
1667  Vector& ri, const Vector& rj, const Vector& rk,
1668  Real distance_f, Real scale_f)
1669 {
1670  BigReal distance = distance_f;
1671  BigReal scale = scale_f;
1672  Vector r_jk = rj - rk;
1673  BigReal r2 = r_jk.length2();
1674  if (r2 < 1e-10 || 100. < r2) { // same low tolerance as used in CHARMM
1675  iout << iWARN << "Large/small distance between lonepair reference atoms: ("
1676  << rj << ") (" << rk << ")\n" << endi;
1677  }
1678  ri = rj + (scale + distance*r_jk.rlength())*r_jk;
1679 }
1680 
1695 void HomePatch::reposition_relative_lonepair(
1696  Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
1697  Real distance, Real angle, Real dihedral)
1698 {
1699  if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
1700  iout << iWARN << "Large distance between lonepair reference atoms: ("
1701  << rj << ") (" << rk << ") (" << rl << ")\n" << endi;
1702  }
1703  BigReal r, t, p, cst, snt, csp, snp, invlen;
1704  Vector v, w, a, b, c;
1705 
1706  if (distance >= 0) {
1707  v = rk;
1708  r = distance;
1709  }
1710  else {
1711  v = 0.5*(rk + rl);
1712  r = -distance;
1713  }
1714 
1715  t = angle;
1716  p = dihedral;
1717  cst = cos(t);
1718  snt = sin(t);
1719  csp = cos(p);
1720  snp = sin(p);
1721  a = v - rj;
1722  b = rl - v;
1723  invlen = a.rlength();
1724  a *= invlen;
1725  c = cross(b, a);
1726  invlen = c.rlength();
1727  c *= invlen;
1728  b = cross(a, c);
1729  w.x = r*cst;
1730  w.y = r*snt*csp;
1731  w.z = r*snt*snp;
1732  ri.x = rj.x + w.x*a.x + w.y*b.x + w.z*c.x;
1733  ri.y = rj.y + w.x*a.y + w.y*b.y + w.z*c.y;
1734  ri.z = rj.z + w.x*a.z + w.y*b.z + w.z*c.z;
1735 }
1736 
1737 void HomePatch::reposition_alchpair(Vector& ri, Vector& rj, Mass& Mi, Mass& Mj) {
1738  Vector ri_old, rj_old;
1739  Mass mi, mj;
1740  ri_old.x = ri.x;
1741  ri_old.y = ri.y;
1742  ri_old.z = ri.z;
1743  rj_old.x = rj.x;
1744  rj_old.y = rj.y;
1745  rj_old.z = rj.z;
1746 
1747  mi = Mi; mj = Mj;
1748  ri.x = (mi * ri_old.x + mj * rj_old.x)/(mi + mj);
1749  ri.y = (mi * ri_old.y + mj * rj_old.y)/(mi + mj);
1750  ri.z = (mi * ri_old.z + mj * rj_old.z)/(mi + mj);
1751  rj.x = ri.x;
1752  rj.y = ri.y;
1753  rj.z = ri.z;
1754 }
1755 
1756 void HomePatch::reposition_all_lonepairs(void) {
1757  // ASSERT: simParams->lonepairs == TRUE
1758  for (int i=0; i < numAtoms; i++) {
1759  if (atom[i].mass < 0.01) {
1760  // found a lone pair
1761  AtomID aid = atom[i].id; // global atom ID of lp
1762  Lphost *lph = Node::Object()->molecule->get_lphost(aid); // its lphost
1763  if (lph == NULL) {
1764  char errmsg[512];
1765  sprintf(errmsg, "reposition lone pairs: "
1766  "no Lphost exists for LP %d\n", aid);
1767  NAMD_die(errmsg);
1768  }
1769  LocalID j = AtomMap::Object()->localID(lph->atom2);
1770  LocalID k = AtomMap::Object()->localID(lph->atom3);
1771  LocalID l = AtomMap::Object()->localID(lph->atom4);
1772  if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
1773  char errmsg[512];
1774  sprintf(errmsg, "reposition lone pairs: "
1775  "LP %d has some Lphost atom off patch\n", aid);
1776  NAMD_die(errmsg);
1777  }
1778  // reposition this lone pair
1779  if (lph->numhosts == 2) {
1780  reposition_colinear_lonepair(atom[i].position, atom[j.index].position,
1781  atom[k.index].position, lph->distance, lph->angle);
1782  }
1783  else if (lph->numhosts == 3) {
1784  reposition_relative_lonepair(atom[i].position, atom[j.index].position,
1785  atom[k.index].position, atom[l.index].position,
1786  lph->distance, lph->angle, lph->dihedral);
1787  }
1788  }
1789  }
1790 }
1791 
1792 void HomePatch::reposition_all_alchpairs(void) {
1793  Molecule *mol = Node::Object()->molecule;
1794  int numFepInitial = mol->numFepInitial;
1795  int alchPair_id;
1796  for (int i = 0; i < numAtoms; i++) {
1797  if (atom[i].partition == 4 ) {
1798  alchPair_id = atom[i].id + numFepInitial; //global id of the pair atom.
1799  LocalID j = AtomMap::Object()->localID(alchPair_id);
1800  reposition_alchpair(atom[i].position, atom[j.index].position, atom[i].mass, atom[j.index].mass);
1801  }
1802  }
1803 }
1804 
1805 void HomePatch::swm4_omrepos(Vector *ref, Vector *pos, Vector *vel,
1806  BigReal invdt) {
1807  // Reposition lonepair (Om) particle of Drude SWM4 water.
1808  // Same comments apply as to tip4_omrepos(), but the ordering of atoms
1809  // is different: O, D, LP, H1, H2.
1810  pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
1811  // Now, adjust velocity of particle to get it to appropriate place
1812  // during next integration "drift-step"
1813  if (invdt != 0) {
1814  vel[2] = (pos[2] - ref[2]) * invdt;
1815  }
1816  // No virial correction needed since lonepair is massless
1817 }
1818 
1819 void HomePatch::tip4_omrepos(Vector* ref, Vector* pos, Vector* vel, BigReal invdt) {
1820  /* Reposition the om particle of a tip4p water
1821  * A little geometry shows that the appropriate position is given by
1822  * R_O + (1 / 2 r_ohc) * ( 0.5 (R_H1 + R_H2) - R_O )
1823  * Here r_om is the distance from the oxygen to Om site, and r_ohc
1824  * is the altitude from the oxygen to the hydrogen center of mass
1825  * Those quantities are precalculated upon initialization of HomePatch
1826  *
1827  * Ordering of TIP4P atoms: O, H1, H2, LP.
1828  */
1829 
1830  //printf("rom/rohc are %f %f and invdt is %f\n", r_om, r_ohc, invdt);
1831  //printf("Other positions are: \n 0: %f %f %f\n 1: %f %f %f\n 2: %f %f %f\n", pos[0].x, pos[0].y, pos[0].z, pos[1].x, pos[1].y, pos[1].z, pos[2].x, pos[2].y, pos[2].z);
1832  pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc);
1833  //printf("New position for lp is %f %f %f\n", pos[3].x, pos[3].y, pos[3].z);
1834 
1835  // Now, adjust the velocity of the particle to get it to the appropriate place
1836  if (invdt != 0) {
1837  vel[3] = (pos[3] - ref[3]) * invdt;
1838  }
1839 
1840  // No virial correction needed, since this is a massless particle
1841  return;
1842 }
1843 
1844 void HomePatch::redistrib_lonepair_forces(const int ftag, Tensor *virial) {
1845  // ASSERT: simParams->lonepairs == TRUE
1846  ForceList *f_mod = f;
1847  for (int i = 0; i < numAtoms; i++) {
1848  if (atom[i].mass < 0.01) {
1849  // found a lone pair
1850  AtomID aid = atom[i].id; // global atom ID of lp
1851  Lphost *lph = Node::Object()->molecule->get_lphost(aid); // its lphost
1852  if (lph == NULL) {
1853  char errmsg[512];
1854  sprintf(errmsg, "redistrib lone pair forces: "
1855  "no Lphost exists for LP %d\n", aid);
1856  NAMD_die(errmsg);
1857  }
1858  LocalID j = AtomMap::Object()->localID(lph->atom2);
1859  LocalID k = AtomMap::Object()->localID(lph->atom3);
1860  LocalID l = AtomMap::Object()->localID(lph->atom4);
1861  if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
1862  char errmsg[512];
1863  sprintf(errmsg, "redistrib lone pair forces: "
1864  "LP %d has some Lphost atom off patch\n", aid);
1865  NAMD_die(errmsg);
1866  }
1867  // redistribute forces from this lone pair
1868  if (lph->numhosts == 2) {
1869  redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
1870  f_mod[ftag][k.index], atom[i].position, atom[j.index].position,
1871  atom[k.index].position, lph->distance, lph->angle);
1872  }
1873  else if (lph->numhosts == 3) {
1874  int midpt = (lph->distance < 0);
1875  redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
1876  f_mod[ftag][k.index], f_mod[ftag][l.index],
1877  atom[i].position, atom[j.index].position,
1878  atom[k.index].position, atom[l.index].position, virial, midpt);
1879  }
1880  }
1881  }
1882 }
1883 
1884 void HomePatch::redistrib_alchpair_forces(const int ftag) { //Virial unchanged
1885  SimParameters *simParams = Node::Object()->simParameters;
1886  Molecule *mol = Node::Object()->molecule;
1887  ForceList *f_mod = f;
1888  int numFepInitial = mol->numFepInitial;
1889  BigReal lambda = simParams->alchLambda;
1890  int alchPair_id;
1891  for (int i = 0; i < numAtoms; i++) {
1892  if (atom[i].partition == 4 ) {
1893  alchPair_id = atom[i].id + numFepInitial; //global id of the pair atom.
1894  LocalID j = AtomMap::Object()->localID(alchPair_id);
1895  redistrib_ap_force(f_mod[ftag][i],f_mod[ftag][j.index]);
1896  }
1897  }
1898 }
1899 
1900 void HomePatch::redistrib_swm4_forces(const int ftag, Tensor *virial) {
1901  // Loop over the patch's atoms and apply the appropriate corrections
1902  // to get all forces off of lone pairs
1903  ForceList *f_mod = f;
1904  for (int i = 0; i < numAtoms; i++) {
1905  if (atom[i].mass < 0.01) {
1906  // found lonepair
1907  redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
1908  f_mod[ftag][i+2], f_mod[ftag][i],
1909  atom[i-2].position, atom[i+1].position,
1910  atom[i+2].position, atom[i].position, virial);
1911  }
1912  }
1913 }
1914 
1915 void HomePatch::redistrib_tip4p_forces(const int ftag, Tensor* virial) {
1916  // Loop over the patch's atoms and apply the appropriate corrections
1917  // to get all forces off of lone pairs
1918  // Atom ordering: O H1 H2 LP
1919 
1920  ForceList *f_mod =f;
1921  for (int i=0; i<numAtoms; i++) {
1922  if (atom[i].mass < 0.01) {
1923  // found lonepair
1924  redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
1925  f_mod[ftag][i-1], f_mod[ftag][i],
1926  atom[i-3].position, atom[i-2].position,
1927  atom[i-1].position, atom[i].position, virial);
1928  }
1929  }
1930 }
1931 
1932 
1934  FullAtom * __restrict atom_arr,
1935  const Force * __restrict force_arr,
1936  const BigReal dt,
1937  int num_atoms
1938  ) {
1939  SimParameters *simParams = Node::Object()->simParameters;
1940 
1941  if ( simParams->fixedAtomsOn ) {
1942  for ( int i = 0; i < num_atoms; ++i ) {
1943  if ( atom_arr[i].atomFixed ) {
1944  atom_arr[i].velocity = 0;
1945  } else {
1946  BigReal dt_mass = dt * atom_arr[i].recipMass; // dt/mass
1947  atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1948  atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1949  atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1950  }
1951  }
1952  } else {
1953  for ( int i = 0; i < num_atoms; ++i ) {
1954  BigReal dt_mass = dt * atom_arr[i].recipMass; // dt/mass
1955  atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1956  atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1957  atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1958  }
1959  }
1960 }
1961 
1963  FullAtom * __restrict atom_arr,
1964  const Force * __restrict force_arr1,
1965  const Force * __restrict force_arr2,
1966  const Force * __restrict force_arr3,
1967  const BigReal dt1,
1968  const BigReal dt2,
1969  const BigReal dt3,
1970  int num_atoms
1971  ) {
1972  SimParameters *simParams = Node::Object()->simParameters;
1973 
1974  if ( simParams->fixedAtomsOn ) {
1975  for ( int i = 0; i < num_atoms; ++i ) {
1976  if ( atom_arr[i].atomFixed ) {
1977  atom_arr[i].velocity = 0;
1978  } else {
1979  BigReal rmass = atom_arr[i].recipMass; // 1/mass
1980  atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1981  + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1982  atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1983  + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1984  atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1985  + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1986  }
1987  }
1988  } else {
1989  for ( int i = 0; i < num_atoms; ++i ) {
1990  BigReal rmass = atom_arr[i].recipMass; // 1/mass
1991  atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1992  + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1993  atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1994  + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1995  atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1996  + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1997  }
1998  }
1999 }
2000 
2002  FullAtom * __restrict atom_arr,
2003  const BigReal dt,
2004  int num_atoms
2005  ) {
2006  SimParameters *simParams = Node::Object()->simParameters;
2007  if ( simParams->fixedAtomsOn ) {
2008  for ( int i = 0; i < num_atoms; ++i ) {
2009  if ( ! atom_arr[i].atomFixed ) {
2010  atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2011  atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2012  atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2013  }
2014  }
2015  } else {
2016  for ( int i = 0; i < num_atoms; ++i ) {
2017  atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2018  atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2019  atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2020  }
2021  }
2022 }
2023 
2024 int HomePatch::hardWallDrude(const BigReal timestep, Tensor *virial,
2025  SubmitReduction *ppreduction)
2026 {
2027  Molecule *mol = Node::Object()->molecule;
2028  SimParameters *simParams = Node::Object()->simParameters;
2029  const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
2030  const int fixedAtomsOn = simParams->fixedAtomsOn;
2031  const BigReal dt = timestep / TIMEFACTOR;
2032  const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2033  int i, ia, ib, j;
2034  int dieOnError = simParams->rigidDie;
2035  Tensor wc; // constraint virial
2036  BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
2037  int nslabs;
2038 
2039  // start data for hard wall boundary between drude and its host atom
2040  // static int Count=0;
2041  int Idx;
2042  double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
2043  Vector v_ab, vb_1, vp_1, vb_2, vp_2, new_vel_a, new_vel_b, new_pos_a, new_pos_b, *new_pos, *new_vel;
2044  double dot_v_r_1, dot_v_r_2;
2045  double vb_cm, dr_a, dr_b;
2046  // end data for hard wall boundary between drude and its host atom
2047 
2048  // start calculation of hard wall boundary between drude and its host atom
2049  if (simParams->drudeHardWallOn) {
2050  if (ppreduction) {
2051  nslabs = simParams->pressureProfileSlabs;
2052  idz = nslabs/lattice.c().z;
2053  zmin = lattice.origin().z - 0.5*lattice.c().z;
2054  }
2055 
2056  r_wall = simParams->drudeBondLen;
2057  r_wall_SQ = r_wall*r_wall;
2058  // Count++;
2059  for (i=1; i<numAtoms; i++) {
2060  if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) { // drude particle
2061  ia = i-1;
2062  ib = i;
2063 
2064  v_ab = atom[ib].position - atom[ia].position;
2065  rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
2066 
2067  if (rab_SQ > r_wall_SQ) { // to impose the hard wall constraint
2068  rab = sqrt(rab_SQ);
2069  if ( (rab > (2.0*r_wall)) && dieOnError ) { // unexpected situation
2070  iout << iERROR << "HardWallDrude> "
2071  << "The drude is too far away from atom "
2072  << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
2073  return -1; // triggers early exit
2074  }
2075 
2076  v_ab.x /= rab;
2077  v_ab.y /= rab;
2078  v_ab.z /= rab;
2079 
2080  if ( fixedAtomsOn && atom[ia].atomFixed ) { // the heavy atom is fixed
2081  if (atom[ib].atomFixed) { // the drude is fixed too
2082  continue;
2083  }
2084  else { // only the heavy atom is fixed
2085  dot_v_r_2 = atom[ib].velocity.x*v_ab.x
2086  + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
2087  vb_2 = dot_v_r_2 * v_ab;
2088  vp_2 = atom[ib].velocity - vb_2;
2089 
2090  dr = rab - r_wall;
2091  if(dot_v_r_2 == 0.0) {
2092  delta_T = maxtime;
2093  }
2094  else {
2095  delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
2096  if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
2097  }
2098 
2099  dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
2100 
2101  vb_2 = dot_v_r_2 * v_ab;
2102 
2103  new_vel_a = atom[ia].velocity;
2104  new_vel_b = vp_2 + vb_2;
2105 
2106  dr_b = -dr + delta_T*dot_v_r_2; // L = L_0 + dT *v_new, v was flipped
2107 
2108  new_pos_a = atom[ia].position;
2109  new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
2110  }
2111  }
2112  else {
2113  mass_a = atom[ia].mass;
2114  mass_b = atom[ib].mass;
2115  mass_sum = mass_a+mass_b;
2116 
2117  dot_v_r_1 = atom[ia].velocity.x*v_ab.x
2118  + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
2119  vb_1 = dot_v_r_1 * v_ab;
2120  vp_1 = atom[ia].velocity - vb_1;
2121 
2122  dot_v_r_2 = atom[ib].velocity.x*v_ab.x
2123  + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
2124  vb_2 = dot_v_r_2 * v_ab;
2125  vp_2 = atom[ib].velocity - vb_2;
2126 
2127  vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
2128 
2129  dot_v_r_1 -= vb_cm;
2130  dot_v_r_2 -= vb_cm;
2131 
2132  dr = rab - r_wall;
2133 
2134  if(dot_v_r_2 == dot_v_r_1) {
2135  delta_T = maxtime;
2136  }
2137  else {
2138  delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1); // the time since the collision occurs
2139  if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
2140  }
2141 
2142  // the relative velocity between ia and ib. Drawn according to T_Drude
2143  v_Bond = sqrt(kbt/mass_b);
2144 
2145  // reflect the velocity along bond vector and scale down
2146  dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
2147  dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
2148 
2149  dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
2150  dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
2151 
2152  new_pos_a = atom[ia].position + dr_a*v_ab; // correct the position
2153  new_pos_b = atom[ib].position + dr_b*v_ab;
2154  // atom[ia].position += (dr_a*v_ab); // correct the position
2155  // atom[ib].position += (dr_b*v_ab);
2156 
2157  dot_v_r_1 += vb_cm;
2158  dot_v_r_2 += vb_cm;
2159 
2160  vb_1 = dot_v_r_1 * v_ab;
2161  vb_2 = dot_v_r_2 * v_ab;
2162 
2163  new_vel_a = vp_1 + vb_1;
2164  new_vel_b = vp_2 + vb_2;
2165  }
2166 
2167  int ppoffset, partition;
2168  if ( invdt == 0 ) {
2169  atom[ia].position = new_pos_a;
2170  atom[ib].position = new_pos_b;
2171  }
2172  else if ( virial == 0 ) {
2173  atom[ia].velocity = new_vel_a;
2174  atom[ib].velocity = new_vel_b;
2175  }
2176  else {
2177  for ( j = 0; j < 2; j++ ) {
2178  if (j==0) { // atom ia, heavy atom
2179  Idx = ia;
2180  new_pos = &new_pos_a;
2181  new_vel = &new_vel_a;
2182  }
2183  else if (j==1) { // atom ib, drude
2184  Idx = ib;
2185  new_pos = &new_pos_b;
2186  new_vel = &new_vel_b;
2187  }
2188  Force df = (*new_vel - atom[Idx].velocity) *
2189  ( atom[Idx].mass * invdt );
2190  Tensor vir = outer(df, atom[Idx].position);
2191  wc += vir;
2192  atom[Idx].velocity = *new_vel;
2193  atom[Idx].position = *new_pos;
2194 
2195  if (ppreduction) {
2196  if (!j) {
2197  BigReal z = new_pos->z;
2198  int partition = atom[Idx].partition;
2199  int slab = (int)floor((z-zmin)*idz);
2200  if (slab < 0) slab += nslabs;
2201  else if (slab >= nslabs) slab -= nslabs;
2202  ppoffset = 3*(slab + nslabs*partition);
2203  }
2204  ppreduction->item(ppoffset ) += vir.xx;
2205  ppreduction->item(ppoffset+1) += vir.yy;
2206  ppreduction->item(ppoffset+2) += vir.zz;
2207  }
2208 
2209  }
2210  }
2211  }
2212  }
2213  }
2214 
2215  // if ( (Count>10000) && (Count%10==0) ) {
2216  // v_ab = atom[1].position - atom[0].position;
2217  // rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
2218  // iout << "DBG_R: " << Count << " " << sqrt(rab_SQ) << "\n" << endi;
2219  // }
2220 
2221  }
2222 
2223  // end calculation of hard wall boundary between drude and its host atom
2224 
2225  if ( dt && virial ) *virial += wc;
2226 
2227  return 0;
2228 }
2229 
2231 
2232  SimParameters *simParams = Node::Object()->simParameters;
2233  const int fixedAtomsOn = simParams->fixedAtomsOn;
2234  const int useSettle = simParams->useSettle;
2235 
2236  // Re-size to containt numAtoms elements
2237  velNew.resize(numAtoms);
2238  posNew.resize(numAtoms);
2239 
2240  // Size of a hydrogen group for water
2241  int wathgsize = 3;
2242  int watmodel = simParams->watmodel;
2243  if (watmodel == WAT_TIP4) wathgsize = 4;
2244  else if (watmodel == WAT_SWM4) wathgsize = 5;
2245 
2246  // Initialize the settle algorithm with water parameters
2247  // settle1() assumes all waters are identical,
2248  // and will generate bad results if they are not.
2249  // XXX this will move to Molecule::build_atom_status when that
2250  // version is debugged
2251  if ( ! settle_initialized ) {
2252  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2253  // find a water
2254  if (atom[ig].rigidBondLength > 0) {
2255  int oatm;
2256  if (simParams->watmodel == WAT_SWM4) {
2257  oatm = ig+3; // skip over Drude and Lonepair
2258  //printf("ig=%d mass_ig=%g oatm=%d mass_oatm=%g\n",
2259  // ig, atom[ig].mass, oatm, atom[oatm].mass);
2260  }
2261  else {
2262  oatm = ig+1;
2263  // Avoid using the Om site to set this by mistake
2264  if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2265  oatm += 1;
2266  }
2267  }
2268 
2269  // initialize settle water parameters
2270  settle1init(atom[ig].mass, atom[oatm].mass,
2271  atom[ig].rigidBondLength,
2272  atom[oatm].rigidBondLength,
2273  settle_mOrmT, settle_mHrmT, settle_ra,
2274  settle_rb, settle_rc, settle_rra);
2275  settle_initialized = 1;
2276  break; // done with init
2277  }
2278  }
2279  }
2280 
2281  Vector ref[10];
2282  BigReal rmass[10];
2283  BigReal dsq[10];
2284  int fixed[10];
2285  int ial[10];
2286  int ibl[10];
2287 
2288  int numSettle = 0;
2289  int numRattle = 0;
2290  int posRattleParam = 0;
2291 
2292  settleList.clear();
2293  rattleList.clear();
2294  noconstList.clear();
2295  rattleParam.clear();
2296 
2297  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2298  int hgs = atom[ig].hydrogenGroupSize;
2299  if ( hgs == 1 ) {
2300  // only one atom in group
2301  noconstList.push_back(ig);
2302  continue;
2303  }
2304  int anyfixed = 0;
2305  for (int i = 0; i < hgs; ++i ) {
2306  ref[i] = atom[ig+i].position;
2307  rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2308  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2309  if ( fixed[i] ) {
2310  anyfixed = 1;
2311  rmass[i] = 0.;
2312  }
2313  }
2314  int icnt = 0;
2315  BigReal tmp = atom[ig].rigidBondLength;
2316  if (tmp > 0.0) { // for water
2317  if (hgs != wathgsize) {
2318  char errmsg[256];
2319  sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
2320  "but the specified water model requires %d atoms.\n",
2321  atom[ig].id+1, hgs, wathgsize);
2322  NAMD_die(errmsg);
2323  }
2324  // Use SETTLE for water unless some of the water atoms are fixed,
2325  if (useSettle && !anyfixed) {
2326  // Store to Settle -list
2327  settleList.push_back(ig);
2328  continue;
2329  }
2330  if ( !(fixed[1] && fixed[2]) ) {
2331  dsq[icnt] = tmp * tmp;
2332  ial[icnt] = 1;
2333  ibl[icnt] = 2;
2334  ++icnt;
2335  }
2336  } // if (tmp > 0.0)
2337  for (int i = 1; i < hgs; ++i ) { // normal bonds to mother atom
2338  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2339  if ( !(fixed[0] && fixed[i]) ) {
2340  dsq[icnt] = tmp * tmp;
2341  ial[icnt] = 0;
2342  ibl[icnt] = i;
2343  ++icnt;
2344  }
2345  }
2346  }
2347  if ( icnt == 0 ) {
2348  // no constraints
2349  noconstList.push_back(ig);
2350  continue;
2351  }
2352  // Store to Rattle -list
2353  RattleList rattleListElem;
2354  rattleListElem.ig = ig;
2355  rattleListElem.icnt = icnt;
2356  rattleList.push_back(rattleListElem);
2357  for (int i = 0; i < icnt; ++i ) {
2358  int a = ial[i];
2359  int b = ibl[i];
2360  RattleParam rattleParamElem;
2361  rattleParamElem.ia = a;
2362  rattleParamElem.ib = b;
2363  rattleParamElem.dsq = dsq[i];
2364  rattleParamElem.rma = rmass[a];
2365  rattleParamElem.rmb = rmass[b];
2366  rattleParam.push_back(rattleParamElem);
2367  }
2368  }
2369 
2370 }
2371 
2372 void HomePatch::addRattleForce(const BigReal invdt, Tensor& wc) {
2373  for (int ig = 0; ig < numAtoms; ++ig ) {
2374  Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
2375  Tensor vir = outer(df, atom[ig].position);
2376  wc += vir;
2377  f[Results::normal][ig] += df;
2378  atom[ig].velocity = velNew[ig];
2379  }
2380 }
2381 
2382 int HomePatch::rattle1(const BigReal timestep, Tensor *virial,
2383  SubmitReduction *ppreduction) {
2384 
2385  SimParameters *simParams = Node::Object()->simParameters;
2386  if (simParams->watmodel != WAT_TIP3 || ppreduction) {
2387  // Call old rattle1 -method instead
2388  return rattle1old(timestep, virial, ppreduction);
2389  }
2390 
2391  if (!rattleListValid) {
2392  buildRattleList();
2393  rattleListValid = true;
2394  }
2395 
2396  const int fixedAtomsOn = simParams->fixedAtomsOn;
2397  const int useSettle = simParams->useSettle;
2398  const BigReal dt = timestep / TIMEFACTOR;
2399  const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2400  const BigReal tol2 = 2.0 * simParams->rigidTol;
2401  int maxiter = simParams->rigidIter;
2402  int dieOnError = simParams->rigidDie;
2403 
2404  Vector ref[10]; // reference position
2405  Vector pos[10]; // new position
2406  Vector vel[10]; // new velocity
2407 
2408  // Manual un-roll
2409  int n = (settleList.size()/2)*2;
2410  for (int j=0;j < n;j+=2) {
2411  int ig;
2412  ig = settleList[j];
2413  for (int i = 0; i < 3; ++i ) {
2414  ref[i] = atom[ig+i].position;
2415  pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2416  }
2417  ig = settleList[j+1];
2418  for (int i = 0; i < 3; ++i ) {
2419  ref[i+3] = atom[ig+i].position;
2420  pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
2421  }
2422  settle1_SIMD<2>(ref, pos,
2423  settle_mOrmT, settle_mHrmT, settle_ra,
2424  settle_rb, settle_rc, settle_rra);
2425 
2426  ig = settleList[j];
2427  for (int i = 0; i < 3; ++i ) {
2428  velNew[ig+i] = (pos[i] - ref[i])*invdt;
2429  posNew[ig+i] = pos[i];
2430  }
2431  ig = settleList[j+1];
2432  for (int i = 0; i < 3; ++i ) {
2433  velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
2434  posNew[ig+i] = pos[i+3];
2435  }
2436 
2437  }
2438 
2439  if (settleList.size() % 2) {
2440  int ig = settleList[settleList.size()-1];
2441  for (int i = 0; i < 3; ++i ) {
2442  ref[i] = atom[ig+i].position;
2443  pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2444  }
2445  settle1_SIMD<1>(ref, pos,
2446  settle_mOrmT, settle_mHrmT, settle_ra,
2447  settle_rb, settle_rc, settle_rra);
2448  for (int i = 0; i < 3; ++i ) {
2449  velNew[ig+i] = (pos[i] - ref[i])*invdt;
2450  posNew[ig+i] = pos[i];
2451  }
2452  }
2453 
2454  int posParam = 0;
2455  for (int j=0;j < rattleList.size();++j) {
2456 
2457  BigReal refx[10];
2458  BigReal refy[10];
2459  BigReal refz[10];
2460 
2461  BigReal posx[10];
2462  BigReal posy[10];
2463  BigReal posz[10];
2464 
2465  int ig = rattleList[j].ig;
2466  int icnt = rattleList[j].icnt;
2467  int hgs = atom[ig].hydrogenGroupSize;
2468  for (int i = 0; i < hgs; ++i ) {
2469  ref[i] = atom[ig+i].position;
2470  pos[i] = atom[ig+i].position;
2471  if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
2472  pos[i] += atom[ig+i].velocity * dt;
2473  }
2474  refx[i] = ref[i].x;
2475  refy[i] = ref[i].y;
2476  refz[i] = ref[i].z;
2477  posx[i] = pos[i].x;
2478  posy[i] = pos[i].y;
2479  posz[i] = pos[i].z;
2480  }
2481 
2482  bool done;
2483  bool consFailure;
2484  if (icnt == 1) {
2485  rattlePair<1>(&rattleParam[posParam],
2486  refx, refy, refz,
2487  posx, posy, posz,
2488  consFailure);
2489  done = true;
2490  } else {
2491  rattleN(icnt, &rattleParam[posParam],
2492  refx, refy, refz,
2493  posx, posy, posz,
2494  tol2, maxiter,
2495  done, consFailure);
2496  }
2497 
2498  // Advance position in rattleParam
2499  posParam += icnt;
2500 
2501  for (int i = 0; i < hgs; ++i ) {
2502  pos[i].x = posx[i];
2503  pos[i].y = posy[i];
2504  pos[i].z = posz[i];
2505  }
2506 
2507  for (int i = 0; i < hgs; ++i ) {
2508  velNew[ig+i] = (pos[i] - ref[i])*invdt;
2509  posNew[ig+i] = pos[i];
2510  }
2511 
2512  if ( consFailure ) {
2513  if ( dieOnError ) {
2514  iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
2515  << (atom[ig].id + 1) << "!\n" << endi;
2516  return -1; // triggers early exit
2517  } else {
2518  iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
2519  << (atom[ig].id + 1) << "!\n" << endi;
2520  }
2521  } else if ( ! done ) {
2522  if ( dieOnError ) {
2523  iout << iERROR << "Exceeded RATTLE iteration limit for atom "
2524  << (atom[ig].id + 1) << "!\n" << endi;
2525  return -1; // triggers early exit
2526  } else {
2527  iout << iWARN << "Exceeded RATTLE iteration limit for atom "
2528  << (atom[ig].id + 1) << "!\n" << endi;
2529  }
2530  }
2531  }
2532  // Finally, we have to go through atoms that are not involved in rattle just so that we have
2533  // their positions and velocities up-to-date in posNew and velNew
2534  for (int j=0;j < noconstList.size();++j) {
2535  int ig = noconstList[j];
2536  int hgs = atom[ig].hydrogenGroupSize;
2537  for (int i = 0; i < hgs; ++i ) {
2538  velNew[ig+i] = atom[ig+i].velocity;
2539  posNew[ig+i] = atom[ig+i].position;
2540  }
2541  }
2542 
2543  if ( invdt == 0 ) {
2544  for (int ig = 0; ig < numAtoms; ++ig )
2545  atom[ig].position = posNew[ig];
2546  } else if ( virial == 0 ) {
2547  for (int ig = 0; ig < numAtoms; ++ig )
2548  atom[ig].velocity = velNew[ig];
2549  } else {
2550  Tensor wc; // constraint virial
2551  addRattleForce(invdt, wc);
2552  *virial += wc;
2553  }
2554 
2555  return 0;
2556 }
2557 
2558 // RATTLE algorithm from Allen & Tildesley
2559 int HomePatch::rattle1old(const BigReal timestep, Tensor *virial,
2560  SubmitReduction *ppreduction)
2561 {
2562  Molecule *mol = Node::Object()->molecule;
2563  SimParameters *simParams = Node::Object()->simParameters;
2564  const int fixedAtomsOn = simParams->fixedAtomsOn;
2565  const int useSettle = simParams->useSettle;
2566  const BigReal dt = timestep / TIMEFACTOR;
2567  const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2568  BigReal tol2 = 2.0 * simParams->rigidTol;
2569  int maxiter = simParams->rigidIter;
2570  int dieOnError = simParams->rigidDie;
2571  int i, iter;
2572  BigReal dsq[10], tmp;
2573  int ial[10], ibl[10];
2574  Vector ref[10]; // reference position
2575  Vector refab[10]; // reference vector
2576  Vector pos[10]; // new position
2577  Vector vel[10]; // new velocity
2578  Vector netdp[10]; // total momentum change from constraint
2579  BigReal rmass[10]; // 1 / mass
2580  int fixed[10]; // is atom fixed?
2581  Tensor wc; // constraint virial
2582  BigReal idz, zmin;
2583  int nslabs;
2584 
2585  // Initialize the settle algorithm with water parameters
2586  // settle1() assumes all waters are identical,
2587  // and will generate bad results if they are not.
2588  // XXX this will move to Molecule::build_atom_status when that
2589  // version is debugged
2590  if ( ! settle_initialized ) {
2591  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2592  // find a water
2593  if (atom[ig].rigidBondLength > 0) {
2594  int oatm;
2595  if (simParams->watmodel == WAT_SWM4) {
2596  oatm = ig+3; // skip over Drude and Lonepair
2597  //printf("ig=%d mass_ig=%g oatm=%d mass_oatm=%g\n",
2598  // ig, atom[ig].mass, oatm, atom[oatm].mass);
2599  }
2600  else {
2601  oatm = ig+1;
2602  // Avoid using the Om site to set this by mistake
2603  if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2604  oatm += 1;
2605  }
2606  }
2607 
2608  // initialize settle water parameters
2609  settle1init(atom[ig].mass, atom[oatm].mass,
2610  atom[ig].rigidBondLength,
2611  atom[oatm].rigidBondLength,
2612  settle_mOrmT, settle_mHrmT, settle_ra,
2613  settle_rb, settle_rc, settle_rra);
2614  settle_initialized = 1;
2615  break; // done with init
2616  }
2617  }
2618  }
2619 
2620  if (ppreduction) {
2621  nslabs = simParams->pressureProfileSlabs;
2622  idz = nslabs/lattice.c().z;
2623  zmin = lattice.origin().z - 0.5*lattice.c().z;
2624  }
2625 
2626  // Size of a hydrogen group for water
2627  int wathgsize = 3;
2628  int watmodel = simParams->watmodel;
2629  if (watmodel == WAT_TIP4) wathgsize = 4;
2630  else if (watmodel == WAT_SWM4) wathgsize = 5;
2631 
2632  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2633  int hgs = atom[ig].hydrogenGroupSize;
2634  if ( hgs == 1 ) continue; // only one atom in group
2635  // cache data in local arrays and integrate positions normally
2636  int anyfixed = 0;
2637  for ( i = 0; i < hgs; ++i ) {
2638  ref[i] = atom[ig+i].position;
2639  pos[i] = atom[ig+i].position;
2640  vel[i] = atom[ig+i].velocity;
2641  rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2642  //printf("rmass of %i is %f\n", ig+i, rmass[i]);
2643  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2644  //printf("fixed status of %i is %i\n", i, fixed[i]);
2645  // undo addVelocityToPosition to get proper reference coordinates
2646  if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
2647  }
2648  int icnt = 0;
2649  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
2650  if (hgs != wathgsize) {
2651  char errmsg[256];
2652  sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
2653  "but the specified water model requires %d atoms.\n",
2654  atom[ig].id+1, hgs, wathgsize);
2655  NAMD_die(errmsg);
2656  }
2657  // Use SETTLE for water unless some of the water atoms are fixed,
2658  if (useSettle && !anyfixed) {
2659  if (simParams->watmodel == WAT_SWM4) {
2660  // SWM4 ordering: O D LP H1 H2
2661  // do swap(O,LP) and call settle with subarray O H1 H2
2662  // swap back after we return
2663  Vector lp_ref = ref[2];
2664  Vector lp_pos = pos[2];
2665  Vector lp_vel = vel[2];
2666  ref[2] = ref[0];
2667  pos[2] = pos[0];
2668  vel[2] = vel[0];
2669  settle1(ref+2, pos+2, vel+2, invdt,
2670  settle_mOrmT, settle_mHrmT, settle_ra,
2671  settle_rb, settle_rc, settle_rra);
2672  ref[0] = ref[2];
2673  pos[0] = pos[2];
2674  vel[0] = vel[2];
2675  ref[2] = lp_ref;
2676  pos[2] = lp_pos;
2677  vel[2] = lp_vel;
2678  // determine for LP updated pos and vel
2679  swm4_omrepos(ref, pos, vel, invdt);
2680  }
2681  else {
2682  settle1(ref, pos, vel, invdt,
2683  settle_mOrmT, settle_mHrmT, settle_ra,
2684  settle_rb, settle_rc, settle_rra);
2685  if (simParams->watmodel == WAT_TIP4) {
2686  tip4_omrepos(ref, pos, vel, invdt);
2687  }
2688  }
2689 
2690  // which slab the hydrogen group will belong to
2691  // for pprofile calculations.
2692  int ppoffset, partition;
2693  if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
2694  atom[ig+i].position = pos[i];
2695  } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
2696  atom[ig+i].velocity = vel[i];
2697  } else for ( i = 0; i < wathgsize; ++i ) {
2698  Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
2699  Tensor vir = outer(df, ref[i]);
2700  wc += vir;
2701  f[Results::normal][ig+i] += df;
2702  atom[ig+i].velocity = vel[i];
2703  if (ppreduction) {
2704  // put all the atoms from a water in the same slab. Atom 0
2705  // should be the parent atom.
2706  if (!i) {
2707  BigReal z = pos[i].z;
2708  partition = atom[ig].partition;
2709  int slab = (int)floor((z-zmin)*idz);
2710  if (slab < 0) slab += nslabs;
2711  else if (slab >= nslabs) slab -= nslabs;
2712  ppoffset = 3*(slab + nslabs*partition);
2713  }
2714  ppreduction->item(ppoffset ) += vir.xx;
2715  ppreduction->item(ppoffset+1) += vir.yy;
2716  ppreduction->item(ppoffset+2) += vir.zz;
2717  }
2718  }
2719  continue;
2720  }
2721  if ( !(fixed[1] && fixed[2]) ) {
2722  dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2723  }
2724  }
2725  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
2726  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2727  if ( !(fixed[0] && fixed[i]) ) {
2728  dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
2729  }
2730  }
2731  }
2732  if ( icnt == 0 ) continue; // no constraints
2733  for ( i = 0; i < icnt; ++i ) {
2734  refab[i] = ref[ial[i]] - ref[ibl[i]];
2735  }
2736  for ( i = 0; i < hgs; ++i ) {
2737  netdp[i] = 0.;
2738  }
2739  int done;
2740  int consFailure;
2741  for ( iter = 0; iter < maxiter; ++iter ) {
2742 //if (iter > 0) CkPrintf("iteration %d\n", iter);
2743  done = 1;
2744  consFailure = 0;
2745  for ( i = 0; i < icnt; ++i ) {
2746  int a = ial[i]; int b = ibl[i];
2747  Vector pab = pos[a] - pos[b];
2748  BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
2749  BigReal rabsq = dsq[i];
2750  BigReal diffsq = rabsq - pabsq;
2751  if ( fabs(diffsq) > (rabsq * tol2) ) {
2752  Vector &rab = refab[i];
2753  BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
2754  if ( rpab < ( rabsq * 1.0e-6 ) ) {
2755  done = 0;
2756  consFailure = 1;
2757  continue;
2758  }
2759  BigReal rma = rmass[a];
2760  BigReal rmb = rmass[b];
2761  BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
2762  Vector dp = rab * gab;
2763  pos[a] += rma * dp;
2764  pos[b] -= rmb * dp;
2765  if ( invdt != 0. ) {
2766  dp *= invdt;
2767  netdp[a] += dp;
2768  netdp[b] -= dp;
2769  }
2770  done = 0;
2771  }
2772  }
2773  if ( done ) break;
2774  }
2775 
2776  if ( consFailure ) {
2777  if ( dieOnError ) {
2778  iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
2779  << (atom[ig].id + 1) << "!\n" << endi;
2780  return -1; // triggers early exit
2781  } else {
2782  iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
2783  << (atom[ig].id + 1) << "!\n" << endi;
2784  }
2785  } else if ( ! done ) {
2786  if ( dieOnError ) {
2787  iout << iERROR << "Exceeded RATTLE iteration limit for atom "
2788  << (atom[ig].id + 1) << "!\n" << endi;
2789  return -1; // triggers early exit
2790  } else {
2791  iout << iWARN << "Exceeded RATTLE iteration limit for atom "
2792  << (atom[ig].id + 1) << "!\n" << endi;
2793  }
2794  }
2795 
2796  // store data back to patch
2797  int ppoffset, partition;
2798  if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
2799  atom[ig+i].position = pos[i];
2800  } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
2801  atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2802  } else for ( i = 0; i < hgs; ++i ) {
2803  Force df = netdp[i] * invdt;
2804  Tensor vir = outer(df, ref[i]);
2805  wc += vir;
2806  f[Results::normal][ig+i] += df;
2807  atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2808  if (ppreduction) {
2809  if (!i) {
2810  BigReal z = pos[i].z;
2811  int partition = atom[ig].partition;
2812  int slab = (int)floor((z-zmin)*idz);
2813  if (slab < 0) slab += nslabs;
2814  else if (slab >= nslabs) slab -= nslabs;
2815  ppoffset = 3*(slab + nslabs*partition);
2816  }
2817  ppreduction->item(ppoffset ) += vir.xx;
2818  ppreduction->item(ppoffset+1) += vir.yy;
2819  ppreduction->item(ppoffset+2) += vir.zz;
2820  }
2821  }
2822  }
2823  if ( dt && virial ) *virial += wc;
2824 
2825  return 0;
2826 }
2827 
2828 // RATTLE algorithm from Allen & Tildesley
2829 void HomePatch::rattle2(const BigReal timestep, Tensor *virial)
2830 {
2831  Molecule *mol = Node::Object()->molecule;
2832  SimParameters *simParams = Node::Object()->simParameters;
2833  const int fixedAtomsOn = simParams->fixedAtomsOn;
2834  const int useSettle = simParams->useSettle;
2835  const BigReal dt = timestep / TIMEFACTOR;
2836  Tensor wc; // constraint virial
2837  BigReal tol = simParams->rigidTol;
2838  int maxiter = simParams->rigidIter;
2839  int dieOnError = simParams->rigidDie;
2840  int i, iter;
2841  BigReal dsqi[10], tmp;
2842  int ial[10], ibl[10];
2843  Vector ref[10]; // reference position
2844  Vector refab[10]; // reference vector
2845  Vector vel[10]; // new velocity
2846  BigReal rmass[10]; // 1 / mass
2847  BigReal redmass[10]; // reduced mass
2848  int fixed[10]; // is atom fixed?
2849 
2850  // Size of a hydrogen group for water
2851  int wathgsize = 3;
2852  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
2853  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
2854 
2855  // CkPrintf("In rattle2!\n");
2856  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2857  // CkPrintf("ig=%d\n",ig);
2858  int hgs = atom[ig].hydrogenGroupSize;
2859  if ( hgs == 1 ) continue; // only one atom in group
2860  // cache data in local arrays and integrate positions normally
2861  int anyfixed = 0;
2862  for ( i = 0; i < hgs; ++i ) {
2863  ref[i] = atom[ig+i].position;
2864  vel[i] = atom[ig+i].velocity;
2865  rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
2866  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2867  if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
2868  }
2869  int icnt = 0;
2870  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
2871  if (hgs != wathgsize) {
2872  NAMD_bug("Hydrogen group error caught in rattle2().");
2873  }
2874  // Use SETTLE for water unless some of the water atoms are fixed,
2875  if (useSettle && !anyfixed) {
2876  if (simParams->watmodel == WAT_SWM4) {
2877  // SWM4 ordering: O D LP H1 H2
2878  // do swap(O,LP) and call settle with subarray O H1 H2
2879  // swap back after we return
2880  Vector lp_ref = ref[2];
2881  // Vector lp_vel = vel[2];
2882  ref[2] = ref[0];
2883  vel[2] = vel[0];
2884  settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
2885  ref[0] = ref[2];
2886  vel[0] = vel[2];
2887  ref[2] = lp_ref;
2888  // vel[2] = vel[0]; // set LP vel to O vel
2889  }
2890  else {
2891  settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
2892  if (simParams->watmodel == WAT_TIP4) {
2893  vel[3] = vel[0];
2894  }
2895  }
2896  for (i=0; i<hgs; i++) {
2897  atom[ig+i].velocity = vel[i];
2898  }
2899  continue;
2900  }
2901  if ( !(fixed[1] && fixed[2]) ) {
2902  redmass[icnt] = 1. / (rmass[1] + rmass[2]);
2903  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2904  }
2905  }
2906  // CkPrintf("Loop 2\n");
2907  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
2908  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2909  if ( !(fixed[0] && fixed[i]) ) {
2910  redmass[icnt] = 1. / (rmass[0] + rmass[i]);
2911  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
2912  ibl[icnt] = i; ++icnt;
2913  }
2914  }
2915  }
2916  if ( icnt == 0 ) continue; // no constraints
2917  // CkPrintf("Loop 3\n");
2918  for ( i = 0; i < icnt; ++i ) {
2919  refab[i] = ref[ial[i]] - ref[ibl[i]];
2920  }
2921  // CkPrintf("Loop 4\n");
2922  int done;
2923  for ( iter = 0; iter < maxiter; ++iter ) {
2924  done = 1;
2925  for ( i = 0; i < icnt; ++i ) {
2926  int a = ial[i]; int b = ibl[i];
2927  Vector vab = vel[a] - vel[b];
2928  Vector &rab = refab[i];
2929  BigReal rabsqi = dsqi[i];
2930  BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
2931  if ( (fabs(rvab) * dt * rabsqi) > tol ) {
2932  Vector dp = rab * (-rvab * redmass[i] * rabsqi);
2933  wc += outer(dp,rab);
2934  vel[a] += rmass[a] * dp;
2935  vel[b] -= rmass[b] * dp;
2936  done = 0;
2937  }
2938  }
2939  if ( done ) break;
2940  //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
2941  }
2942  if ( ! done ) {
2943  if ( dieOnError ) {
2944  NAMD_die("Exceeded maximum number of iterations in rattle2().");
2945  } else {
2946  iout << iWARN <<
2947  "Exceeded maximum number of iterations in rattle2().\n" << endi;
2948  }
2949  }
2950  // store data back to patch
2951  for ( i = 0; i < hgs; ++i ) {
2952  atom[ig+i].velocity = vel[i];
2953  }
2954  }
2955  // CkPrintf("Leaving rattle2!\n");
2956  // check that there isn't a constant needed here!
2957  *virial += wc / ( 0.5 * dt );
2958 
2959 }
2960 
2961 
2962 // Adjust gradients for minimizer
2963 void HomePatch::minimize_rattle2(const BigReal timestep, Tensor *virial, bool forces)
2964 {
2965  Molecule *mol = Node::Object()->molecule;
2966  SimParameters *simParams = Node::Object()->simParameters;
2967  Force *f1 = f[Results::normal].begin();
2968  const int fixedAtomsOn = simParams->fixedAtomsOn;
2969  const int useSettle = simParams->useSettle;
2970  const BigReal dt = timestep / TIMEFACTOR;
2971  Tensor wc; // constraint virial
2972  BigReal tol = simParams->rigidTol;
2973  int maxiter = simParams->rigidIter;
2974  int dieOnError = simParams->rigidDie;
2975  int i, iter;
2976  BigReal dsqi[10], tmp;
2977  int ial[10], ibl[10];
2978  Vector ref[10]; // reference position
2979  Vector refab[10]; // reference vector
2980  Vector vel[10]; // new velocity
2981  BigReal rmass[10]; // 1 / mass
2982  BigReal redmass[10]; // reduced mass
2983  int fixed[10]; // is atom fixed?
2984 
2985  // Size of a hydrogen group for water
2986  int wathgsize = 3;
2987  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
2988  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
2989 
2990  // CkPrintf("In rattle2!\n");
2991  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2992  // CkPrintf("ig=%d\n",ig);
2993  int hgs = atom[ig].hydrogenGroupSize;
2994  if ( hgs == 1 ) continue; // only one atom in group
2995  // cache data in local arrays and integrate positions normally
2996  int anyfixed = 0;
2997  for ( i = 0; i < hgs; ++i ) {
2998  ref[i] = atom[ig+i].position;
2999  vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
3000  rmass[i] = 1.0;
3001  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
3002  if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
3003  }
3004  int icnt = 0;
3005  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
3006  if (hgs != wathgsize) {
3007  NAMD_bug("Hydrogen group error caught in rattle2().");
3008  }
3009  // Use SETTLE for water unless some of the water atoms are fixed,
3010  if (useSettle && !anyfixed) {
3011  if (simParams->watmodel == WAT_SWM4) {
3012  // SWM4 ordering: O D LP H1 H2
3013  // do swap(O,LP) and call settle with subarray O H1 H2
3014  // swap back after we return
3015  Vector lp_ref = ref[2];
3016  // Vector lp_vel = vel[2];
3017  ref[2] = ref[0];
3018  vel[2] = vel[0];
3019  settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
3020  ref[0] = ref[2];
3021  vel[0] = vel[2];
3022  ref[2] = lp_ref;
3023  // vel[2] = vel[0]; // set LP vel to O vel
3024  }
3025  else {
3026  settle2(1.0, 1.0, ref, vel, dt, virial);
3027  if (simParams->watmodel == WAT_TIP4) {
3028  vel[3] = vel[0];
3029  }
3030  }
3031  for (i=0; i<hgs; i++) {
3032  ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3033  }
3034  continue;
3035  }
3036  if ( !(fixed[1] && fixed[2]) ) {
3037  redmass[icnt] = 1. / (rmass[1] + rmass[2]);
3038  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3039  }
3040  }
3041  // CkPrintf("Loop 2\n");
3042  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
3043  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3044  if ( !(fixed[0] && fixed[i]) ) {
3045  redmass[icnt] = 1. / (rmass[0] + rmass[i]);
3046  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
3047  ibl[icnt] = i; ++icnt;
3048  }
3049  }
3050  }
3051  if ( icnt == 0 ) continue; // no constraints
3052  // CkPrintf("Loop 3\n");
3053  for ( i = 0; i < icnt; ++i ) {
3054  refab[i] = ref[ial[i]] - ref[ibl[i]];
3055  }
3056  // CkPrintf("Loop 4\n");
3057  int done;
3058  for ( iter = 0; iter < maxiter; ++iter ) {
3059  done = 1;
3060  for ( i = 0; i < icnt; ++i ) {
3061  int a = ial[i]; int b = ibl[i];
3062  Vector vab = vel[a] - vel[b];
3063  Vector &rab = refab[i];
3064  BigReal rabsqi = dsqi[i];
3065  BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
3066  if ( (fabs(rvab) * dt * rabsqi) > tol ) {
3067  Vector dp = rab * (-rvab * redmass[i] * rabsqi);
3068  wc += outer(dp,rab);
3069  vel[a] += rmass[a] * dp;
3070  vel[b] -= rmass[b] * dp;
3071  done = 0;
3072  }
3073  }
3074  if ( done ) break;
3075  //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
3076  }
3077  if ( ! done ) {
3078  if ( dieOnError ) {
3079  NAMD_die("Exceeded maximum number of iterations in rattle2().");
3080  } else {
3081  iout << iWARN <<
3082  "Exceeded maximum number of iterations in rattle2().\n" << endi;
3083  }
3084  }
3085  // store data back to patch
3086  for ( i = 0; i < hgs; ++i ) {
3087  ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3088  }
3089  }
3090  // CkPrintf("Leaving rattle2!\n");
3091  // check that there isn't a constant needed here!
3092  *virial += wc / ( 0.5 * dt );
3093 
3094 }
3095 
3096 
3097 // BEGIN LA
3099 {
3100  DebugM(2, "loweAndersenVelocities\n");
3101  Molecule *mol = Node::Object()->molecule;
3102  SimParameters *simParams = Node::Object()->simParameters;
3103  v.resize(numAtoms);
3104  for (int i = 0; i < numAtoms; ++i) {
3105  //v[i] = p[i];
3106  // co-opt CompAtom structure to pass velocity and mass information
3107  v[i].position = atom[i].velocity;
3108  v[i].charge = atom[i].mass;
3109  }
3110  DebugM(2, "loweAndersenVelocities\n");
3111 }
3112 
3114 {
3115  DebugM(2, "loweAndersenFinish\n");
3116  v.resize(0);
3117 }
3118 // END LA
3119 
3120 //LCPO
3122  Molecule *mol = Node::Object()->molecule;
3123  const int *lcpoParamType = mol->getLcpoParamType();
3124 
3125  lcpoType.resize(numAtoms);
3126  for (int i = 0; i < numAtoms; i++) {
3127  lcpoType[i] = lcpoParamType[pExt[i].id];
3128  }
3129 }
3130 
3131 //set intrinsic radii of atom when doMigration
3133  intRad.resize(numAtoms*2);
3134  intRad.setall(0);
3135  Molecule *mol = Node::Object()->molecule;
3136  SimParameters *simParams = Node::Object()->simParameters;
3137  Real offset = simParams->coulomb_radius_offset;
3138  for (int i = 0; i < numAtoms; i++) {
3139  Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
3140  Real screen = MassToScreen(atom[i].mass);//same
3141  intRad[2*i+0] = rad - offset;//r0
3142  intRad[2*i+1] = screen*(rad - offset);//s0
3143  }
3144 }
3145 
3146 //compute born radius after phase 1, before phase 2
3148 
3149  SimParameters *simParams = Node::Object()->simParameters;
3150  BigReal alphaMax = simParams->alpha_max;
3151  BigReal delta = simParams->gbis_delta;
3152  BigReal beta = simParams->gbis_beta;
3153  BigReal gamma = simParams->gbis_gamma;
3154  BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
3155 
3156  BigReal rhoi;
3157  BigReal rhoi0;
3158  //calculate bornRad from psiSum
3159  for (int i = 0; i < numAtoms; i++) {
3160  rhoi0 = intRad[2*i];
3161  rhoi = rhoi0+coulomb_radius_offset;
3162  psiFin[i] += psiSum[i];
3163  psiFin[i] *= rhoi0;
3164  bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
3165  bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
3166 #ifdef PRINT_COMP
3167  CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
3168 #endif
3169  }
3170 
3171  gbisP2Ready();
3172 }
3173 
3174 //compute dHdrPrefix after phase 2, before phase 3
3176 
3177  SimParameters *simParams = Node::Object()->simParameters;
3178  BigReal delta = simParams->gbis_delta;
3179  BigReal beta = simParams->gbis_beta;
3180  BigReal gamma = simParams->gbis_gamma;
3181  BigReal epsilon_s = simParams->solvent_dielectric;
3182  BigReal epsilon_p = simParams->dielectric;
3183  BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
3184  BigReal epsilon_p_i = 1/simParams->dielectric;
3185  BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
3186  BigReal kappa = simParams->kappa;
3187  BigReal fij, expkappa, Dij, dEdai, dedasum;
3188  BigReal rhoi, rhoi0, psii, nbetapsi;
3189  BigReal gammapsi2, tanhi, daidr;
3190  for (int i = 0; i < numAtoms; i++) {
3191  //add diagonal dEda term
3192  dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
3193  fij = bornRad[i];//inf
3194  expkappa = exp(-kappa*fij);//0
3195  Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
3196  //calculate dHij prefix
3197  dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
3198  *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
3199  dHdrPrefix[i] += dEdai;
3200  dedasum = dHdrPrefix[i];
3201 
3202  rhoi0 = intRad[2*i];
3203  rhoi = rhoi0+coulomb_radius_offset;
3204  psii = psiFin[i];
3205  nbetapsi = -beta*psii;
3206  gammapsi2 = gamma*psii*psii;
3207  tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
3208  daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
3209  * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
3210  dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
3211 #ifdef PRINT_COMP
3212  CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
3213 #endif
3214  }
3215  gbisP3Ready();
3216 }
3217 
3218 //send born radius to proxies to begin phase 2
3220  if (proxy.size() > 0) {
3221  CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3222  for (int i = 0; i < proxy.size(); i++) {
3223  int node = proxy[i];
3225  msg->patch = patchID;
3226  msg->origPe = CkMyPe();
3227  memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
3228  msg->destPe = node;
3229  int seq = flags.sequence;
3231  SET_PRIORITY(msg,seq,priority);
3232  cp[node].recvData(msg);
3233  }
3234  }
3236 }
3237 
3238 //send dHdrPrefix to proxies to begin phase 3
3240  if (proxy.size() > 0) {
3241  CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3242  //only nonzero message should be sent for doFullElec
3243  int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
3244  for (int i = 0; i < proxy.size(); i++) {
3245  int node = proxy[i];
3247  msg->patch = patchID;
3248  msg->dHdrPrefixLen = msgAtoms;
3249  msg->origPe = CkMyPe();
3250  memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
3251  msg->destPe = node;
3252  int seq = flags.sequence;
3254  SET_PRIORITY(msg,seq,priority);
3255  cp[node].recvData(msg);
3256  }
3257  }
3259 }
3260 
3261 //receive proxy results from phase 1
3263  ++numGBISP1Arrived;
3264  for ( int i = 0; i < msg->psiSumLen; ++i ) {
3265  psiFin[i] += msg->psiSum[i];
3266  }
3267  delete msg;
3268 
3269  if (flags.doNonbonded) {
3270  //awaken if phase 1 done
3271  if (phase1BoxClosedCalled == true &&
3272  numGBISP1Arrived==proxy.size() ) {
3273  sequencer->awaken();
3274  }
3275  } else {
3276  //awaken if all phases done on noWork step
3277  if (boxesOpen == 0 &&
3278  numGBISP1Arrived == proxy.size() &&
3279  numGBISP2Arrived == proxy.size() &&
3280  numGBISP3Arrived == proxy.size()) {
3281  sequencer->awaken();
3282  }
3283  }
3284 }
3285 
3286 //receive proxy results from phase 2
3288  ++numGBISP2Arrived;
3289  //accumulate dEda
3290  for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
3291  dHdrPrefix[i] += msg->dEdaSum[i];
3292  }
3293  delete msg;
3294 
3295  if (flags.doNonbonded) {
3296  //awaken if phase 2 done
3297  if (phase2BoxClosedCalled == true &&
3298  numGBISP2Arrived==proxy.size() ) {
3299  sequencer->awaken();
3300  }
3301  } else {
3302  //awaken if all phases done on noWork step
3303  if (boxesOpen == 0 &&
3304  numGBISP1Arrived == proxy.size() &&
3305  numGBISP2Arrived == proxy.size() &&
3306  numGBISP3Arrived == proxy.size()) {
3307  sequencer->awaken();
3308  }
3309  }
3310 }
3311 
3312 // MOLLY algorithm part 1
3314 {
3315  Molecule *mol = Node::Object()->molecule;
3316  SimParameters *simParams = Node::Object()->simParameters;
3317  BigReal tol = simParams->mollyTol;
3318  int maxiter = simParams->mollyIter;
3319  int i, iter;
3320  HGArrayBigReal dsq;
3321  BigReal tmp;
3322  HGArrayInt ial, ibl;
3323  HGArrayVector ref; // reference position
3324  HGArrayVector refab; // reference vector
3325  HGArrayBigReal rmass; // 1 / mass
3326  BigReal *lambda; // Lagrange multipliers
3327  CompAtom *avg; // averaged position
3328  int numLambdas = 0;
3329  HGArrayInt fixed; // is atom fixed?
3330 
3331  // iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
3332  p_avg.resize(numAtoms);
3333  for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
3334 
3335  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3336  int hgs = atom[ig].hydrogenGroupSize;
3337  if ( hgs == 1 ) continue; // only one atom in group
3338  for ( i = 0; i < hgs; ++i ) {
3339  ref[i] = atom[ig+i].position;
3340  rmass[i] = 1. / atom[ig+i].mass;
3341  fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
3342  if ( fixed[i] ) rmass[i] = 0.;
3343  }
3344  avg = &(p_avg[ig]);
3345  int icnt = 0;
3346 
3347  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
3348  if ( hgs != 3 ) {
3349  NAMD_die("Hydrogen group error caught in mollyAverage(). It's a bug!\n");
3350  }
3351  if ( !(fixed[1] && fixed[2]) ) {
3352  dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3353  }
3354  }
3355  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
3356  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3357  if ( !(fixed[0] && fixed[i]) ) {
3358  dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3359  }
3360  }
3361  }
3362  if ( icnt == 0 ) continue; // no constraints
3363  numLambdas += icnt;
3364  molly_lambda.resize(numLambdas);
3365  lambda = &(molly_lambda[numLambdas - icnt]);
3366  for ( i = 0; i < icnt; ++i ) {
3367  refab[i] = ref[ial[i]] - ref[ibl[i]];
3368  }
3369  // iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
3370  iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
3371  if ( iter == maxiter ) {
3372  iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
3373  }
3374  }
3375 
3376  // for ( i=0; i<numAtoms; ++i ) {
3377  // if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
3378  // iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
3379  // << p[i].position << " to " << p_avg[i].position << "\n" << endi;
3380  // }
3381  // }
3382 
3383 }
3384 
3385 
3386 // MOLLY algorithm part 2
3388 {
3389  Molecule *mol = Node::Object()->molecule;
3390  SimParameters *simParams = Node::Object()->simParameters;
3391  Tensor wc; // constraint virial
3392  int i;
3393  HGArrayInt ial, ibl;
3394  HGArrayVector ref; // reference position
3395  CompAtom *avg; // averaged position
3396  HGArrayVector refab; // reference vector
3397  HGArrayVector force; // new force
3398  HGArrayBigReal rmass; // 1 / mass
3399  BigReal *lambda; // Lagrange multipliers
3400  int numLambdas = 0;
3401  HGArrayInt fixed; // is atom fixed?
3402 
3403  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3404  int hgs = atom[ig].hydrogenGroupSize;
3405  if (hgs == 1 ) continue; // only one atom in group
3406  for ( i = 0; i < hgs; ++i ) {
3407  ref[i] = atom[ig+i].position;
3408  force[i] = f[Results::slow][ig+i];
3409  rmass[i] = 1. / atom[ig+i].mass;
3410  fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
3411  if ( fixed[i] ) rmass[i] = 0.;
3412  }
3413  int icnt = 0;
3414  // c-ji I'm only going to mollify water for now
3415  if ( atom[ig].rigidBondLength > 0 ) { // for water
3416  if ( hgs != 3 ) {
3417  NAMD_die("Hydrogen group error caught in mollyMollify(). It's a bug!\n");
3418  }
3419  if ( !(fixed[1] && fixed[2]) ) {
3420  ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3421  }
3422  }
3423  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
3424  if ( atom[ig+i].rigidBondLength > 0 ) {
3425  if ( !(fixed[0] && fixed[i]) ) {
3426  ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3427  }
3428  }
3429  }
3430 
3431  if ( icnt == 0 ) continue; // no constraints
3432  lambda = &(molly_lambda[numLambdas]);
3433  numLambdas += icnt;
3434  for ( i = 0; i < icnt; ++i ) {
3435  refab[i] = ref[ial[i]] - ref[ibl[i]];
3436  }
3437  avg = &(p_avg[ig]);
3438  mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
3439  // store data back to patch
3440  for ( i = 0; i < hgs; ++i ) {
3441  wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
3442  f[Results::slow][ig+i] = force[i];
3443  }
3444  }
3445  // check that there isn't a constant needed here!
3446  *virial += wc;
3447  p_avg.resize(0);
3448 }
3449 
3451  checkpoint_atom.copy(atom);
3452  checkpoint_lattice = lattice;
3453 
3454  // DMK - Atom Separation (water vs. non-water)
3455  #if NAMD_SeparateWaters != 0
3456  checkpoint_numWaterAtoms = numWaterAtoms;
3457  #endif
3458 }
3459 
3460 void HomePatch::revert(void) {
3461  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3462 
3463  atom.copy(checkpoint_atom);
3464  numAtoms = atom.size();
3465  lattice = checkpoint_lattice;
3466 
3467  doAtomUpdate = true;
3468  rattleListValid = false;
3469 
3470  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3471 
3472  // DMK - Atom Separation (water vs. non-water)
3473  #if NAMD_SeparateWaters != 0
3474  numWaterAtoms = checkpoint_numWaterAtoms;
3475  #endif
3476 }
3477 
3478 void HomePatch::exchangeCheckpoint(int scriptTask, int &bpc) { // initiating replica
3479  SimParameters *simParams = Node::Object()->simParameters;
3480  checkpoint_task = scriptTask;
3481  const int remote = simParams->scriptIntArg1;
3482  const char *key = simParams->scriptStringArg1;
3483  PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
3484 }
3485 
3486 void HomePatch::recvCheckpointReq(int task, const char *key, int replica, int pe) { // responding replica
3487  if ( task == SCRIPT_CHECKPOINT_FREE ) {
3488  if ( ! checkpoints.count(key) ) {
3489  NAMD_die("Unable to free checkpoint, requested key was never stored.");
3490  }
3491  delete checkpoints[key];
3492  checkpoints.erase(key);
3493  PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
3494  return;
3495  }
3496  CheckpointAtomsMsg *msg;
3497  if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
3498  if ( ! checkpoints.count(key) ) {
3499  NAMD_die("Unable to load checkpoint, requested key was never stored.");
3500  }
3501  checkpoint_t &cp = *checkpoints[key];
3502  msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
3503  msg->lattice = cp.lattice;
3505  msg->numAtoms = cp.numAtoms;
3506  memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
3507  } else {
3508  msg = new (0,1,0) CheckpointAtomsMsg;
3509  }
3510  msg->pid = patchID;
3511  msg->replica = CmiMyPartition();
3512  msg->pe = CkMyPe();
3513  PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
3514 }
3515 
3516 void HomePatch::recvCheckpointLoad(CheckpointAtomsMsg *msg) { // initiating replica
3517  SimParameters *simParams = Node::Object()->simParameters;
3518  const int remote = simParams->scriptIntArg1;
3519  const char *key = simParams->scriptStringArg1;
3521  NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
3522  }
3523  if ( msg->replica != remote ) {
3524  NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
3525  }
3527  CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
3528  strcpy(newmsg->key,key);
3529  newmsg->lattice = lattice;
3530  newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
3531  newmsg->pid = patchID;
3532  newmsg->pe = CkMyPe();
3533  newmsg->replica = CmiMyPartition();
3534  newmsg->numAtoms = numAtoms;
3535  memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
3536  PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
3537  }
3539  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3540  lattice = msg->lattice;
3542  numAtoms = msg->numAtoms;
3543  atom.resize(numAtoms);
3544  memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
3545  doAtomUpdate = true;
3546  rattleListValid = false;
3547  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3548  }
3551  }
3552  delete msg;
3553 }
3554 
3555 void HomePatch::recvCheckpointStore(CheckpointAtomsMsg *msg) { // responding replica
3556  if ( ! checkpoints.count(msg->key) ) {
3557  checkpoints[msg->key] = new checkpoint_t;
3558  }
3559  checkpoint_t &cp = *checkpoints[msg->key];
3560  cp.lattice = msg->lattice;
3562  cp.numAtoms = msg->numAtoms;
3563  cp.atoms.resize(cp.numAtoms);
3564  memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
3566  delete msg;
3567 }
3568 
3569 void HomePatch::recvCheckpointAck() { // initiating replica
3570  CkpvAccess(_qd)->process();
3571 }
3572 
3573 
3574 void HomePatch::exchangeAtoms(int scriptTask) {
3575  SimParameters *simParams = Node::Object()->simParameters;
3576  // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
3577  if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
3578  exchange_dst = (int) simParams->scriptArg1;
3579  // create and save outgoing message
3580  exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
3584  memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
3585  if ( exchange_req >= 0 ) {
3587  }
3588  }
3589  if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
3590  exchange_src = (int) simParams->scriptArg2;
3592  }
3593 }
3594 
3596  exchange_req = req;
3597  if ( exchange_msg ) {
3598  // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
3600  exchange_msg = 0;
3601  exchange_req = -1;
3602  CkpvAccess(_qd)->process();
3603  }
3604 }
3605 
3607  // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
3608  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3609  lattice = msg->lattice;
3610  numAtoms = msg->numAtoms;
3611  atom.resize(numAtoms);
3612  memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
3613  delete msg;
3614  CkpvAccess(_qd)->process();
3615  doAtomUpdate = true;
3616  rattleListValid = false;
3617  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3618 }
3619 
3620 void HomePatch::submitLoadStats(int timestep)
3621 {
3622  LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
3623 }
3624 
3625 
3626 //
3627 // XXX operates on CompAtom, not FullAtom
3628 //
3630 {
3631 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3632  char dpcbuf[32];
3633  sprintf(dpcbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::DO_PAIRLIST_CHECK], this->getPatchID());
3634  NAMD_EVENT_START_EX(1, NamdProfileEvent::DO_PAIRLIST_CHECK, dpcbuf);
3635 #endif
3636 
3637  SimParameters *simParams = Node::Object()->simParameters;
3638 
3639  if ( numAtoms == 0 || ! flags.usePairlists ) {
3640  flags.pairlistTolerance = 0.;
3641  flags.maxAtomMovement = 99999.;
3642 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3643  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
3644 #endif
3645  return;
3646  }
3647 
3648  int i; int n = numAtoms;
3649  CompAtom *p_i = p.begin();
3650 
3651  if ( flags.savePairlists ) {
3652  flags.pairlistTolerance = doPairlistCheck_newTolerance;
3653  flags.maxAtomMovement = 0.;
3654  doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
3655  doPairlistCheck_lattice = lattice;
3656  doPairlistCheck_positions.resize(numAtoms);
3657  CompAtom *psave_i = doPairlistCheck_positions.begin();
3658  for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
3659 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3660  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
3661 #endif
3662  return;
3663  }
3664 
3665  Lattice &lattice_old = doPairlistCheck_lattice;
3666  Position center_cur = lattice.unscale(center);
3667  Position center_old = lattice_old.unscale(center);
3668  Vector center_delta = center_cur - center_old;
3669 
3670  // find max deviation to corner (any neighbor shares a corner)
3671  BigReal max_cd = 0.;
3672  for ( i=0; i<2; ++i ) {
3673  for ( int j=0; j<2; ++j ) {
3674  for ( int k=0; k<2; ++k ) {
3675  ScaledPosition corner( i ? min.x : max.x ,
3676  j ? min.y : max.y ,
3677  k ? min.z : max.z );
3678  Vector corner_delta =
3679  lattice.unscale(corner) - lattice_old.unscale(corner);
3680  corner_delta -= center_delta;
3681  BigReal cd = corner_delta.length2();
3682  if ( cd > max_cd ) max_cd = cd;
3683  }
3684  }
3685  }
3686  max_cd = sqrt(max_cd);
3687 
3688  // find max deviation of atoms relative to center
3689  BigReal max_pd = 0.;
3690  CompAtom *p_old_i = doPairlistCheck_positions.begin();
3691  for ( i=0; i<n; ++i ) {
3692  Vector p_delta = p_i[i].position - p_old_i[i].position;
3693  p_delta -= center_delta;
3694  BigReal pd = p_delta.length2();
3695  if ( pd > max_pd ) max_pd = pd;
3696  }
3697  max_pd = sqrt(max_pd);
3698 
3699  BigReal max_tol = max_pd + max_cd;
3700 
3701  flags.maxAtomMovement = max_tol;
3702 
3703  // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
3704 
3705  if ( max_tol > ( (1. - simParams->pairlistTrigger) *
3706  doPairlistCheck_newTolerance ) ) {
3707  doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
3708  }
3709 
3710  if ( max_tol > doPairlistCheck_newTolerance ) {
3711  doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
3712  }
3713 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3714  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
3715 #endif
3716 }
3717 
3719 {
3720  if ( ! flags.doNonbonded ) return;
3721 
3722  SimParameters *simParams = Node::Object()->simParameters;
3723  BigReal hgcut = 0.5 * simParams->hgroupCutoff; hgcut *= hgcut;
3724  BigReal maxrad2 = 0.;
3725 
3726  FullAtomList::iterator p_i = atom.begin();
3727  FullAtomList::iterator p_e = atom.end();
3728 
3729  while ( p_i != p_e ) {
3730  const int hgs = p_i->hydrogenGroupSize;
3731  if ( ! hgs ) break; // avoid infinite loop on bug
3732  int ngs = hgs;
3733  if ( ngs > 5 ) ngs = 5; // XXX why? limit to at most 5 atoms per group
3734  BigReal x = p_i->position.x;
3735  BigReal y = p_i->position.y;
3736  BigReal z = p_i->position.z;
3737  int i;
3738  for ( i = 1; i < ngs; ++i ) { // limit spatial extent
3739  p_i[i].nonbondedGroupSize = 0;
3740  BigReal dx = p_i[i].position.x - x;
3741  BigReal dy = p_i[i].position.y - y;
3742  BigReal dz = p_i[i].position.z - z;
3743  BigReal r2 = dx * dx + dy * dy + dz * dz;
3744  if ( r2 > hgcut ) break;
3745  else if ( r2 > maxrad2 ) maxrad2 = r2;
3746  }
3747  p_i->nonbondedGroupSize = i;
3748  for ( ; i < hgs; ++i ) {
3749  p_i[i].nonbondedGroupSize = 1;
3750  }
3751  p_i += hgs;
3752  }
3753 
3754  if ( p_i != p_e ) {
3755  NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
3756  }
3757 
3758  flags.maxGroupRadius = sqrt(maxrad2);
3759 
3760 }
3761 
3763 {
3764  SimParameters *simParams = Node::Object()->simParameters;
3765 
3766  BigReal sysdima = lattice.a_r().unit() * lattice.a();
3767  BigReal sysdimb = lattice.b_r().unit() * lattice.b();
3768  BigReal sysdimc = lattice.c_r().unit() * lattice.c();
3769 
3770  BigReal minSize = simParams->patchDimension - simParams->margin;
3771 
3772  if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
3773  ( bAwayDist*sysdimb < minSize*0.9999 ) ||
3774  ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
3775 
3776  NAMD_die("Periodic cell has become too small for original patch grid!\n"
3777  "Possible solutions are to restart from a recent checkpoint,\n"
3778  "increase margin, or disable useFlexibleCell for liquid simulation.");
3779  }
3780 
3781  BigReal cutoff = simParams->cutoff;
3782 
3783  BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
3784  BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
3785  BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
3786 
3787  if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
3788  NAMD_die("Periodic cell has become too small for original patch grid!\n"
3789  "There are probably many margin violations already reported.\n"
3790  "Possible solutions are to restart from a recent checkpoint,\n"
3791  "increase margin, or disable useFlexibleCell for liquid simulation.");
3792  }
3793 
3794  BigReal minx = min.x - margina;
3795  BigReal miny = min.y - marginb;
3796  BigReal minz = min.z - marginc;
3797  BigReal maxx = max.x + margina;
3798  BigReal maxy = max.y + marginb;
3799  BigReal maxz = max.z + marginc;
3800 
3801  int xdev, ydev, zdev;
3802  int problemCount = 0;
3803 
3804  FullAtomList::iterator p_i = atom.begin();
3805  FullAtomList::iterator p_e = atom.end();
3806  for ( ; p_i != p_e; ++p_i ) {
3807 
3808  ScaledPosition s = lattice.scale(p_i->position);
3809 
3810  // check if atom is within bounds
3811  if (s.x < minx) xdev = 0;
3812  else if (maxx <= s.x) xdev = 2;
3813  else xdev = 1;
3814 
3815  if (s.y < miny) ydev = 0;
3816  else if (maxy <= s.y) ydev = 2;
3817  else ydev = 1;
3818 
3819  if (s.z < minz) zdev = 0;
3820  else if (maxz <= s.z) zdev = 2;
3821  else zdev = 1;
3822 
3823  if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
3824  ++problemCount;
3825  }
3826 
3827  }
3828 
3829  marginViolations = problemCount;
3830  // if ( problemCount ) {
3831  // iout << iERROR <<
3832  // "Found " << problemCount << " margin violations!\n" << endi;
3833  // }
3834 
3835 }
3836 
3837 
3838 void
3840 {
3841  int i;
3842 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3843  NAMD_EVENT_START(1, NamdProfileEvent::ATOM_MIGRATIONS);
3844 #endif
3845  for (i=0; i<numNeighbors; i++) {
3846  realInfo[i].mList.resize(0);
3847  }
3848 
3849  // Purge the AtomMap
3850  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3851 
3852  // Determine atoms that need to migrate
3853 
3854  BigReal minx = min.x;
3855  BigReal miny = min.y;
3856  BigReal minz = min.z;
3857  BigReal maxx = max.x;
3858  BigReal maxy = max.y;
3859  BigReal maxz = max.z;
3860 
3861  int xdev, ydev, zdev;
3862  int delnum = 0;
3863 
3864  FullAtomList::iterator atom_i = atom.begin();
3865  FullAtomList::iterator atom_e = atom.end();
3866 
3867  // DMK - Atom Separation (water vs. non-water)
3868  #if NAMD_SeparateWaters != 0
3869  FullAtomList::iterator atom_first = atom_i;
3870  int numLostWaterAtoms = 0;
3871  #endif
3872 
3873  while ( atom_i != atom_e ) {
3874  // Even though this code iterates through all atoms successively
3875  // it moves entire hydrogen/migration groups as follows:
3876  // Only the parent atom of the hydrogen/migration group has
3877  // nonzero migrationGroupSize. Values determined for xdev,ydev,zdev
3878  // will persist through the remaining group members so that each
3879  // following atom will again be added to the same mList.
3880  if ( atom_i->migrationGroupSize ) {
3881  Position pos = atom_i->position;
3882  if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
3883  // If there are multiple hydrogen groups in a migration group
3884  // (e.g. for supporting lone pairs)
3885  // the following code takes the average position (midpoint)
3886  // of their parents.
3887  int mgs = atom_i->migrationGroupSize;
3888  int c = 1;
3889  for ( int j=atom_i->hydrogenGroupSize; j<mgs;
3890  j+=(atom_i+j)->hydrogenGroupSize ) {
3891  pos += (atom_i+j)->position;
3892  ++c;
3893  }
3894  pos *= 1./c;
3895  // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
3896  }
3897 
3898  // Scaling the position below transforms space within patch from
3899  // what could have been a rotated parallelepiped into
3900  // orthogonal coordinates, where we can use minmax comparison
3901  // to detect which of our nearest neighbors this
3902  // parent atom might have entered.
3903  ScaledPosition s = lattice.scale(pos);
3904 
3905  // check if atom is within bounds
3906  if (s.x < minx) xdev = 0;
3907  else if (maxx <= s.x) xdev = 2;
3908  else xdev = 1;
3909 
3910  if (s.y < miny) ydev = 0;
3911  else if (maxy <= s.y) ydev = 2;
3912  else ydev = 1;
3913 
3914  if (s.z < minz) zdev = 0;
3915  else if (maxz <= s.z) zdev = 2;
3916  else zdev = 1;
3917 
3918  }
3919 
3920  if (mInfo[xdev][ydev][zdev]) { // process atom for migration
3921  // Don't migrate if destination is myself
3922 
3923  // See if we have a migration list already
3924  MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
3925  DebugM(3,"Migrating atom " << atom_i->id << " from patch "
3926  << patchID << " with position " << atom_i->position << "\n");
3927  mCur.add(*atom_i);
3928 
3929  ++delnum;
3930 
3931 
3932  // DMK - Atom Separation (water vs. non-water)
3933  #if NAMD_SeparateWaters != 0
3934  // Check to see if this atom is part of a water molecule. If
3935  // so, numWaterAtoms needs to be adjusted to refect the lost of
3936  // this atom.
3937  // NOTE: The atom separation code assumes that if the oxygen
3938  // atom of the hydrogen group making up the water molecule is
3939  // migrated to another HomePatch, the hydrogens will also
3940  // move!!!
3941  int atomIndex = atom_i - atom_first;
3942  if (atomIndex < numWaterAtoms)
3943  numLostWaterAtoms++;
3944  #endif
3945 
3946 
3947  } else {
3948  // By keeping track of delnum total being deleted from FullAtomList
3949  // the else clause allows us to fill holes as we visit each atom.
3950 
3951  if ( delnum ) { *(atom_i-delnum) = *atom_i; }
3952 
3953  }
3954 
3955  ++atom_i;
3956  }
3957 
3958  // DMK - Atom Separation (water vs. non-water)
3959  #if NAMD_SeparateWaters != 0
3960  numWaterAtoms -= numLostWaterAtoms;
3961  #endif
3962 
3963 
3964  int delpos = numAtoms - delnum;
3965  DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
3966  atom.del(delpos,delnum);
3967 
3968  numAtoms = atom.size();
3969 
3970  PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
3971 
3972  // signal depositMigration() that we are inMigration mode
3973  inMigration = true;
3974 
3975  // Drain the migration message buffer
3976  for (i=0; i<numMlBuf; i++) {
3977  DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
3979  }
3980  numMlBuf = 0;
3981 
3982 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
3983  NAMD_EVENT_STOP(1, NamdProfileEvent::ATOM_MIGRATIONS);
3984 #endif
3985 
3986  if (!allMigrationIn) {
3987  DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
3988  migrationSuspended = true;
3989  sequencer->suspend();
3990  migrationSuspended = false;
3991  }
3992  allMigrationIn = false;
3993 
3994  inMigration = false;
3995  marginViolations = 0;
3996 }
3997 
3998 void
4000 {
4001 
4002  if (!inMigration) { // We have to buffer changes due to migration
4003  // until our patch is in migration mode
4004  msgbuf[numMlBuf++] = msg;
4005  return;
4006  }
4007 
4008 
4009  // DMK - Atom Separation (water vs. non-water)
4010  #if NAMD_SeparateWaters != 0
4011 
4012 
4013  // Merge the incoming list of atoms with the current list of
4014  // atoms. Note that mergeSeparatedAtomList() will apply any
4015  // required transformations to the incoming atoms as it is
4016  // separating them.
4017  mergeAtomList(msg->migrationList);
4018 
4019 
4020  #else
4021 
4022  // Merge the incoming list of atoms with the current list of
4023  // atoms. Apply transformations to the incoming atoms as they are
4024  // added to this patch's list.
4025  {
4026  MigrationList &migrationList = msg->migrationList;
4028  Transform mother_transform;
4029  for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
4030  DebugM(1,"Migrating atom " << mi->id << " to patch "
4031  << patchID << " with position " << mi->position << "\n");
4032  if ( mi->migrationGroupSize ) {
4033  if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
4034  Position pos = mi->position;
4035  int mgs = mi->migrationGroupSize;
4036  int c = 1;
4037  for ( int j=mi->hydrogenGroupSize; j<mgs;
4038  j+=(mi+j)->hydrogenGroupSize ) {
4039  pos += (mi+j)->position;
4040  ++c;
4041  }
4042  pos *= 1./c;
4043  // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
4044  mother_transform = mi->transform;
4045  pos = lattice.nearest(pos,center,&mother_transform);
4046  mi->position = lattice.reverse_transform(mi->position,mi->transform);
4047  mi->position = lattice.apply_transform(mi->position,mother_transform);
4048  mi->transform = mother_transform;
4049  } else {
4050  mi->position = lattice.nearest(mi->position,center,&(mi->transform));
4051  mother_transform = mi->transform;
4052  }
4053  } else {
4054  mi->position = lattice.reverse_transform(mi->position,mi->transform);
4055  mi->position = lattice.apply_transform(mi->position,mother_transform);
4056  mi->transform = mother_transform;
4057  }
4058  atom.add(*mi);
4059  }
4060  }
4061 
4062 
4063  #endif // if (NAMD_SeparateWaters != 0)
4064 
4065 
4066  numAtoms = atom.size();
4067  delete msg;
4068 
4069  DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
4070  if (!--patchMigrationCounter) {
4071  DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
4072  allMigrationIn = true;
4073  patchMigrationCounter = numNeighbors;
4074  if (migrationSuspended) {
4075  DebugM(3,"patch "<<patchID<<" is being awakened\n");
4076  migrationSuspended = false;
4077  sequencer->awaken();
4078  return;
4079  }
4080  else {
4081  DebugM(3,"patch "<<patchID<<" was not suspended\n");
4082  }
4083  }
4084 }
4085 
4086 
4087 
4088 // DMK - Atom Separation (water vs. non-water)
4089 #if NAMD_SeparateWaters != 0
4090 
4091 // This function will separate waters from non-waters in the current
4092 // atom list (regardless of whether or not the atom list is has been
4093 // sorted yet or not).
4094 void HomePatch::separateAtoms() {
4095  SimParameters *simParams = Node::Object()->simParameters;
4096 
4097  // Basic Idea: Iterate through all the atoms in the current list
4098  // of atoms. Pack the waters in the current atoms list and move
4099  // the non-waters to the scratch list. Once the atoms have all
4100  // been separated, copy the non-waters to the end of the waters.
4101  // NOTE: This code does not assume that the atoms list has been
4102  // separated in any manner.
4103 
4104  // NOTE: Sanity check - Doesn't look like the default constructor actually
4105  // adds any atoms but does set numAtoms. ???
4106  if (atom.size() < 0) return; // Nothing to do.
4107 
4108  // Resize the scratch FullAtomList (tempAtom)
4109  tempAtom.resize(numAtoms); // NOTE: Worst case: all non-water
4110 
4111  // Define size of a water hydrogen group
4112  int wathgsize = 3;
4113  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
4114  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
4115 
4116  // Iterate through all the atoms
4117  int i = 0;
4118  int waterIndex = 0;
4119  int nonWaterIndex = 0;
4120  while (i < numAtoms) {
4121 
4122  FullAtom &atom_i = atom[i];
4123  Mass mass = atom_i.mass;
4124  int hgs = atom_i.hydrogenGroupSize;
4125  // Check to see if this hydrogen group is a water molecule
4126  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4127 
4128  // Move this hydrogen group up in the current atom list
4129  if (waterIndex != i) {
4130  atom[waterIndex ] = atom[i ]; // Oxygen
4131  atom[waterIndex + 1] = atom[i + 1]; // Hydrogen
4132  atom[waterIndex + 2] = atom[i + 2]; // Hydrogen
4133  if (wathgsize > 3) atom[waterIndex + 3] = atom[i + 3]; // lonepair
4134  if (wathgsize > 4) atom[waterIndex + 4] = atom[i + 4]; // drude
4135  // actual Drude water is arranged: O D LP H H
4136  }
4137 
4138  waterIndex += wathgsize;
4139  i += wathgsize;
4140 
4141  } else {
4142 
4143  // Move this hydrogen group into non-water (scratch) atom list
4144  for (int j = 0; j < hgs; j++)
4145  tempAtom[nonWaterIndex + j] = atom[i + j];
4146 
4147  nonWaterIndex += hgs;
4148  i += hgs;
4149  }
4150 
4151  } // end iterating through atoms
4152 
4153  // Iterate through the non-water (scratch) atom list, adding the
4154  // atoms to the end of the atom list.
4155  // NOTE: This could be done with a straight memcpy if the internal
4156  // data structures of ResizeArray could be accessed directly.
4157  // Or, perhaps add a member to ResizeArray that can add a consecutive
4158  // list of elements starting at a particular index (would be cleaner).
4159  for (i = 0; i < nonWaterIndex; i++)
4160  atom[waterIndex + i] = tempAtom[i];
4161 
4162  // Set numWaterAtoms
4163  numWaterAtoms = waterIndex;
4164 }
4165 
4166 
4167 // This function will merge the given list of atoms (not assumed to
4168 // be separated) with the current list of atoms (already assumed
4169 // to be separated).
4170 // NOTE: This function applies the transformations to the incoming
4171 // atoms as it is separating them.
4172 void HomePatch::mergeAtomList(FullAtomList &al) {
4173  SimParameters *simParams = Node::Object()->simParameters;
4174 
4175  // Sanity check
4176  if (al.size() <= 0) return; // Nothing to do
4177 
4178  const int orig_atomSize = atom.size();
4179  const int orig_alSize = al.size();
4180 
4181  // Resize the atom list (will eventually hold contents of both lists)
4182  atom.resize(orig_atomSize + orig_alSize); // NOTE: Will have contents of both
4183 
4184 
4185  #if 0 // version where non-waters are moved to scratch first
4186 
4187 
4188  // Basic Idea: The current list is separated already so copy the
4189  // non-water atoms out of it into the scratch atom array. Then
4190  // separate the incoming/given list (al), adding the waters to the
4191  // end of the waters in atom list and non-waters to the end of the
4192  // scratch list. At this point, all waters are in atom list and all
4193  // non-waters are in the scratch list so just copy the scratch list
4194  // to the end of the atom list.
4195  // NOTE: If al is already separated and the number of waters in it
4196  // is know, could simply move the non-waters in atoms back by that
4197  // amount and directly copy the waters in al into the created gap
4198  // and the non-waters in al to the end. Leave this as an
4199  // optimization for later since I'm not sure if this would actually
4200  // do better as the combining code (for combining migration
4201  // messages) would also have to merge the contents of the atom lists
4202  // they carry. Generally speaking, there is probably a faster way
4203  // to do this, but this will get it working.
4204 
4205  // Copy all the non-waters in the current atom list into the
4206  // scratch atom list.
4207  const int orig_atom_numNonWaters = orig_atomSize - numWaterAtoms;
4208  tempAtom.resize(orig_atom_numNonWaters + al.size()); // NOTE: Worst case
4209  for (int i = 0; i < orig_atom_numNonWaters; i++)
4210  tempAtom[i] = atom[numWaterAtoms + i];
4211 
4212  // Separate the contents of the given atom list (applying the
4213  // transforms as needed)
4214  int atom_waterIndex = numWaterAtoms;
4215  int atom_nonWaterIndex = orig_atom_numNonWaters;
4216  int i = 0;
4217  while (i < orig_alSize) {
4218 
4219  FullAtom &atom_i = al[i];
4220  int hgs = atom_i.hydrogenGroupSize;
4221  if ( hgs != atom_i.migrationGroupSize ) {
4222  NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
4223  }
4224  Mass mass = atom_i.mass;
4225 
4226  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4227 
4228  // Apply the transforms
4229 
4230  // Oxygen (@ +0)
4231  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4232  Transform mother_transform = al[i].transform;
4233 
4234  // Hydrogen (@ +1)
4235  al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
4236  al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
4237  al[i+1].transform = mother_transform;
4238 
4239  // Hydrogen (@ +2)
4240  al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
4241  al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
4242  al[i+2].transform = mother_transform;
4243 
4244  // Add to the end of the waters in the current list of atoms
4245  atom[atom_waterIndex ] = al[i ];
4246  atom[atom_waterIndex + 1] = al[i + 1];
4247  atom[atom_waterIndex + 2] = al[i + 2];
4248 
4249  atom_waterIndex += 3;
4250  i += 3;
4251 
4252  } else {
4253 
4254  // Apply the transforms
4255 
4256  // Non-Hydrogen (@ +0)
4257  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4258  Transform mother_transform = al[i].transform;
4259 
4260  // Hydrogens (@ +1 -> +(hgs-1))
4261  for (int j = 1; j < hgs; j++) {
4262  al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
4263  al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
4264  al[i+j].transform = mother_transform;
4265  }
4266 
4267  // Add to the end of the non-waters (scratch) atom list
4268  for (int j = 0; j < hgs; j++)
4269  tempAtom[atom_nonWaterIndex + j] = al[i + j];
4270 
4271  atom_nonWaterIndex += hgs;
4272  i += hgs;
4273  }
4274 
4275  } // end while iterating through given atom list
4276 
4277  // Copy all the non-waters to the end of the current atom list
4278  for (int i = 0; i < atom_nonWaterIndex; i++)
4279  atom[atom_waterIndex + i] = tempAtom[i];
4280 
4281  // Set numWaterAtoms and numAtoms
4282  numWaterAtoms = atom_waterIndex;
4283  numAtoms = atom.size();
4284 
4285 
4286  #else
4287 
4288 
4289  // Basic Idea: Count the number of water atoms in the incoming atom
4290  // list then move the non-waters back in the current atom list to
4291  // make room for the incoming waters. Once there is room in the
4292  // current list, separate the incoming list as the atoms are being
4293  // added to the current list.
4294  // NOTE: Since the incoming atom list is likely to be small,
4295  // iterating over its hydrogen groups twice should not be too bad.
4296  // NOTE: This code assumes the current list is already separated,
4297  // the incoming list may not be separated, and the transforms are
4298  // applied to the incoming atoms as the separation occurs.
4299 
4300  // size of a water hydrogen group
4301  int wathgsize = 3;
4302  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
4303  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
4304 
4305  // Count the number of waters in the given atom list
4306  int al_numWaterAtoms = 0;
4307  int i = 0;
4308  while (i < orig_alSize) {
4309 
4310  FullAtom &atom_i = al[i];
4311  int hgs = atom_i.hydrogenGroupSize;
4312  Mass mass = atom_i.mass;
4313 
4314  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4315  al_numWaterAtoms += wathgsize;
4316  }
4317 
4318  i += hgs;
4319  }
4320 
4321  // Move all of the non-waters in the current atom list back (to a
4322  // higher index) by the number of waters in the given list.
4323  if (al_numWaterAtoms > 0) {
4324  for (i = orig_atomSize - 1; i >= numWaterAtoms; i--) {
4325  atom[i + al_numWaterAtoms] = atom[i];
4326  }
4327  }
4328 
4329  // Separte the atoms in the given atom list. Perform the
4330  // transformations on them and then add them to the appropriate
4331  // location in the current atom list.
4332  int atom_waterIndex = numWaterAtoms;
4333  int atom_nonWaterIndex = orig_atomSize + al_numWaterAtoms;
4334  i = 0;
4335  while (i < orig_alSize) {
4336 
4337  FullAtom &atom_i = al[i];
4338  int hgs = atom_i.hydrogenGroupSize;
4339  if ( hgs != atom_i.migrationGroupSize ) {
4340  NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
4341  }
4342  Mass mass = atom_i.mass;
4343 
4344  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4345 
4346  // Apply the transforms
4347 
4348  // Oxygen (@ +0)
4349  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4350  Transform mother_transform = al[i].transform;
4351 
4352  // Hydrogen (@ +1)
4353  al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
4354  al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
4355  al[i+1].transform = mother_transform;
4356 
4357  // Hydrogen (@ +2)
4358  al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
4359  al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
4360  al[i+2].transform = mother_transform;
4361 
4362  // Add to the end of the waters in the current list of atoms
4363  atom[atom_waterIndex ] = al[i ];
4364  atom[atom_waterIndex + 1] = al[i + 1];
4365  atom[atom_waterIndex + 2] = al[i + 2];
4366 
4367  if (wathgsize > 3) atom[atom_waterIndex + 3] = al[i + 3];
4368 
4369  atom_waterIndex += wathgsize;
4370  i += wathgsize;
4371 
4372  } else {
4373 
4374  // Apply the transforms
4375 
4376  // Non-Hydrogen (@ +0)
4377  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4378  Transform mother_transform = al[i].transform;
4379 
4380  // Hydrogens (@ +1 -> +(hgs-1))
4381  for (int j = 1; j < hgs; j++) {
4382  al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
4383  al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
4384  al[i+j].transform = mother_transform;
4385  }
4386 
4387  // Add to the end of the non-waters (scratch) atom list
4388  for (int j = 0; j < hgs; j++)
4389  atom[atom_nonWaterIndex + j] = al[i + j];
4390 
4391  atom_nonWaterIndex += hgs;
4392  i += hgs;
4393  }
4394 
4395  } // end while iterating through given atom list
4396 
4397  // Set numWaterAtoms and numAtoms
4398  numWaterAtoms = atom_waterIndex;
4399  numAtoms = atom_nonWaterIndex;
4400 
4401  #endif
4402 }
4403 
4404 #endif
4405 
4406 
4407 
4408 inline void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx,
4409  HGArrayBigReal &b)
4410 {
4411  int i,ii=-1,ip,j;
4412  double sum;
4413 
4414  for (i=0;i<n;i++) {
4415  ip=indx[i];
4416  sum=b[ip];
4417  b[ip]=b[i];
4418  if (ii >= 0)
4419  for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
4420  else if (sum) ii=i;
4421  b[i]=sum;
4422  }
4423  for (i=n-1;i>=0;i--) {
4424  sum=b[i];
4425  for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
4426  b[i]=sum/a[i][i];
4427  }
4428 }
4429 
4430 
4431 inline void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
4432 {
4433 
4434  int i,imax,j,k;
4435  double big,dum,sum,temp;
4436  HGArrayBigReal vv;
4437  *d=1.0;
4438  for (i=0;i<n;i++) {
4439  big=0.0;
4440  for (j=0;j<n;j++)
4441  if ((temp=fabs(a[i][j])) > big) big=temp;
4442  if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
4443  vv[i]=1.0/big;
4444  }
4445  for (j=0;j<n;j++) {
4446  for (i=0;i<j;i++) {
4447  sum=a[i][j];
4448  for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
4449  a[i][j]=sum;
4450  }
4451  big=0.0;
4452  for (i=j;i<n;i++) {
4453  sum=a[i][j];
4454  for (k=0;k<j;k++)
4455  sum -= a[i][k]*a[k][j];
4456  a[i][j]=sum;
4457  if ( (dum=vv[i]*fabs(sum)) >= big) {
4458  big=dum;
4459  imax=i;
4460  }
4461  }
4462  if (j != imax) {
4463  for (k=0;k<n;k++) {
4464  dum=a[imax][k];
4465  a[imax][k]=a[j][k];
4466  a[j][k]=dum;
4467  }
4468  *d = -(*d);
4469  vv[imax]=vv[j];
4470  }
4471  indx[j]=imax;
4472  if (a[j][j] == 0.0) a[j][j]=TINY;
4473  if (j != n-1) {
4474  dum=1.0/(a[j][j]);
4475  for (i=j+1;i<n;i++) a[i][j] *= dum;
4476  }
4477  }
4478 }
4479 
4480 
4481 inline void G_q(const HGArrayVector &refab, HGMatrixVector &gqij,
4482  const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl) {
4483  int i;
4484  // step through the rows of the matrix
4485  for(i=0;i<m;i++) {
4486  gqij[i][ial[i]]=2.0*refab[i];
4487  gqij[i][ibl[i]]=-gqij[i][ial[i]];
4488  }
4489 };
4490 
4491 
4492 // c-ji code for MOLLY 7-31-99
4493 int average(CompAtom *qtilde,const HGArrayVector &q,BigReal *lambda,const int n,const int m, const HGArrayBigReal &imass, const HGArrayBigReal &length2, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab, const BigReal tolf, const int ntrial) {
4494  // input: n = length of hydrogen group to be averaged (shaked)
4495  // q[n] = vector array of original positions
4496  // m = number of constraints
4497  // imass[n] = inverse mass for each atom
4498  // length2[m] = square of reference bond length for each constraint
4499  // ial[m] = atom a in each constraint
4500  // ibl[m] = atom b in each constraint
4501  // refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
4502  // tolf = function error tolerance for Newton's iteration
4503  // ntrial = max number of Newton's iterations
4504  // output: lambda[m] = double array of lagrange multipliers (used by mollify)
4505  // qtilde[n] = vector array of averaged (shaked) positions
4506 
4507  int k,k1,i,j;
4508  BigReal errx,errf,d,tolx;
4509 
4510  HGArrayInt indx;
4511  HGArrayBigReal p;
4512  HGArrayBigReal fvec;
4513  HGMatrixBigReal fjac;
4514  HGArrayVector avgab;
4515  HGMatrixVector grhs;
4516  HGMatrixVector auxrhs;
4517  HGMatrixVector glhs;
4518 
4519  // iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
4520  tolx=tolf;
4521 
4522  // initialize lambda, globalGrhs
4523 
4524  for (i=0;i<m;i++) {
4525  lambda[i]=0.0;
4526  }
4527 
4528  // define grhs, auxrhs for all iterations
4529  // grhs= g_x(q)
4530  //
4531  G_q(refab,grhs,n,m,ial,ibl);
4532  for (k=1;k<=ntrial;k++) {
4533  // usrfun(qtilde,q0,lambda,fvec,fjac,n,water);
4534  HGArrayBigReal gij;
4535  // this used to be main loop of usrfun
4536  // compute qtilde given q0, lambda, IMASSes
4537  {
4538  BigReal multiplier;
4539  HGArrayVector tmp;
4540  for (i=0;i<m;i++) {
4541  multiplier = lambda[i];
4542  // auxrhs = M^{-1}grhs^{T}
4543  for (j=0;j<n;j++) {
4544  auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
4545  }
4546  }
4547  for (j=0;j<n;j++) {
4548  // tmp[j]=0.0;
4549  for (i=0;i<m;i++) {
4550  tmp[j]+=auxrhs[i][j];
4551  }
4552  }
4553 
4554  for (j=0;j<n;j++) {
4555  qtilde[j].position=q[j]+tmp[j];
4556  }
4557  // delete [] tmp;
4558  }
4559 
4560  for ( i = 0; i < m; i++ ) {
4561  avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
4562  }
4563 
4564  // iout<<iINFO << "Calling Jac" << std::endl<<endi;
4565  // Jac(qtilde, q0, fjac,n,water);
4566  {
4567  // Vector glhs[3*n+3];
4568 
4569  HGMatrixVector grhs2;
4570 
4571  G_q(avgab,glhs,n,m,ial,ibl);
4572 #ifdef DEBUG0
4573  iout<<iINFO << "G_q:" << std::endl<<endi;
4574  for (i=0;i<m;i++) {
4575  iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
4576  }
4577 #endif
4578  // G_q(refab,grhs2,m,ial,ibl);
4579  // update with the masses
4580  for (j=0; j<n; j++) { // number of atoms
4581  for (i=0; i<m;i++) { // number of constraints
4582  grhs2[i][j] = grhs[i][j]*imass[j];
4583  }
4584  }
4585 
4586  // G_q(qtilde) * M^-1 G_q'(q0) =
4587  // G_q(qtilde) * grhs'
4588  for (i=0;i<m;i++) { // number of constraints
4589  for (j=0;j<m;j++) { // number of constraints
4590  fjac[i][j] = 0;
4591  for (k1=0;k1<n;k1++) {
4592  fjac[i][j] += glhs[i][k1]*grhs2[j][k1];
4593  }
4594  }
4595  }
4596 #ifdef DEBUG0
4597  iout<<iINFO << "glhs" <<endi;
4598  for(i=0;i<9;i++) {
4599  iout<<iINFO << glhs[i] << ","<<endi;
4600  }
4601  iout<<iINFO << std::endl<<endi;
4602  for(i=0;i<9;i++) {
4603  iout<<iINFO << grhs2[i] << ","<<endi;
4604  }
4605  iout<<iINFO << std::endl<<endi;
4606 #endif
4607  // delete[] grhs2;
4608  }
4609  // end of Jac calculation
4610 #ifdef DEBUG0
4611  iout<<iINFO << "Jac" << std::endl<<endi;
4612  for (i=0;i<m;i++)
4613  for (j=0;j<m;j++)
4614  iout<<iINFO << fjac[i][j] << " "<<endi;
4615  iout<< std::endl<<endi;
4616 #endif
4617  // calculate constraints in gij for n constraints this being a water
4618  // G(qtilde, gij, n, water);
4619  for (i=0;i<m;i++) {
4620  gij[i]=avgab[i]*avgab[i]-length2[i];
4621  }
4622 #ifdef DEBUG0
4623  iout<<iINFO << "G" << std::endl<<endi;
4624  iout<<iINFO << "( "<<endi;
4625  for(i=0;i<m-1;i++) {
4626  iout<<iINFO << gij[i] << ", "<<endi;
4627  }
4628  iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
4629 #endif
4630  // fill the return vector
4631  for(i=0;i<m;i++) {
4632  fvec[i] = gij[i];
4633  }
4634  // free up the constraints
4635  // delete[] gij;
4636  // continue Newton's iteration
4637  errf=0.0;
4638  for (i=0;i<m;i++) errf += fabs(fvec[i]);
4639 #ifdef DEBUG0
4640  iout<<iINFO << "errf: " << errf << std::endl<<endi;
4641 #endif
4642  if (errf <= tolf) {
4643  break;
4644  }
4645  for (i=0;i<m;i++) p[i] = -fvec[i];
4646  // iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
4647  ludcmp(fjac,m,indx,&d);
4648  lubksb(fjac,m,indx,p);
4649 
4650  errx=0.0;
4651  for (i=0;i<m;i++) {
4652  errx += fabs(p[i]);
4653  }
4654  for (i=0;i<m;i++)
4655  lambda[i] += p[i];
4656 
4657 #ifdef DEBUG0
4658  iout<<iINFO << "lambda:" << lambda[0]
4659  << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
4660  iout<<iINFO << "errx: " << errx << std::endl<<endi;
4661 #endif
4662  if (errx <= tolx) break;
4663 #ifdef DEBUG0
4664  iout<<iINFO << "Qtilde:" << std::endl<<endi;
4665  iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi;
4666 #endif
4667  }
4668 #ifdef DEBUG
4669  iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
4670 #endif
4671 
4672  return k; //
4673 }
4674 
4675 void mollify(CompAtom *qtilde,const HGArrayVector &q0,const BigReal *lambda, HGArrayVector &force,const int n, const int m, const HGArrayBigReal &imass,const HGArrayInt &ial,const HGArrayInt &ibl,const HGArrayVector &refab) {
4676  int i,j,k;
4677  BigReal d;
4678  HGMatrixBigReal fjac;
4679  Vector zero(0.0,0.0,0.0);
4680 
4681  HGArrayVector tmpforce;
4682  HGArrayVector tmpforce2;
4683  HGArrayVector y;
4684  HGMatrixVector grhs;
4685  HGMatrixVector glhs;
4686  HGArrayBigReal aux;
4687  HGArrayInt indx;
4688 
4689  for(i=0;i<n;i++) {
4690  tmpforce[i]=imass[i]*force[i];
4691  }
4692 
4693  HGMatrixVector grhs2;
4694  HGArrayVector avgab;
4695 
4696  for ( i = 0; i < m; i++ ) {
4697  avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
4698  }
4699 
4700  G_q(avgab,glhs,n,m,ial,ibl);
4701  G_q(refab,grhs,n,m,ial,ibl);
4702  // update with the masses
4703  for (j=0; j<n; j++) { // number of atoms
4704  for (i=0; i<m;i++) { // number of constraints
4705  grhs2[i][j] = grhs[i][j]*imass[j];
4706  }
4707  }
4708 
4709  // G_q(qtilde) * M^-1 G_q'(q0) =
4710  // G_q(qtilde) * grhs'
4711  for (i=0;i<m;i++) { // number of constraints
4712  for (j=0;j<m;j++) { // number of constraints
4713  fjac[j][i] = 0;
4714  for (k=0;k<n;k++) {
4715  fjac[j][i] += glhs[i][k]*grhs2[j][k];
4716  }
4717  }
4718  }
4719 
4720  // aux=gqij*tmpforce
4721  // globalGrhs::computeGlobalGrhs(q0,n,water);
4722  // G_q(refab,grhs,m,ial,ibl);
4723  for(i=0;i<m;i++) {
4724  aux[i]=0.0;
4725  for(j=0;j<n;j++) {
4726  aux[i]+=grhs[i][j]*tmpforce[j];
4727  }
4728  }
4729 
4730  ludcmp(fjac,m,indx,&d);
4731  lubksb(fjac,m,indx,aux);
4732 
4733  for(j=0;j<n;j++) {
4734  y[j] = zero;
4735  for(i=0;i<m;i++) {
4736  y[j] += aux[i]*glhs[i][j];
4737  }
4738  }
4739  for(i=0;i<n;i++) {
4740  y[i]=force[i]-y[i];
4741  }
4742 
4743  // gqq12*y
4744  for(i=0;i<n;i++) {
4745  tmpforce2[i]=imass[i]*y[i];
4746  }
4747 
4748  // here we assume that tmpforce is initialized to zero.
4749  for (i=0;i<n;i++) {
4750  tmpforce[i]=zero;
4751  }
4752 
4753  for (j=0;j<m;j++) {
4754  Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
4755  tmpforce[ial[j]] += tmpf;
4756  tmpforce[ibl[j]] -= tmpf;
4757  }
4758  // c-ji the other bug for 2 constraint water was this line (2-4-99)
4759  // for(i=0;i<m;i++) {
4760  for(i=0;i<n;i++) {
4761  force[i]=tmpforce[i]+y[i];
4762  }
4763 
4764 }
static Node * Object()
Definition: Node.h:86
void depositMigration(MigrateAtomsMsg *)
Definition: HomePatch.C:3999
void recvCheckpointLoad(CheckpointAtomsMsg *msg)
Definition: HomePatch.C:3516
template void settle1_SIMD< 1 >(const Vector *ref, Vector *pos, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
BigReal hgroupCutoff
char scriptStringArg1[128]
void sendProxies()
Definition: HomePatch.C:468
#define NAMD_EVENT_STOP(eon, id)
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Data * clientOpen(int count=1)
Definition: OwnerBox.h:58
void addVelocityToPosition(FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms)
Definition: HomePatch.C:2001
int ib
Definition: Settle.h:57
BigReal scriptArg2
int * lcpoTypeList
Definition: ProxyMgr.h:112
float q
Definition: NamdTypes.h:154
int numNodesWithPatches(void)
Definition: PatchMap.h:61
void runSequencer(void)
Definition: HomePatch.C:269
void minimize_rattle2(const BigReal, Tensor *virial, bool forces=false)
Definition: HomePatch.C:2963
RealList intRad
Definition: Patch.h:155
BigReal pairlistTrigger
void registerProxy(RegisterProxyMsg *)
Definition: HomePatch.C:402
Vector a_r() const
Definition: Lattice.h:268
BigReal solvent_dielectric
void copy(ResizeArray< Elem > &ra)
Definition: ResizeArray.h:59
int sortOrder
Definition: NamdTypes.h:87
OwnerBox< Patch, Results > forceBox
Definition: Patch.h:239
static ProxyMgr * Object()
Definition: ProxyMgr.h:394
int marginViolations
Definition: HomePatch.h:333
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
#define BOLTZMANN
Definition: common.h:47
int32 atom2
Definition: structures.h:121
SortedArray< LSSSubsDat > & lssSubs(ComputeQMMgr *mgr)
Definition: ComputeQM.C:595
#define TINY
Definition: HomePatch.C:51
int AtomID
Definition: NamdTypes.h:29
int flLen[Results::maxNumForces]
Definition: ProxyMgr.h:179
const int * get_qmAtmIndx()
Definition: Molecule.h:823
Lattice & lattice
Definition: Patch.h:126
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
void del(int index, int num=1)
Definition: ResizeArray.h:104
static PatchMap * Object()
Definition: PatchMap.h:27
void clientRemove()
Definition: OwnerBox.h:88
void sendProxies(int pid, int *list, int n)
Definition: ProxyMgr.C:600
void registerIDsFullAtom(const FullAtom *begin, const FullAtom *end)
Definition: AtomMap.C:50
BigReal max_c(int pid) const
Definition: PatchMap.h:96
int find(const Elem &e) const
Definition: ResizeArray.h:137
CompAtom * velocityPtrEnd
Definition: Patch.h:202
int exchange_dst
Definition: HomePatch.h:441
Definition: Vector.h:64
int32 numhosts
Either 2 or 3 host atoms, depending on LP type.
Definition: structures.h:124
BigReal min_a(int pid) const
Definition: PatchMap.h:91
#define NAMD_SeparateWaters
Definition: common.h:173
SimParameters * simParameters
Definition: Node.h:178
int get_numQMAtoms()
Definition: Molecule.h:825
Vector c_r() const
Definition: Lattice.h:270
void sendNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *)
Definition: ProxyMgr.C:1160
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2382
float z
Definition: NamdTypes.h:154
int index_a(int pid) const
Definition: PatchMap.h:86
MigrationList migrationList
void addForceToMomentum(FullAtom *__restrict atom_arr, const Force *__restrict force_arr, const BigReal dt, int num_atoms)
Definition: HomePatch.C:1933
void rattle2(const BigReal, Tensor *virial)
Definition: HomePatch.C:2829
static __thread float4 * forces
int savePairlists
Definition: PatchTypes.h:39
void recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg)
Definition: HomePatch.C:630
float Real
Definition: common.h:109
#define COULOMB
Definition: common.h:46
BigReal & item(int i)
Definition: ReductionMgr.h:312
void gbisComputeAfterP2()
Definition: HomePatch.C:3175
#define DebugM(x, y)
Definition: Debug.h:59
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
void sendExchangeReq(int pid, int src)
Definition: PatchMgr.C:386
Real dihedral
Definition: structures.h:127
BigReal gbis_gamma
BigReal z
Definition: Vector.h:66
Real distance
Definition: structures.h:125
int32 atom3
Definition: structures.h:122
float x
Definition: NamdTypes.h:154
CompAtom * avgPositionPtrEnd
Definition: Patch.h:198
void receiveResults(ProxyResultVarsizeMsg *msg)
Definition: HomePatch.C:796
RealList dHdrPrefix
Definition: Patch.h:159
char const *const NamdProfileEventStr[]
int usePairlists
Definition: PatchTypes.h:38
Position position
Definition: NamdTypes.h:53
BigReal dsq
Definition: Settle.h:58
ExchangeAtomsMsg * exchange_msg
Definition: HomePatch.h:444
void positionsReady(int n=0)
Definition: Patch.C:403
#define PROXY_DATA_PRIORITY
Definition: Priorities.h:40
CompAtomList v
Definition: Patch.h:149
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
virtual void boxClosed(int)
Definition: HomePatch.C:334
int inMigration
Definition: HomePatch.h:492
GBRealList psiFin
Definition: Patch.h:157
if(ComputeNonbondedUtil::goMethod==2)
void doGroupSizeCheck()
Definition: HomePatch.C:3718
BigReal rigidTol
BigReal HGMatrixBigReal[MAXHGS][MAXHGS]
Definition: HomePatch.C:65
CompAtom * avgPositionList
Definition: ProxyMgr.h:104
void loweAndersenVelocities()
Definition: HomePatch.C:3098
#define iout
Definition: InfoStream.h:87
int doLoweAndersen
Definition: PatchTypes.h:26
int pressureProfileSlabs
int numaway_b(void) const
Definition: PatchMap.h:69
CompAtom * velocityList
Definition: ProxyMgr.h:107
BigReal length(void) const
Definition: Vector.h:169
void sortAtomsForCUDA(int *order, const FullAtom *atoms, int nfree, int n)
Definition: SortAtoms.C:116
void recvCheckpointReq(int task, const char *key, int replica, int pe)
Definition: HomePatch.C:3486
Vector origin() const
Definition: Lattice.h:262
CudaAtom * cudaAtomPtr
Definition: Patch.h:205
BigReal alchLambda
void exchangeCheckpoint(int scriptTask, int &bpc)
Definition: HomePatch.C:3478
Position nearest(Position data, ScaledPosition ref) const
Definition: Lattice.h:90
AtomMapper * atomMapper
Definition: Patch.h:152
template void settle1_SIMD< 2 >(const Vector *ref, Vector *pos, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
FullAtom * atoms
Definition: PatchMgr.h:89
GBRealList dEdaSum
Definition: Patch.h:160
void gbisP3Ready()
Definition: HomePatch.C:3239
Vector b_r() const
Definition: Lattice.h:269
BigReal min_b(int pid) const
Definition: PatchMap.h:93
BigReal coulomb_radius_offset
bool rattleListValid
Definition: HomePatch.h:385
BigReal mollyTol
void positionsReady(int doMigration=0)
Definition: HomePatch.C:925
void recvCheckpointAck()
Definition: HomePatch.C:3569
int settle2(BigReal mO, BigReal mH, const Vector *pos, Vector *vel, BigReal dt, Tensor *virial)
Definition: Settle.C:705
std::vector< RattleList > rattleList
Definition: HomePatch.h:381
Definition: Patch.h:35
void patchLoad(PatchID id, int nAtoms, int timestep)
Flags flags
Definition: Patch.h:127
void revert(void)
Definition: HomePatch.C:3460
Charge charge
Definition: NamdTypes.h:54
int32 atom4
Definition: structures.h:123
void sendCheckpointAck(int pid, int dst, int dstpe)
Definition: PatchMgr.C:358
static ProxyNodeAwareSpanningTreeMsg * getANewMsg(PatchID pid, NodeID nid, proxyTreeNode *tree, int size)
Definition: ProxyMgr.C:197
void reorder(Elem *a, int n)
Definition: Random.h:220
NodeIDList tree
Definition: ProxyMgr.h:265
BigReal alpha_max
int periodic_a(void) const
Definition: PatchMap.h:73
int newVdWType
Definition: ComputeQM.h:33
const Real * get_qmAtomGroup() const
Definition: Molecule.h:819
#define PRIORITY_SIZE
Definition: Priorities.h:13
Real r_ohc
Definition: Molecule.h:467
Lattice lattice
Definition: PatchMgr.h:109
void qmSwapAtoms()
Definition: HomePatch.C:866
BigReal HGArrayBigReal[MAXHGS]
Definition: HomePatch.C:63
BigReal rlength(void)
Definition: Vector.h:177
Definition: Random.h:37
int const * getLcpoParamType()
Definition: Molecule.h:478
void awaken(void)
Definition: Sequencer.h:29
#define NAMD_EVENT_START(eon, id)
int boxesOpen
Definition: Patch.h:243
IntList lcpoType
Definition: Patch.h:164
Lphost * get_lphost(int atomid) const
Definition: Molecule.h:1081
static int compDistance(const void *a, const void *b)
Definition: HomePatch.C:456
void unregisterProxy(UnregisterProxyMsg *)
Definition: HomePatch.C:416
GBRealList psiSum
Definition: Patch.h:156
GBReal * dEdaSum
Definition: ProxyMgr.h:51
CudaAtom * cudaAtomList
Definition: ProxyMgr.h:123
void NAMD_bug(const char *err_msg)
Definition: common.C:123
void gbisP2Ready()
Definition: Patch.C:570
ResizeArray< FullAtom > atoms
Definition: HomePatch.h:433
template void rattlePair< 1 >(const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, bool &consFailure)
int oneAwayNeighbors(int pid, PatchID *neighbor_ids=0)
Definition: PatchMap.C:532
void mollify(CompAtom *qtilde, const HGArrayVector &q0, const BigReal *lambda, HGArrayVector &force, const int n, const int m, const HGArrayBigReal &imass, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab)
Definition: HomePatch.C:4675
zVector HGMatrixVector[MAXHGS][MAXHGS]
Definition: HomePatch.C:66
int index
Definition: NamdTypes.h:195
int doFullElectrostatics
Definition: PatchTypes.h:23
iterator end(void)
Definition: ResizeArray.h:37
std::map< std::string, checkpoint_t * > checkpoints
Definition: HomePatch.h:435
#define WAT_SWM4
Definition: common.h:192
void submitLoadStats(int timestep)
Definition: HomePatch.C:3620
CompAtomList p_avg
Definition: Patch.h:147
void mollyMollify(Tensor *virial)
Definition: HomePatch.C:3387
Bool staticAtomAssignment
int numFepInitial
Definition: Molecule.h:607
int rattle1old(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2559
int plExtLen
Definition: ProxyMgr.h:121
gridSize z
~HomePatch()
Definition: HomePatch.C:321
Force * f[maxNumForces]
Definition: PatchTypes.h:67
void buildRattleList()
Definition: HomePatch.C:2230
NodeID destNodeID
Definition: Migration.h:21
int index_b(int pid) const
Definition: PatchMap.h:87
BigReal drudeBondLen
void clientClose(int count=1)
Definition: OwnerBox.h:62
CompAtomList p
Definition: Patch.h:146
void recvSpanningTree(int *t, int n)
Definition: HomePatch.C:636
LocalID localID(AtomID id)
Definition: AtomMap.h:74
BigReal gbis_beta
Real * get_qmAtmChrg()
Definition: Molecule.h:822
static Sync * Object()
Definition: Sync.h:50
int numAtoms
Definition: Patch.h:144
void setall(const Elem &elem)
Definition: ResizeArray.h:90
void run(void)
Definition: Sequencer.C:126
void replaceForces(ExtForce *f)
Definition: HomePatch.C:1331
std::vector< int > settleList
Definition: HomePatch.h:380
BigReal scriptArg1
BigReal x
Definition: Vector.h:66
void sendMigrationMsgs(PatchID, MigrationInfo *, int)
Definition: PatchMgr.C:175
int sequence
Definition: PatchTypes.h:18
int PatchID
Definition: NamdTypes.h:182
void gbisComputeAfterP1()
Definition: HomePatch.C:3147
int doNonbonded
Definition: PatchTypes.h:22
void NAMD_die(const char *err_msg)
Definition: common.C:83
void sendNodeAwareSpanningTree()
Definition: HomePatch.C:631
Force force
Definition: NamdTypes.h:202
void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx, HGArrayBigReal &b)
Definition: HomePatch.C:4408
static LdbCoordinator * Object()
void setLcpoType()
Definition: HomePatch.C:3121
CompAtom * avgPositionPtrBegin
Definition: Patch.h:197
PatchID pid
Definition: NamdTypes.h:194
int exchange_req
Definition: HomePatch.h:443
int periodic_b(void) const
Definition: PatchMap.h:74
void doPairlistCheck()
Definition: HomePatch.C:3629
int berendsenPressure_count
Definition: Sequencer.h:104
static AtomMap * Object()
Definition: AtomMap.h:36
void gbisP3Ready()
Definition: Patch.C:587
void doAtomMigration()
Definition: HomePatch.C:3839
MigrateAtomsMsg * msgbuf[PatchMap::MaxOneAway]
Definition: HomePatch.h:494
int hardWallDrude(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2024
BigReal maxAtomMovement
Definition: PatchTypes.h:41
Tensor outer(const Vector &v1, const Vector &v2)
Definition: Tensor.h:241
void PatchReady(void)
Definition: Sync.C:150
void saveForce(const int ftag=Results::normal)
Definition: HomePatch.C:1336
void setGBISIntrinsicRadii()
Definition: HomePatch.C:3132
std::vector< Vector > posNew
Definition: HomePatch.h:389
BigReal xx
Definition: Tensor.h:17
PatchID patch
Definition: ProxyMgr.h:97
int berendsenPressure_count
Definition: PatchMgr.h:87
int exchange_src
Definition: HomePatch.h:442
BigReal max_b(int pid) const
Definition: PatchMap.h:94
Flags flags
Definition: ProxyMgr.h:98
void gbisP2Ready()
Definition: HomePatch.C:3219
int index_c(int pid) const
Definition: PatchMap.h:88
Real angle
Definition: structures.h:126
#define WAT_TIP3
Definition: common.h:190
void sendSpanningTree()
Definition: HomePatch.C:659
int add(const Elem &elem)
Definition: ResizeArray.h:97
int migrationGroupSize
Definition: NamdTypes.h:117
void checkpoint(void)
Definition: HomePatch.C:3450
BigReal zz
Definition: Tensor.h:19
std::vector< RattleParam > rattleParam
Definition: HomePatch.h:382
void sendProxyData(ProxyDataMsg *, int, int *)
Definition: ProxyMgr.C:1563
BigReal max_a(int pid) const
Definition: PatchMap.h:92
void sendSpanningTree(ProxySpanningTreeMsg *)
Definition: ProxyMgr.C:1155
Real * intRadList
Definition: ProxyMgr.h:110
float y
Definition: NamdTypes.h:154
PatchID getPatchID()
Definition: Patch.h:114
void suspend(void)
Definition: Sequencer.C:136
void sendCheckpointLoad(CheckpointAtomsMsg *msg, int dst, int dstpe)
Definition: PatchMgr.C:308
BigReal length2(void) const
Definition: Vector.h:173
CompAtom * positionList
Definition: ProxyMgr.h:102
int HGArrayInt[MAXHGS]
Definition: HomePatch.C:62
Real newCharge
Definition: ComputeQM.h:34
#define simParams
Definition: Output.C:127
int settle1(const Vector *ref, Vector *pos, Vector *vel, BigReal invdt, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
optimized settle1 algorithm, reuses water properties as much as possible
Definition: Settle.C:55
short vdwType
Definition: NamdTypes.h:55
RealList bornRad
Definition: Patch.h:158
void printOut(char *tag)
Definition: ProxyMgr.C:218
void resize(int i)
Definition: ResizeArray.h:84
void swap(ResizeArray< Elem > &ra)
Definition: ResizeArray.h:64
int node(int pid) const
Definition: PatchMap.h:114
void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
Definition: HomePatch.C:4431
void sendProxyList(int pid, int *plist, int size)
Definition: ProxyMgr.C:1979
#define NAMD_EVENT_START_EX(eon, id, str)
Position apply_transform(Position data, const Transform &t) const
Definition: Lattice.h:132
void sendProxyAll(ProxyDataMsg *, int, int *)
Definition: ProxyMgr.C:1677
Position unscale(ScaledPosition s) const
Definition: Lattice.h:77
MigrationList mList
Definition: Migration.h:22
void clientAdd()
Definition: OwnerBox.h:74
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ cudaTextureObject_t const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t float const PatchPairRecord *__restrict__ const int *__restrict__ atomIndex
ForceList * forceList[Results::maxNumForces]
Definition: ProxyMgr.h:168
void G_q(const HGArrayVector &refab, HGMatrixVector &gqij, const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl)
Definition: HomePatch.C:4481
BigReal rmb
Definition: Settle.h:60
NodeID node
Definition: ProxyMgr.h:166
const PatchID patchID
Definition: Patch.h:143
Position reverse_transform(Position data, const Transform &t) const
Definition: Lattice.h:138
Definition: Tensor.h:15
void mollyAverage()
Definition: HomePatch.C:3313
OwnerBox< Patch, GBReal > dEdaSumBox
Definition: Patch.h:229
int pid(int aIndex, int bIndex, int cIndex)
Definition: PatchMap.inl:27
BigReal y
Definition: Vector.h:66
static float MassToRadius(Mass mi)
Definition: ComputeGBIS.inl:55
Vector b() const
Definition: Lattice.h:253
int flLen[Results::maxNumForces]
Definition: ProxyMgr.h:233
#define MAXHGS
Definition: HomePatch.C:52
Real * dHdrPrefix
Definition: ProxyMgr.h:59
int doLCPO
Definition: PatchTypes.h:29
BigReal rma
Definition: Settle.h:59
Mass mass
Definition: NamdTypes.h:108
static float MassToScreen(Mass mi)
BigReal yy
Definition: Tensor.h:18
void recvExchangeReq(int req)
Definition: HomePatch.C:3595
void settle1init(BigReal pmO, BigReal pmH, BigReal hhdist, BigReal ohdist, BigReal &mOrmT, BigReal &mHrmT, BigReal &ra, BigReal &rb, BigReal &rc, BigReal &rra)
initialize cached water properties
Definition: Settle.C:40
int ia
Definition: Settle.h:56
Elem * find(const Elem &elem)
CompAtomExt * positionExtList
Definition: ProxyMgr.h:122
void buildSpanningTree(void)
Definition: HomePatch.C:674
BigReal pairlistDist
#define TIMEFACTOR
Definition: common.h:48
int checkpoint_task
Definition: HomePatch.h:428
ScaledPosition scale(Position p) const
Definition: Lattice.h:83
float Mass
Definition: ComputeGBIS.inl:20
BigReal patchDimension
int numPatchesOnNode(int node)
Definition: PatchMap.h:60
gridSize y
#define WAT_TIP4
Definition: common.h:191
BigReal pairlistTolerance
Definition: PatchTypes.h:40
BigReal gbis_delta
std::vector< Vector > velNew
Definition: HomePatch.h:388
BigReal pairlistShrink
int numaway_c(void) const
Definition: PatchMap.h:70
void receiveResult(ProxyGBISP1ResultMsg *msg)
Definition: HomePatch.C:3262
void rattleN(const int icnt, const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, const BigReal tol2, const int maxiter, bool &done, bool &consFailure)
Definition: Settle.C:591
BigReal dielectric
Lattice lattice
Definition: PatchMgr.h:82
int doGBIS
Definition: PatchTypes.h:28
int size(void) const
Definition: ResizeArray.h:127
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
BigReal min_c(int pid) const
Definition: PatchMap.h:95
infostream & endi(infostream &s)
Definition: InfoStream.C:38
#define SET_PRIORITY(MSG, SEQ, PRIO)
Definition: Priorities.h:18
void unregisterIDsFullAtom(const FullAtom *begin, const FullAtom *end)
Definition: AtomMap.C:100
zVector HGArrayVector[MAXHGS]
Definition: HomePatch.C:64
BigReal pairlistGrow
int newID
Definition: ComputeQM.h:32
void addForceToMomentum3(FullAtom *__restrict atom_arr, const Force *__restrict force_arr1, const Force *__restrict force_arr2, const Force *__restrict force_arr3, const BigReal dt1, const BigReal dt2, const BigReal dt3, int num_atoms)
Definition: HomePatch.C:1962
int avgPlLen
Definition: ProxyMgr.h:103
void registerPatch(int patchID, int numPes, int *pes)
Definition: ProxyMgr.C:1840
std::vector< int > noconstList
Definition: HomePatch.h:383
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
#define GB2_PROXY_DATA_PRIORITY
Definition: Priorities.h:58
gridSize x
int proxySpanDim
Definition: ProxyMgr.C:48
FullAtom * atoms
Definition: PatchMgr.h:112
int periodic_c(void) const
Definition: PatchMap.h:75
void loweAndersenFinish()
Definition: HomePatch.C:3113
int isOpen()
Definition: OwnerBox.h:51
int numaway_a(void) const
Definition: PatchMap.h:68
void addRattleForce(const BigReal invdt, Tensor &wc)
Definition: HomePatch.C:2372
CompAtomExtList pExt
Definition: Patch.h:174
static __thread int num_atoms
int NodeID
Definition: NamdTypes.h:184
void sendCheckpointStore(CheckpointAtomsMsg *msg, int dst, int dstpe)
Definition: PatchMgr.C:333
Molecule * molecule
Definition: Node.h:176
BigReal maxGroupRadius
Definition: PatchTypes.h:42
Vector a() const
Definition: Lattice.h:252
int numMlBuf
Definition: HomePatch.h:493
void sendCheckpointReq(int pid, int remote, const char *key, int task)
Definition: PatchMgr.C:270
void useSequencer(Sequencer *sequencerPtr)
Definition: HomePatch.C:265
void doMarginCheck()
Definition: HomePatch.C:3762
static PatchMgr * Object()
Definition: PatchMgr.h:152
OwnerBox< Patch, GBReal > psiSumBox
Definition: Patch.h:225
#define GB3_PROXY_DATA_PRIORITY
Definition: Priorities.h:66
Real r_om
Definition: Molecule.h:466
void sendExchangeMsg(ExchangeAtomsMsg *msg, int dst, int dstpe)
Definition: PatchMgr.C:417
int doMolly
Definition: PatchTypes.h:24
void recvCheckpointStore(CheckpointAtomsMsg *msg)
Definition: HomePatch.C:3555
Vector unit(void) const
Definition: Vector.h:182
Vector c() const
Definition: Lattice.h:254
void recvExchangeMsg(ExchangeAtomsMsg *msg)
Definition: HomePatch.C:3606
double BigReal
Definition: common.h:114
int average(CompAtom *qtilde, const HGArrayVector &q, BigReal *lambda, const int n, const int m, const HGArrayBigReal &imass, const HGArrayBigReal &length2, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab, const BigReal tolf, const int ntrial)
Definition: HomePatch.C:4493
#define PATCH_PRIORITY(PID)
Definition: Priorities.h:25
CompAtom * velocityPtrBegin
Definition: Patch.h:201
int proxySendSpanning
Definition: ProxyMgr.C:45
iterator begin(void)
Definition: ResizeArray.h:36
BigReal drudeTemp
void exchangeAtoms(int scriptTask)
Definition: HomePatch.C:3574