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)
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_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_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 #if ! defined(NAMD_CUDA)
1067 #if defined(NAMD_AVXTILES)
1068  if (!simParams->useAVXTiles)
1069 #endif
1070  doGroupSizeCheck();
1071 #endif
1072 
1073  // Copy information needed by computes and proxys to Patch::p.
1074  // Resize only if atoms were migrated
1075  if (doMigration) {
1076  p.resize(numAtoms);
1077  pExt.resize(numAtoms);
1078  }
1079  CompAtom *p_i = p.begin();
1080  CompAtomExt *pExt_i = pExt.begin();
1081  FullAtom *a_i = atom.begin();
1082  int i; int n = numAtoms;
1083  for ( i=0; i<n; ++i ) {
1084  p_i[i] = a_i[i];
1085  pExt_i[i] = a_i[i];
1086  }
1087 
1088  // Measure atom movement to test pairlist validity
1089  doPairlistCheck();
1090 
1091  if (flags.doMolly) mollyAverage();
1092  // BEGIN LA
1094  // END LA
1095 
1096  if (flags.doGBIS) {
1097  //reset for next time step
1098  numGBISP1Arrived = 0;
1099  phase1BoxClosedCalled = false;
1100  numGBISP2Arrived = 0;
1101  phase2BoxClosedCalled = false;
1102  numGBISP3Arrived = 0;
1103  phase3BoxClosedCalled = false;
1104  if (doMigration || isNewProxyAdded)
1106  }
1107 
1108  if (flags.doLCPO) {
1109  if (doMigration || isNewProxyAdded) {
1110  setLcpoType();
1111  }
1112  }
1113 
1114  // Must Add Proxy Changes when migration completed!
1116  int *pids = NULL;
1117  int pidsPreAllocated = 1;
1118  int npid;
1119  if (proxySendSpanning == 0) {
1120  npid = proxy.size();
1121  pids = new int[npid];
1122  pidsPreAllocated = 0;
1123  int *pidi = pids;
1124  int *pide = pids + proxy.size();
1125  int patchNodesLast =
1126  ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
1127  for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
1128  {
1129  if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
1130  *(--pide) = *pli;
1131  } else {
1132  *(pidi++) = *pli;
1133  }
1134  }
1135  }
1136  else {
1137 #ifdef NODEAWARE_PROXY_SPANNINGTREE
1138  #ifdef USE_NODEPATCHMGR
1139  npid = numNodeChild;
1140  pids = nodeChildren;
1141  #else
1142  npid = nChild;
1143  pids = child;
1144  #endif
1145 #else
1146  npid = nChild;
1147  pidsPreAllocated = 0;
1148  pids = new int[proxySpanDim];
1149  for (int i=0; i<nChild; i++) pids[i] = child[i];
1150 #endif
1151  }
1152  if (npid) { //have proxies
1153  int seq = flags.sequence;
1154  int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
1155  //begin to prepare proxy msg and send it
1156  int pdMsgPLLen = p.size();
1157  int pdMsgAvgPLLen = 0;
1158  if(flags.doMolly) {
1159  pdMsgAvgPLLen = p_avg.size();
1160  }
1161  // BEGIN LA
1162  int pdMsgVLLen = 0;
1163  if (flags.doLoweAndersen) {
1164  pdMsgVLLen = v.size();
1165  }
1166  // END LA
1167 
1168  int intRadLen = 0;
1169  if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
1170  intRadLen = numAtoms * 2;
1171  }
1172 
1173  //LCPO
1174  int lcpoTypeLen = 0;
1175  if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
1176  lcpoTypeLen = numAtoms;
1177  }
1178 
1179  int pdMsgPLExtLen = 0;
1180  if(doMigration || isNewProxyAdded) {
1181  pdMsgPLExtLen = pExt.size();
1182  }
1183 
1184  int cudaAtomLen = 0;
1185 #if defined(NAMD_CUDA)
1186  cudaAtomLen = numAtoms;
1187 #elif defined(NAMD_AVXTILES)
1188  if (simParams->useAVXTiles)
1189  cudaAtomLen = (numAtoms + 15) & (~15);
1190 #elif defined(NAMD_MIC)
1191  #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
1192  cudaAtomLen = numAtoms;
1193  #else
1194  cudaAtomLen = (numAtoms + 15) & (~15);
1195  #endif
1196 #endif
1197 
1198  ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
1199  lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
1200 
1201  SET_PRIORITY(nmsg,seq,priority);
1202  nmsg->patch = patchID;
1203  nmsg->flags = flags;
1204  nmsg->plLen = pdMsgPLLen;
1205  //copying data to the newly created msg
1206  memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
1207  nmsg->avgPlLen = pdMsgAvgPLLen;
1208  if(flags.doMolly) {
1209  memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
1210  }
1211  // BEGIN LA
1212  nmsg->vlLen = pdMsgVLLen;
1213  if (flags.doLoweAndersen) {
1214  memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
1215  }
1216  // END LA
1217 
1218  if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
1219  for (int i = 0; i < numAtoms * 2; i++) {
1220  nmsg->intRadList[i] = intRad[i];
1221  }
1222  }
1223 
1224  if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
1225  for (int i = 0; i < numAtoms; i++) {
1226  nmsg->lcpoTypeList[i] = lcpoType[i];
1227  }
1228  }
1229 
1230  nmsg->plExtLen = pdMsgPLExtLen;
1231  if(doMigration || isNewProxyAdded){
1232  memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
1233  }
1234 
1235 // DMK
1236 #if defined(NAMD_CUDA) || defined(NAMD_MIC) || defined(NAMD_AVXTILES)
1237  #ifdef NAMD_AVXTILES
1238  if (simParams->useAVXTiles)
1239  #endif
1240  memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
1241 #endif
1242 
1243 #if NAMD_SeparateWaters != 0
1244  //DMK - Atom Separation (water vs. non-water)
1245  nmsg->numWaterAtoms = numWaterAtoms;
1246 #endif
1247 
1248 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
1249  nmsg->isFromImmMsgCall = 0;
1250 #endif
1251 
1252  #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
1253  DebugFileTrace *dft = DebugFileTrace::Object();
1254  dft->openTrace();
1255  dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
1256  for(int i=0; i<npid; i++) {
1257  dft->writeTrace("%d ", pids[i]);
1258  }
1259  dft->writeTrace("\n");
1260  dft->closeTrace();
1261  #endif
1262 
1263 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1264  if (isProxyChanged || localphs == NULL)
1265  {
1266 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
1267  //CmiAssert(isProxyChanged);
1268  if (nphs) {
1269  for (int i=0; i<nphs; i++) {
1270  CmiDestoryPersistent(localphs[i]);
1271  }
1272  delete [] localphs;
1273  }
1274  localphs = new PersistentHandle[npid];
1275  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;
1276  for (int i=0; i<npid; i++) {
1277 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
1278  if (proxySendSpanning)
1279  localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
1280  else
1281 #endif
1282  localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
1283  }
1284  nphs = npid;
1285  }
1286  CmiAssert(nphs == npid && localphs != NULL);
1287  CmiUsePersistentHandle(localphs, nphs);
1288 #endif
1289  if(doMigration || isNewProxyAdded) {
1290  ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
1291  }else{
1292  ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
1293  }
1294 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1295  CmiUsePersistentHandle(NULL, 0);
1296 #endif
1297  isNewProxyAdded = 0;
1298  }
1299  isProxyChanged = 0;
1300  if(!pidsPreAllocated) delete [] pids;
1301  DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
1302 
1303 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
1304  positionPtrBegin = p.begin();
1305  positionPtrEnd = p.end();
1306 #endif
1307 
1308  if(flags.doMolly) {
1311  }
1312  // BEGIN LA
1313  if (flags.doLoweAndersen) {
1314  velocityPtrBegin = v.begin();
1315  velocityPtrEnd = v.end();
1316  }
1317  // END LA
1318 
1319  Patch::positionsReady(doMigration);
1320 
1321  patchMapRead = 1;
1322 
1323  // gzheng
1324  Sync::Object()->PatchReady();
1325 
1326  NAMD_EVENT_STOP(1, NamdProfileEvent::POSITIONS_READY);
1327 
1328 }
1329 
1331 {
1332  replacementForces = f;
1333 }
1334 
1335 void HomePatch::saveForce(const int ftag)
1336 {
1337  f_saved[ftag].resize(numAtoms);
1338  for ( int i = 0; i < numAtoms; ++i )
1339  {
1340  f_saved[ftag][i] = f[ftag][i];
1341  }
1342 }
1343 
1344 
1345 #undef DEBUG_REDISTRIB_FORCE
1346 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
1347 //#define DEBUG_REDISTRIB_FORCE
1348 
1363 void HomePatch::redistrib_colinear_lp_force(
1364  Vector& fi, Vector& fj, Vector& fk,
1365  const Vector& ri, const Vector& rj, const Vector& rk,
1366  Real distance_f, Real scale_f) {
1367  BigReal distance = distance_f;
1368  BigReal scale = scale_f;
1369  Vector r_jk = rj - rk;
1370  // TODO: Add a check for small distances?
1371  BigReal r_jk_rlength = r_jk.rlength();
1372  distance *= r_jk_rlength;
1373  BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
1374  fj += (1. + scale + distance)*fi - r_jk*fdot;
1375  fk -= (scale + distance)*fi + r_jk*fdot;
1376 }
1377 
1403 #define FIX_FOR_WATER
1404 void HomePatch::redistrib_relative_lp_force(
1405  Vector& fi, Vector& fj, Vector& fk, Vector& fl,
1406  const Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
1407  Tensor *virial, int midpt) {
1408 #ifdef DEBUG_REDISTRIB_FORCE
1409  Vector foldnet, toldnet; // old net force, old net torque
1410  foldnet = fi + fj + fk + fl;
1411  toldnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
1412 #endif
1413  Vector fja(0), fka(0), fla(0);
1414 
1415  Vector r = ri - rj;
1416  BigReal invr2 = r.rlength();
1417  invr2 *= invr2;
1418  BigReal fdot = (fi*r) * invr2;
1419  Vector fr = r * fdot;
1420 
1421  fja += fr;
1422 
1423  Vector s, t;
1424  if (midpt) {
1425  s = rj - 0.5*(rk + rl);
1426  t = 0.5*(rk - rl);
1427  }
1428  else {
1429  s = rj - rk;
1430  t = rk - rl;
1431  }
1432  BigReal invs2 = s.rlength();
1433  invs2 *= invs2;
1434 
1435  Vector p = cross(r,s);
1436 #if !defined(FIX_FOR_WATER)
1437  BigReal invp = p.rlength();
1438 #else
1439  BigReal p2 = p.length2(); // fix division by zero above
1440 #endif
1441 
1442  Vector q = cross(s,t);
1443  BigReal invq = q.rlength();
1444 
1445 #if !defined(FIX_FOR_WATER)
1446  BigReal fpdot = (fi*p) * invp;
1447  Vector fp = p * fpdot;
1448  Vector ft = fi - fr - fp;
1449 #else
1450  BigReal fpdot;
1451  Vector fp, ft;
1452  if (p2 < 1e-6) { // vector is near zero, assume no fp contribution to force
1453  fpdot = 0;
1454  fp = 0;
1455  ft = fi - fr;
1456  }
1457  else {
1458  fpdot = (fi*p) / sqrt(p2);
1459  fp = p * fpdot;
1460  ft = fi - fr - fp;
1461  }
1462 #endif
1463 
1464  fja += ft;
1465  Vector v = cross(r,ft); // torque
1466  ft = cross(s,v) * invs2;
1467  fja -= ft;
1468 
1469  if (midpt) {
1470  fka += 0.5 * ft;
1471  fla += 0.5 * ft;
1472  }
1473  else {
1474  fka += ft;
1475  }
1476 
1477  BigReal srdot = (s*r) * invs2;
1478  Vector rr = r - s*srdot;
1479  BigReal rrdot = rr.length();
1480  BigReal stdot = (s*t) * invs2;
1481  Vector tt = t - s*stdot;
1482  BigReal invtt = tt.rlength();
1483  BigReal fact = rrdot*fpdot*invtt*invq;
1484  Vector fq = q * fact;
1485 
1486  fla += fq;
1487  fja += fp*(1+srdot) + fq*stdot;
1488 
1489  ft = fq*(1+stdot) + fp*srdot;
1490 
1491  if (midpt) {
1492  fka += -0.5*ft;
1493  fla += -0.5*ft;
1494  }
1495  else {
1496  fka -= ft;
1497  }
1498 
1499  if (virial) {
1500  Tensor va = outer(fja,rj);
1501  va += outer(fka,rk);
1502  va += outer(fla,rl);
1503  va -= outer(fi,ri);
1504  *virial += va;
1505  }
1506 
1507  fi = 0; // lone pair has zero force
1508  fj += fja;
1509  fk += fka;
1510  fl += fla;
1511 
1512 #ifdef DEBUG_REDISTRIB_FORCE
1513 #define TOL_REDISTRIB 1e-4
1514  Vector fnewnet, tnewnet; // new net force, new net torque
1515  fnewnet = fi + fj + fk + fl;
1516  tnewnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
1517  Vector fdiff = fnewnet - foldnet;
1518  Vector tdiff = tnewnet - toldnet;
1519  if (fdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1520  printf("Error: force redistribution for water exceeded tolerance: "
1521  "fdiff=(%f, %f, %f)\n", fdiff.x, fdiff.y, fdiff.z);
1522  }
1523  if (tdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1524  printf("Error: torque redistribution for water exceeded tolerance: "
1525  "tdiff=(%f, %f, %f)\n", tdiff.x, tdiff.y, tdiff.z);
1526  }
1527 #endif
1528 }
1529 
1530 void HomePatch::redistrib_ap_force(Vector& fi, Vector& fj) {
1531  // final state 'atom' force must transfer to initial state atom
1532  // in single topology FEP.
1533  Vector fi_old = fi;
1534  Vector fj_old = fj;
1535  fi = fi_old + fj_old;
1536  fj = fi_old + fj_old;
1537 }
1538 
1539 /* Redistribute forces from the massless lonepair charge particle onto
1540  * the other atoms of the water.
1541  *
1542  * This is done using the same algorithm as charmm uses for TIP4P lonepairs.
1543  *
1544  * Pass by reference the forces (O H1 H2 LP) to be modified,
1545  * pass by constant reference the corresponding positions,
1546  * and a pointer to virial.
1547  */
1548 void HomePatch::redistrib_lp_water_force(
1549  Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
1550  const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
1551  const Vector& p_lp, Tensor *virial) {
1552 
1553 #ifdef DEBUG_REDISTRIB_FORCE
1554  // Debug information to check against results at end
1555 
1556  // total force and torque relative to origin
1557  Vector totforce, tottorque;
1558  totforce = f_ox + f_h1 + f_h2 + f_lp;
1559  tottorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
1560  //printf("Torque without LP is %f/%f/%f\n",
1561  // tottorque.x, tottorque.y, tottorque.z);
1562  tottorque += cross(f_lp, p_lp);
1563  //printf("Torque with LP is %f/%f/%f\n",
1564  // tottorque.x, tottorque.y, tottorque.z);
1565 #endif
1566 
1567  // accumulate force adjustments
1568  Vector fad_ox(0), fad_h(0);
1569 
1570  // Calculate the radial component of the force and add it to the oxygen
1571  Vector r_ox_lp = p_lp - p_ox;
1572  BigReal invlen2_r_ox_lp = r_ox_lp.rlength();
1573  invlen2_r_ox_lp *= invlen2_r_ox_lp;
1574  BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
1575  Vector f_rad = r_ox_lp * rad_factor;
1576 
1577  fad_ox += f_rad;
1578 
1579  // Calculate the angular component
1580  Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
1581  Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5; // half of r_h2_h1
1582 
1583  // deviation from collinearity of charge site
1584  //Vector r_oop = cross(r_ox_lp, r_hcom_ox);
1585  //
1586  // vector out of o-h-h plane
1587  //Vector r_perp = cross(r_hcom_ox, r_h2_h1_2);
1588 
1589  // Here we assume that Ox/Lp/Hcom are linear
1590  // If you want to correct for deviations, this is the place
1591 
1592 // printf("Deviation from linearity for ox %i: %f/%f/%f\n", oxind, r_oop.x, r_oop.y, r_oop.z);
1593 
1594  Vector f_ang = f_lp - f_rad; // leave the angular component
1595 
1596  // now split this component onto the other atoms
1597  BigReal len_r_ox_lp = r_ox_lp.length();
1598  BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
1599  BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
1600  BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
1601 
1602  fad_ox += (f_ang * oxcomp);
1603  fad_h += (f_ang * hydcomp); // adjustment for both hydrogens
1604 
1605  // Add virial contributions
1606  if (virial) {
1607  Tensor vir = outer(fad_ox, p_ox);
1608  vir += outer(fad_h, p_h1);
1609  vir += outer(fad_h, p_h2);
1610  vir -= outer(f_lp, p_lp);
1611  *virial += vir;
1612  }
1613 
1614  //Vector zerovec(0.0, 0.0, 0.0);
1615  f_lp = 0;
1616  f_ox += fad_ox;
1617  f_h1 += fad_h;
1618  f_h2 += fad_h;
1619 
1620 #ifdef DEBUG_REDISTRIB_FORCE
1621  // Check that the total force and torque come out right
1622  Vector newforce, newtorque;
1623  newforce = f_ox + f_h1 + f_h2;
1624  newtorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
1625  Vector fdiff = newforce - totforce;
1626  Vector tdiff = newtorque - tottorque;
1627  BigReal error = fdiff.length();
1628  if (error > 0.0001) {
1629  printf("Error: Force redistribution for water "
1630  "exceeded force tolerance: error=%f\n", error);
1631  }
1632 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1633  printf("Error in net force: %f\n", error);
1634 #endif
1635 
1636  error = tdiff.length();
1637  if (error > 0.0001) {
1638  printf("Error: Force redistribution for water "
1639  "exceeded torque tolerance: error=%f\n", error);
1640  }
1641 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1642  printf("Error in net torque: %f\n", error);
1643 #endif
1644 #endif /* DEBUG */
1645 }
1646 
1665 void HomePatch::reposition_colinear_lonepair(
1666  Vector& ri, const Vector& rj, const Vector& rk,
1667  Real distance_f, Real scale_f)
1668 {
1669  BigReal distance = distance_f;
1670  BigReal scale = scale_f;
1671  Vector r_jk = rj - rk;
1672  BigReal r2 = r_jk.length2();
1673  if (r2 < 1e-10 || 100. < r2) { // same low tolerance as used in CHARMM
1674  iout << iWARN << "Large/small distance between lonepair reference atoms: ("
1675  << rj << ") (" << rk << ")\n" << endi;
1676  }
1677  ri = rj + (scale + distance*r_jk.rlength())*r_jk;
1678 }
1679 
1694 void HomePatch::reposition_relative_lonepair(
1695  Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
1696  Real distance, Real angle, Real dihedral)
1697 {
1698  if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
1699  iout << iWARN << "Large distance between lonepair reference atoms: ("
1700  << rj << ") (" << rk << ") (" << rl << ")\n" << endi;
1701  }
1702  BigReal r, t, p, cst, snt, csp, snp, invlen;
1703  Vector v, w, a, b, c;
1704 
1705  if (distance >= 0) {
1706  v = rk;
1707  r = distance;
1708  }
1709  else {
1710  v = 0.5*(rk + rl);
1711  r = -distance;
1712  }
1713 
1714  t = angle;
1715  p = dihedral;
1716  cst = cos(t);
1717  snt = sin(t);
1718  csp = cos(p);
1719  snp = sin(p);
1720  a = v - rj;
1721  b = rl - v;
1722  invlen = a.rlength();
1723  a *= invlen;
1724  c = cross(b, a);
1725  invlen = c.rlength();
1726  c *= invlen;
1727  b = cross(a, c);
1728  w.x = r*cst;
1729  w.y = r*snt*csp;
1730  w.z = r*snt*snp;
1731  ri.x = rj.x + w.x*a.x + w.y*b.x + w.z*c.x;
1732  ri.y = rj.y + w.x*a.y + w.y*b.y + w.z*c.y;
1733  ri.z = rj.z + w.x*a.z + w.y*b.z + w.z*c.z;
1734 }
1735 
1736 void HomePatch::reposition_alchpair(Vector& ri, Vector& rj, Mass& Mi, Mass& Mj) {
1737  Vector ri_old, rj_old;
1738  Mass mi, mj;
1739  ri_old.x = ri.x;
1740  ri_old.y = ri.y;
1741  ri_old.z = ri.z;
1742  rj_old.x = rj.x;
1743  rj_old.y = rj.y;
1744  rj_old.z = rj.z;
1745 
1746  mi = Mi; mj = Mj;
1747  ri.x = (mi * ri_old.x + mj * rj_old.x)/(mi + mj);
1748  ri.y = (mi * ri_old.y + mj * rj_old.y)/(mi + mj);
1749  ri.z = (mi * ri_old.z + mj * rj_old.z)/(mi + mj);
1750  rj.x = ri.x;
1751  rj.y = ri.y;
1752  rj.z = ri.z;
1753 }
1754 
1755 void HomePatch::reposition_all_lonepairs(void) {
1756  // ASSERT: simParams->lonepairs == TRUE
1757  for (int i=0; i < numAtoms; i++) {
1758  if (atom[i].mass < 0.01) {
1759  // found a lone pair
1760  AtomID aid = atom[i].id; // global atom ID of lp
1761  Lphost *lph = Node::Object()->molecule->get_lphost(aid); // its lphost
1762  if (lph == NULL) {
1763  char errmsg[512];
1764  sprintf(errmsg, "reposition lone pairs: "
1765  "no Lphost exists for LP %d\n", aid);
1766  NAMD_die(errmsg);
1767  }
1768  LocalID j = AtomMap::Object()->localID(lph->atom2);
1769  LocalID k = AtomMap::Object()->localID(lph->atom3);
1770  LocalID l = AtomMap::Object()->localID(lph->atom4);
1771  if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
1772  char errmsg[512];
1773  sprintf(errmsg, "reposition lone pairs: "
1774  "LP %d has some Lphost atom off patch\n", aid);
1775  NAMD_die(errmsg);
1776  }
1777  // reposition this lone pair
1778  if (lph->numhosts == 2) {
1779  reposition_colinear_lonepair(atom[i].position, atom[j.index].position,
1780  atom[k.index].position, lph->distance, lph->angle);
1781  }
1782  else if (lph->numhosts == 3) {
1783  reposition_relative_lonepair(atom[i].position, atom[j.index].position,
1784  atom[k.index].position, atom[l.index].position,
1785  lph->distance, lph->angle, lph->dihedral);
1786  }
1787  }
1788  }
1789 }
1790 
1791 void HomePatch::reposition_all_alchpairs(void) {
1792  Molecule *mol = Node::Object()->molecule;
1793  int numFepInitial = mol->numFepInitial;
1794  int alchPair_id;
1795  for (int i = 0; i < numAtoms; i++) {
1796  if (atom[i].partition == 4 ) {
1797  alchPair_id = atom[i].id + numFepInitial; //global id of the pair atom.
1798  LocalID j = AtomMap::Object()->localID(alchPair_id);
1799  reposition_alchpair(atom[i].position, atom[j.index].position, atom[i].mass, atom[j.index].mass);
1800  }
1801  }
1802 }
1803 
1804 void HomePatch::swm4_omrepos(Vector *ref, Vector *pos, Vector *vel,
1805  BigReal invdt) {
1806  // Reposition lonepair (Om) particle of Drude SWM4 water.
1807  // Same comments apply as to tip4_omrepos(), but the ordering of atoms
1808  // is different: O, D, LP, H1, H2.
1809  pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
1810  // Now, adjust velocity of particle to get it to appropriate place
1811  // during next integration "drift-step"
1812  if (invdt != 0) {
1813  vel[2] = (pos[2] - ref[2]) * invdt;
1814  }
1815  // No virial correction needed since lonepair is massless
1816 }
1817 
1818 void HomePatch::tip4_omrepos(Vector* ref, Vector* pos, Vector* vel, BigReal invdt) {
1819  /* Reposition the om particle of a tip4p water
1820  * A little geometry shows that the appropriate position is given by
1821  * R_O + (1 / 2 r_ohc) * ( 0.5 (R_H1 + R_H2) - R_O )
1822  * Here r_om is the distance from the oxygen to Om site, and r_ohc
1823  * is the altitude from the oxygen to the hydrogen center of mass
1824  * Those quantities are precalculated upon initialization of HomePatch
1825  *
1826  * Ordering of TIP4P atoms: O, H1, H2, LP.
1827  */
1828 
1829  //printf("rom/rohc are %f %f and invdt is %f\n", r_om, r_ohc, invdt);
1830  //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);
1831  pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc);
1832  //printf("New position for lp is %f %f %f\n", pos[3].x, pos[3].y, pos[3].z);
1833 
1834  // Now, adjust the velocity of the particle to get it to the appropriate place
1835  if (invdt != 0) {
1836  vel[3] = (pos[3] - ref[3]) * invdt;
1837  }
1838 
1839  // No virial correction needed, since this is a massless particle
1840  return;
1841 }
1842 
1843 void HomePatch::redistrib_lonepair_forces(const int ftag, Tensor *virial) {
1844  // ASSERT: simParams->lonepairs == TRUE
1845  ForceList *f_mod = f;
1846  for (int i = 0; i < numAtoms; i++) {
1847  if (atom[i].mass < 0.01) {
1848  // found a lone pair
1849  AtomID aid = atom[i].id; // global atom ID of lp
1850  Lphost *lph = Node::Object()->molecule->get_lphost(aid); // its lphost
1851  if (lph == NULL) {
1852  char errmsg[512];
1853  sprintf(errmsg, "redistrib lone pair forces: "
1854  "no Lphost exists for LP %d\n", aid);
1855  NAMD_die(errmsg);
1856  }
1857  LocalID j = AtomMap::Object()->localID(lph->atom2);
1858  LocalID k = AtomMap::Object()->localID(lph->atom3);
1859  LocalID l = AtomMap::Object()->localID(lph->atom4);
1860  if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
1861  char errmsg[512];
1862  sprintf(errmsg, "redistrib lone pair forces: "
1863  "LP %d has some Lphost atom off patch\n", aid);
1864  NAMD_die(errmsg);
1865  }
1866  // redistribute forces from this lone pair
1867  if (lph->numhosts == 2) {
1868  redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
1869  f_mod[ftag][k.index], atom[i].position, atom[j.index].position,
1870  atom[k.index].position, lph->distance, lph->angle);
1871  }
1872  else if (lph->numhosts == 3) {
1873  int midpt = (lph->distance < 0);
1874  redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
1875  f_mod[ftag][k.index], f_mod[ftag][l.index],
1876  atom[i].position, atom[j.index].position,
1877  atom[k.index].position, atom[l.index].position, virial, midpt);
1878  }
1879  }
1880  }
1881 }
1882 
1883 void HomePatch::redistrib_alchpair_forces(const int ftag) { //Virial unchanged
1884  SimParameters *simParams = Node::Object()->simParameters;
1885  Molecule *mol = Node::Object()->molecule;
1886  ForceList *f_mod = f;
1887  int numFepInitial = mol->numFepInitial;
1888  BigReal lambda = simParams->alchLambda;
1889  int alchPair_id;
1890  for (int i = 0; i < numAtoms; i++) {
1891  if (atom[i].partition == 4 ) {
1892  alchPair_id = atom[i].id + numFepInitial; //global id of the pair atom.
1893  LocalID j = AtomMap::Object()->localID(alchPair_id);
1894  redistrib_ap_force(f_mod[ftag][i],f_mod[ftag][j.index]);
1895  }
1896  }
1897 }
1898 
1899 void HomePatch::redistrib_swm4_forces(const int ftag, Tensor *virial) {
1900  // Loop over the patch's atoms and apply the appropriate corrections
1901  // to get all forces off of lone pairs
1902  ForceList *f_mod = f;
1903  for (int i = 0; i < numAtoms; i++) {
1904  if (atom[i].mass < 0.01) {
1905  // found lonepair
1906  redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
1907  f_mod[ftag][i+2], f_mod[ftag][i],
1908  atom[i-2].position, atom[i+1].position,
1909  atom[i+2].position, atom[i].position, virial);
1910  }
1911  }
1912 }
1913 
1914 void HomePatch::redistrib_tip4p_forces(const int ftag, Tensor* virial) {
1915  // Loop over the patch's atoms and apply the appropriate corrections
1916  // to get all forces off of lone pairs
1917  // Atom ordering: O H1 H2 LP
1918 
1919  ForceList *f_mod =f;
1920  for (int i=0; i<numAtoms; i++) {
1921  if (atom[i].mass < 0.01) {
1922  // found lonepair
1923  redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
1924  f_mod[ftag][i-1], f_mod[ftag][i],
1925  atom[i-3].position, atom[i-2].position,
1926  atom[i-1].position, atom[i].position, virial);
1927  }
1928  }
1929 }
1930 
1931 
1933  FullAtom * __restrict atom_arr,
1934  const Force * __restrict force_arr,
1935  const BigReal dt,
1936  int num_atoms
1937  ) {
1938  SimParameters *simParams = Node::Object()->simParameters;
1939 
1940  if ( simParams->fixedAtomsOn ) {
1941  for ( int i = 0; i < num_atoms; ++i ) {
1942  if ( atom_arr[i].atomFixed ) {
1943  atom_arr[i].velocity = 0;
1944  } else {
1945  BigReal dt_mass = dt * atom_arr[i].recipMass; // dt/mass
1946  atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1947  atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1948  atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1949  }
1950  }
1951  } else {
1952  for ( int i = 0; i < num_atoms; ++i ) {
1953  BigReal dt_mass = dt * atom_arr[i].recipMass; // dt/mass
1954  atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1955  atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1956  atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1957  }
1958  }
1959 }
1960 
1962  FullAtom * __restrict atom_arr,
1963  const Force * __restrict force_arr1,
1964  const Force * __restrict force_arr2,
1965  const Force * __restrict force_arr3,
1966  const BigReal dt1,
1967  const BigReal dt2,
1968  const BigReal dt3,
1969  int num_atoms
1970  ) {
1971  SimParameters *simParams = Node::Object()->simParameters;
1972 
1973  if ( simParams->fixedAtomsOn ) {
1974  for ( int i = 0; i < num_atoms; ++i ) {
1975  if ( atom_arr[i].atomFixed ) {
1976  atom_arr[i].velocity = 0;
1977  } else {
1978  BigReal rmass = atom_arr[i].recipMass; // 1/mass
1979  atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1980  + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1981  atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1982  + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1983  atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1984  + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1985  }
1986  }
1987  } else {
1988  for ( int i = 0; i < num_atoms; ++i ) {
1989  BigReal rmass = atom_arr[i].recipMass; // 1/mass
1990  atom_arr[i].velocity.x += (force_arr1[i].x*dt1
1991  + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1992  atom_arr[i].velocity.y += (force_arr1[i].y*dt1
1993  + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1994  atom_arr[i].velocity.z += (force_arr1[i].z*dt1
1995  + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1996  }
1997  }
1998 }
1999 
2001  FullAtom * __restrict atom_arr,
2002  const BigReal dt,
2003  int num_atoms
2004  ) {
2005  SimParameters *simParams = Node::Object()->simParameters;
2006  if ( simParams->fixedAtomsOn ) {
2007  for ( int i = 0; i < num_atoms; ++i ) {
2008  if ( ! atom_arr[i].atomFixed ) {
2009  atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2010  atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2011  atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2012  }
2013  }
2014  } else {
2015  for ( int i = 0; i < num_atoms; ++i ) {
2016  atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
2017  atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
2018  atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
2019  }
2020  }
2021 }
2022 
2023 int HomePatch::hardWallDrude(const BigReal timestep, Tensor *virial,
2024  SubmitReduction *ppreduction)
2025 {
2026  Molecule *mol = Node::Object()->molecule;
2027  SimParameters *simParams = Node::Object()->simParameters;
2028  const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
2029  const int fixedAtomsOn = simParams->fixedAtomsOn;
2030  const BigReal dt = timestep / TIMEFACTOR;
2031  const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2032  int i, ia, ib, j;
2033  int dieOnError = simParams->rigidDie;
2034  Tensor wc; // constraint virial
2035  BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
2036  int nslabs;
2037 
2038  // start data for hard wall boundary between drude and its host atom
2039  // static int Count=0;
2040  int Idx;
2041  double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
2042  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;
2043  double dot_v_r_1, dot_v_r_2;
2044  double vb_cm, dr_a, dr_b;
2045  // end data for hard wall boundary between drude and its host atom
2046 
2047  // start calculation of hard wall boundary between drude and its host atom
2048  if (simParams->drudeHardWallOn) {
2049  if (ppreduction) {
2050  nslabs = simParams->pressureProfileSlabs;
2051  idz = nslabs/lattice.c().z;
2052  zmin = lattice.origin().z - 0.5*lattice.c().z;
2053  }
2054 
2055  r_wall = simParams->drudeBondLen;
2056  r_wall_SQ = r_wall*r_wall;
2057  // Count++;
2058  for (i=1; i<numAtoms; i++) {
2059  if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) { // drude particle
2060  ia = i-1;
2061  ib = i;
2062 
2063  v_ab = atom[ib].position - atom[ia].position;
2064  rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
2065 
2066  if (rab_SQ > r_wall_SQ) { // to impose the hard wall constraint
2067  rab = sqrt(rab_SQ);
2068  if ( (rab > (2.0*r_wall)) && dieOnError ) { // unexpected situation
2069  iout << iERROR << "HardWallDrude> "
2070  << "The drude is too far away from atom "
2071  << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
2072  return -1; // triggers early exit
2073  }
2074 
2075  v_ab.x /= rab;
2076  v_ab.y /= rab;
2077  v_ab.z /= rab;
2078 
2079  if ( fixedAtomsOn && atom[ia].atomFixed ) { // the heavy atom is fixed
2080  if (atom[ib].atomFixed) { // the drude is fixed too
2081  continue;
2082  }
2083  else { // only the heavy atom is fixed
2084  dot_v_r_2 = atom[ib].velocity.x*v_ab.x
2085  + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
2086  vb_2 = dot_v_r_2 * v_ab;
2087  vp_2 = atom[ib].velocity - vb_2;
2088 
2089  dr = rab - r_wall;
2090  if(dot_v_r_2 == 0.0) {
2091  delta_T = maxtime;
2092  }
2093  else {
2094  delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
2095  if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
2096  }
2097 
2098  dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
2099 
2100  vb_2 = dot_v_r_2 * v_ab;
2101 
2102  new_vel_a = atom[ia].velocity;
2103  new_vel_b = vp_2 + vb_2;
2104 
2105  dr_b = -dr + delta_T*dot_v_r_2; // L = L_0 + dT *v_new, v was flipped
2106 
2107  new_pos_a = atom[ia].position;
2108  new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
2109  }
2110  }
2111  else {
2112  mass_a = atom[ia].mass;
2113  mass_b = atom[ib].mass;
2114  mass_sum = mass_a+mass_b;
2115 
2116  dot_v_r_1 = atom[ia].velocity.x*v_ab.x
2117  + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
2118  vb_1 = dot_v_r_1 * v_ab;
2119  vp_1 = atom[ia].velocity - vb_1;
2120 
2121  dot_v_r_2 = atom[ib].velocity.x*v_ab.x
2122  + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
2123  vb_2 = dot_v_r_2 * v_ab;
2124  vp_2 = atom[ib].velocity - vb_2;
2125 
2126  vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
2127 
2128  dot_v_r_1 -= vb_cm;
2129  dot_v_r_2 -= vb_cm;
2130 
2131  dr = rab - r_wall;
2132 
2133  if(dot_v_r_2 == dot_v_r_1) {
2134  delta_T = maxtime;
2135  }
2136  else {
2137  delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1); // the time since the collision occurs
2138  if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
2139  }
2140 
2141  // the relative velocity between ia and ib. Drawn according to T_Drude
2142  v_Bond = sqrt(kbt/mass_b);
2143 
2144  // reflect the velocity along bond vector and scale down
2145  dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
2146  dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
2147 
2148  dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
2149  dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
2150 
2151  new_pos_a = atom[ia].position + dr_a*v_ab; // correct the position
2152  new_pos_b = atom[ib].position + dr_b*v_ab;
2153  // atom[ia].position += (dr_a*v_ab); // correct the position
2154  // atom[ib].position += (dr_b*v_ab);
2155 
2156  dot_v_r_1 += vb_cm;
2157  dot_v_r_2 += vb_cm;
2158 
2159  vb_1 = dot_v_r_1 * v_ab;
2160  vb_2 = dot_v_r_2 * v_ab;
2161 
2162  new_vel_a = vp_1 + vb_1;
2163  new_vel_b = vp_2 + vb_2;
2164  }
2165 
2166  int ppoffset, partition;
2167  if ( invdt == 0 ) {
2168  atom[ia].position = new_pos_a;
2169  atom[ib].position = new_pos_b;
2170  }
2171  else if ( virial == 0 ) {
2172  atom[ia].velocity = new_vel_a;
2173  atom[ib].velocity = new_vel_b;
2174  }
2175  else {
2176  for ( j = 0; j < 2; j++ ) {
2177  if (j==0) { // atom ia, heavy atom
2178  Idx = ia;
2179  new_pos = &new_pos_a;
2180  new_vel = &new_vel_a;
2181  }
2182  else if (j==1) { // atom ib, drude
2183  Idx = ib;
2184  new_pos = &new_pos_b;
2185  new_vel = &new_vel_b;
2186  }
2187  Force df = (*new_vel - atom[Idx].velocity) *
2188  ( atom[Idx].mass * invdt );
2189  Tensor vir = outer(df, atom[Idx].position);
2190  wc += vir;
2191  atom[Idx].velocity = *new_vel;
2192  atom[Idx].position = *new_pos;
2193 
2194  if (ppreduction) {
2195  if (!j) {
2196  BigReal z = new_pos->z;
2197  int partition = atom[Idx].partition;
2198  int slab = (int)floor((z-zmin)*idz);
2199  if (slab < 0) slab += nslabs;
2200  else if (slab >= nslabs) slab -= nslabs;
2201  ppoffset = 3*(slab + nslabs*partition);
2202  }
2203  ppreduction->item(ppoffset ) += vir.xx;
2204  ppreduction->item(ppoffset+1) += vir.yy;
2205  ppreduction->item(ppoffset+2) += vir.zz;
2206  }
2207 
2208  }
2209  }
2210  }
2211  }
2212  }
2213 
2214  // if ( (Count>10000) && (Count%10==0) ) {
2215  // v_ab = atom[1].position - atom[0].position;
2216  // rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
2217  // iout << "DBG_R: " << Count << " " << sqrt(rab_SQ) << "\n" << endi;
2218  // }
2219 
2220  }
2221 
2222  // end calculation of hard wall boundary between drude and its host atom
2223 
2224  if ( dt && virial ) *virial += wc;
2225 
2226  return 0;
2227 }
2228 
2230 
2231  SimParameters *simParams = Node::Object()->simParameters;
2232  const int fixedAtomsOn = simParams->fixedAtomsOn;
2233  const int useSettle = simParams->useSettle;
2234 
2235  // Re-size to containt numAtoms elements
2236  velNew.resize(numAtoms);
2237  posNew.resize(numAtoms);
2238 
2239  // Size of a hydrogen group for water
2240  int wathgsize = 3;
2241  int watmodel = simParams->watmodel;
2242  if (watmodel == WAT_TIP4) wathgsize = 4;
2243  else if (watmodel == WAT_SWM4) wathgsize = 5;
2244 
2245  // Initialize the settle algorithm with water parameters
2246  // settle1() assumes all waters are identical,
2247  // and will generate bad results if they are not.
2248  // XXX this will move to Molecule::build_atom_status when that
2249  // version is debugged
2250  if ( ! settle_initialized ) {
2251  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2252  // find a water
2253  if (atom[ig].rigidBondLength > 0) {
2254  int oatm;
2255  if (simParams->watmodel == WAT_SWM4) {
2256  oatm = ig+3; // skip over Drude and Lonepair
2257  //printf("ig=%d mass_ig=%g oatm=%d mass_oatm=%g\n",
2258  // ig, atom[ig].mass, oatm, atom[oatm].mass);
2259  }
2260  else {
2261  oatm = ig+1;
2262  // Avoid using the Om site to set this by mistake
2263  if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2264  oatm += 1;
2265  }
2266  }
2267 
2268  // initialize settle water parameters
2269  settle1init(atom[ig].mass, atom[oatm].mass,
2270  atom[ig].rigidBondLength,
2271  atom[oatm].rigidBondLength,
2272  settle_mOrmT, settle_mHrmT, settle_ra,
2273  settle_rb, settle_rc, settle_rra);
2274  settle_initialized = 1;
2275  break; // done with init
2276  }
2277  }
2278  }
2279 
2280  Vector ref[10];
2281  BigReal rmass[10];
2282  BigReal dsq[10];
2283  int fixed[10];
2284  int ial[10];
2285  int ibl[10];
2286 
2287  int numSettle = 0;
2288  int numRattle = 0;
2289  int posRattleParam = 0;
2290 
2291  settleList.clear();
2292  rattleList.clear();
2293  noconstList.clear();
2294  rattleParam.clear();
2295 
2296  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2297  int hgs = atom[ig].hydrogenGroupSize;
2298  if ( hgs == 1 ) {
2299  // only one atom in group
2300  noconstList.push_back(ig);
2301  continue;
2302  }
2303  int anyfixed = 0;
2304  for (int i = 0; i < hgs; ++i ) {
2305  ref[i] = atom[ig+i].position;
2306  rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2307  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2308  if ( fixed[i] ) {
2309  anyfixed = 1;
2310  rmass[i] = 0.;
2311  }
2312  }
2313  int icnt = 0;
2314  BigReal tmp = atom[ig].rigidBondLength;
2315  if (tmp > 0.0) { // for water
2316  if (hgs != wathgsize) {
2317  char errmsg[256];
2318  sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
2319  "but the specified water model requires %d atoms.\n",
2320  atom[ig].id+1, hgs, wathgsize);
2321  NAMD_die(errmsg);
2322  }
2323  // Use SETTLE for water unless some of the water atoms are fixed,
2324  if (useSettle && !anyfixed) {
2325  // Store to Settle -list
2326  settleList.push_back(ig);
2327  continue;
2328  }
2329  if ( !(fixed[1] && fixed[2]) ) {
2330  dsq[icnt] = tmp * tmp;
2331  ial[icnt] = 1;
2332  ibl[icnt] = 2;
2333  ++icnt;
2334  }
2335  } // if (tmp > 0.0)
2336  for (int i = 1; i < hgs; ++i ) { // normal bonds to mother atom
2337  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2338  if ( !(fixed[0] && fixed[i]) ) {
2339  dsq[icnt] = tmp * tmp;
2340  ial[icnt] = 0;
2341  ibl[icnt] = i;
2342  ++icnt;
2343  }
2344  }
2345  }
2346  if ( icnt == 0 ) {
2347  // no constraints
2348  noconstList.push_back(ig);
2349  continue;
2350  }
2351  // Store to Rattle -list
2352  RattleList rattleListElem;
2353  rattleListElem.ig = ig;
2354  rattleListElem.icnt = icnt;
2355  rattleList.push_back(rattleListElem);
2356  for (int i = 0; i < icnt; ++i ) {
2357  int a = ial[i];
2358  int b = ibl[i];
2359  RattleParam rattleParamElem;
2360  rattleParamElem.ia = a;
2361  rattleParamElem.ib = b;
2362  rattleParamElem.dsq = dsq[i];
2363  rattleParamElem.rma = rmass[a];
2364  rattleParamElem.rmb = rmass[b];
2365  rattleParam.push_back(rattleParamElem);
2366  }
2367  }
2368 
2369 }
2370 
2371 void HomePatch::addRattleForce(const BigReal invdt, Tensor& wc) {
2372  for (int ig = 0; ig < numAtoms; ++ig ) {
2373  Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
2374  Tensor vir = outer(df, atom[ig].position);
2375  wc += vir;
2376  f[Results::normal][ig] += df;
2377  atom[ig].velocity = velNew[ig];
2378  }
2379 }
2380 
2381 int HomePatch::rattle1(const BigReal timestep, Tensor *virial,
2382  SubmitReduction *ppreduction) {
2383 
2384  SimParameters *simParams = Node::Object()->simParameters;
2385  if (simParams->watmodel != WAT_TIP3 || ppreduction) {
2386  // Call old rattle1 -method instead
2387  return rattle1old(timestep, virial, ppreduction);
2388  }
2389 
2390  if (!rattleListValid) {
2391  buildRattleList();
2392  rattleListValid = true;
2393  }
2394 
2395  const int fixedAtomsOn = simParams->fixedAtomsOn;
2396  const int useSettle = simParams->useSettle;
2397  const BigReal dt = timestep / TIMEFACTOR;
2398  const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2399  const BigReal tol2 = 2.0 * simParams->rigidTol;
2400  int maxiter = simParams->rigidIter;
2401  int dieOnError = simParams->rigidDie;
2402 
2403  Vector ref[10]; // reference position
2404  Vector pos[10]; // new position
2405  Vector vel[10]; // new velocity
2406 
2407  // Manual un-roll
2408  int n = (settleList.size()/2)*2;
2409  for (int j=0;j < n;j+=2) {
2410  int ig;
2411  ig = settleList[j];
2412  for (int i = 0; i < 3; ++i ) {
2413  ref[i] = atom[ig+i].position;
2414  pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2415  }
2416  ig = settleList[j+1];
2417  for (int i = 0; i < 3; ++i ) {
2418  ref[i+3] = atom[ig+i].position;
2419  pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
2420  }
2421  settle1_SIMD<2>(ref, pos,
2422  settle_mOrmT, settle_mHrmT, settle_ra,
2423  settle_rb, settle_rc, settle_rra);
2424 
2425  ig = settleList[j];
2426  for (int i = 0; i < 3; ++i ) {
2427  velNew[ig+i] = (pos[i] - ref[i])*invdt;
2428  posNew[ig+i] = pos[i];
2429  }
2430  ig = settleList[j+1];
2431  for (int i = 0; i < 3; ++i ) {
2432  velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
2433  posNew[ig+i] = pos[i+3];
2434  }
2435 
2436  }
2437 
2438  if (settleList.size() % 2) {
2439  int ig = settleList[settleList.size()-1];
2440  for (int i = 0; i < 3; ++i ) {
2441  ref[i] = atom[ig+i].position;
2442  pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2443  }
2444  settle1_SIMD<1>(ref, pos,
2445  settle_mOrmT, settle_mHrmT, settle_ra,
2446  settle_rb, settle_rc, settle_rra);
2447  for (int i = 0; i < 3; ++i ) {
2448  velNew[ig+i] = (pos[i] - ref[i])*invdt;
2449  posNew[ig+i] = pos[i];
2450  }
2451  }
2452 
2453  int posParam = 0;
2454  for (int j=0;j < rattleList.size();++j) {
2455 
2456  BigReal refx[10];
2457  BigReal refy[10];
2458  BigReal refz[10];
2459 
2460  BigReal posx[10];
2461  BigReal posy[10];
2462  BigReal posz[10];
2463 
2464  int ig = rattleList[j].ig;
2465  int icnt = rattleList[j].icnt;
2466  int hgs = atom[ig].hydrogenGroupSize;
2467  for (int i = 0; i < hgs; ++i ) {
2468  ref[i] = atom[ig+i].position;
2469  pos[i] = atom[ig+i].position;
2470  if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
2471  pos[i] += atom[ig+i].velocity * dt;
2472  }
2473  refx[i] = ref[i].x;
2474  refy[i] = ref[i].y;
2475  refz[i] = ref[i].z;
2476  posx[i] = pos[i].x;
2477  posy[i] = pos[i].y;
2478  posz[i] = pos[i].z;
2479  }
2480 
2481  bool done;
2482  bool consFailure;
2483  if (icnt == 1) {
2484  rattlePair<1>(&rattleParam[posParam],
2485  refx, refy, refz,
2486  posx, posy, posz,
2487  consFailure);
2488  done = true;
2489  } else {
2490  rattleN(icnt, &rattleParam[posParam],
2491  refx, refy, refz,
2492  posx, posy, posz,
2493  tol2, maxiter,
2494  done, consFailure);
2495  }
2496 
2497  // Advance position in rattleParam
2498  posParam += icnt;
2499 
2500  for (int i = 0; i < hgs; ++i ) {
2501  pos[i].x = posx[i];
2502  pos[i].y = posy[i];
2503  pos[i].z = posz[i];
2504  }
2505 
2506  for (int i = 0; i < hgs; ++i ) {
2507  velNew[ig+i] = (pos[i] - ref[i])*invdt;
2508  posNew[ig+i] = pos[i];
2509  }
2510 
2511  if ( consFailure ) {
2512  if ( dieOnError ) {
2513  iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
2514  << (atom[ig].id + 1) << "!\n" << endi;
2515  return -1; // triggers early exit
2516  } else {
2517  iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
2518  << (atom[ig].id + 1) << "!\n" << endi;
2519  }
2520  } else if ( ! done ) {
2521  if ( dieOnError ) {
2522  iout << iERROR << "Exceeded RATTLE iteration limit for atom "
2523  << (atom[ig].id + 1) << "!\n" << endi;
2524  return -1; // triggers early exit
2525  } else {
2526  iout << iWARN << "Exceeded RATTLE iteration limit for atom "
2527  << (atom[ig].id + 1) << "!\n" << endi;
2528  }
2529  }
2530  }
2531  // Finally, we have to go through atoms that are not involved in rattle just so that we have
2532  // their positions and velocities up-to-date in posNew and velNew
2533  for (int j=0;j < noconstList.size();++j) {
2534  int ig = noconstList[j];
2535  int hgs = atom[ig].hydrogenGroupSize;
2536  for (int i = 0; i < hgs; ++i ) {
2537  velNew[ig+i] = atom[ig+i].velocity;
2538  posNew[ig+i] = atom[ig+i].position;
2539  }
2540  }
2541 
2542  if ( invdt == 0 ) {
2543  for (int ig = 0; ig < numAtoms; ++ig )
2544  atom[ig].position = posNew[ig];
2545  } else if ( virial == 0 ) {
2546  for (int ig = 0; ig < numAtoms; ++ig )
2547  atom[ig].velocity = velNew[ig];
2548  } else {
2549  Tensor wc; // constraint virial
2550  addRattleForce(invdt, wc);
2551  *virial += wc;
2552  }
2553 
2554  return 0;
2555 }
2556 
2557 // RATTLE algorithm from Allen & Tildesley
2558 int HomePatch::rattle1old(const BigReal timestep, Tensor *virial,
2559  SubmitReduction *ppreduction)
2560 {
2561  Molecule *mol = Node::Object()->molecule;
2562  SimParameters *simParams = Node::Object()->simParameters;
2563  const int fixedAtomsOn = simParams->fixedAtomsOn;
2564  const int useSettle = simParams->useSettle;
2565  const BigReal dt = timestep / TIMEFACTOR;
2566  const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2567  BigReal tol2 = 2.0 * simParams->rigidTol;
2568  int maxiter = simParams->rigidIter;
2569  int dieOnError = simParams->rigidDie;
2570  int i, iter;
2571  BigReal dsq[10], tmp;
2572  int ial[10], ibl[10];
2573  Vector ref[10]; // reference position
2574  Vector refab[10]; // reference vector
2575  Vector pos[10]; // new position
2576  Vector vel[10]; // new velocity
2577  Vector netdp[10]; // total momentum change from constraint
2578  BigReal rmass[10]; // 1 / mass
2579  int fixed[10]; // is atom fixed?
2580  Tensor wc; // constraint virial
2581  BigReal idz, zmin;
2582  int nslabs;
2583 
2584  // Initialize the settle algorithm with water parameters
2585  // settle1() assumes all waters are identical,
2586  // and will generate bad results if they are not.
2587  // XXX this will move to Molecule::build_atom_status when that
2588  // version is debugged
2589  if ( ! settle_initialized ) {
2590  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2591  // find a water
2592  if (atom[ig].rigidBondLength > 0) {
2593  int oatm;
2594  if (simParams->watmodel == WAT_SWM4) {
2595  oatm = ig+3; // skip over Drude and Lonepair
2596  //printf("ig=%d mass_ig=%g oatm=%d mass_oatm=%g\n",
2597  // ig, atom[ig].mass, oatm, atom[oatm].mass);
2598  }
2599  else {
2600  oatm = ig+1;
2601  // Avoid using the Om site to set this by mistake
2602  if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2603  oatm += 1;
2604  }
2605  }
2606 
2607  // initialize settle water parameters
2608  settle1init(atom[ig].mass, atom[oatm].mass,
2609  atom[ig].rigidBondLength,
2610  atom[oatm].rigidBondLength,
2611  settle_mOrmT, settle_mHrmT, settle_ra,
2612  settle_rb, settle_rc, settle_rra);
2613  settle_initialized = 1;
2614  break; // done with init
2615  }
2616  }
2617  }
2618 
2619  if (ppreduction) {
2620  nslabs = simParams->pressureProfileSlabs;
2621  idz = nslabs/lattice.c().z;
2622  zmin = lattice.origin().z - 0.5*lattice.c().z;
2623  }
2624 
2625  // Size of a hydrogen group for water
2626  int wathgsize = 3;
2627  int watmodel = simParams->watmodel;
2628  if (watmodel == WAT_TIP4) wathgsize = 4;
2629  else if (watmodel == WAT_SWM4) wathgsize = 5;
2630 
2631  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2632  int hgs = atom[ig].hydrogenGroupSize;
2633  if ( hgs == 1 ) continue; // only one atom in group
2634  // cache data in local arrays and integrate positions normally
2635  int anyfixed = 0;
2636  for ( i = 0; i < hgs; ++i ) {
2637  ref[i] = atom[ig+i].position;
2638  pos[i] = atom[ig+i].position;
2639  vel[i] = atom[ig+i].velocity;
2640  rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2641  //printf("rmass of %i is %f\n", ig+i, rmass[i]);
2642  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2643  //printf("fixed status of %i is %i\n", i, fixed[i]);
2644  // undo addVelocityToPosition to get proper reference coordinates
2645  if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
2646  }
2647  int icnt = 0;
2648  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
2649  if (hgs != wathgsize) {
2650  char errmsg[256];
2651  sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
2652  "but the specified water model requires %d atoms.\n",
2653  atom[ig].id+1, hgs, wathgsize);
2654  NAMD_die(errmsg);
2655  }
2656  // Use SETTLE for water unless some of the water atoms are fixed,
2657  if (useSettle && !anyfixed) {
2658  if (simParams->watmodel == WAT_SWM4) {
2659  // SWM4 ordering: O D LP H1 H2
2660  // do swap(O,LP) and call settle with subarray O H1 H2
2661  // swap back after we return
2662  Vector lp_ref = ref[2];
2663  Vector lp_pos = pos[2];
2664  Vector lp_vel = vel[2];
2665  ref[2] = ref[0];
2666  pos[2] = pos[0];
2667  vel[2] = vel[0];
2668  settle1(ref+2, pos+2, vel+2, invdt,
2669  settle_mOrmT, settle_mHrmT, settle_ra,
2670  settle_rb, settle_rc, settle_rra);
2671  ref[0] = ref[2];
2672  pos[0] = pos[2];
2673  vel[0] = vel[2];
2674  ref[2] = lp_ref;
2675  pos[2] = lp_pos;
2676  vel[2] = lp_vel;
2677  // determine for LP updated pos and vel
2678  swm4_omrepos(ref, pos, vel, invdt);
2679  }
2680  else {
2681  settle1(ref, pos, vel, invdt,
2682  settle_mOrmT, settle_mHrmT, settle_ra,
2683  settle_rb, settle_rc, settle_rra);
2684  if (simParams->watmodel == WAT_TIP4) {
2685  tip4_omrepos(ref, pos, vel, invdt);
2686  }
2687  }
2688 
2689  // which slab the hydrogen group will belong to
2690  // for pprofile calculations.
2691  int ppoffset, partition;
2692  if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
2693  atom[ig+i].position = pos[i];
2694  } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
2695  atom[ig+i].velocity = vel[i];
2696  } else for ( i = 0; i < wathgsize; ++i ) {
2697  Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
2698  Tensor vir = outer(df, ref[i]);
2699  wc += vir;
2700  f[Results::normal][ig+i] += df;
2701  atom[ig+i].velocity = vel[i];
2702  if (ppreduction) {
2703  // put all the atoms from a water in the same slab. Atom 0
2704  // should be the parent atom.
2705  if (!i) {
2706  BigReal z = pos[i].z;
2707  partition = atom[ig].partition;
2708  int slab = (int)floor((z-zmin)*idz);
2709  if (slab < 0) slab += nslabs;
2710  else if (slab >= nslabs) slab -= nslabs;
2711  ppoffset = 3*(slab + nslabs*partition);
2712  }
2713  ppreduction->item(ppoffset ) += vir.xx;
2714  ppreduction->item(ppoffset+1) += vir.yy;
2715  ppreduction->item(ppoffset+2) += vir.zz;
2716  }
2717  }
2718  continue;
2719  }
2720  if ( !(fixed[1] && fixed[2]) ) {
2721  dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2722  }
2723  }
2724  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
2725  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2726  if ( !(fixed[0] && fixed[i]) ) {
2727  dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
2728  }
2729  }
2730  }
2731  if ( icnt == 0 ) continue; // no constraints
2732  for ( i = 0; i < icnt; ++i ) {
2733  refab[i] = ref[ial[i]] - ref[ibl[i]];
2734  }
2735  for ( i = 0; i < hgs; ++i ) {
2736  netdp[i] = 0.;
2737  }
2738  int done;
2739  int consFailure;
2740  for ( iter = 0; iter < maxiter; ++iter ) {
2741 //if (iter > 0) CkPrintf("iteration %d\n", iter);
2742  done = 1;
2743  consFailure = 0;
2744  for ( i = 0; i < icnt; ++i ) {
2745  int a = ial[i]; int b = ibl[i];
2746  Vector pab = pos[a] - pos[b];
2747  BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
2748  BigReal rabsq = dsq[i];
2749  BigReal diffsq = rabsq - pabsq;
2750  if ( fabs(diffsq) > (rabsq * tol2) ) {
2751  Vector &rab = refab[i];
2752  BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
2753  if ( rpab < ( rabsq * 1.0e-6 ) ) {
2754  done = 0;
2755  consFailure = 1;
2756  continue;
2757  }
2758  BigReal rma = rmass[a];
2759  BigReal rmb = rmass[b];
2760  BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
2761  Vector dp = rab * gab;
2762  pos[a] += rma * dp;
2763  pos[b] -= rmb * dp;
2764  if ( invdt != 0. ) {
2765  dp *= invdt;
2766  netdp[a] += dp;
2767  netdp[b] -= dp;
2768  }
2769  done = 0;
2770  }
2771  }
2772  if ( done ) break;
2773  }
2774 
2775  if ( consFailure ) {
2776  if ( dieOnError ) {
2777  iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
2778  << (atom[ig].id + 1) << "!\n" << endi;
2779  return -1; // triggers early exit
2780  } else {
2781  iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
2782  << (atom[ig].id + 1) << "!\n" << endi;
2783  }
2784  } else if ( ! done ) {
2785  if ( dieOnError ) {
2786  iout << iERROR << "Exceeded RATTLE iteration limit for atom "
2787  << (atom[ig].id + 1) << "!\n" << endi;
2788  return -1; // triggers early exit
2789  } else {
2790  iout << iWARN << "Exceeded RATTLE iteration limit for atom "
2791  << (atom[ig].id + 1) << "!\n" << endi;
2792  }
2793  }
2794 
2795  // store data back to patch
2796  int ppoffset, partition;
2797  if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
2798  atom[ig+i].position = pos[i];
2799  } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
2800  atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2801  } else for ( i = 0; i < hgs; ++i ) {
2802  Force df = netdp[i] * invdt;
2803  Tensor vir = outer(df, ref[i]);
2804  wc += vir;
2805  f[Results::normal][ig+i] += df;
2806  atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2807  if (ppreduction) {
2808  if (!i) {
2809  BigReal z = pos[i].z;
2810  int partition = atom[ig].partition;
2811  int slab = (int)floor((z-zmin)*idz);
2812  if (slab < 0) slab += nslabs;
2813  else if (slab >= nslabs) slab -= nslabs;
2814  ppoffset = 3*(slab + nslabs*partition);
2815  }
2816  ppreduction->item(ppoffset ) += vir.xx;
2817  ppreduction->item(ppoffset+1) += vir.yy;
2818  ppreduction->item(ppoffset+2) += vir.zz;
2819  }
2820  }
2821  }
2822  if ( dt && virial ) *virial += wc;
2823 
2824  return 0;
2825 }
2826 
2827 // RATTLE algorithm from Allen & Tildesley
2828 void HomePatch::rattle2(const BigReal timestep, Tensor *virial)
2829 {
2830  Molecule *mol = Node::Object()->molecule;
2831  SimParameters *simParams = Node::Object()->simParameters;
2832  const int fixedAtomsOn = simParams->fixedAtomsOn;
2833  const int useSettle = simParams->useSettle;
2834  const BigReal dt = timestep / TIMEFACTOR;
2835  Tensor wc; // constraint virial
2836  BigReal tol = simParams->rigidTol;
2837  int maxiter = simParams->rigidIter;
2838  int dieOnError = simParams->rigidDie;
2839  int i, iter;
2840  BigReal dsqi[10], tmp;
2841  int ial[10], ibl[10];
2842  Vector ref[10]; // reference position
2843  Vector refab[10]; // reference vector
2844  Vector vel[10]; // new velocity
2845  BigReal rmass[10]; // 1 / mass
2846  BigReal redmass[10]; // reduced mass
2847  int fixed[10]; // is atom fixed?
2848 
2849  // Size of a hydrogen group for water
2850  int wathgsize = 3;
2851  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
2852  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
2853 
2854  // CkPrintf("In rattle2!\n");
2855  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2856  // CkPrintf("ig=%d\n",ig);
2857  int hgs = atom[ig].hydrogenGroupSize;
2858  if ( hgs == 1 ) continue; // only one atom in group
2859  // cache data in local arrays and integrate positions normally
2860  int anyfixed = 0;
2861  for ( i = 0; i < hgs; ++i ) {
2862  ref[i] = atom[ig+i].position;
2863  vel[i] = atom[ig+i].velocity;
2864  rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
2865  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2866  if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
2867  }
2868  int icnt = 0;
2869  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
2870  if (hgs != wathgsize) {
2871  NAMD_bug("Hydrogen group error caught in rattle2().");
2872  }
2873  // Use SETTLE for water unless some of the water atoms are fixed,
2874  if (useSettle && !anyfixed) {
2875  if (simParams->watmodel == WAT_SWM4) {
2876  // SWM4 ordering: O D LP H1 H2
2877  // do swap(O,LP) and call settle with subarray O H1 H2
2878  // swap back after we return
2879  Vector lp_ref = ref[2];
2880  // Vector lp_vel = vel[2];
2881  ref[2] = ref[0];
2882  vel[2] = vel[0];
2883  settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
2884  ref[0] = ref[2];
2885  vel[0] = vel[2];
2886  ref[2] = lp_ref;
2887  // vel[2] = vel[0]; // set LP vel to O vel
2888  }
2889  else {
2890  settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
2891  if (simParams->watmodel == WAT_TIP4) {
2892  vel[3] = vel[0];
2893  }
2894  }
2895  for (i=0; i<hgs; i++) {
2896  atom[ig+i].velocity = vel[i];
2897  }
2898  continue;
2899  }
2900  if ( !(fixed[1] && fixed[2]) ) {
2901  redmass[icnt] = 1. / (rmass[1] + rmass[2]);
2902  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
2903  }
2904  }
2905  // CkPrintf("Loop 2\n");
2906  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
2907  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2908  if ( !(fixed[0] && fixed[i]) ) {
2909  redmass[icnt] = 1. / (rmass[0] + rmass[i]);
2910  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
2911  ibl[icnt] = i; ++icnt;
2912  }
2913  }
2914  }
2915  if ( icnt == 0 ) continue; // no constraints
2916  // CkPrintf("Loop 3\n");
2917  for ( i = 0; i < icnt; ++i ) {
2918  refab[i] = ref[ial[i]] - ref[ibl[i]];
2919  }
2920  // CkPrintf("Loop 4\n");
2921  int done;
2922  for ( iter = 0; iter < maxiter; ++iter ) {
2923  done = 1;
2924  for ( i = 0; i < icnt; ++i ) {
2925  int a = ial[i]; int b = ibl[i];
2926  Vector vab = vel[a] - vel[b];
2927  Vector &rab = refab[i];
2928  BigReal rabsqi = dsqi[i];
2929  BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
2930  if ( (fabs(rvab) * dt * rabsqi) > tol ) {
2931  Vector dp = rab * (-rvab * redmass[i] * rabsqi);
2932  wc += outer(dp,rab);
2933  vel[a] += rmass[a] * dp;
2934  vel[b] -= rmass[b] * dp;
2935  done = 0;
2936  }
2937  }
2938  if ( done ) break;
2939  //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
2940  }
2941  if ( ! done ) {
2942  if ( dieOnError ) {
2943  NAMD_die("Exceeded maximum number of iterations in rattle2().");
2944  } else {
2945  iout << iWARN <<
2946  "Exceeded maximum number of iterations in rattle2().\n" << endi;
2947  }
2948  }
2949  // store data back to patch
2950  for ( i = 0; i < hgs; ++i ) {
2951  atom[ig+i].velocity = vel[i];
2952  }
2953  }
2954  // CkPrintf("Leaving rattle2!\n");
2955  // check that there isn't a constant needed here!
2956  *virial += wc / ( 0.5 * dt );
2957 
2958 }
2959 
2960 
2961 // Adjust gradients for minimizer
2962 void HomePatch::minimize_rattle2(const BigReal timestep, Tensor *virial, bool forces)
2963 {
2964  Molecule *mol = Node::Object()->molecule;
2965  SimParameters *simParams = Node::Object()->simParameters;
2966  Force *f1 = f[Results::normal].begin();
2967  const int fixedAtomsOn = simParams->fixedAtomsOn;
2968  const int useSettle = simParams->useSettle;
2969  const BigReal dt = timestep / TIMEFACTOR;
2970  Tensor wc; // constraint virial
2971  BigReal tol = simParams->rigidTol;
2972  int maxiter = simParams->rigidIter;
2973  int dieOnError = simParams->rigidDie;
2974  int i, iter;
2975  BigReal dsqi[10], tmp;
2976  int ial[10], ibl[10];
2977  Vector ref[10]; // reference position
2978  Vector refab[10]; // reference vector
2979  Vector vel[10]; // new velocity
2980  BigReal rmass[10]; // 1 / mass
2981  BigReal redmass[10]; // reduced mass
2982  int fixed[10]; // is atom fixed?
2983 
2984  // Size of a hydrogen group for water
2985  int wathgsize = 3;
2986  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
2987  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
2988 
2989  // CkPrintf("In rattle2!\n");
2990  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2991  // CkPrintf("ig=%d\n",ig);
2992  int hgs = atom[ig].hydrogenGroupSize;
2993  if ( hgs == 1 ) continue; // only one atom in group
2994  // cache data in local arrays and integrate positions normally
2995  int anyfixed = 0;
2996  for ( i = 0; i < hgs; ++i ) {
2997  ref[i] = atom[ig+i].position;
2998  vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
2999  rmass[i] = 1.0;
3000  fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
3001  if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
3002  }
3003  int icnt = 0;
3004  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
3005  if (hgs != wathgsize) {
3006  NAMD_bug("Hydrogen group error caught in rattle2().");
3007  }
3008  // Use SETTLE for water unless some of the water atoms are fixed,
3009  if (useSettle && !anyfixed) {
3010  if (simParams->watmodel == WAT_SWM4) {
3011  // SWM4 ordering: O D LP H1 H2
3012  // do swap(O,LP) and call settle with subarray O H1 H2
3013  // swap back after we return
3014  Vector lp_ref = ref[2];
3015  // Vector lp_vel = vel[2];
3016  ref[2] = ref[0];
3017  vel[2] = vel[0];
3018  settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
3019  ref[0] = ref[2];
3020  vel[0] = vel[2];
3021  ref[2] = lp_ref;
3022  // vel[2] = vel[0]; // set LP vel to O vel
3023  }
3024  else {
3025  settle2(1.0, 1.0, ref, vel, dt, virial);
3026  if (simParams->watmodel == WAT_TIP4) {
3027  vel[3] = vel[0];
3028  }
3029  }
3030  for (i=0; i<hgs; i++) {
3031  ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3032  }
3033  continue;
3034  }
3035  if ( !(fixed[1] && fixed[2]) ) {
3036  redmass[icnt] = 1. / (rmass[1] + rmass[2]);
3037  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3038  }
3039  }
3040  // CkPrintf("Loop 2\n");
3041  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
3042  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3043  if ( !(fixed[0] && fixed[i]) ) {
3044  redmass[icnt] = 1. / (rmass[0] + rmass[i]);
3045  dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
3046  ibl[icnt] = i; ++icnt;
3047  }
3048  }
3049  }
3050  if ( icnt == 0 ) continue; // no constraints
3051  // CkPrintf("Loop 3\n");
3052  for ( i = 0; i < icnt; ++i ) {
3053  refab[i] = ref[ial[i]] - ref[ibl[i]];
3054  }
3055  // CkPrintf("Loop 4\n");
3056  int done;
3057  for ( iter = 0; iter < maxiter; ++iter ) {
3058  done = 1;
3059  for ( i = 0; i < icnt; ++i ) {
3060  int a = ial[i]; int b = ibl[i];
3061  Vector vab = vel[a] - vel[b];
3062  Vector &rab = refab[i];
3063  BigReal rabsqi = dsqi[i];
3064  BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
3065  if ( (fabs(rvab) * dt * rabsqi) > tol ) {
3066  Vector dp = rab * (-rvab * redmass[i] * rabsqi);
3067  wc += outer(dp,rab);
3068  vel[a] += rmass[a] * dp;
3069  vel[b] -= rmass[b] * dp;
3070  done = 0;
3071  }
3072  }
3073  if ( done ) break;
3074  //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
3075  }
3076  if ( ! done ) {
3077  if ( dieOnError ) {
3078  NAMD_die("Exceeded maximum number of iterations in rattle2().");
3079  } else {
3080  iout << iWARN <<
3081  "Exceeded maximum number of iterations in rattle2().\n" << endi;
3082  }
3083  }
3084  // store data back to patch
3085  for ( i = 0; i < hgs; ++i ) {
3086  ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
3087  }
3088  }
3089  // CkPrintf("Leaving rattle2!\n");
3090  // check that there isn't a constant needed here!
3091  *virial += wc / ( 0.5 * dt );
3092 
3093 }
3094 
3095 
3096 // BEGIN LA
3098 {
3099  DebugM(2, "loweAndersenVelocities\n");
3100  Molecule *mol = Node::Object()->molecule;
3101  SimParameters *simParams = Node::Object()->simParameters;
3102  v.resize(numAtoms);
3103  for (int i = 0; i < numAtoms; ++i) {
3104  //v[i] = p[i];
3105  // co-opt CompAtom structure to pass velocity and mass information
3106  v[i].position = atom[i].velocity;
3107  v[i].charge = atom[i].mass;
3108  }
3109  DebugM(2, "loweAndersenVelocities\n");
3110 }
3111 
3113 {
3114  DebugM(2, "loweAndersenFinish\n");
3115  v.resize(0);
3116 }
3117 // END LA
3118 
3119 //LCPO
3121  Molecule *mol = Node::Object()->molecule;
3122  const int *lcpoParamType = mol->getLcpoParamType();
3123 
3124  lcpoType.resize(numAtoms);
3125  for (int i = 0; i < numAtoms; i++) {
3126  lcpoType[i] = lcpoParamType[pExt[i].id];
3127  }
3128 }
3129 
3130 //set intrinsic radii of atom when doMigration
3132  intRad.resize(numAtoms*2);
3133  intRad.setall(0);
3134  Molecule *mol = Node::Object()->molecule;
3135  SimParameters *simParams = Node::Object()->simParameters;
3136  Real offset = simParams->coulomb_radius_offset;
3137  for (int i = 0; i < numAtoms; i++) {
3138  Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
3139  Real screen = MassToScreen(atom[i].mass);//same
3140  intRad[2*i+0] = rad - offset;//r0
3141  intRad[2*i+1] = screen*(rad - offset);//s0
3142  }
3143 }
3144 
3145 //compute born radius after phase 1, before phase 2
3147 
3148  SimParameters *simParams = Node::Object()->simParameters;
3149  BigReal alphaMax = simParams->alpha_max;
3150  BigReal delta = simParams->gbis_delta;
3151  BigReal beta = simParams->gbis_beta;
3152  BigReal gamma = simParams->gbis_gamma;
3153  BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
3154 
3155  BigReal rhoi;
3156  BigReal rhoi0;
3157  //calculate bornRad from psiSum
3158  for (int i = 0; i < numAtoms; i++) {
3159  rhoi0 = intRad[2*i];
3160  rhoi = rhoi0+coulomb_radius_offset;
3161  psiFin[i] += psiSum[i];
3162  psiFin[i] *= rhoi0;
3163  bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
3164  bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
3165 #ifdef PRINT_COMP
3166  CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
3167 #endif
3168  }
3169 
3170  gbisP2Ready();
3171 }
3172 
3173 //compute dHdrPrefix after phase 2, before phase 3
3175 
3176  SimParameters *simParams = Node::Object()->simParameters;
3177  BigReal delta = simParams->gbis_delta;
3178  BigReal beta = simParams->gbis_beta;
3179  BigReal gamma = simParams->gbis_gamma;
3180  BigReal epsilon_s = simParams->solvent_dielectric;
3181  BigReal epsilon_p = simParams->dielectric;
3182  BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
3183  BigReal epsilon_p_i = 1/simParams->dielectric;
3184  BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
3185  BigReal kappa = simParams->kappa;
3186  BigReal fij, expkappa, Dij, dEdai, dedasum;
3187  BigReal rhoi, rhoi0, psii, nbetapsi;
3188  BigReal gammapsi2, tanhi, daidr;
3189  for (int i = 0; i < numAtoms; i++) {
3190  //add diagonal dEda term
3191  dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
3192  fij = bornRad[i];//inf
3193  expkappa = exp(-kappa*fij);//0
3194  Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
3195  //calculate dHij prefix
3196  dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
3197  *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
3198  dHdrPrefix[i] += dEdai;
3199  dedasum = dHdrPrefix[i];
3200 
3201  rhoi0 = intRad[2*i];
3202  rhoi = rhoi0+coulomb_radius_offset;
3203  psii = psiFin[i];
3204  nbetapsi = -beta*psii;
3205  gammapsi2 = gamma*psii*psii;
3206  tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
3207  daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
3208  * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
3209  dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
3210 #ifdef PRINT_COMP
3211  CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
3212 #endif
3213  }
3214  gbisP3Ready();
3215 }
3216 
3217 //send born radius to proxies to begin phase 2
3219  if (proxy.size() > 0) {
3220  CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3221  for (int i = 0; i < proxy.size(); i++) {
3222  int node = proxy[i];
3224  msg->patch = patchID;
3225  msg->origPe = CkMyPe();
3226  memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
3227  msg->destPe = node;
3228  int seq = flags.sequence;
3230  SET_PRIORITY(msg,seq,priority);
3231  cp[node].recvData(msg);
3232  }
3233  }
3235 }
3236 
3237 //send dHdrPrefix to proxies to begin phase 3
3239  if (proxy.size() > 0) {
3240  CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3241  //only nonzero message should be sent for doFullElec
3242  int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
3243  for (int i = 0; i < proxy.size(); i++) {
3244  int node = proxy[i];
3246  msg->patch = patchID;
3247  msg->dHdrPrefixLen = msgAtoms;
3248  msg->origPe = CkMyPe();
3249  memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
3250  msg->destPe = node;
3251  int seq = flags.sequence;
3253  SET_PRIORITY(msg,seq,priority);
3254  cp[node].recvData(msg);
3255  }
3256  }
3258 }
3259 
3260 //receive proxy results from phase 1
3262  ++numGBISP1Arrived;
3263  for ( int i = 0; i < msg->psiSumLen; ++i ) {
3264  psiFin[i] += msg->psiSum[i];
3265  }
3266  delete msg;
3267 
3268  if (flags.doNonbonded) {
3269  //awaken if phase 1 done
3270  if (phase1BoxClosedCalled == true &&
3271  numGBISP1Arrived==proxy.size() ) {
3272  sequencer->awaken();
3273  }
3274  } else {
3275  //awaken if all phases done on noWork step
3276  if (boxesOpen == 0 &&
3277  numGBISP1Arrived == proxy.size() &&
3278  numGBISP2Arrived == proxy.size() &&
3279  numGBISP3Arrived == proxy.size()) {
3280  sequencer->awaken();
3281  }
3282  }
3283 }
3284 
3285 //receive proxy results from phase 2
3287  ++numGBISP2Arrived;
3288  //accumulate dEda
3289  for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
3290  dHdrPrefix[i] += msg->dEdaSum[i];
3291  }
3292  delete msg;
3293 
3294  if (flags.doNonbonded) {
3295  //awaken if phase 2 done
3296  if (phase2BoxClosedCalled == true &&
3297  numGBISP2Arrived==proxy.size() ) {
3298  sequencer->awaken();
3299  }
3300  } else {
3301  //awaken if all phases done on noWork step
3302  if (boxesOpen == 0 &&
3303  numGBISP1Arrived == proxy.size() &&
3304  numGBISP2Arrived == proxy.size() &&
3305  numGBISP3Arrived == proxy.size()) {
3306  sequencer->awaken();
3307  }
3308  }
3309 }
3310 
3311 // MOLLY algorithm part 1
3313 {
3314  Molecule *mol = Node::Object()->molecule;
3315  SimParameters *simParams = Node::Object()->simParameters;
3316  BigReal tol = simParams->mollyTol;
3317  int maxiter = simParams->mollyIter;
3318  int i, iter;
3319  HGArrayBigReal dsq;
3320  BigReal tmp;
3321  HGArrayInt ial, ibl;
3322  HGArrayVector ref; // reference position
3323  HGArrayVector refab; // reference vector
3324  HGArrayBigReal rmass; // 1 / mass
3325  BigReal *lambda; // Lagrange multipliers
3326  CompAtom *avg; // averaged position
3327  int numLambdas = 0;
3328  HGArrayInt fixed; // is atom fixed?
3329 
3330  // iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
3331  p_avg.resize(numAtoms);
3332  for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
3333 
3334  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3335  int hgs = atom[ig].hydrogenGroupSize;
3336  if ( hgs == 1 ) continue; // only one atom in group
3337  for ( i = 0; i < hgs; ++i ) {
3338  ref[i] = atom[ig+i].position;
3339  rmass[i] = 1. / atom[ig+i].mass;
3340  fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
3341  if ( fixed[i] ) rmass[i] = 0.;
3342  }
3343  avg = &(p_avg[ig]);
3344  int icnt = 0;
3345 
3346  if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) { // for water
3347  if ( hgs != 3 ) {
3348  NAMD_die("Hydrogen group error caught in mollyAverage(). It's a bug!\n");
3349  }
3350  if ( !(fixed[1] && fixed[2]) ) {
3351  dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3352  }
3353  }
3354  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
3355  if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3356  if ( !(fixed[0] && fixed[i]) ) {
3357  dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3358  }
3359  }
3360  }
3361  if ( icnt == 0 ) continue; // no constraints
3362  numLambdas += icnt;
3363  molly_lambda.resize(numLambdas);
3364  lambda = &(molly_lambda[numLambdas - icnt]);
3365  for ( i = 0; i < icnt; ++i ) {
3366  refab[i] = ref[ial[i]] - ref[ibl[i]];
3367  }
3368  // iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
3369  iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
3370  if ( iter == maxiter ) {
3371  iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
3372  }
3373  }
3374 
3375  // for ( i=0; i<numAtoms; ++i ) {
3376  // if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
3377  // iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
3378  // << p[i].position << " to " << p_avg[i].position << "\n" << endi;
3379  // }
3380  // }
3381 
3382 }
3383 
3384 
3385 // MOLLY algorithm part 2
3387 {
3388  Molecule *mol = Node::Object()->molecule;
3389  SimParameters *simParams = Node::Object()->simParameters;
3390  Tensor wc; // constraint virial
3391  int i;
3392  HGArrayInt ial, ibl;
3393  HGArrayVector ref; // reference position
3394  CompAtom *avg; // averaged position
3395  HGArrayVector refab; // reference vector
3396  HGArrayVector force; // new force
3397  HGArrayBigReal rmass; // 1 / mass
3398  BigReal *lambda; // Lagrange multipliers
3399  int numLambdas = 0;
3400  HGArrayInt fixed; // is atom fixed?
3401 
3402  for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3403  int hgs = atom[ig].hydrogenGroupSize;
3404  if (hgs == 1 ) continue; // only one atom in group
3405  for ( i = 0; i < hgs; ++i ) {
3406  ref[i] = atom[ig+i].position;
3407  force[i] = f[Results::slow][ig+i];
3408  rmass[i] = 1. / atom[ig+i].mass;
3409  fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
3410  if ( fixed[i] ) rmass[i] = 0.;
3411  }
3412  int icnt = 0;
3413  // c-ji I'm only going to mollify water for now
3414  if ( atom[ig].rigidBondLength > 0 ) { // for water
3415  if ( hgs != 3 ) {
3416  NAMD_die("Hydrogen group error caught in mollyMollify(). It's a bug!\n");
3417  }
3418  if ( !(fixed[1] && fixed[2]) ) {
3419  ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
3420  }
3421  }
3422  for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
3423  if ( atom[ig+i].rigidBondLength > 0 ) {
3424  if ( !(fixed[0] && fixed[i]) ) {
3425  ial[icnt] = 0; ibl[icnt] = i; ++icnt;
3426  }
3427  }
3428  }
3429 
3430  if ( icnt == 0 ) continue; // no constraints
3431  lambda = &(molly_lambda[numLambdas]);
3432  numLambdas += icnt;
3433  for ( i = 0; i < icnt; ++i ) {
3434  refab[i] = ref[ial[i]] - ref[ibl[i]];
3435  }
3436  avg = &(p_avg[ig]);
3437  mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
3438  // store data back to patch
3439  for ( i = 0; i < hgs; ++i ) {
3440  wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
3441  f[Results::slow][ig+i] = force[i];
3442  }
3443  }
3444  // check that there isn't a constant needed here!
3445  *virial += wc;
3446  p_avg.resize(0);
3447 }
3448 
3450  checkpoint_atom.copy(atom);
3451  checkpoint_lattice = lattice;
3452 
3453  // DMK - Atom Separation (water vs. non-water)
3454  #if NAMD_SeparateWaters != 0
3455  checkpoint_numWaterAtoms = numWaterAtoms;
3456  #endif
3457 }
3458 
3459 void HomePatch::revert(void) {
3460  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3461 
3462  atom.copy(checkpoint_atom);
3463  numAtoms = atom.size();
3464  lattice = checkpoint_lattice;
3465 
3466  doAtomUpdate = true;
3467  rattleListValid = false;
3468 
3469  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3470 
3471  // DMK - Atom Separation (water vs. non-water)
3472  #if NAMD_SeparateWaters != 0
3473  numWaterAtoms = checkpoint_numWaterAtoms;
3474  #endif
3475 }
3476 
3477 void HomePatch::exchangeCheckpoint(int scriptTask, int &bpc) { // initiating replica
3478  SimParameters *simParams = Node::Object()->simParameters;
3479  checkpoint_task = scriptTask;
3480  const int remote = simParams->scriptIntArg1;
3481  const char *key = simParams->scriptStringArg1;
3482  PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
3483 }
3484 
3485 void HomePatch::recvCheckpointReq(int task, const char *key, int replica, int pe) { // responding replica
3486  if ( task == SCRIPT_CHECKPOINT_FREE ) {
3487  if ( ! checkpoints.count(key) ) {
3488  NAMD_die("Unable to free checkpoint, requested key was never stored.");
3489  }
3490  delete checkpoints[key];
3491  checkpoints.erase(key);
3492  PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
3493  return;
3494  }
3495  CheckpointAtomsMsg *msg;
3496  if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
3497  if ( ! checkpoints.count(key) ) {
3498  NAMD_die("Unable to load checkpoint, requested key was never stored.");
3499  }
3500  checkpoint_t &cp = *checkpoints[key];
3501  msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
3502  msg->lattice = cp.lattice;
3504  msg->numAtoms = cp.numAtoms;
3505  memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
3506  } else {
3507  msg = new (0,1,0) CheckpointAtomsMsg;
3508  }
3509  msg->pid = patchID;
3510  msg->replica = CmiMyPartition();
3511  msg->pe = CkMyPe();
3512  PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
3513 }
3514 
3515 void HomePatch::recvCheckpointLoad(CheckpointAtomsMsg *msg) { // initiating replica
3516  SimParameters *simParams = Node::Object()->simParameters;
3517  const int remote = simParams->scriptIntArg1;
3518  const char *key = simParams->scriptStringArg1;
3520  NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
3521  }
3522  if ( msg->replica != remote ) {
3523  NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
3524  }
3526  CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
3527  strcpy(newmsg->key,key);
3528  newmsg->lattice = lattice;
3529  newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
3530  newmsg->pid = patchID;
3531  newmsg->pe = CkMyPe();
3532  newmsg->replica = CmiMyPartition();
3533  newmsg->numAtoms = numAtoms;
3534  memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
3535  PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
3536  }
3538  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3539  lattice = msg->lattice;
3541  numAtoms = msg->numAtoms;
3542  atom.resize(numAtoms);
3543  memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
3544  doAtomUpdate = true;
3545  rattleListValid = false;
3546  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3547  }
3550  }
3551  delete msg;
3552 }
3553 
3554 void HomePatch::recvCheckpointStore(CheckpointAtomsMsg *msg) { // responding replica
3555  if ( ! checkpoints.count(msg->key) ) {
3556  checkpoints[msg->key] = new checkpoint_t;
3557  }
3558  checkpoint_t &cp = *checkpoints[msg->key];
3559  cp.lattice = msg->lattice;
3561  cp.numAtoms = msg->numAtoms;
3562  cp.atoms.resize(cp.numAtoms);
3563  memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
3565  delete msg;
3566 }
3567 
3568 void HomePatch::recvCheckpointAck() { // initiating replica
3569  CkpvAccess(_qd)->process();
3570 }
3571 
3572 
3573 void HomePatch::exchangeAtoms(int scriptTask) {
3574  SimParameters *simParams = Node::Object()->simParameters;
3575  // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
3576  if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
3577  exchange_dst = (int) simParams->scriptArg1;
3578  // create and save outgoing message
3579  exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
3583  memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
3584  if ( exchange_req >= 0 ) {
3586  }
3587  }
3588  if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
3589  exchange_src = (int) simParams->scriptArg2;
3591  }
3592 }
3593 
3595  exchange_req = req;
3596  if ( exchange_msg ) {
3597  // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
3599  exchange_msg = 0;
3600  exchange_req = -1;
3601  CkpvAccess(_qd)->process();
3602  }
3603 }
3604 
3606  // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
3607  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3608  lattice = msg->lattice;
3609  numAtoms = msg->numAtoms;
3610  atom.resize(numAtoms);
3611  memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
3612  delete msg;
3613  CkpvAccess(_qd)->process();
3614  doAtomUpdate = true;
3615  rattleListValid = false;
3616  if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3617 }
3618 
3619 void HomePatch::submitLoadStats(int timestep)
3620 {
3621  LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
3622 }
3623 
3624 
3625 //
3626 // XXX operates on CompAtom, not FullAtom
3627 //
3629 {
3630 #if 0
3631 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_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 #endif
3637 
3638  SimParameters *simParams = Node::Object()->simParameters;
3639 
3640  if ( numAtoms == 0 || ! flags.usePairlists ) {
3641  flags.pairlistTolerance = 0.;
3642  flags.maxAtomMovement = 99999.;
3643 #if 0
3644  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
3645 #endif
3646  return;
3647  }
3648 
3649  int i; int n = numAtoms;
3650  CompAtom *p_i = p.begin();
3651 
3652  if ( flags.savePairlists ) {
3653  flags.pairlistTolerance = doPairlistCheck_newTolerance;
3654  flags.maxAtomMovement = 0.;
3655  doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
3656  doPairlistCheck_lattice = lattice;
3657  doPairlistCheck_positions.resize(numAtoms);
3658  CompAtom *psave_i = doPairlistCheck_positions.begin();
3659  for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
3660 #if 0
3661  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
3662 #endif
3663  return;
3664  }
3665 
3666  Lattice &lattice_old = doPairlistCheck_lattice;
3667  Position center_cur = lattice.unscale(center);
3668  Position center_old = lattice_old.unscale(center);
3669  Vector center_delta = center_cur - center_old;
3670 
3671  // find max deviation to corner (any neighbor shares a corner)
3672  BigReal max_cd = 0.;
3673  for ( i=0; i<2; ++i ) {
3674  for ( int j=0; j<2; ++j ) {
3675  for ( int k=0; k<2; ++k ) {
3676  ScaledPosition corner( i ? min.x : max.x ,
3677  j ? min.y : max.y ,
3678  k ? min.z : max.z );
3679  Vector corner_delta =
3680  lattice.unscale(corner) - lattice_old.unscale(corner);
3681  corner_delta -= center_delta;
3682  BigReal cd = corner_delta.length2();
3683  if ( cd > max_cd ) max_cd = cd;
3684  }
3685  }
3686  }
3687  max_cd = sqrt(max_cd);
3688 
3689  // find max deviation of atoms relative to center
3690  BigReal max_pd = 0.;
3691  CompAtom *p_old_i = doPairlistCheck_positions.begin();
3692  for ( i=0; i<n; ++i ) {
3693  Vector p_delta = p_i[i].position - p_old_i[i].position;
3694  p_delta -= center_delta;
3695  BigReal pd = p_delta.length2();
3696  if ( pd > max_pd ) max_pd = pd;
3697  }
3698  max_pd = sqrt(max_pd);
3699 
3700  BigReal max_tol = max_pd + max_cd;
3701 
3702  flags.maxAtomMovement = max_tol;
3703 
3704  // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
3705 
3706  if ( max_tol > ( (1. - simParams->pairlistTrigger) *
3707  doPairlistCheck_newTolerance ) ) {
3708  doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
3709  }
3710 
3711  if ( max_tol > doPairlistCheck_newTolerance ) {
3712  doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
3713  }
3714 #if 0
3715  NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
3716 #endif
3717 }
3718 
3720 {
3721  if ( ! flags.doNonbonded ) return;
3722 
3723  SimParameters *simParams = Node::Object()->simParameters;
3724  BigReal hgcut = 0.5 * simParams->hgroupCutoff; hgcut *= hgcut;
3725  BigReal maxrad2 = 0.;
3726 
3727  FullAtomList::iterator p_i = atom.begin();
3728  FullAtomList::iterator p_e = atom.end();
3729 
3730  while ( p_i != p_e ) {
3731  const int hgs = p_i->hydrogenGroupSize;
3732  if ( ! hgs ) break; // avoid infinite loop on bug
3733  int ngs = hgs;
3734  if ( ngs > 5 ) ngs = 5; // XXX why? limit to at most 5 atoms per group
3735  BigReal x = p_i->position.x;
3736  BigReal y = p_i->position.y;
3737  BigReal z = p_i->position.z;
3738  int i;
3739  for ( i = 1; i < ngs; ++i ) { // limit spatial extent
3740  p_i[i].nonbondedGroupSize = 0;
3741  BigReal dx = p_i[i].position.x - x;
3742  BigReal dy = p_i[i].position.y - y;
3743  BigReal dz = p_i[i].position.z - z;
3744  BigReal r2 = dx * dx + dy * dy + dz * dz;
3745  if ( r2 > hgcut ) break;
3746  else if ( r2 > maxrad2 ) maxrad2 = r2;
3747  }
3748  p_i->nonbondedGroupSize = i;
3749  for ( ; i < hgs; ++i ) {
3750  p_i[i].nonbondedGroupSize = 1;
3751  }
3752  p_i += hgs;
3753  }
3754 
3755  if ( p_i != p_e ) {
3756  NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
3757  }
3758 
3759  flags.maxGroupRadius = sqrt(maxrad2);
3760 
3761 }
3762 
3764 {
3765  SimParameters *simParams = Node::Object()->simParameters;
3766 
3767  BigReal sysdima = lattice.a_r().unit() * lattice.a();
3768  BigReal sysdimb = lattice.b_r().unit() * lattice.b();
3769  BigReal sysdimc = lattice.c_r().unit() * lattice.c();
3770 
3771  BigReal minSize = simParams->patchDimension - simParams->margin;
3772 
3773  if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
3774  ( bAwayDist*sysdimb < minSize*0.9999 ) ||
3775  ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
3776 
3777  NAMD_die("Periodic cell has become too small for original patch grid!\n"
3778  "Possible solutions are to restart from a recent checkpoint,\n"
3779  "increase margin, or disable useFlexibleCell for liquid simulation.");
3780  }
3781 
3782  BigReal cutoff = simParams->cutoff;
3783 
3784  BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
3785  BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
3786  BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
3787 
3788  if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
3789  NAMD_die("Periodic cell has become too small for original patch grid!\n"
3790  "There are probably many margin violations already reported.\n"
3791  "Possible solutions are to restart from a recent checkpoint,\n"
3792  "increase margin, or disable useFlexibleCell for liquid simulation.");
3793  }
3794 
3795  BigReal minx = min.x - margina;
3796  BigReal miny = min.y - marginb;
3797  BigReal minz = min.z - marginc;
3798  BigReal maxx = max.x + margina;
3799  BigReal maxy = max.y + marginb;
3800  BigReal maxz = max.z + marginc;
3801 
3802  int xdev, ydev, zdev;
3803  int problemCount = 0;
3804 
3805  FullAtomList::iterator p_i = atom.begin();
3806  FullAtomList::iterator p_e = atom.end();
3807  for ( ; p_i != p_e; ++p_i ) {
3808 
3809  ScaledPosition s = lattice.scale(p_i->position);
3810 
3811  // check if atom is within bounds
3812  if (s.x < minx) xdev = 0;
3813  else if (maxx <= s.x) xdev = 2;
3814  else xdev = 1;
3815 
3816  if (s.y < miny) ydev = 0;
3817  else if (maxy <= s.y) ydev = 2;
3818  else ydev = 1;
3819 
3820  if (s.z < minz) zdev = 0;
3821  else if (maxz <= s.z) zdev = 2;
3822  else zdev = 1;
3823 
3824  if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
3825  ++problemCount;
3826  }
3827 
3828  }
3829 
3830  marginViolations = problemCount;
3831  // if ( problemCount ) {
3832  // iout << iERROR <<
3833  // "Found " << problemCount << " margin violations!\n" << endi;
3834  // }
3835 
3836 }
3837 
3838 
3839 void
3841 {
3842  int i;
3843 
3844  for (i=0; i<numNeighbors; i++) {
3845  realInfo[i].mList.resize(0);
3846  }
3847 
3848  // Purge the AtomMap
3849  atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3850 
3851  // Determine atoms that need to migrate
3852 
3853  BigReal minx = min.x;
3854  BigReal miny = min.y;
3855  BigReal minz = min.z;
3856  BigReal maxx = max.x;
3857  BigReal maxy = max.y;
3858  BigReal maxz = max.z;
3859 
3860  int xdev, ydev, zdev;
3861  int delnum = 0;
3862 
3863  FullAtomList::iterator atom_i = atom.begin();
3864  FullAtomList::iterator atom_e = atom.end();
3865 
3866  // DMK - Atom Separation (water vs. non-water)
3867  #if NAMD_SeparateWaters != 0
3868  FullAtomList::iterator atom_first = atom_i;
3869  int numLostWaterAtoms = 0;
3870  #endif
3871 
3872  while ( atom_i != atom_e ) {
3873  // Even though this code iterates through all atoms successively
3874  // it moves entire hydrogen/migration groups as follows:
3875  // Only the parent atom of the hydrogen/migration group has
3876  // nonzero migrationGroupSize. Values determined for xdev,ydev,zdev
3877  // will persist through the remaining group members so that each
3878  // following atom will again be added to the same mList.
3879  if ( atom_i->migrationGroupSize ) {
3880  Position pos = atom_i->position;
3881  if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
3882  // If there are multiple hydrogen groups in a migration group
3883  // (e.g. for supporting lone pairs)
3884  // the following code takes the average position (midpoint)
3885  // of their parents.
3886  int mgs = atom_i->migrationGroupSize;
3887  int c = 1;
3888  for ( int j=atom_i->hydrogenGroupSize; j<mgs;
3889  j+=(atom_i+j)->hydrogenGroupSize ) {
3890  pos += (atom_i+j)->position;
3891  ++c;
3892  }
3893  pos *= 1./c;
3894  // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
3895  }
3896 
3897  // Scaling the position below transforms space within patch from
3898  // what could have been a rotated parallelepiped into
3899  // orthogonal coordinates, where we can use minmax comparison
3900  // to detect which of our nearest neighbors this
3901  // parent atom might have entered.
3902  ScaledPosition s = lattice.scale(pos);
3903 
3904  // check if atom is within bounds
3905  if (s.x < minx) xdev = 0;
3906  else if (maxx <= s.x) xdev = 2;
3907  else xdev = 1;
3908 
3909  if (s.y < miny) ydev = 0;
3910  else if (maxy <= s.y) ydev = 2;
3911  else ydev = 1;
3912 
3913  if (s.z < minz) zdev = 0;
3914  else if (maxz <= s.z) zdev = 2;
3915  else zdev = 1;
3916 
3917  }
3918 
3919  if (mInfo[xdev][ydev][zdev]) { // process atom for migration
3920  // Don't migrate if destination is myself
3921 
3922  // See if we have a migration list already
3923  MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
3924  DebugM(3,"Migrating atom " << atom_i->id << " from patch "
3925  << patchID << " with position " << atom_i->position << "\n");
3926  mCur.add(*atom_i);
3927 
3928  ++delnum;
3929 
3930 
3931  // DMK - Atom Separation (water vs. non-water)
3932  #if NAMD_SeparateWaters != 0
3933  // Check to see if this atom is part of a water molecule. If
3934  // so, numWaterAtoms needs to be adjusted to refect the lost of
3935  // this atom.
3936  // NOTE: The atom separation code assumes that if the oxygen
3937  // atom of the hydrogen group making up the water molecule is
3938  // migrated to another HomePatch, the hydrogens will also
3939  // move!!!
3940  int atomIndex = atom_i - atom_first;
3941  if (atomIndex < numWaterAtoms)
3942  numLostWaterAtoms++;
3943  #endif
3944 
3945 
3946  } else {
3947  // By keeping track of delnum total being deleted from FullAtomList
3948  // the else clause allows us to fill holes as we visit each atom.
3949 
3950  if ( delnum ) { *(atom_i-delnum) = *atom_i; }
3951 
3952  }
3953 
3954  ++atom_i;
3955  }
3956 
3957  // DMK - Atom Separation (water vs. non-water)
3958  #if NAMD_SeparateWaters != 0
3959  numWaterAtoms -= numLostWaterAtoms;
3960  #endif
3961 
3962 
3963  int delpos = numAtoms - delnum;
3964  DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
3965  atom.del(delpos,delnum);
3966 
3967  numAtoms = atom.size();
3968 
3969  PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
3970 
3971  // signal depositMigration() that we are inMigration mode
3972  inMigration = true;
3973 
3974  // Drain the migration message buffer
3975  for (i=0; i<numMlBuf; i++) {
3976  DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
3978  }
3979  numMlBuf = 0;
3980 
3981  NAMD_EVENT_STOP(1, NamdProfileEvent::ATOM_MIGRATIONS);
3982 
3983  if (!allMigrationIn) {
3984  DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
3985  migrationSuspended = true;
3986  sequencer->suspend();
3987  migrationSuspended = false;
3988  }
3989  allMigrationIn = false;
3990 
3991  inMigration = false;
3992  marginViolations = 0;
3993 }
3994 
3995 void
3997 {
3998 
3999  if (!inMigration) { // We have to buffer changes due to migration
4000  // until our patch is in migration mode
4001  msgbuf[numMlBuf++] = msg;
4002  return;
4003  }
4004 
4005 
4006  // DMK - Atom Separation (water vs. non-water)
4007  #if NAMD_SeparateWaters != 0
4008 
4009 
4010  // Merge the incoming list of atoms with the current list of
4011  // atoms. Note that mergeSeparatedAtomList() will apply any
4012  // required transformations to the incoming atoms as it is
4013  // separating them.
4014  mergeAtomList(msg->migrationList);
4015 
4016 
4017  #else
4018 
4019  // Merge the incoming list of atoms with the current list of
4020  // atoms. Apply transformations to the incoming atoms as they are
4021  // added to this patch's list.
4022  {
4023  MigrationList &migrationList = msg->migrationList;
4025  Transform mother_transform;
4026  for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
4027  DebugM(1,"Migrating atom " << mi->id << " to patch "
4028  << patchID << " with position " << mi->position << "\n");
4029  if ( mi->migrationGroupSize ) {
4030  if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
4031  Position pos = mi->position;
4032  int mgs = mi->migrationGroupSize;
4033  int c = 1;
4034  for ( int j=mi->hydrogenGroupSize; j<mgs;
4035  j+=(mi+j)->hydrogenGroupSize ) {
4036  pos += (mi+j)->position;
4037  ++c;
4038  }
4039  pos *= 1./c;
4040  // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
4041  mother_transform = mi->transform;
4042  pos = lattice.nearest(pos,center,&mother_transform);
4043  mi->position = lattice.reverse_transform(mi->position,mi->transform);
4044  mi->position = lattice.apply_transform(mi->position,mother_transform);
4045  mi->transform = mother_transform;
4046  } else {
4047  mi->position = lattice.nearest(mi->position,center,&(mi->transform));
4048  mother_transform = mi->transform;
4049  }
4050  } else {
4051  mi->position = lattice.reverse_transform(mi->position,mi->transform);
4052  mi->position = lattice.apply_transform(mi->position,mother_transform);
4053  mi->transform = mother_transform;
4054  }
4055  atom.add(*mi);
4056  }
4057  }
4058 
4059 
4060  #endif // if (NAMD_SeparateWaters != 0)
4061 
4062 
4063  numAtoms = atom.size();
4064  delete msg;
4065 
4066  DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
4067  if (!--patchMigrationCounter) {
4068  DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
4069  allMigrationIn = true;
4070  patchMigrationCounter = numNeighbors;
4071  if (migrationSuspended) {
4072  DebugM(3,"patch "<<patchID<<" is being awakened\n");
4073  migrationSuspended = false;
4074  sequencer->awaken();
4075  return;
4076  }
4077  else {
4078  DebugM(3,"patch "<<patchID<<" was not suspended\n");
4079  }
4080  }
4081 }
4082 
4083 
4084 
4085 // DMK - Atom Separation (water vs. non-water)
4086 #if NAMD_SeparateWaters != 0
4087 
4088 // This function will separate waters from non-waters in the current
4089 // atom list (regardless of whether or not the atom list is has been
4090 // sorted yet or not).
4091 void HomePatch::separateAtoms() {
4092  SimParameters *simParams = Node::Object()->simParameters;
4093 
4094  // Basic Idea: Iterate through all the atoms in the current list
4095  // of atoms. Pack the waters in the current atoms list and move
4096  // the non-waters to the scratch list. Once the atoms have all
4097  // been separated, copy the non-waters to the end of the waters.
4098  // NOTE: This code does not assume that the atoms list has been
4099  // separated in any manner.
4100 
4101  // NOTE: Sanity check - Doesn't look like the default constructor actually
4102  // adds any atoms but does set numAtoms. ???
4103  if (atom.size() < 0) return; // Nothing to do.
4104 
4105  // Resize the scratch FullAtomList (tempAtom)
4106  tempAtom.resize(numAtoms); // NOTE: Worst case: all non-water
4107 
4108  // Define size of a water hydrogen group
4109  int wathgsize = 3;
4110  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
4111  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
4112 
4113  // Iterate through all the atoms
4114  int i = 0;
4115  int waterIndex = 0;
4116  int nonWaterIndex = 0;
4117  while (i < numAtoms) {
4118 
4119  FullAtom &atom_i = atom[i];
4120  Mass mass = atom_i.mass;
4121  int hgs = atom_i.hydrogenGroupSize;
4122  // Check to see if this hydrogen group is a water molecule
4123  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4124 
4125  // Move this hydrogen group up in the current atom list
4126  if (waterIndex != i) {
4127  atom[waterIndex ] = atom[i ]; // Oxygen
4128  atom[waterIndex + 1] = atom[i + 1]; // Hydrogen
4129  atom[waterIndex + 2] = atom[i + 2]; // Hydrogen
4130  if (wathgsize > 3) atom[waterIndex + 3] = atom[i + 3]; // lonepair
4131  if (wathgsize > 4) atom[waterIndex + 4] = atom[i + 4]; // drude
4132  // actual Drude water is arranged: O D LP H H
4133  }
4134 
4135  waterIndex += wathgsize;
4136  i += wathgsize;
4137 
4138  } else {
4139 
4140  // Move this hydrogen group into non-water (scratch) atom list
4141  for (int j = 0; j < hgs; j++)
4142  tempAtom[nonWaterIndex + j] = atom[i + j];
4143 
4144  nonWaterIndex += hgs;
4145  i += hgs;
4146  }
4147 
4148  } // end iterating through atoms
4149 
4150  // Iterate through the non-water (scratch) atom list, adding the
4151  // atoms to the end of the atom list.
4152  // NOTE: This could be done with a straight memcpy if the internal
4153  // data structures of ResizeArray could be accessed directly.
4154  // Or, perhaps add a member to ResizeArray that can add a consecutive
4155  // list of elements starting at a particular index (would be cleaner).
4156  for (i = 0; i < nonWaterIndex; i++)
4157  atom[waterIndex + i] = tempAtom[i];
4158 
4159  // Set numWaterAtoms
4160  numWaterAtoms = waterIndex;
4161 }
4162 
4163 
4164 // This function will merge the given list of atoms (not assumed to
4165 // be separated) with the current list of atoms (already assumed
4166 // to be separated).
4167 // NOTE: This function applies the transformations to the incoming
4168 // atoms as it is separating them.
4169 void HomePatch::mergeAtomList(FullAtomList &al) {
4170  SimParameters *simParams = Node::Object()->simParameters;
4171 
4172  // Sanity check
4173  if (al.size() <= 0) return; // Nothing to do
4174 
4175  const int orig_atomSize = atom.size();
4176  const int orig_alSize = al.size();
4177 
4178  // Resize the atom list (will eventually hold contents of both lists)
4179  atom.resize(orig_atomSize + orig_alSize); // NOTE: Will have contents of both
4180 
4181 
4182  #if 0 // version where non-waters are moved to scratch first
4183 
4184 
4185  // Basic Idea: The current list is separated already so copy the
4186  // non-water atoms out of it into the scratch atom array. Then
4187  // separate the incoming/given list (al), adding the waters to the
4188  // end of the waters in atom list and non-waters to the end of the
4189  // scratch list. At this point, all waters are in atom list and all
4190  // non-waters are in the scratch list so just copy the scratch list
4191  // to the end of the atom list.
4192  // NOTE: If al is already separated and the number of waters in it
4193  // is know, could simply move the non-waters in atoms back by that
4194  // amount and directly copy the waters in al into the created gap
4195  // and the non-waters in al to the end. Leave this as an
4196  // optimization for later since I'm not sure if this would actually
4197  // do better as the combining code (for combining migration
4198  // messages) would also have to merge the contents of the atom lists
4199  // they carry. Generally speaking, there is probably a faster way
4200  // to do this, but this will get it working.
4201 
4202  // Copy all the non-waters in the current atom list into the
4203  // scratch atom list.
4204  const int orig_atom_numNonWaters = orig_atomSize - numWaterAtoms;
4205  tempAtom.resize(orig_atom_numNonWaters + al.size()); // NOTE: Worst case
4206  for (int i = 0; i < orig_atom_numNonWaters; i++)
4207  tempAtom[i] = atom[numWaterAtoms + i];
4208 
4209  // Separate the contents of the given atom list (applying the
4210  // transforms as needed)
4211  int atom_waterIndex = numWaterAtoms;
4212  int atom_nonWaterIndex = orig_atom_numNonWaters;
4213  int i = 0;
4214  while (i < orig_alSize) {
4215 
4216  FullAtom &atom_i = al[i];
4217  int hgs = atom_i.hydrogenGroupSize;
4218  if ( hgs != atom_i.migrationGroupSize ) {
4219  NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
4220  }
4221  Mass mass = atom_i.mass;
4222 
4223  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4224 
4225  // Apply the transforms
4226 
4227  // Oxygen (@ +0)
4228  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4229  Transform mother_transform = al[i].transform;
4230 
4231  // Hydrogen (@ +1)
4232  al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
4233  al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
4234  al[i+1].transform = mother_transform;
4235 
4236  // Hydrogen (@ +2)
4237  al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
4238  al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
4239  al[i+2].transform = mother_transform;
4240 
4241  // Add to the end of the waters in the current list of atoms
4242  atom[atom_waterIndex ] = al[i ];
4243  atom[atom_waterIndex + 1] = al[i + 1];
4244  atom[atom_waterIndex + 2] = al[i + 2];
4245 
4246  atom_waterIndex += 3;
4247  i += 3;
4248 
4249  } else {
4250 
4251  // Apply the transforms
4252 
4253  // Non-Hydrogen (@ +0)
4254  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4255  Transform mother_transform = al[i].transform;
4256 
4257  // Hydrogens (@ +1 -> +(hgs-1))
4258  for (int j = 1; j < hgs; j++) {
4259  al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
4260  al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
4261  al[i+j].transform = mother_transform;
4262  }
4263 
4264  // Add to the end of the non-waters (scratch) atom list
4265  for (int j = 0; j < hgs; j++)
4266  tempAtom[atom_nonWaterIndex + j] = al[i + j];
4267 
4268  atom_nonWaterIndex += hgs;
4269  i += hgs;
4270  }
4271 
4272  } // end while iterating through given atom list
4273 
4274  // Copy all the non-waters to the end of the current atom list
4275  for (int i = 0; i < atom_nonWaterIndex; i++)
4276  atom[atom_waterIndex + i] = tempAtom[i];
4277 
4278  // Set numWaterAtoms and numAtoms
4279  numWaterAtoms = atom_waterIndex;
4280  numAtoms = atom.size();
4281 
4282 
4283  #else
4284 
4285 
4286  // Basic Idea: Count the number of water atoms in the incoming atom
4287  // list then move the non-waters back in the current atom list to
4288  // make room for the incoming waters. Once there is room in the
4289  // current list, separate the incoming list as the atoms are being
4290  // added to the current list.
4291  // NOTE: Since the incoming atom list is likely to be small,
4292  // iterating over its hydrogen groups twice should not be too bad.
4293  // NOTE: This code assumes the current list is already separated,
4294  // the incoming list may not be separated, and the transforms are
4295  // applied to the incoming atoms as the separation occurs.
4296 
4297  // size of a water hydrogen group
4298  int wathgsize = 3;
4299  if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
4300  else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
4301 
4302  // Count the number of waters in the given atom list
4303  int al_numWaterAtoms = 0;
4304  int i = 0;
4305  while (i < orig_alSize) {
4306 
4307  FullAtom &atom_i = al[i];
4308  int hgs = atom_i.hydrogenGroupSize;
4309  Mass mass = atom_i.mass;
4310 
4311  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4312  al_numWaterAtoms += wathgsize;
4313  }
4314 
4315  i += hgs;
4316  }
4317 
4318  // Move all of the non-waters in the current atom list back (to a
4319  // higher index) by the number of waters in the given list.
4320  if (al_numWaterAtoms > 0) {
4321  for (i = orig_atomSize - 1; i >= numWaterAtoms; i--) {
4322  atom[i + al_numWaterAtoms] = atom[i];
4323  }
4324  }
4325 
4326  // Separte the atoms in the given atom list. Perform the
4327  // transformations on them and then add them to the appropriate
4328  // location in the current atom list.
4329  int atom_waterIndex = numWaterAtoms;
4330  int atom_nonWaterIndex = orig_atomSize + al_numWaterAtoms;
4331  i = 0;
4332  while (i < orig_alSize) {
4333 
4334  FullAtom &atom_i = al[i];
4335  int hgs = atom_i.hydrogenGroupSize;
4336  if ( hgs != atom_i.migrationGroupSize ) {
4337  NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
4338  }
4339  Mass mass = atom_i.mass;
4340 
4341  if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
4342 
4343  // Apply the transforms
4344 
4345  // Oxygen (@ +0)
4346  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4347  Transform mother_transform = al[i].transform;
4348 
4349  // Hydrogen (@ +1)
4350  al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
4351  al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
4352  al[i+1].transform = mother_transform;
4353 
4354  // Hydrogen (@ +2)
4355  al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
4356  al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
4357  al[i+2].transform = mother_transform;
4358 
4359  // Add to the end of the waters in the current list of atoms
4360  atom[atom_waterIndex ] = al[i ];
4361  atom[atom_waterIndex + 1] = al[i + 1];
4362  atom[atom_waterIndex + 2] = al[i + 2];
4363 
4364  if (wathgsize > 3) atom[atom_waterIndex + 3] = al[i + 3];
4365 
4366  atom_waterIndex += wathgsize;
4367  i += wathgsize;
4368 
4369  } else {
4370 
4371  // Apply the transforms
4372 
4373  // Non-Hydrogen (@ +0)
4374  al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
4375  Transform mother_transform = al[i].transform;
4376 
4377  // Hydrogens (@ +1 -> +(hgs-1))
4378  for (int j = 1; j < hgs; j++) {
4379  al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
4380  al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
4381  al[i+j].transform = mother_transform;
4382  }
4383 
4384  // Add to the end of the non-waters (scratch) atom list
4385  for (int j = 0; j < hgs; j++)
4386  atom[atom_nonWaterIndex + j] = al[i + j];
4387 
4388  atom_nonWaterIndex += hgs;
4389  i += hgs;
4390  }
4391 
4392  } // end while iterating through given atom list
4393 
4394  // Set numWaterAtoms and numAtoms
4395  numWaterAtoms = atom_waterIndex;
4396  numAtoms = atom_nonWaterIndex;
4397 
4398  #endif
4399 }
4400 
4401 #endif
4402 
4403 
4404 
4405 inline void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx,
4406  HGArrayBigReal &b)
4407 {
4408  int i,ii=-1,ip,j;
4409  double sum;
4410 
4411  for (i=0;i<n;i++) {
4412  ip=indx[i];
4413  sum=b[ip];
4414  b[ip]=b[i];
4415  if (ii >= 0)
4416  for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
4417  else if (sum) ii=i;
4418  b[i]=sum;
4419  }
4420  for (i=n-1;i>=0;i--) {
4421  sum=b[i];
4422  for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
4423  b[i]=sum/a[i][i];
4424  }
4425 }
4426 
4427 
4428 inline void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
4429 {
4430 
4431  int i,imax,j,k;
4432  double big,dum,sum,temp;
4433  HGArrayBigReal vv;
4434  *d=1.0;
4435  for (i=0;i<n;i++) {
4436  big=0.0;
4437  for (j=0;j<n;j++)
4438  if ((temp=fabs(a[i][j])) > big) big=temp;
4439  if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
4440  vv[i]=1.0/big;
4441  }
4442  for (j=0;j<n;j++) {
4443  for (i=0;i<j;i++) {
4444  sum=a[i][j];
4445  for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
4446  a[i][j]=sum;
4447  }
4448  big=0.0;
4449  for (i=j;i<n;i++) {
4450  sum=a[i][j];
4451  for (k=0;k<j;k++)
4452  sum -= a[i][k]*a[k][j];
4453  a[i][j]=sum;
4454  if ( (dum=vv[i]*fabs(sum)) >= big) {
4455  big=dum;
4456  imax=i;
4457  }
4458  }
4459  if (j != imax) {
4460  for (k=0;k<n;k++) {
4461  dum=a[imax][k];
4462  a[imax][k]=a[j][k];
4463  a[j][k]=dum;
4464  }
4465  *d = -(*d);
4466  vv[imax]=vv[j];
4467  }
4468  indx[j]=imax;
4469  if (a[j][j] == 0.0) a[j][j]=TINY;
4470  if (j != n-1) {
4471  dum=1.0/(a[j][j]);
4472  for (i=j+1;i<n;i++) a[i][j] *= dum;
4473  }
4474  }
4475 }
4476 
4477 
4478 inline void G_q(const HGArrayVector &refab, HGMatrixVector &gqij,
4479  const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl) {
4480  int i;
4481  // step through the rows of the matrix
4482  for(i=0;i<m;i++) {
4483  gqij[i][ial[i]]=2.0*refab[i];
4484  gqij[i][ibl[i]]=-gqij[i][ial[i]];
4485  }
4486 };
4487 
4488 
4489 // c-ji code for MOLLY 7-31-99
4490 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) {
4491  // input: n = length of hydrogen group to be averaged (shaked)
4492  // q[n] = vector array of original positions
4493  // m = number of constraints
4494  // imass[n] = inverse mass for each atom
4495  // length2[m] = square of reference bond length for each constraint
4496  // ial[m] = atom a in each constraint
4497  // ibl[m] = atom b in each constraint
4498  // refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
4499  // tolf = function error tolerance for Newton's iteration
4500  // ntrial = max number of Newton's iterations
4501  // output: lambda[m] = double array of lagrange multipliers (used by mollify)
4502  // qtilde[n] = vector array of averaged (shaked) positions
4503 
4504  int k,k1,i,j;
4505  BigReal errx,errf,d,tolx;
4506 
4507  HGArrayInt indx;
4508  HGArrayBigReal p;
4509  HGArrayBigReal fvec;
4510  HGMatrixBigReal fjac;
4511  HGArrayVector avgab;
4512  HGMatrixVector grhs;
4513  HGMatrixVector auxrhs;
4514  HGMatrixVector glhs;
4515 
4516  // iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
4517  tolx=tolf;
4518 
4519  // initialize lambda, globalGrhs
4520 
4521  for (i=0;i<m;i++) {
4522  lambda[i]=0.0;
4523  }
4524 
4525  // define grhs, auxrhs for all iterations
4526  // grhs= g_x(q)
4527  //
4528  G_q(refab,grhs,n,m,ial,ibl);
4529  for (k=1;k<=ntrial;k++) {
4530  // usrfun(qtilde,q0,lambda,fvec,fjac,n,water);
4531  HGArrayBigReal gij;
4532  // this used to be main loop of usrfun
4533  // compute qtilde given q0, lambda, IMASSes
4534  {
4535  BigReal multiplier;
4536  HGArrayVector tmp;
4537  for (i=0;i<m;i++) {
4538  multiplier = lambda[i];
4539  // auxrhs = M^{-1}grhs^{T}
4540  for (j=0;j<n;j++) {
4541  auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
4542  }
4543  }
4544  for (j=0;j<n;j++) {
4545  // tmp[j]=0.0;
4546  for (i=0;i<m;i++) {
4547  tmp[j]+=auxrhs[i][j];
4548  }
4549  }
4550 
4551  for (j=0;j<n;j++) {
4552  qtilde[j].position=q[j]+tmp[j];
4553  }
4554  // delete [] tmp;
4555  }
4556 
4557  for ( i = 0; i < m; i++ ) {
4558  avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
4559  }
4560 
4561  // iout<<iINFO << "Calling Jac" << std::endl<<endi;
4562  // Jac(qtilde, q0, fjac,n,water);
4563  {
4564  // Vector glhs[3*n+3];
4565 
4566  HGMatrixVector grhs2;
4567 
4568  G_q(avgab,glhs,n,m,ial,ibl);
4569 #ifdef DEBUG0
4570  iout<<iINFO << "G_q:" << std::endl<<endi;
4571  for (i=0;i<m;i++) {
4572  iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
4573  }
4574 #endif
4575  // G_q(refab,grhs2,m,ial,ibl);
4576  // update with the masses
4577  for (j=0; j<n; j++) { // number of atoms
4578  for (i=0; i<m;i++) { // number of constraints
4579  grhs2[i][j] = grhs[i][j]*imass[j];
4580  }
4581  }
4582 
4583  // G_q(qtilde) * M^-1 G_q'(q0) =
4584  // G_q(qtilde) * grhs'
4585  for (i=0;i<m;i++) { // number of constraints
4586  for (j=0;j<m;j++) { // number of constraints
4587  fjac[i][j] = 0;
4588  for (k1=0;k1<n;k1++) {
4589  fjac[i][j] += glhs[i][k1]*grhs2[j][k1];
4590  }
4591  }
4592  }
4593 #ifdef DEBUG0
4594  iout<<iINFO << "glhs" <<endi;
4595  for(i=0;i<9;i++) {
4596  iout<<iINFO << glhs[i] << ","<<endi;
4597  }
4598  iout<<iINFO << std::endl<<endi;
4599  for(i=0;i<9;i++) {
4600  iout<<iINFO << grhs2[i] << ","<<endi;
4601  }
4602  iout<<iINFO << std::endl<<endi;
4603 #endif
4604  // delete[] grhs2;
4605  }
4606  // end of Jac calculation
4607 #ifdef DEBUG0
4608  iout<<iINFO << "Jac" << std::endl<<endi;
4609  for (i=0;i<m;i++)
4610  for (j=0;j<m;j++)
4611  iout<<iINFO << fjac[i][j] << " "<<endi;
4612  iout<< std::endl<<endi;
4613 #endif
4614  // calculate constraints in gij for n constraints this being a water
4615  // G(qtilde, gij, n, water);
4616  for (i=0;i<m;i++) {
4617  gij[i]=avgab[i]*avgab[i]-length2[i];
4618  }
4619 #ifdef DEBUG0
4620  iout<<iINFO << "G" << std::endl<<endi;
4621  iout<<iINFO << "( "<<endi;
4622  for(i=0;i<m-1;i++) {
4623  iout<<iINFO << gij[i] << ", "<<endi;
4624  }
4625  iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
4626 #endif
4627  // fill the return vector
4628  for(i=0;i<m;i++) {
4629  fvec[i] = gij[i];
4630  }
4631  // free up the constraints
4632  // delete[] gij;
4633  // continue Newton's iteration
4634  errf=0.0;
4635  for (i=0;i<m;i++) errf += fabs(fvec[i]);
4636 #ifdef DEBUG0
4637  iout<<iINFO << "errf: " << errf << std::endl<<endi;
4638 #endif
4639  if (errf <= tolf) {
4640  break;
4641  }
4642  for (i=0;i<m;i++) p[i] = -fvec[i];
4643  // iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
4644  ludcmp(fjac,m,indx,&d);
4645  lubksb(fjac,m,indx,p);
4646 
4647  errx=0.0;
4648  for (i=0;i<m;i++) {
4649  errx += fabs(p[i]);
4650  }
4651  for (i=0;i<m;i++)
4652  lambda[i] += p[i];
4653 
4654 #ifdef DEBUG0
4655  iout<<iINFO << "lambda:" << lambda[0]
4656  << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
4657  iout<<iINFO << "errx: " << errx << std::endl<<endi;
4658 #endif
4659  if (errx <= tolx) break;
4660 #ifdef DEBUG0
4661  iout<<iINFO << "Qtilde:" << std::endl<<endi;
4662  iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi;
4663 #endif
4664  }
4665 #ifdef DEBUG
4666  iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
4667 #endif
4668 
4669  return k; //
4670 }
4671 
4672 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) {
4673  int i,j,k;
4674  BigReal d;
4675  HGMatrixBigReal fjac;
4676  Vector zero(0.0,0.0,0.0);
4677 
4678  HGArrayVector tmpforce;
4679  HGArrayVector tmpforce2;
4680  HGArrayVector y;
4681  HGMatrixVector grhs;
4682  HGMatrixVector glhs;
4683  HGArrayBigReal aux;
4684  HGArrayInt indx;
4685 
4686  for(i=0;i<n;i++) {
4687  tmpforce[i]=imass[i]*force[i];
4688  }
4689 
4690  HGMatrixVector grhs2;
4691  HGArrayVector avgab;
4692 
4693  for ( i = 0; i < m; i++ ) {
4694  avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
4695  }
4696 
4697  G_q(avgab,glhs,n,m,ial,ibl);
4698  G_q(refab,grhs,n,m,ial,ibl);
4699  // update with the masses
4700  for (j=0; j<n; j++) { // number of atoms
4701  for (i=0; i<m;i++) { // number of constraints
4702  grhs2[i][j] = grhs[i][j]*imass[j];
4703  }
4704  }
4705 
4706  // G_q(qtilde) * M^-1 G_q'(q0) =
4707  // G_q(qtilde) * grhs'
4708  for (i=0;i<m;i++) { // number of constraints
4709  for (j=0;j<m;j++) { // number of constraints
4710  fjac[j][i] = 0;
4711  for (k=0;k<n;k++) {
4712  fjac[j][i] += glhs[i][k]*grhs2[j][k];
4713  }
4714  }
4715  }
4716 
4717  // aux=gqij*tmpforce
4718  // globalGrhs::computeGlobalGrhs(q0,n,water);
4719  // G_q(refab,grhs,m,ial,ibl);
4720  for(i=0;i<m;i++) {
4721  aux[i]=0.0;
4722  for(j=0;j<n;j++) {
4723  aux[i]+=grhs[i][j]*tmpforce[j];
4724  }
4725  }
4726 
4727  ludcmp(fjac,m,indx,&d);
4728  lubksb(fjac,m,indx,aux);
4729 
4730  for(j=0;j<n;j++) {
4731  y[j] = zero;
4732  for(i=0;i<m;i++) {
4733  y[j] += aux[i]*glhs[i][j];
4734  }
4735  }
4736  for(i=0;i<n;i++) {
4737  y[i]=force[i]-y[i];
4738  }
4739 
4740  // gqq12*y
4741  for(i=0;i<n;i++) {
4742  tmpforce2[i]=imass[i]*y[i];
4743  }
4744 
4745  // here we assume that tmpforce is initialized to zero.
4746  for (i=0;i<n;i++) {
4747  tmpforce[i]=zero;
4748  }
4749 
4750  for (j=0;j<m;j++) {
4751  Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
4752  tmpforce[ial[j]] += tmpf;
4753  tmpforce[ibl[j]] -= tmpf;
4754  }
4755  // c-ji the other bug for 2 constraint water was this line (2-4-99)
4756  // for(i=0;i<m;i++) {
4757  for(i=0;i<n;i++) {
4758  force[i]=tmpforce[i]+y[i];
4759  }
4760 
4761 }
static Node * Object()
Definition: Node.h:86
void depositMigration(MigrateAtomsMsg *)
Definition: HomePatch.C:3996
void recvCheckpointLoad(CheckpointAtomsMsg *msg)
Definition: HomePatch.C:3515
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:2000
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:2962
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:41
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:2381
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:1932
void rattle2(const BigReal, Tensor *virial)
Definition: HomePatch.C:2828
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:3174
#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:402
#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:3719
BigReal rigidTol
BigReal HGMatrixBigReal[MAXHGS][MAXHGS]
Definition: HomePatch.C:65
CompAtom * avgPositionList
Definition: ProxyMgr.h:104
void loweAndersenVelocities()
Definition: HomePatch.C:3097
#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:119
void recvCheckpointReq(int task, const char *key, int replica, int pe)
Definition: HomePatch.C:3485
Vector origin() const
Definition: Lattice.h:262
CudaAtom * cudaAtomPtr
Definition: Patch.h:205
BigReal alchLambda
void exchangeCheckpoint(int scriptTask, int &bpc)
Definition: HomePatch.C:3477
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:3238
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:3568
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:3459
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
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
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t cudaTextureObject_t const int const float const PatchPairRecord *__restrict__ const int *__restrict__ atomIndex
void NAMD_bug(const char *err_msg)
Definition: common.C:123
void gbisP2Ready()
Definition: Patch.C:569
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:4672
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:3619
CompAtomList p_avg
Definition: Patch.h:147
void mollyMollify(Tensor *virial)
Definition: HomePatch.C:3386
Bool staticAtomAssignment
int numFepInitial
Definition: Molecule.h:607
int rattle1old(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2558
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:2229
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:1330
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:3146
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:4405
static LdbCoordinator * Object()
void setLcpoType()
Definition: HomePatch.C:3120
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:3628
int berendsenPressure_count
Definition: Sequencer.h:104
static AtomMap * Object()
Definition: AtomMap.h:36
void gbisP3Ready()
Definition: Patch.C:586
void doAtomMigration()
Definition: HomePatch.C:3840
MigrateAtomsMsg * msgbuf[PatchMap::MaxOneAway]
Definition: HomePatch.h:494
int hardWallDrude(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2023
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:1335
void setGBISIntrinsicRadii()
Definition: HomePatch.C:3131
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:3218
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:3449
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:4428
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
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:4478
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:3312
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:3594
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:3261
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:1961
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:3112
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:2371
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:3763
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:3554
Vector unit(void) const
Definition: Vector.h:182
Vector c() const
Definition: Lattice.h:254
void recvExchangeMsg(ExchangeAtomsMsg *msg)
Definition: HomePatch.C:3605
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:4490
#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:3573