Update default behavior and documentation for Drude
[namd.git] / src / HomePatch.C
1
2 /**
3 ***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
4 ***  The Board of Trustees of the University of Illinois.
5 ***  All rights reserved.
6 **/
7
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 //#define PRINT_COMP
49 #define TINY 1.0e-20;
50 #define MAXHGS 10
51 #define MIN_DEBUG_LEVEL 2
52 //#define DEBUGM
53 #include "Debug.h"
54
55 #include <vector>
56 #include <algorithm>
57 using namespace std;
58
59 typedef int HGArrayInt[MAXHGS];
60 typedef BigReal HGArrayBigReal[MAXHGS];
61 typedef zVector HGArrayVector[MAXHGS];
62 typedef BigReal HGMatrixBigReal[MAXHGS][MAXHGS];
63 typedef zVector HGMatrixVector[MAXHGS][MAXHGS];
64
65 #include "ComputeNonbondedMICKernel.h"
66
67 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);
68
69 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);
70
71 #define MASS_EPSILON  (1.0e-35)  //a very small floating point number
72
73
74 // DMK - Atom Separation (water vs. non-water)
75 #if NAMD_SeparateWaters != 0
76
77 // Macro to test if a hydrogen group represents a water molecule.
78 // NOTE: This test is the same test in Molecule.C for setting the
79 //   OxygenAtom flag in status.
80 // hgtype should be the number of atoms in a water hydrogen group
81 // It must now be set based on simulation parameters because we might
82 // be using tip4p
83
84 // DJH: This will give false positive for full Drude model,
85 //      e.g. O D H is not water but has hgs==3
86 #define IS_HYDROGEN_GROUP_WATER(hgs, mass)                 \
87   ((hgs >= 3) && ((mass >= 14.0) && (mass <= 18.0)))
88
89 #endif
90
91
92 HomePatch::HomePatch(PatchID pd, FullAtomList &al) : Patch(pd)
93 // DMK - Atom Separation (water vs. non-water)
94 #if NAMD_SeparateWaters != 0
95   ,tempAtom()
96 #endif
97
98   atom.swap(al);
99   settle_initialized = 0;
100
101   doAtomUpdate = true;
102   rattleListValid = false;
103
104   exchange_msg = 0;
105   exchange_req = -1;
106
107   numGBISP1Arrived = 0;
108   numGBISP2Arrived = 0;
109   numGBISP3Arrived = 0;
110   phase1BoxClosedCalled = false;
111   phase2BoxClosedCalled = false;
112   phase3BoxClosedCalled = false;
113
114   min.x = PatchMap::Object()->min_a(patchID);
115   min.y = PatchMap::Object()->min_b(patchID);
116   min.z = PatchMap::Object()->min_c(patchID);
117   max.x = PatchMap::Object()->max_a(patchID);
118   max.y = PatchMap::Object()->max_b(patchID);
119   max.z = PatchMap::Object()->max_c(patchID);
120   center = 0.5*(min+max);
121
122   int aAway = PatchMap::Object()->numaway_a();
123   if ( PatchMap::Object()->periodic_a() ||
124        PatchMap::Object()->gridsize_a() > aAway + 1 ) {
125     aAwayDist = (max.x - min.x) * aAway;
126   } else {
127     aAwayDist = Node::Object()->simParameters->patchDimension;
128   }
129   int bAway = PatchMap::Object()->numaway_b();
130   if ( PatchMap::Object()->periodic_b() ||
131        PatchMap::Object()->gridsize_b() > bAway + 1 ) {
132     bAwayDist = (max.y - min.y) * bAway;
133   } else {
134     bAwayDist = Node::Object()->simParameters->patchDimension;
135   }
136   int cAway = PatchMap::Object()->numaway_c();
137   if ( PatchMap::Object()->periodic_c() ||
138        PatchMap::Object()->gridsize_c() > cAway + 1 ) {
139     cAwayDist = (max.z - min.z) * cAway;
140   } else {
141     cAwayDist = Node::Object()->simParameters->patchDimension;
142   }
143
144   migrationSuspended = false;
145   allMigrationIn = false;
146   marginViolations = 0;
147   patchMapRead = 0; // We delay read of PatchMap data
148                     // to make sure it is really valid
149   inMigration = false;
150   numMlBuf = 0;
151   flags.sequence = -1;
152   flags.maxForceUsed = -1;
153
154   numAtoms = atom.size();
155   replacementForces = 0;
156
157   SimParameters *simParams = Node::Object()->simParameters;
158   doPairlistCheck_newTolerance = 
159         0.5 * ( simParams->pairlistDist - simParams->cutoff );
160
161
162   numFixedAtoms = 0;
163   if ( simParams->fixedAtomsOn ) {
164     for ( int i = 0; i < numAtoms; ++i ) {
165       numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
166     }
167   }
168
169 #ifdef NODEAWARE_PROXY_SPANNINGTREE
170   ptnTree.resize(0);
171   /*children = NULL;
172   numChild = 0;*/
173 #else
174   child =  new int[proxySpanDim];
175   nChild = 0;   // number of proxy spanning tree children
176 #endif
177
178 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
179   nphs = 0;
180   localphs = NULL;
181   isProxyChanged = 0;
182 #endif
183
184
185   // DMK - Atom Separation (water vs. non-water)
186   #if NAMD_SeparateWaters != 0
187
188     // Create the scratch memory for separating atoms
189     tempAtom.resize(numAtoms);
190     numWaterAtoms = -1;
191
192     // Separate the current list of atoms
193     separateAtoms();
194
195   #endif
196     
197   // Handle unusual water models here
198   if (simParams->watmodel == WAT_TIP4) init_tip4();
199   else if (simParams->watmodel == WAT_SWM4) init_swm4();
200
201   isNewProxyAdded = 0;
202 }
203
204 void HomePatch::write_tip4_props() {
205   printf("Writing r_om and r_ohc: %f | %f\n", r_om, r_ohc);
206 }
207
208 void HomePatch::init_tip4() {
209   // initialize the distances needed for the tip4p water model
210   Molecule *mol = Node::Object()->molecule;
211   r_om = mol->r_om;
212   r_ohc = mol->r_ohc;
213 }
214
215
216 void ::HomePatch::init_swm4() {
217   // initialize the distances needed for the SWM4 water model
218   Molecule *mol = Node::Object()->molecule;
219   r_om = mol->r_om;
220   r_ohc = mol->r_ohc;
221 }
222
223
224 void HomePatch::reinitAtoms(FullAtomList &al) {
225   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
226
227   atom.swap(al);
228   numAtoms = atom.size();
229
230   doAtomUpdate = true;
231   rattleListValid = false;
232
233   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
234
235   // DMK - Atom Separation (water vs. non-water)
236   #if NAMD_SeparateWaters != 0
237
238     // Reset the numWaterAtoms value
239     numWaterAtoms = -1;
240
241     // Separate the atoms
242     separateAtoms();
243
244   #endif
245 }
246
247 // Bind a Sequencer to this HomePatch
248 void HomePatch::useSequencer(Sequencer *sequencerPtr)
249 { sequencer=sequencerPtr; }
250
251 // start simulation over this Patch of atoms
252 void HomePatch::runSequencer(void)
253 { sequencer->run(); }
254
255 void HomePatch::readPatchMap() {
256   // iout << "Patch " << patchID << " has " << proxy.size() << " proxies.\n" << endi;
257   PatchMap *p = PatchMap::Object();
258   PatchID nnPatchID[PatchMap::MaxOneAway];
259
260   patchMigrationCounter = numNeighbors 
261     = PatchMap::Object()->oneAwayNeighbors(patchID, nnPatchID);
262   DebugM( 1, "NumNeighbors for pid " <<patchID<<" is "<< numNeighbors << "\n");
263   int n;
264   for (n=0; n<numNeighbors; n++) {
265     realInfo[n].destNodeID = p->node(realInfo[n].destPatchID = nnPatchID[n]);
266      DebugM( 1, " nnPatchID=" <<nnPatchID[n]<<" nnNodeID="<< realInfo[n].destNodeID<< "\n");
267     realInfo[n].mList.resize(0);
268   }
269
270   // Make mapping from the 3x3x3 cube of pointers to real migration info
271   for (int i=0; i<3; i++)
272     for (int j=0; j<3; j++)
273       for (int k=0; k<3; k++)
274       {
275         int pid =  p->pid(p->index_a(patchID)+i-1, 
276             p->index_b(patchID)+j-1, p->index_c(patchID)+k-1);
277         if (pid < 0) {
278            DebugM(5, "ERROR, for patchID " << patchID <<" I got neigh pid = " << pid << "\n");
279         }
280         if (pid == patchID && ! (
281                 ( (i-1) && p->periodic_a() ) ||
282                 ( (j-1) && p->periodic_b() ) ||
283                 ( (k-1) && p->periodic_c() ) )) {
284           mInfo[i][j][k] = NULL;
285         }
286         else {
287           // Does not work as expected for periodic with only two patches.
288           // Also need to check which image we want, but OK for now.  -JCP
289           for (n = 0; n<numNeighbors; n++) {
290             if (pid == realInfo[n].destPatchID) {
291               mInfo[i][j][k] = &realInfo[n];
292               break;
293             }
294           }
295           if (n == numNeighbors) { // disaster! 
296             DebugM(4,"BAD News, I could not find PID " << pid << "\n");
297           }
298         }
299       }
300
301   DebugM(1,"Patch("<<patchID<<") # of neighbors = " << numNeighbors << "\n");
302 }
303
304 HomePatch::~HomePatch()
305 {
306     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
307 #ifdef NODEAWARE_PROXY_SPANNINGTREE
308     ptnTree.resize(0);
309     #ifdef USE_NODEPATCHMGR
310     delete [] nodeChildren;    
311     #endif
312 #endif
313   delete [] child;
314 }
315
316
317 void HomePatch::boxClosed(int box) {
318   // begin gbis
319   if (box == 5) {// end of phase 1
320     phase1BoxClosedCalled = true;
321     if (!psiSumBox.isOpen() && numGBISP1Arrived == proxy.size()) {
322       if (flags.doGBIS && flags.doNonbonded) {
323         sequencer->awaken();
324       }
325     } else {
326       //need to wait until proxies arrive before awakening
327     }
328   } else if (box == 6) {// intRad
329     //do nothing
330   } else if (box == 7) {// bornRad
331     //do nothing
332   } else if (box == 8) {// end of phase 2
333     phase2BoxClosedCalled = true;
334     //if no proxies, AfterP1 can't be called from receive
335     //so it will be called from here
336     if (!dEdaSumBox.isOpen() && numGBISP2Arrived == proxy.size()) {
337       if (flags.doGBIS && flags.doNonbonded) {
338         sequencer->awaken();
339       }
340     } else {
341       //need to wait until proxies arrive before awakening
342     }
343   } else if (box == 9) {
344     //do nothing
345   } else if (box == 10) {
346     //lcpoType Box closed: do nothing
347   } else {
348     //do nothing
349   }
350   // end gbis
351
352   if ( ! --boxesOpen ) {
353     if ( replacementForces ) {
354       for ( int i = 0; i < numAtoms; ++i ) {
355         if ( replacementForces[i].replace ) {
356           for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
357           f[Results::normal][i] = replacementForces[i].force;
358         }
359       }
360       replacementForces = 0;
361     }
362     DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
363         << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
364     // only awaken suspended threads.  Then say it is suspended.
365
366     phase3BoxClosedCalled = true;
367     if (flags.doGBIS) {
368       if (flags.doNonbonded) {
369         sequencer->awaken();
370       } else {
371         if (numGBISP1Arrived == proxy.size() &&
372           numGBISP2Arrived == proxy.size() &&
373           numGBISP3Arrived == proxy.size()) {
374           sequencer->awaken();//all boxes closed and all proxies arrived
375         }
376       }
377     } else {//non-gbis awaken
378       sequencer->awaken();
379     }
380   } else {
381     DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
382   }
383 }
384
385 void HomePatch::registerProxy(RegisterProxyMsg *msg) {
386   DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
387   proxy.add(msg->node);
388   forceBox.clientAdd();
389
390   isNewProxyAdded = 1;
391 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
392   isProxyChanged = 1;
393 #endif
394
395   Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
396   delete msg;
397 }
398
399 void HomePatch::unregisterProxy(UnregisterProxyMsg *msg) {
400 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
401   isProxyChanged = 1;
402 #endif
403   int n = msg->node;
404   NodeID *pe = proxy.begin();
405   for ( ; *pe != n; ++pe );
406   forceBox.clientRemove();
407   proxy.del(pe - proxy.begin());
408   delete msg;
409 }
410
411 #if USE_TOPOMAP && USE_SPANNING_TREE
412
413 int HomePatch::findSubroots(int dim, int* subroots, int psize, int* pidscopy){
414   int nChild = 0;
415   int cones[6][proxySpanDim*proxySpanDim+proxySpanDim];
416   int conesizes[6] = {0,0,0,0,0,0};
417   int conecounters[6] = {0,0,0,0,0,0};
418   int childcounter = 0;
419   nChild = (psize>proxySpanDim)?proxySpanDim:psize;
420   TopoManager tmgr;
421   for(int i=0;i<psize;i++){
422     int cone = tmgr.getConeNumberForRank(pidscopy[i]);
423     cones[cone][conesizes[cone]++] = pidscopy[i];
424   }
425
426   while(childcounter<nChild){
427     for(int i=0;i<6;i++){
428       if(conecounters[i]<conesizes[i]){
429         subroots[childcounter++] = cones[i][conecounters[i]++];
430       }
431     }
432   }
433   for(int i=nChild;i<proxySpanDim;i++)
434     subroots[i] = -1;
435   return nChild;
436 }
437 #endif // USE_TOPOMAP 
438
439 static int compDistance(const void *a, const void *b)
440 {
441   int d1 = abs(*(int *)a - CkMyPe());
442   int d2 = abs(*(int *)b - CkMyPe());
443   if (d1 < d2) 
444     return -1;
445   else if (d1 == d2) 
446     return 0;
447   else 
448     return 1;
449 }
450
451 void HomePatch::sendProxies()
452 {
453 #if USE_NODEPATCHMGR
454         CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
455         NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
456         npm->sendProxyList(patchID, proxy.begin(), proxy.size());
457 #else
458         ProxyMgr::Object()->sendProxies(patchID, proxy.begin(), proxy.size());
459 #endif
460 }
461
462 #ifdef NODEAWARE_PROXY_SPANNINGTREE
463 void HomePatch::buildNodeAwareSpanningTree(void){
464 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
465         DebugFileTrace *dft = DebugFileTrace::Object();
466         dft->openTrace();
467         dft->writeTrace("HomePatch[%d] has %d proxy on proc[%d] node[%d]\n", patchID, proxy.size(), CkMyPe(), CkMyNode());
468         dft->writeTrace("Proxies are: ");
469         for(int i=0; i<proxy.size(); i++) dft->writeTrace("%d(%d), ", proxy[i], CkNodeOf(proxy[i]));
470         dft->writeTrace("\n");
471         dft->closeTrace();
472 #endif
473  
474     //build the naive spanning tree for this home patch
475     if(! proxy.size()) {
476         //this case will not happen in practice.
477         //In debugging state where spanning tree is enforced, then this could happen
478         //Chao Mei        
479        return;
480     }
481     ProxyMgr::buildSinglePatchNodeAwareSpanningTree(patchID, proxy, ptnTree);
482     //optimize on the naive spanning tree
483
484     //setup the children
485     setupChildrenFromProxySpanningTree();
486     //send down to children
487     sendNodeAwareSpanningTree();
488 }
489
490 void HomePatch::setupChildrenFromProxySpanningTree(){
491 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
492     isProxyChanged = 1;
493 #endif
494     if(ptnTree.size()==0) {
495         nChild = 0;
496         delete [] child;
497         child = NULL;
498         #ifdef USE_NODEPATCHMGR
499         numNodeChild = 0;
500         delete [] nodeChildren;
501         nodeChildren = NULL;        
502         #endif
503         return;
504     }
505     proxyTreeNode *rootnode = &ptnTree.item(0);
506     CmiAssert(rootnode->peIDs[0] == CkMyPe());
507     //set up children
508     //1. add external children (the first proc inside the proxy tree node)    
509     //2. add internal children (with threshold that would enable spanning    
510     int internalChild = rootnode->numPes-1;
511     int externalChild = ptnTree.size()-1;
512     externalChild = (proxySpanDim>externalChild)?externalChild:proxySpanDim;
513     int internalSlots = proxySpanDim-externalChild;
514     if(internalChild>0){
515       if(internalSlots==0) {
516          //at least having one internal child
517         internalChild = 1;
518       }else{
519         internalChild = (internalSlots>internalChild)?internalChild:internalSlots;
520       }
521     }
522     
523     nChild = externalChild+internalChild;
524     CmiAssert(nChild>0);
525
526     //exclude the root node        
527     delete [] child;
528     child = new int[nChild];    
529
530     for(int i=0; i<externalChild; i++) {
531         child[i] = ptnTree.item(i+1).peIDs[0];
532     }
533     for(int i=externalChild, j=1; i<nChild; i++, j++) {
534         child[i] = rootnode->peIDs[j];
535     }
536
537 #ifdef USE_NODEPATCHMGR
538     //only register the cores that have proxy patches. The HomePach's core
539     //doesn't need to be registered.
540     CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
541     NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
542     if(rootnode->numPes==1){
543         npm->registerPatch(patchID, 0, NULL);        
544     }
545     else{
546         npm->registerPatch(patchID, rootnode->numPes-1, &rootnode->peIDs[1]);        
547     }
548
549     //set up childrens in terms of node ids
550     numNodeChild = externalChild;
551     if(internalChild) numNodeChild++;
552     delete [] nodeChildren;
553     nodeChildren = new int[numNodeChild];    
554     for(int i=0; i<externalChild; i++) {
555         nodeChildren[i] = ptnTree.item(i+1).nodeID;        
556     }
557     //the last entry always stores this node id if there are proxies
558     //on other cores of the same node for this patch.
559     if(internalChild)
560         nodeChildren[numNodeChild-1] = rootnode->nodeID;
561 #endif
562     
563 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
564     DebugFileTrace *dft = DebugFileTrace::Object();
565     dft->openTrace();
566     dft->writeTrace("HomePatch[%d] has %d children: ", patchID, nChild);
567     for(int i=0; i<nChild; i++)
568         dft->writeTrace("%d ", child[i]);
569     dft->writeTrace("\n");
570     dft->closeTrace();
571 #endif   
572 }
573 #endif
574
575 #ifdef NODEAWARE_PROXY_SPANNINGTREE
576 //This is not an entry method, but takes an argument of message type
577 void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){
578     //set up the whole tree ptnTree
579     int treesize = msg->numNodesWithProxies;    
580     ptnTree.resize(treesize);    
581     int *pAllPes = msg->allPes;
582     for(int i=0; i<treesize; i++) {
583         proxyTreeNode *oneNode = &ptnTree.item(i);
584         delete [] oneNode->peIDs;
585         oneNode->numPes = msg->numPesOfNode[i];
586         oneNode->nodeID = CkNodeOf(*pAllPes);
587         oneNode->peIDs = new int[oneNode->numPes];
588         for(int j=0; j<oneNode->numPes; j++) {
589             oneNode->peIDs[j] = *pAllPes;
590             pAllPes++;
591         }
592     }
593     //setup children
594     setupChildrenFromProxySpanningTree();
595     //send down to children
596     sendNodeAwareSpanningTree();
597 }
598
599 void HomePatch::sendNodeAwareSpanningTree(){
600     if(ptnTree.size()==0) return;    
601     ProxyNodeAwareSpanningTreeMsg *msg = 
602         ProxyNodeAwareSpanningTreeMsg::getANewMsg(patchID, CkMyPe(), ptnTree.begin(), ptnTree.size());
603    
604     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
605     msg->printOut("HP::sendST");
606     #endif
607
608     CmiAssert(CkMyPe() == msg->allPes[0]);
609     ProxyMgr::Object()->sendNodeAwareSpanningTree(msg);
610
611 }
612 #else
613 void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){}
614 void HomePatch::sendNodeAwareSpanningTree(){}
615 #endif
616
617 #ifndef NODEAWARE_PROXY_SPANNINGTREE
618 // recv a spanning tree from load balancer
619 void HomePatch::recvSpanningTree(int *t, int n)
620 {
621   int i;
622   nChild=0;
623   tree.resize(n);
624   for (i=0; i<n; i++) {
625     tree[i] = t[i];
626   }
627
628   for (i=1; i<=proxySpanDim; i++) {
629     if (tree.size() <= i) break;
630     child[i-1] = tree[i];
631     nChild++;
632   }
633
634 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
635   isProxyChanged = 1;
636 #endif
637
638   // send down to children
639   sendSpanningTree();
640 }
641
642 void HomePatch::sendSpanningTree()
643 {
644   if (tree.size() == 0) return;
645   ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
646   msg->patch = patchID;
647   msg->node = CkMyPe();
648   msg->tree.copy(tree);  // copy data for thread safety
649   ProxyMgr::Object()->sendSpanningTree(msg);  
650 }
651 #else
652 void HomePatch::recvSpanningTree(int *t, int n){}
653 void HomePatch::sendSpanningTree(){}
654 #endif
655
656 #ifndef NODEAWARE_PROXY_SPANNINGTREE
657 void HomePatch::buildSpanningTree(void)
658 {
659   nChild = 0;
660   int psize = proxy.size();
661   if (psize == 0) return;
662   NodeIDList oldtree;  oldtree.swap(tree);
663   int oldsize = oldtree.size();
664   tree.resize(psize + 1);
665   tree.setall(-1);
666   tree[0] = CkMyPe();
667   int s=1, e=psize+1;
668   NodeIDList::iterator pli;
669   int patchNodesLast =
670     ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
671   int nNonPatch = 0;
672 #if 1
673     // try to put it to the same old tree
674   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
675   {
676     int oldindex = oldtree.find(*pli);
677     if (oldindex != -1 && oldindex < psize) {
678       tree[oldindex] = *pli;
679     }
680   }
681   s=1; e=psize;
682   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
683   {
684     if (tree.find(*pli) != -1) continue;    // already assigned
685     if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
686       while (tree[e] != -1) { e--; if (e==-1) e = psize; }
687       tree[e] = *pli;
688     } else {
689       while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
690       tree[s] = *pli;
691       nNonPatch++;
692     }
693   }
694 #if 1
695   if (oldsize==0 && nNonPatch) {
696     // first time, sort by distance
697     qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
698   }
699 #endif
700
701   //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
702
703 #if USE_TOPOMAP && USE_SPANNING_TREE
704   
705   //Right now only works for spanning trees with two levels
706   int *treecopy = new int [psize];
707   int subroots[proxySpanDim];
708   int subsizes[proxySpanDim];
709   int subtrees[proxySpanDim][proxySpanDim];
710   int idxes[proxySpanDim];
711   int i = 0;
712
713   for(i=0;i<proxySpanDim;i++){
714     subsizes[i] = 0;
715     idxes[i] = i;
716   }
717   
718   for(i=0;i<psize;i++){
719     treecopy[i] = tree[i+1];
720   }
721   
722   TopoManager tmgr;
723   tmgr.sortRanksByHops(treecopy,nNonPatch);
724   tmgr.sortRanksByHops(treecopy+nNonPatch,
725                                                 psize-nNonPatch);  
726   
727   /* build tree and subtrees */
728   nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
729   delete [] treecopy;
730   
731   for(int i=1;i<psize+1;i++){
732     int isSubroot=0;
733     for(int j=0;j<nChild;j++)
734       if(tree[i]==subroots[j]){
735         isSubroot=1;
736         break;
737       }
738     if(isSubroot) continue;
739     
740     int bAdded = 0;
741     tmgr.sortIndexByHops(tree[i], subroots,
742                                                   idxes, proxySpanDim);
743     for(int j=0;j<proxySpanDim;j++){
744       if(subsizes[idxes[j]]<proxySpanDim){
745         subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
746         bAdded = 1; 
747         break;
748       }
749     }
750     if( psize > proxySpanDim && ! bAdded ) {
751       NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
752     }
753   }
754
755 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
756   
757   for (int i=1; i<=proxySpanDim; i++) {
758     if (tree.size() <= i) break;
759     child[i-1] = tree[i];
760     nChild++;
761   }
762 #endif
763 #endif
764   
765 #if 0
766   // for debugging
767   CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
768   for (int i=0; i<psize+1; i++) {
769     CkPrintf("%d ", tree[i]);
770   }
771   CkPrintf("\n");
772 #endif
773   // send to children nodes
774   sendSpanningTree();
775 }
776 #endif
777
778
779 void HomePatch::receiveResults(ProxyResultVarsizeMsg *msg) {
780
781     numGBISP3Arrived++;
782     DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
783     int n = msg->node;
784     Results *r = forceBox.clientOpen();
785
786     char *iszeroPtr = msg->isZero;
787     Force *msgFPtr = msg->forceArr;
788
789     for ( int k = 0; k < Results::maxNumForces; ++k )
790     {
791       Force *rfPtr = r->f[k];
792       for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
793           if((*iszeroPtr)!=1) {
794               *rfPtr += *msgFPtr;
795               msgFPtr++;
796           }
797       }      
798     }
799     forceBox.clientClose();
800     delete msg;
801 }
802
803 void HomePatch::receiveResults(ProxyResultMsg *msg) {
804   numGBISP3Arrived++;
805   DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
806   int n = msg->node;
807   Results *r = forceBox.clientOpen();
808   for ( int k = 0; k < Results::maxNumForces; ++k )
809   {
810     Force *f = r->f[k];
811     register ForceList::iterator f_i, f_e;
812     f_i = msg->forceList[k]->begin();
813     f_e = msg->forceList[k]->end();
814     for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
815   }
816   forceBox.clientClose();
817   delete msg;
818 }
819
820 void HomePatch::receiveResults(ProxyCombinedResultRawMsg* msg)
821 {
822     numGBISP3Arrived++;
823   DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodeSize<<")\n");
824     Results *r = forceBox.clientOpen(msg->nodeSize);
825       register char* isNonZero = msg->isForceNonZero;
826       register Force* f_i = msg->forceArr;
827       for ( int k = 0; k < Results::maxNumForces; ++k )
828       {
829         Force *f = r->f[k];
830                 int nf = msg->flLen[k];
831 #ifdef ARCH_POWERPC
832 #pragma disjoint (*f_i, *f)
833 #endif
834         for (int count = 0; count < nf; count++) {
835           if(*isNonZero){
836                         f[count].x += f_i->x;
837                         f[count].y += f_i->y;
838                         f[count].z += f_i->z;
839                         f_i++;
840           }
841           isNonZero++;
842         }
843       }
844     forceBox.clientClose(msg->nodeSize);
845
846     delete msg;
847 }
848
849 void HomePatch::qmSwapAtoms()
850 {
851     // This is used for LSS in QM/MM simulations.
852     // Changes atom labels so that we effectively exchange solvent
853     // molecules between classical and quantum modes.
854     
855     SimParameters *simParams = Node::Object()->simParameters;
856     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
857     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
858     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
859     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
860     
861     ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
862             CkpvAccess(BOCclass_group).computeQMMgr) ;
863     
864     FullAtom *a_i = atom.begin();
865     
866     for (int i=0; i<numAtoms; ++i ) { 
867         
868         LSSSubsDat *subP = lssSubs(mgrP).find( LSSSubsDat(a_i[i].id) ) ;
869         
870         if ( subP != NULL ) {
871             a_i[i].id = subP->newID ;
872             a_i[i].vdwType = subP->newVdWType ;
873             
874             // If we are swappign a classical atom with a QM one, the charge
875             // will need extra handling.
876             if (qmAtomGroup[subP->newID] > 0 && simParams->PMEOn) {
877                 // We make sure that the last atom charge calculated for the 
878                 // QM atom being transfered here is available for PME
879                 // in the next step.
880                 
881                 // Loops over all QM atoms (in all QM groups) comparing their 
882                 // global indices (the sequential atom ID from NAMD).
883                 for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
884                     
885                     if (qmAtmIndx[qmIter] == subP->newID) {
886                         qmAtmChrg[qmIter] = subP->newCharge;
887                         break;
888                     }
889                     
890                 }
891                 
892                 // For QM atoms, the charge in the full atom structure is zero.
893                 // Electrostatic interactions between QM atoms and their 
894                 // environment is handled in ComputeQM.
895                 a_i[i].charge = 0;
896             }
897             else {
898                 // If we are swappign a QM atom with a Classical one, only the charge
899                 // in the full atomstructure needs updating, since it used to be zero.
900                 a_i[i].charge = subP->newCharge ;
901             }
902         }
903     }
904     
905     return;
906 }
907
908 void HomePatch::positionsReady(int doMigration)
909 {    
910   SimParameters *simParams = Node::Object()->simParameters;
911
912   flags.sequence++;
913
914   if (!patchMapRead) { readPatchMap(); }
915       
916   if (numNeighbors && ! simParams->staticAtomAssignment) {
917     if (doMigration) {
918       rattleListValid = false;
919       doAtomMigration();
920     } else {
921       doMarginCheck();
922     }
923   }
924   
925   if (doMigration && simParams->qmLSSOn)
926       qmSwapAtoms();
927
928 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
929   if ( doMigration ) {
930     int n = numAtoms;
931     FullAtom *a_i = atom.begin();
932     #if defined(NAMD_CUDA) || (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
933     int *ao = new int[n];
934     int nfree;
935     if ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces ) {
936       int k = 0;
937       int k2 = n;
938       for ( int j=0; j<n; ++j ) {
939         // put fixed atoms at end
940         if ( a_i[j].atomFixed ) ao[--k2] = j;
941         else ao[k++] = j;
942       }
943       nfree = k;
944     } else {
945       nfree = n;
946       for ( int j=0; j<n; ++j ) {
947         ao[j] = j;
948       }
949     }
950
951     sortAtomsForCUDA(ao,a_i,nfree,n);
952   
953     for ( int i=0; i<n; ++i ) { 
954       a_i[i].sortOrder = ao[i];
955     }
956     delete [] ao;
957     #else
958       for (int i = 0; i < n; ++i) {
959         a_i[i].sortOrder = i;
960       }
961     #endif
962   }
963
964   { 
965     const double charge_scaling = sqrt(COULOMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1);
966     const Vector ucenter = lattice.unscale(center);
967     const BigReal ucenter_x = ucenter.x;
968     const BigReal ucenter_y = ucenter.y;
969     const BigReal ucenter_z = ucenter.z;
970     const int n = numAtoms;
971     #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
972       int n_16 = n;
973       n_16 = (n + 15) & (~15);
974       cudaAtomList.resize(n_16);
975       CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
976       mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
977       mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
978       mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
979       mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
980     #else
981     cudaAtomList.resize(n);
982     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
983     #endif
984     const FullAtom *a = atom.begin();
985     for ( int k=0; k<n; ++k ) {
986       #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
987         int j = a[k].sortOrder;
988         atom_x[k] = a[j].position.x - ucenter_x;
989         atom_y[k] = a[j].position.y - ucenter_y;
990         atom_z[k] = a[j].position.z - ucenter_z;
991         atom_q[k] = charge_scaling * a[j].charge;
992       #else
993       int j = a[k].sortOrder;
994       ac[k].x = a[j].position.x - ucenter_x;
995       ac[k].y = a[j].position.y - ucenter_y;
996       ac[k].z = a[j].position.z - ucenter_z;
997       ac[k].q = charge_scaling * a[j].charge;
998       #endif
999     }
1000   }
1001 #else
1002   doMigration = doMigration && numNeighbors;
1003 #endif
1004   doMigration = doMigration || ! patchMapRead;
1005
1006   doMigration = doMigration || doAtomUpdate;
1007   doAtomUpdate = false;
1008
1009   // Workaround for oversize groups
1010   doGroupSizeCheck();
1011
1012   // Copy information needed by computes and proxys to Patch::p.
1013   p.resize(numAtoms);
1014   CompAtom *p_i = p.begin();
1015   pExt.resize(numAtoms);
1016   CompAtomExt *pExt_i = pExt.begin();
1017   FullAtom *a_i = atom.begin();
1018   int i; int n = numAtoms;
1019   for ( i=0; i<n; ++i ) { 
1020     p_i[i] = a_i[i]; 
1021     pExt_i[i] = a_i[i];
1022   }
1023
1024   // Measure atom movement to test pairlist validity
1025   doPairlistCheck();
1026
1027   if (flags.doMolly) mollyAverage();
1028   // BEGIN LA
1029   if (flags.doLoweAndersen) loweAndersenVelocities();
1030   // END LA
1031
1032     if (flags.doGBIS) {
1033       //reset for next time step
1034       numGBISP1Arrived = 0;
1035       phase1BoxClosedCalled = false;
1036       numGBISP2Arrived = 0;
1037       phase2BoxClosedCalled = false;
1038       numGBISP3Arrived = 0;
1039       phase3BoxClosedCalled = false;
1040       if (doMigration || isNewProxyAdded)
1041         setGBISIntrinsicRadii();
1042     }
1043
1044   if (flags.doLCPO) {
1045     if (doMigration || isNewProxyAdded) {
1046       setLcpoType();
1047     }
1048   }
1049
1050   // Must Add Proxy Changes when migration completed!
1051   NodeIDList::iterator pli;
1052   int *pids = NULL;
1053   int pidsPreAllocated = 1;
1054   int npid;
1055   if (proxySendSpanning == 0) {
1056     npid = proxy.size();
1057     pids = new int[npid];
1058     pidsPreAllocated = 0;
1059     int *pidi = pids;
1060     int *pide = pids + proxy.size();
1061     int patchNodesLast =
1062       ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
1063     for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
1064     {
1065       if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
1066         *(--pide) = *pli;
1067       } else {
1068         *(pidi++) = *pli;
1069       }
1070     }
1071   }
1072   else {
1073 #ifdef NODEAWARE_PROXY_SPANNINGTREE
1074     #ifdef USE_NODEPATCHMGR
1075     npid = numNodeChild;
1076     pids = nodeChildren;
1077     #else
1078     npid = nChild;
1079     pids = child;
1080     #endif
1081 #else
1082     npid = nChild;
1083     pidsPreAllocated = 0;
1084     pids = new int[proxySpanDim];
1085     for (int i=0; i<nChild; i++) pids[i] = child[i];
1086 #endif
1087   }
1088   if (npid) { //have proxies
1089     int seq = flags.sequence;
1090     int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
1091     //begin to prepare proxy msg and send it
1092     int pdMsgPLLen = p.size();
1093     int pdMsgAvgPLLen = 0;
1094     if(flags.doMolly) {
1095         pdMsgAvgPLLen = p_avg.size();
1096     }
1097     // BEGIN LA
1098     int pdMsgVLLen = 0;
1099     if (flags.doLoweAndersen) {
1100         pdMsgVLLen = v.size();
1101     }
1102     // END LA
1103
1104     int intRadLen = 0;
1105     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
1106             intRadLen = numAtoms * 2;
1107     }
1108
1109     //LCPO
1110     int lcpoTypeLen = 0;
1111     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
1112             lcpoTypeLen = numAtoms;
1113     }
1114
1115     int pdMsgPLExtLen = 0;
1116     if(doMigration || isNewProxyAdded) {
1117         pdMsgPLExtLen = pExt.size();
1118     }
1119
1120     int cudaAtomLen = 0;
1121 #ifdef NAMD_CUDA
1122     cudaAtomLen = numAtoms;
1123 #endif
1124
1125     #ifdef NAMD_MIC
1126       #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
1127         cudaAtomLen = numAtoms;
1128       #else
1129         cudaAtomLen = (numAtoms + 15) & (~15);
1130       #endif
1131     #endif
1132
1133     ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
1134       lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
1135
1136     SET_PRIORITY(nmsg,seq,priority);
1137     nmsg->patch = patchID;
1138     nmsg->flags = flags;
1139     nmsg->plLen = pdMsgPLLen;                
1140     //copying data to the newly created msg
1141     memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
1142     nmsg->avgPlLen = pdMsgAvgPLLen;        
1143     if(flags.doMolly) {
1144         memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
1145     }
1146     // BEGIN LA
1147     nmsg->vlLen = pdMsgVLLen;
1148     if (flags.doLoweAndersen) {
1149         memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
1150     }
1151     // END LA
1152
1153     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
1154       for (int i = 0; i < numAtoms * 2; i++) {
1155         nmsg->intRadList[i] = intRad[i];
1156       }
1157     }
1158
1159     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
1160       for (int i = 0; i < numAtoms; i++) {
1161         nmsg->lcpoTypeList[i] = lcpoType[i];
1162       }
1163     }
1164
1165     nmsg->plExtLen = pdMsgPLExtLen;
1166     if(doMigration || isNewProxyAdded){     
1167         memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
1168     }
1169
1170 // DMK
1171 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
1172     memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
1173 #endif
1174     
1175 #if NAMD_SeparateWaters != 0
1176     //DMK - Atom Separation (water vs. non-water)
1177     nmsg->numWaterAtoms = numWaterAtoms;
1178 #endif
1179
1180 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
1181     nmsg->isFromImmMsgCall = 0;
1182 #endif
1183     
1184     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
1185     DebugFileTrace *dft = DebugFileTrace::Object();
1186     dft->openTrace();
1187     dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
1188     for(int i=0; i<npid; i++) {
1189         dft->writeTrace("%d ", pids[i]);
1190     }
1191     dft->writeTrace("\n");
1192     dft->closeTrace();
1193     #endif
1194
1195 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1196     if (isProxyChanged || localphs == NULL)
1197     {
1198 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
1199      //CmiAssert(isProxyChanged);
1200      if (nphs) {
1201        for (int i=0; i<nphs; i++) {
1202          CmiDestoryPersistent(localphs[i]);
1203        }
1204        delete [] localphs;
1205      }
1206      localphs = new PersistentHandle[npid];
1207      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;
1208      for (int i=0; i<npid; i++) {
1209 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
1210        if (proxySendSpanning)
1211            localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
1212        else
1213 #endif
1214        localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
1215      }
1216      nphs = npid;
1217     }
1218     CmiAssert(nphs == npid && localphs != NULL);
1219     CmiUsePersistentHandle(localphs, nphs);
1220 #endif
1221     if(doMigration || isNewProxyAdded) {
1222         ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
1223     }else{
1224         ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
1225     }
1226 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
1227     CmiUsePersistentHandle(NULL, 0);
1228 #endif
1229     isNewProxyAdded = 0;
1230   }
1231   isProxyChanged = 0;
1232   if(!pidsPreAllocated) delete [] pids;
1233   DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
1234
1235 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
1236   positionPtrBegin = p.begin();
1237   positionPtrEnd = p.end();
1238 #endif
1239
1240   if(flags.doMolly) {
1241       avgPositionPtrBegin = p_avg.begin();
1242       avgPositionPtrEnd = p_avg.end();
1243   }
1244   // BEGIN LA
1245   if (flags.doLoweAndersen) {
1246       velocityPtrBegin = v.begin();
1247       velocityPtrEnd = v.end();
1248   }
1249   // END LA
1250
1251   Patch::positionsReady(doMigration);
1252
1253   patchMapRead = 1;
1254
1255   // gzheng
1256   Sync::Object()->PatchReady();
1257 }
1258
1259 void HomePatch::replaceForces(ExtForce *f)
1260 {
1261   replacementForces = f;
1262 }
1263
1264 void HomePatch::saveForce(const int ftag)
1265 {
1266   f_saved[ftag].resize(numAtoms);
1267   for ( int i = 0; i < numAtoms; ++i )
1268   {
1269     f_saved[ftag][i] = f[ftag][i];
1270   }
1271 }
1272
1273
1274 #undef DEBUG_REDISTRIB_FORCE 
1275 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
1276 //#define DEBUG_REDISTRIB_FORCE
1277
1278 /**
1279  * Redistribute forces from lone pair site onto its colinear host atoms.
1280  *
1281  * /param fi the force vector on lonepair i
1282  * /param fj the force vector on the first host j
1283  * /param fk the force vector on the second host k
1284  * /param ri the position vector of the lone pair i
1285  * /param rj the position vector of the first host j
1286  * /param rk the position vector of the second host k
1287  * /param distance the distance to be placed from the reference point along the
1288  *     vector from k to j
1289  * /param scale the reference point is placed using a scaled value of the
1290  *     vector from k to j.
1291  */
1292 void HomePatch::redistrib_colinear_lp_force(
1293     Vector& fi, Vector& fj, Vector& fk,
1294     const Vector& ri, const Vector& rj, const Vector& rk,
1295     Real distance_f, Real scale_f) {
1296   BigReal distance = distance_f;
1297   BigReal scale = scale_f;
1298   Vector r_jk = rj - rk;
1299   // TODO: Add a check for small distances?
1300   BigReal r_jk_rlength = r_jk.rlength();
1301   distance *= r_jk_rlength;
1302   BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
1303   fj += (1. + scale + distance)*fi - r_jk*fdot;
1304   fk -= (scale + distance)*fi + r_jk*fdot;
1305 }
1306
1307 /**
1308  * Redistribute forces from lone pair site onto its host atoms.
1309  *
1310  * Atoms are "labeled" i, j, k, l, where atom i is the lone pair.
1311  * Positions of atoms are ri, rj, rk, rl.
1312  * Forces of atoms are fi, fj, fk, fl; updated forces are returned.
1313  * Accumulate updates to the virial.
1314  *
1315  * The forces on the atoms are updated so that:
1316  *   - the force fi on the lone pair site is 0
1317  *   - the net force fi+fj+fk+fl is conserved
1318  *   - the net torque cross(ri,fi)+cross(rj,fj)+cross(rk,fk)+cross(rl,fl)
1319  *     is conserved
1320  *
1321  * If "midpt" is true (nonzero), then use the midpoint of rk and rl
1322  * (e.g. rk and rl are the hydrogen atoms for water).  Otherwise use
1323  * planes defined by ri,rj,rk and rj,rk,rl.
1324  *
1325  * Having "midpt" set true corresponds in CHARMM to having a negative
1326  * distance parameter in Lphost.
1327  *
1328  * Use FIX_FOR_WATER below to fix the case that occurs when the lone pair
1329  * site for water lies on the perpendicular bisector of rk and rl, making
1330  * the cross(r,s) zero.
1331  */
1332 #define FIX_FOR_WATER
1333 void HomePatch::redistrib_relative_lp_force(
1334     Vector& fi, Vector& fj, Vector& fk, Vector& fl,
1335     const Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
1336     Tensor *virial, int midpt) {
1337 #ifdef DEBUG_REDISTRIB_FORCE
1338   Vector foldnet, toldnet;  // old net force, old net torque
1339   foldnet = fi + fj + fk + fl;
1340   toldnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
1341 #endif
1342   Vector fja(0), fka(0), fla(0);
1343
1344   Vector r = ri - rj;
1345   BigReal invr2 = r.rlength();
1346   invr2 *= invr2;
1347   BigReal fdot = (fi*r) * invr2;
1348   Vector fr = r * fdot;
1349
1350   fja += fr;
1351
1352   Vector s, t;
1353   if (midpt) {
1354     s = rj - 0.5*(rk + rl);
1355     t = 0.5*(rk - rl);
1356   }
1357   else {
1358     s = rj - rk;
1359     t = rk - rl;
1360   }
1361   BigReal invs2 = s.rlength();
1362   invs2 *= invs2;
1363
1364   Vector p = cross(r,s);
1365 #if !defined(FIX_FOR_WATER)
1366   BigReal invp = p.rlength();
1367 #else
1368   BigReal p2 = p.length2();  // fix division by zero above
1369 #endif
1370
1371   Vector q = cross(s,t);
1372   BigReal invq = q.rlength();
1373
1374 #if !defined(FIX_FOR_WATER)
1375   BigReal fpdot = (fi*p) * invp;
1376   Vector fp = p * fpdot;
1377   Vector ft = fi - fr - fp;
1378 #else
1379   BigReal fpdot;
1380   Vector fp, ft;
1381   if (p2 < 1e-6) {  // vector is near zero, assume no fp contribution to force
1382     fpdot = 0;
1383     fp = 0;
1384     ft = fi - fr;
1385   }
1386   else {
1387     fpdot = (fi*p) / sqrt(p2);
1388     fp = p * fpdot;
1389     ft = fi - fr - fp;
1390   }
1391 #endif
1392
1393   fja += ft;
1394   Vector v = cross(r,ft);  // torque
1395   ft = cross(s,v) * invs2;
1396   fja -= ft;
1397
1398   if (midpt) {
1399     fka += 0.5 * ft;
1400     fla += 0.5 * ft;
1401   }
1402   else {
1403     fka += ft;
1404   }
1405
1406   BigReal srdot = (s*r) * invs2;
1407   Vector rr = r - s*srdot;
1408   BigReal rrdot = rr.length();
1409   BigReal stdot = (s*t) * invs2;
1410   Vector tt = t - s*stdot;
1411   BigReal invtt = tt.rlength();
1412   BigReal fact = rrdot*fpdot*invtt*invq;
1413   Vector fq = q * fact;
1414
1415   fla += fq;
1416   fja += fp*(1+srdot) + fq*stdot;
1417
1418   ft = fq*(1+stdot) + fp*srdot;
1419
1420   if (midpt) {
1421     fka += -0.5*ft;
1422     fla += -0.5*ft;
1423   }
1424   else {
1425     fka -= ft;
1426   }
1427
1428   if (virial) {
1429     Tensor va = outer(fja,rj);
1430     va += outer(fka,rk);
1431     va += outer(fla,rl);
1432     va -= outer(fi,ri);
1433     *virial += va;
1434   }
1435
1436   fi = 0;  // lone pair has zero force
1437   fj += fja;
1438   fk += fka;
1439   fl += fla;
1440
1441 #ifdef DEBUG_REDISTRIB_FORCE
1442 #define TOL_REDISTRIB  1e-4
1443   Vector fnewnet, tnewnet;  // new net force, new net torque
1444   fnewnet = fi + fj + fk + fl;
1445   tnewnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
1446   Vector fdiff = fnewnet - foldnet;
1447   Vector tdiff = tnewnet - toldnet;
1448   if (fdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1449     printf("Error:  force redistribution for water exceeded tolerance:  "
1450         "fdiff=(%f, %f, %f)\n", fdiff.x, fdiff.y, fdiff.z);
1451   }
1452   if (tdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
1453     printf("Error:  torque redistribution for water exceeded tolerance:  "
1454         "tdiff=(%f, %f, %f)\n", tdiff.x, tdiff.y, tdiff.z);
1455   }
1456 #endif
1457 }
1458
1459
1460 /* Redistribute forces from the massless lonepair charge particle onto
1461  * the other atoms of the water.
1462  *
1463  * This is done using the same algorithm as charmm uses for TIP4P lonepairs.
1464  *
1465  * Pass by reference the forces (O H1 H2 LP) to be modified,
1466  * pass by constant reference the corresponding positions,
1467  * and a pointer to virial.
1468  */
1469 void HomePatch::redistrib_lp_water_force(
1470     Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
1471     const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
1472     const Vector& p_lp, Tensor *virial) {
1473
1474 #ifdef DEBUG_REDISTRIB_FORCE 
1475   // Debug information to check against results at end
1476
1477   // total force and torque relative to origin
1478   Vector totforce, tottorque;
1479   totforce = f_ox + f_h1 + f_h2 + f_lp;
1480   tottorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
1481   //printf("Torque without LP is %f/%f/%f\n",
1482   //    tottorque.x, tottorque.y, tottorque.z);
1483   tottorque += cross(f_lp, p_lp);
1484   //printf("Torque with LP is %f/%f/%f\n",
1485   //    tottorque.x, tottorque.y, tottorque.z);
1486 #endif
1487
1488   // accumulate force adjustments
1489   Vector fad_ox(0), fad_h(0);
1490
1491   // Calculate the radial component of the force and add it to the oxygen
1492   Vector r_ox_lp = p_lp - p_ox;
1493   BigReal invlen2_r_ox_lp = r_ox_lp.rlength();
1494   invlen2_r_ox_lp *= invlen2_r_ox_lp;
1495   BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
1496   Vector f_rad = r_ox_lp * rad_factor;
1497
1498   fad_ox += f_rad;
1499
1500   // Calculate the angular component
1501   Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
1502   Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5; // half of r_h2_h1
1503
1504   // deviation from collinearity of charge site
1505   //Vector r_oop = cross(r_ox_lp, r_hcom_ox);
1506   //
1507   // vector out of o-h-h plane
1508   //Vector r_perp = cross(r_hcom_ox, r_h2_h1_2);
1509
1510   // Here we assume that Ox/Lp/Hcom are linear
1511   // If you want to correct for deviations, this is the place
1512
1513 //  printf("Deviation from linearity for ox %i: %f/%f/%f\n", oxind, r_oop.x, r_oop.y, r_oop.z);
1514
1515   Vector f_ang = f_lp - f_rad; // leave the angular component
1516
1517   // now split this component onto the other atoms
1518   BigReal len_r_ox_lp = r_ox_lp.length();
1519   BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
1520   BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
1521   BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
1522
1523   fad_ox += (f_ang * oxcomp);
1524   fad_h += (f_ang * hydcomp);  // adjustment for both hydrogens
1525
1526   // Add virial contributions
1527   if (virial) {
1528     Tensor vir = outer(fad_ox, p_ox);
1529     vir += outer(fad_h, p_h1);
1530     vir += outer(fad_h, p_h2);
1531     vir -= outer(f_lp, p_lp);
1532     *virial += vir;
1533   }
1534
1535   //Vector zerovec(0.0, 0.0, 0.0);
1536   f_lp = 0;
1537   f_ox += fad_ox;
1538   f_h1 += fad_h;
1539   f_h2 += fad_h;
1540
1541 #ifdef DEBUG_REDISTRIB_FORCE 
1542   // Check that the total force and torque come out right
1543   Vector newforce, newtorque;
1544   newforce = f_ox + f_h1 + f_h2;
1545   newtorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
1546   Vector fdiff = newforce - totforce;
1547   Vector tdiff = newtorque - tottorque;
1548   BigReal error = fdiff.length();
1549   if (error > 0.0001) {
1550      printf("Error:  Force redistribution for water "
1551          "exceeded force tolerance:  error=%f\n", error);
1552   }
1553 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1554   printf("Error in net force:  %f\n", error);
1555 #endif
1556
1557   error = tdiff.length();
1558   if (error > 0.0001) {
1559      printf("Error:  Force redistribution for water "
1560          "exceeded torque tolerance:  error=%f\n", error);
1561   }
1562 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
1563   printf("Error in net torque:  %f\n", error);
1564 #endif
1565 #endif /* DEBUG */
1566 }
1567
1568 /**
1569  * Reposition lonepair i in a colinear fashion relative to its hosts j and k.
1570  *
1571  * The lonepair is placed at a fixed (possibly negative) distance along the
1572  * vector from k to j relative to a reference point. The reference point is
1573  * computed by multiplying the vector from k to j by a (possibly negative)
1574  * scale factor. For example, a scale value of 0 (the default) sets the
1575  * reference point as rj, while a value of -1 sets the reference point as rk.
1576  * A scale value of -0.5 would use the center of the two hosts.
1577  *
1578  * /param ri the position vector of the lonepair, i
1579  * /param rj the position vector of the first host, j
1580  * /param rk the position vector of the second host, k
1581  * /param distance the distance to be placed from the reference point along the
1582  *     vector from k to j
1583  * /param scale the reference point is placed using a scaled value of the
1584  *     vector from k to j.
1585  */
1586 void HomePatch::reposition_colinear_lonepair(
1587     Vector& ri, const Vector& rj, const Vector& rk,
1588     Real distance_f, Real scale_f)
1589 {
1590   BigReal distance = distance_f;
1591   BigReal scale = scale_f;
1592   Vector r_jk = rj - rk;
1593   BigReal r2 = r_jk.length2();
1594   if (r2 < 1e-10 || 100. < r2) { // same low tolerance as used in CHARMM
1595     iout << iWARN << "Large/small distance between lonepair reference atoms: ("
1596          << rj << ") (" << rk << ")\n" << endi;
1597   }
1598   ri = rj + (scale + distance*r_jk.rlength())*r_jk;
1599 }
1600
1601 /**
1602  * Reposition lonepair i relative to its hosts j, k, and l.
1603  *
1604  * The lonepair is placed in a bond-angle-torsion coordinate frame. An angle
1605  * and torsion of 0 corresponds to a bisector placement, as in TIP4P water.
1606  *
1607  * /param ri the position vector of the lonepair, i
1608  * /param rj the position vector of the first host, j
1609  * /param rk the position vector of the second host, k
1610  * /param rl the position vector of the third host, j
1611  * /param distance the negative distance from j
1612  * /param angle the angle made by i-j-k
1613  * /param dihedral the dihedral made by i-j-k-l
1614  */
1615 void HomePatch::reposition_relative_lonepair(
1616     Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
1617     Real distance, Real angle, Real dihedral)
1618 {
1619   if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
1620     iout << iWARN << "Large distance between lonepair reference atoms: ("
1621       << rj << ") (" << rk << ") (" << rl << ")\n" << endi;
1622   }
1623   BigReal r, t, p, cst, snt, csp, snp, invlen;
1624   Vector v, w, a, b, c;
1625
1626   if (distance >= 0) {
1627     v = rk;
1628     r = distance;
1629   }
1630   else {
1631     v = 0.5*(rk + rl);
1632     r = -distance;
1633   }
1634
1635   t = angle;
1636   p = dihedral;
1637   cst = cos(t);
1638   snt = sin(t);
1639   csp = cos(p);
1640   snp = sin(p);
1641   a = v - rj;
1642   b = rl - v;
1643   invlen = a.rlength();
1644   a *= invlen;
1645   c = cross(b, a);
1646   invlen = c.rlength();
1647   c *= invlen;
1648   b = cross(a, c);
1649   w.x = r*cst;
1650   w.y = r*snt*csp;
1651   w.z = r*snt*snp;
1652   ri.x = rj.x + w.x*a.x + w.y*b.x + w.z*c.x;
1653   ri.y = rj.y + w.x*a.y + w.y*b.y + w.z*c.y;
1654   ri.z = rj.z + w.x*a.z + w.y*b.z + w.z*c.z;
1655 }
1656
1657 void HomePatch::reposition_all_lonepairs(void) {
1658   // ASSERT: simParams->lonepairs == TRUE
1659   for (int i=0;  i < numAtoms;  i++) {
1660     if (atom[i].mass < 0.01) {
1661       // found a lone pair
1662       AtomID aid = atom[i].id;  // global atom ID of lp
1663       Lphost *lph = Node::Object()->molecule->get_lphost(aid);  // its lphost
1664       if (lph == NULL) {
1665         char errmsg[512];
1666         sprintf(errmsg, "reposition lone pairs: "
1667             "no Lphost exists for LP %d\n", aid);
1668         NAMD_die(errmsg);
1669       }
1670       LocalID j = AtomMap::Object()->localID(lph->atom2);
1671       LocalID k = AtomMap::Object()->localID(lph->atom3);
1672       LocalID l = AtomMap::Object()->localID(lph->atom4);
1673       if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
1674         char errmsg[512];
1675         sprintf(errmsg, "reposition lone pairs: "
1676             "LP %d has some Lphost atom off patch\n", aid);
1677         NAMD_die(errmsg);
1678       }
1679       // reposition this lone pair
1680       if (lph->numhosts == 2) {
1681         reposition_colinear_lonepair(atom[i].position, atom[j.index].position,
1682             atom[k.index].position, lph->distance, lph->angle);
1683       }
1684       else if (lph->numhosts == 3) {
1685         reposition_relative_lonepair(atom[i].position, atom[j.index].position,
1686             atom[k.index].position, atom[l.index].position,
1687             lph->distance, lph->angle, lph->dihedral);
1688       }
1689     }
1690   }
1691 }
1692
1693 void HomePatch::swm4_omrepos(Vector *ref, Vector *pos, Vector *vel,
1694     BigReal invdt) {
1695   // Reposition lonepair (Om) particle of Drude SWM4 water.
1696   // Same comments apply as to tip4_omrepos(), but the ordering of atoms
1697   // is different: O, D, LP, H1, H2.
1698   pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
1699   // Now, adjust velocity of particle to get it to appropriate place
1700   // during next integration "drift-step"
1701   if (invdt != 0) {
1702     vel[2] = (pos[2] - ref[2]) * invdt;
1703   }
1704   // No virial correction needed since lonepair is massless
1705 }
1706
1707 void HomePatch::tip4_omrepos(Vector* ref, Vector* pos, Vector* vel, BigReal invdt) {
1708   /* Reposition the om particle of a tip4p water
1709    * A little geometry shows that the appropriate position is given by
1710    * R_O + (1 / 2 r_ohc) * ( 0.5 (R_H1 + R_H2) - R_O ) 
1711    * Here r_om is the distance from the oxygen to Om site, and r_ohc
1712    * is the altitude from the oxygen to the hydrogen center of mass
1713    * Those quantities are precalculated upon initialization of HomePatch
1714    *
1715    * Ordering of TIP4P atoms: O, H1, H2, LP.
1716    */
1717
1718   //printf("rom/rohc are %f %f and invdt is %f\n", r_om, r_ohc, invdt);
1719   //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);
1720   pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc); 
1721   //printf("New position for lp is %f %f %f\n", pos[3].x, pos[3].y, pos[3].z);
1722
1723   // Now, adjust the velocity of the particle to get it to the appropriate place
1724   if (invdt != 0) {
1725     vel[3] = (pos[3] - ref[3]) * invdt;
1726   }
1727
1728   // No virial correction needed, since this is a massless particle
1729   return;
1730 }
1731
1732 void HomePatch::redistrib_lonepair_forces(const int ftag, Tensor *virial) {
1733   // ASSERT: simParams->lonepairs == TRUE
1734   ForceList *f_mod = f;
1735   for (int i = 0;  i < numAtoms;  i++) {
1736     if (atom[i].mass < 0.01) {
1737       // found a lone pair
1738       AtomID aid = atom[i].id;  // global atom ID of lp
1739       Lphost *lph = Node::Object()->molecule->get_lphost(aid);  // its lphost
1740       if (lph == NULL) {
1741         char errmsg[512];
1742         sprintf(errmsg, "redistrib lone pair forces: "
1743             "no Lphost exists for LP %d\n", aid);
1744         NAMD_die(errmsg);
1745       }
1746       LocalID j = AtomMap::Object()->localID(lph->atom2);
1747       LocalID k = AtomMap::Object()->localID(lph->atom3);
1748       LocalID l = AtomMap::Object()->localID(lph->atom4);
1749       if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
1750         char errmsg[512];
1751         sprintf(errmsg, "redistrib lone pair forces: "
1752             "LP %d has some Lphost atom off patch\n", aid);
1753         NAMD_die(errmsg);
1754       }
1755       // redistribute forces from this lone pair
1756       if (lph->numhosts == 2) {
1757         redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
1758             f_mod[ftag][k.index], atom[i].position, atom[j.index].position,
1759             atom[k.index].position, lph->distance, lph->angle);
1760       }
1761       else if (lph->numhosts == 3) {
1762         int midpt = (lph->distance < 0);
1763         redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
1764             f_mod[ftag][k.index], f_mod[ftag][l.index],
1765             atom[i].position, atom[j.index].position,
1766             atom[k.index].position, atom[l.index].position, virial, midpt);
1767       }
1768     }
1769   }
1770 }
1771
1772 void HomePatch::redistrib_swm4_forces(const int ftag, Tensor *virial) {
1773   // Loop over the patch's atoms and apply the appropriate corrections
1774   // to get all forces off of lone pairs
1775   ForceList *f_mod = f;
1776   for (int i = 0;  i < numAtoms;  i++) {
1777     if (atom[i].mass < 0.01) {
1778       // found lonepair
1779       redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
1780           f_mod[ftag][i+2], f_mod[ftag][i],
1781           atom[i-2].position, atom[i+1].position,
1782           atom[i+2].position, atom[i].position, virial);
1783     }
1784   }
1785 }
1786
1787 void HomePatch::redistrib_tip4p_forces(const int ftag, Tensor* virial) {
1788   // Loop over the patch's atoms and apply the appropriate corrections
1789   // to get all forces off of lone pairs
1790   // Atom ordering:  O H1 H2 LP
1791
1792   ForceList *f_mod =f;
1793   for (int i=0; i<numAtoms; i++) {
1794     if (atom[i].mass < 0.01) {
1795       // found lonepair
1796       redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
1797           f_mod[ftag][i-1], f_mod[ftag][i],
1798           atom[i-3].position, atom[i-2].position,
1799           atom[i-1].position, atom[i].position, virial);
1800     }
1801   }
1802 }
1803
1804
1805 void HomePatch::addForceToMomentum(
1806     FullAtom       * __restrict atom_arr,
1807     const Force    * __restrict force_arr,
1808     const BigReal    dt,
1809     int              num_atoms
1810     ) {
1811   SimParameters *simParams = Node::Object()->simParameters;
1812
1813   if ( simParams->fixedAtomsOn ) {
1814     for ( int i = 0; i < num_atoms; ++i ) {
1815       if ( atom_arr[i].atomFixed ) {
1816         atom_arr[i].velocity = 0;
1817       } else {
1818         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
1819         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1820         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1821         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1822       }
1823     }
1824   } else {
1825     for ( int i = 0; i < num_atoms; ++i ) {
1826       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
1827       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
1828       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
1829       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
1830     }
1831   }
1832 }
1833
1834 void HomePatch::addForceToMomentum3(
1835     FullAtom       * __restrict atom_arr,
1836     const Force    * __restrict force_arr1,
1837     const Force    * __restrict force_arr2,
1838     const Force    * __restrict force_arr3,
1839     const BigReal    dt1,
1840     const BigReal    dt2,
1841     const BigReal    dt3,
1842     int              num_atoms
1843     ) {
1844   SimParameters *simParams = Node::Object()->simParameters;
1845
1846   if ( simParams->fixedAtomsOn ) {
1847     for ( int i = 0; i < num_atoms; ++i ) {
1848       if ( atom_arr[i].atomFixed ) {
1849         atom_arr[i].velocity = 0;
1850       } else {
1851         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
1852         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
1853             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1854         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
1855             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1856         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
1857             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1858       }
1859     }
1860   } else {
1861     for ( int i = 0; i < num_atoms; ++i ) {
1862       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
1863       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
1864           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
1865       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
1866           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
1867       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
1868           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
1869     }
1870   }
1871 }
1872
1873 void HomePatch::addVelocityToPosition(
1874     FullAtom       * __restrict atom_arr,
1875     const BigReal    dt,
1876     int              num_atoms
1877     ) {
1878   SimParameters *simParams = Node::Object()->simParameters;
1879   if ( simParams->fixedAtomsOn ) {
1880     for ( int i = 0; i < num_atoms; ++i ) {
1881       if ( ! atom_arr[i].atomFixed ) {
1882         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
1883         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
1884         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
1885       }
1886     }
1887   } else {
1888     for ( int i = 0; i < num_atoms; ++i ) {
1889       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
1890       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
1891       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
1892     }
1893   }
1894 }
1895
1896 int HomePatch::hardWallDrude(const BigReal timestep, Tensor *virial,
1897     SubmitReduction *ppreduction)
1898 {
1899   Molecule *mol = Node::Object()->molecule;
1900   SimParameters *simParams = Node::Object()->simParameters;
1901   const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
1902   const int fixedAtomsOn = simParams->fixedAtomsOn;
1903   const BigReal dt = timestep / TIMEFACTOR;
1904   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
1905   int i, ia, ib, j;
1906   int dieOnError = simParams->rigidDie;
1907   Tensor wc;  // constraint virial
1908   BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
1909   int nslabs;
1910
1911   // start data for hard wall boundary between drude and its host atom
1912   // static int Count=0;
1913   int Idx;
1914   double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
1915   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;
1916   double dot_v_r_1, dot_v_r_2;
1917   double vb_cm, dr_a, dr_b;
1918   // end data for hard wall boundary between drude and its host atom
1919
1920   // start calculation of hard wall boundary between drude and its host atom
1921   if (simParams->drudeHardWallOn) {
1922     if (ppreduction) {
1923       nslabs = simParams->pressureProfileSlabs;
1924       idz = nslabs/lattice.c().z;
1925       zmin = lattice.origin().z - 0.5*lattice.c().z;
1926     }
1927
1928     r_wall = simParams->drudeBondLen;
1929     r_wall_SQ = r_wall*r_wall;
1930     // Count++;
1931     for (i=1; i<numAtoms; i++)  {
1932       if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) { // drude particle
1933         ia = i-1;
1934         ib = i;
1935
1936         v_ab = atom[ib].position - atom[ia].position;
1937         rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
1938
1939         if (rab_SQ > r_wall_SQ) {  // to impose the hard wall constraint
1940           rab = sqrt(rab_SQ);
1941           if ( (rab > (2.0*r_wall)) && dieOnError ) {  // unexpected situation
1942             iout << iERROR << "HardWallDrude> "
1943               << "The drude is too far away from atom "
1944               << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
1945             return -1;  // triggers early exit
1946           }
1947
1948           v_ab.x /= rab;
1949           v_ab.y /= rab;
1950           v_ab.z /= rab;
1951
1952           if ( fixedAtomsOn && atom[ia].atomFixed ) {  // the heavy atom is fixed
1953             if (atom[ib].atomFixed) {  // the drude is fixed too
1954               continue;
1955             }
1956             else {  // only the heavy atom is fixed
1957               dot_v_r_2 = atom[ib].velocity.x*v_ab.x
1958                 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
1959               vb_2 = dot_v_r_2 * v_ab;
1960               vp_2 = atom[ib].velocity - vb_2;
1961
1962               dr = rab - r_wall;
1963               if(dot_v_r_2 == 0.0) {
1964                 delta_T = maxtime;
1965               }
1966               else {
1967                 delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
1968                 if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
1969               }
1970
1971               dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
1972
1973               vb_2 = dot_v_r_2 * v_ab;
1974
1975               new_vel_a = atom[ia].velocity;
1976               new_vel_b = vp_2 + vb_2;
1977
1978               dr_b = -dr + delta_T*dot_v_r_2;  // L = L_0 + dT *v_new, v was flipped
1979
1980               new_pos_a = atom[ia].position;
1981               new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
1982             }
1983           }
1984           else {
1985             mass_a = atom[ia].mass;
1986             mass_b = atom[ib].mass;
1987             mass_sum = mass_a+mass_b;
1988
1989             dot_v_r_1 = atom[ia].velocity.x*v_ab.x
1990               + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
1991             vb_1 = dot_v_r_1 * v_ab;
1992             vp_1 = atom[ia].velocity - vb_1;
1993
1994             dot_v_r_2 = atom[ib].velocity.x*v_ab.x
1995               + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
1996             vb_2 = dot_v_r_2 * v_ab;
1997             vp_2 = atom[ib].velocity - vb_2;
1998
1999             vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
2000
2001             dot_v_r_1 -= vb_cm;
2002             dot_v_r_2 -= vb_cm;
2003
2004             dr = rab - r_wall;
2005
2006             if(dot_v_r_2 == dot_v_r_1) {
2007                 delta_T = maxtime;
2008             }
2009             else {
2010               delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);  // the time since the collision occurs
2011               if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
2012             }
2013             
2014             // the relative velocity between ia and ib. Drawn according to T_Drude
2015             v_Bond = sqrt(kbt/mass_b);
2016
2017             // reflect the velocity along bond vector and scale down
2018             dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
2019             dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
2020
2021             dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
2022             dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
2023
2024             new_pos_a = atom[ia].position + dr_a*v_ab;  // correct the position
2025             new_pos_b = atom[ib].position + dr_b*v_ab;
2026             // atom[ia].position += (dr_a*v_ab);  // correct the position
2027             // atom[ib].position += (dr_b*v_ab);
2028             
2029             dot_v_r_1 += vb_cm;
2030             dot_v_r_2 += vb_cm;
2031
2032             vb_1 = dot_v_r_1 * v_ab;
2033             vb_2 = dot_v_r_2 * v_ab;
2034
2035             new_vel_a = vp_1 + vb_1;
2036             new_vel_b = vp_2 + vb_2;
2037           }
2038
2039           int ppoffset, partition;
2040           if ( invdt == 0 ) {
2041             atom[ia].position = new_pos_a;
2042             atom[ib].position = new_pos_b;
2043           }
2044           else if ( virial == 0 ) {
2045             atom[ia].velocity = new_vel_a;
2046             atom[ib].velocity = new_vel_b;
2047           }
2048           else {
2049             for ( j = 0; j < 2; j++ ) {
2050               if (j==0) {  // atom ia, heavy atom
2051                 Idx = ia;
2052                 new_pos = &new_pos_a;
2053                 new_vel = &new_vel_a;
2054               }
2055               else if (j==1) {  // atom ib, drude
2056                 Idx = ib;
2057                 new_pos = &new_pos_b;
2058                 new_vel = &new_vel_b;
2059               }
2060               Force df = (*new_vel - atom[Idx].velocity) *
2061                 ( atom[Idx].mass * invdt );
2062               Tensor vir = outer(df, atom[Idx].position);
2063               wc += vir;
2064               atom[Idx].velocity = *new_vel;
2065               atom[Idx].position = *new_pos;
2066
2067               if (ppreduction) {
2068                 if (!j) {
2069                   BigReal z = new_pos->z;
2070                   int partition = atom[Idx].partition;
2071                   int slab = (int)floor((z-zmin)*idz);
2072                   if (slab < 0) slab += nslabs;
2073                   else if (slab >= nslabs) slab -= nslabs;
2074                   ppoffset = 3*(slab + nslabs*partition);
2075                 }
2076                 ppreduction->item(ppoffset  ) += vir.xx;
2077                 ppreduction->item(ppoffset+1) += vir.yy;
2078                 ppreduction->item(ppoffset+2) += vir.zz;
2079               }
2080
2081             }
2082           }
2083         }                               
2084       }
2085     }
2086
2087     // if ( (Count>10000) && (Count%10==0) ) {
2088     //   v_ab = atom[1].position - atom[0].position;
2089     //   rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
2090     //   iout << "DBG_R: " << Count << "  " << sqrt(rab_SQ) << "\n" << endi;
2091     // }
2092
2093   }
2094
2095   // end calculation of hard wall boundary between drude and its host atom
2096
2097   if ( dt && virial ) *virial += wc;
2098
2099   return 0;
2100 }
2101
2102 void HomePatch::buildRattleList() {
2103
2104   SimParameters *simParams = Node::Object()->simParameters;
2105   const int fixedAtomsOn = simParams->fixedAtomsOn;
2106   const int useSettle = simParams->useSettle;
2107
2108   // Re-size to containt numAtoms elements
2109   velNew.resize(numAtoms);
2110   posNew.resize(numAtoms);
2111
2112   // Size of a hydrogen group for water
2113   int wathgsize = 3;
2114   int watmodel = simParams->watmodel;
2115   if (watmodel == WAT_TIP4) wathgsize = 4;
2116   else if (watmodel == WAT_SWM4) wathgsize = 5;
2117
2118   // Initialize the settle algorithm with water parameters
2119   // settle1() assumes all waters are identical,
2120   // and will generate bad results if they are not.
2121   // XXX this will move to Molecule::build_atom_status when that 
2122   // version is debugged
2123   if ( ! settle_initialized ) {
2124     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2125       // find a water
2126       if (atom[ig].rigidBondLength > 0) {
2127         int oatm;
2128         if (simParams->watmodel == WAT_SWM4) {
2129           oatm = ig+3;  // skip over Drude and Lonepair
2130           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
2131           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
2132         }
2133         else {
2134           oatm = ig+1;
2135           // Avoid using the Om site to set this by mistake
2136           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2137             oatm += 1;
2138           }
2139         }
2140
2141         // initialize settle water parameters
2142         settle1init(atom[ig].mass, atom[oatm].mass, 
2143                     atom[ig].rigidBondLength,
2144                     atom[oatm].rigidBondLength,
2145                     settle_mOrmT, settle_mHrmT, settle_ra,
2146                     settle_rb, settle_rc, settle_rra);
2147         settle_initialized = 1;
2148         break; // done with init
2149       }
2150     }
2151   }
2152   
2153   Vector ref[10];
2154   BigReal rmass[10];
2155   BigReal dsq[10];
2156   int fixed[10];
2157   int ial[10];
2158   int ibl[10];
2159
2160   int numSettle = 0;
2161   int numRattle = 0;
2162   int posRattleParam = 0;
2163
2164   settleList.clear();
2165   rattleList.clear();
2166   noconstList.clear();
2167   rattleParam.clear();
2168
2169   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2170     int hgs = atom[ig].hydrogenGroupSize;
2171     if ( hgs == 1 ) {
2172       // only one atom in group
2173       noconstList.push_back(ig);
2174       continue;
2175     }
2176     int anyfixed = 0;
2177     for (int i = 0; i < hgs; ++i ) {
2178       ref[i] = atom[ig+i].position;
2179       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2180       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2181       if ( fixed[i] ) {
2182         anyfixed = 1;
2183         rmass[i] = 0.;
2184       }
2185     }
2186     int icnt = 0;
2187     BigReal tmp = atom[ig].rigidBondLength;
2188     if (tmp > 0.0) {  // for water
2189       if (hgs != wathgsize) {
2190         char errmsg[256];
2191         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
2192                          "but the specified water model requires %d atoms.\n",
2193                          atom[ig].id+1, hgs, wathgsize);
2194         NAMD_die(errmsg);
2195       }
2196       // Use SETTLE for water unless some of the water atoms are fixed,
2197       if (useSettle && !anyfixed) {
2198         // Store to Settle -list
2199         settleList.push_back(ig);
2200         continue;
2201       }
2202       if ( !(fixed[1] && fixed[2]) ) {
2203         dsq[icnt] = tmp * tmp;
2204         ial[icnt] = 1;
2205         ibl[icnt] = 2;
2206         ++icnt;
2207       }
2208     } // if (tmp > 0.0)
2209     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
2210       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2211         if ( !(fixed[0] && fixed[i]) ) {
2212           dsq[icnt] = tmp * tmp;
2213           ial[icnt] = 0;
2214           ibl[icnt] = i;
2215           ++icnt;
2216         }
2217       }
2218     }
2219     if ( icnt == 0 ) {
2220       // no constraints
2221       noconstList.push_back(ig);
2222       continue;  
2223     }
2224     // Store to Rattle -list
2225     RattleList rattleListElem;
2226     rattleListElem.ig  = ig;
2227     rattleListElem.icnt = icnt;
2228     rattleList.push_back(rattleListElem);
2229     for (int i = 0; i < icnt; ++i ) {
2230       int a = ial[i];
2231       int b = ibl[i];
2232       RattleParam rattleParamElem;
2233       rattleParamElem.ia = a;
2234       rattleParamElem.ib = b;
2235       rattleParamElem.dsq = dsq[i];
2236       rattleParamElem.rma = rmass[a];
2237       rattleParamElem.rmb = rmass[b];
2238       rattleParam.push_back(rattleParamElem);
2239     }
2240   }
2241
2242 }
2243
2244 void HomePatch::addRattleForce(const BigReal invdt, Tensor& wc) {
2245   for (int ig = 0; ig < numAtoms; ++ig ) {
2246     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
2247     Tensor vir = outer(df, atom[ig].position);
2248     wc += vir;
2249     f[Results::normal][ig] += df;
2250     atom[ig].velocity = velNew[ig];
2251   }
2252 }
2253
2254 int HomePatch::rattle1(const BigReal timestep, Tensor *virial, 
2255   SubmitReduction *ppreduction) {
2256
2257   SimParameters *simParams = Node::Object()->simParameters;
2258   if (simParams->watmodel != WAT_TIP3 || ppreduction) {
2259     // Call old rattle1 -method instead
2260     return rattle1old(timestep, virial, ppreduction);
2261   }
2262
2263   if (!rattleListValid) {
2264     buildRattleList();
2265     rattleListValid = true;
2266   }
2267
2268   const int fixedAtomsOn = simParams->fixedAtomsOn;
2269   const int useSettle = simParams->useSettle;
2270   const BigReal dt = timestep / TIMEFACTOR;
2271   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2272   const BigReal tol2 = 2.0 * simParams->rigidTol;
2273   int maxiter = simParams->rigidIter;
2274   int dieOnError = simParams->rigidDie;
2275
2276   Vector ref[10];  // reference position
2277   Vector pos[10];  // new position
2278   Vector vel[10];  // new velocity
2279
2280   // Manual un-roll
2281   int n = (settleList.size()/2)*2;
2282   for (int j=0;j < n;j+=2) {
2283     int ig;
2284     ig = settleList[j];
2285     for (int i = 0; i < 3; ++i ) {
2286       ref[i] = atom[ig+i].position;
2287       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2288     }
2289     ig = settleList[j+1];
2290     for (int i = 0; i < 3; ++i ) {
2291       ref[i+3] = atom[ig+i].position;
2292       pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
2293     }
2294     settle1_SIMD<2>(ref, pos,
2295       settle_mOrmT, settle_mHrmT, settle_ra,
2296       settle_rb, settle_rc, settle_rra);
2297
2298     ig = settleList[j];
2299     for (int i = 0; i < 3; ++i ) {
2300       velNew[ig+i] = (pos[i] - ref[i])*invdt;
2301       posNew[ig+i] = pos[i];
2302     }
2303     ig = settleList[j+1];
2304     for (int i = 0; i < 3; ++i ) {
2305       velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
2306       posNew[ig+i] = pos[i+3];
2307     }
2308
2309   }
2310
2311   if (settleList.size() % 2) {
2312     int ig = settleList[settleList.size()-1];
2313     for (int i = 0; i < 3; ++i ) {
2314       ref[i] = atom[ig+i].position;
2315       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
2316     }
2317     settle1_SIMD<1>(ref, pos,
2318             settle_mOrmT, settle_mHrmT, settle_ra,
2319             settle_rb, settle_rc, settle_rra);        
2320     for (int i = 0; i < 3; ++i ) {
2321       velNew[ig+i] = (pos[i] - ref[i])*invdt;
2322       posNew[ig+i] = pos[i];
2323     }
2324   }
2325
2326   int posParam = 0;
2327   for (int j=0;j < rattleList.size();++j) {
2328
2329     BigReal refx[10];
2330     BigReal refy[10];
2331     BigReal refz[10];
2332
2333     BigReal posx[10];
2334     BigReal posy[10];
2335     BigReal posz[10];
2336
2337     int ig = rattleList[j].ig;
2338     int icnt = rattleList[j].icnt;
2339     int hgs = atom[ig].hydrogenGroupSize;
2340     for (int i = 0; i < hgs; ++i ) {
2341       ref[i] = atom[ig+i].position;
2342       pos[i] = atom[ig+i].position;
2343       if (!(fixedAtomsOn && atom[ig+i].atomFixed))
2344         pos[i] += atom[ig+i].velocity * dt;
2345       refx[i] = ref[i].x;
2346       refy[i] = ref[i].y;
2347       refz[i] = ref[i].z;
2348       posx[i] = pos[i].x;
2349       posy[i] = pos[i].y;
2350       posz[i] = pos[i].z;
2351     }
2352
2353
2354     bool done;
2355     bool consFailure;
2356     if (icnt == 1) {
2357       rattlePair<1>(&rattleParam[posParam],
2358         refx, refy, refz,
2359         posx, posy, posz,
2360         consFailure);
2361       done = true;
2362     } else {
2363       rattleN(icnt, &rattleParam[posParam],
2364         refx, refy, refz,
2365         posx, posy, posz,
2366         tol2, maxiter,
2367         done, consFailure);
2368     }
2369
2370     // Advance position in rattleParam
2371     posParam += icnt;
2372
2373     for (int i = 0; i < hgs; ++i ) {
2374       pos[i].x = posx[i];
2375       pos[i].y = posy[i];
2376       pos[i].z = posz[i];
2377     }
2378
2379     for (int i = 0; i < hgs; ++i ) {
2380       velNew[ig+i] = (pos[i] - ref[i])*invdt;
2381       posNew[ig+i] = pos[i];
2382     }
2383
2384     if ( consFailure ) {
2385       if ( dieOnError ) {
2386         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
2387         << (atom[ig].id + 1) << "!\n" << endi;
2388         return -1;  // triggers early exit
2389       } else {
2390         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
2391         << (atom[ig].id + 1) << "!\n" << endi;
2392       }
2393     } else if ( ! done ) {
2394       if ( dieOnError ) {
2395         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
2396         << (atom[ig].id + 1) << "!\n" << endi;
2397         return -1;  // triggers early exit
2398       } else {
2399         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
2400         << (atom[ig].id + 1) << "!\n" << endi;
2401       }
2402     }
2403
2404   }
2405
2406   // Finally, we have to go through atoms that are not involved in rattle just so that we have
2407   // their positions and velocities up-to-date in posNew and velNew
2408   for (int j=0;j < noconstList.size();++j) {
2409     int ig = noconstList[j];
2410     int hgs = atom[ig].hydrogenGroupSize;
2411     for (int i = 0; i < hgs; ++i ) {
2412       velNew[ig+i] = atom[ig+i].velocity;
2413       posNew[ig+i] = atom[ig+i].position;
2414     }
2415   }
2416
2417   if ( invdt == 0 ) {
2418     for (int ig = 0; ig < numAtoms; ++ig )
2419       atom[ig].position = posNew[ig];
2420   } else if ( virial == 0 ) {
2421     for (int ig = 0; ig < numAtoms; ++ig )
2422       atom[ig].velocity = velNew[ig];
2423   } else {
2424     Tensor wc;  // constraint virial
2425     addRattleForce(invdt, wc);
2426     *virial += wc;
2427   }
2428
2429   return 0;
2430 }
2431
2432 //  RATTLE algorithm from Allen & Tildesley
2433 int HomePatch::rattle1old(const BigReal timestep, Tensor *virial, 
2434     SubmitReduction *ppreduction)
2435 {
2436   Molecule *mol = Node::Object()->molecule;
2437   SimParameters *simParams = Node::Object()->simParameters;
2438   const int fixedAtomsOn = simParams->fixedAtomsOn;
2439   const int useSettle = simParams->useSettle;
2440   const BigReal dt = timestep / TIMEFACTOR;
2441   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
2442   BigReal tol2 = 2.0 * simParams->rigidTol;
2443   int maxiter = simParams->rigidIter;
2444   int dieOnError = simParams->rigidDie;
2445   int i, iter;
2446   BigReal dsq[10], tmp;
2447   int ial[10], ibl[10];
2448   Vector ref[10];  // reference position
2449   Vector refab[10];  // reference vector
2450   Vector pos[10];  // new position
2451   Vector vel[10];  // new velocity
2452   Vector netdp[10];  // total momentum change from constraint
2453   BigReal rmass[10];  // 1 / mass
2454   int fixed[10];  // is atom fixed?
2455   Tensor wc;  // constraint virial
2456   BigReal idz, zmin;
2457   int nslabs;
2458
2459   // Initialize the settle algorithm with water parameters
2460   // settle1() assumes all waters are identical,
2461   // and will generate bad results if they are not.
2462   // XXX this will move to Molecule::build_atom_status when that 
2463   // version is debugged
2464   if ( ! settle_initialized ) {
2465     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2466       // find a water
2467       if (atom[ig].rigidBondLength > 0) {
2468         int oatm;
2469         if (simParams->watmodel == WAT_SWM4) {
2470           oatm = ig+3;  // skip over Drude and Lonepair
2471           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
2472           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
2473         }
2474         else {
2475           oatm = ig+1;
2476           // Avoid using the Om site to set this by mistake
2477           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
2478             oatm += 1;
2479           }
2480         }
2481
2482         // initialize settle water parameters
2483         settle1init(atom[ig].mass, atom[oatm].mass, 
2484                     atom[ig].rigidBondLength,
2485                     atom[oatm].rigidBondLength,
2486                     settle_mOrmT, settle_mHrmT, settle_ra,
2487                     settle_rb, settle_rc, settle_rra);
2488         settle_initialized = 1;
2489         break; // done with init
2490       }
2491     }
2492   }
2493
2494   if (ppreduction) {
2495     nslabs = simParams->pressureProfileSlabs;
2496     idz = nslabs/lattice.c().z;
2497     zmin = lattice.origin().z - 0.5*lattice.c().z;
2498   }
2499
2500   // Size of a hydrogen group for water
2501   int wathgsize = 3;
2502   int watmodel = simParams->watmodel;
2503   if (watmodel == WAT_TIP4) wathgsize = 4;
2504   else if (watmodel == WAT_SWM4) wathgsize = 5;
2505   
2506   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2507     int hgs = atom[ig].hydrogenGroupSize;
2508     if ( hgs == 1 ) continue;  // only one atom in group
2509     // cache data in local arrays and integrate positions normally
2510     int anyfixed = 0;
2511     for ( i = 0; i < hgs; ++i ) {
2512       ref[i] = atom[ig+i].position;
2513       pos[i] = atom[ig+i].position;
2514       vel[i] = atom[ig+i].velocity;
2515       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
2516       //printf("rmass of %i is %f\n", ig+i, rmass[i]);
2517       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2518       //printf("fixed status of %i is %i\n", i, fixed[i]);
2519       // undo addVelocityToPosition to get proper reference coordinates
2520       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
2521     }
2522     int icnt = 0;
2523     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
2524       if (hgs != wathgsize) {
2525         char errmsg[256];
2526         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
2527                          "but the specified water model requires %d atoms.\n",
2528                          atom[ig].id+1, hgs, wathgsize);
2529         NAMD_die(errmsg);
2530       }
2531       // Use SETTLE for water unless some of the water atoms are fixed,
2532       if (useSettle && !anyfixed) {
2533         if (simParams->watmodel == WAT_SWM4) {
2534           // SWM4 ordering:  O D LP H1 H2
2535           // do swap(O,LP) and call settle with subarray O H1 H2
2536           // swap back after we return
2537           Vector lp_ref = ref[2];
2538           Vector lp_pos = pos[2];
2539           Vector lp_vel = vel[2];
2540           ref[2] = ref[0];
2541           pos[2] = pos[0];
2542           vel[2] = vel[0];
2543           settle1(ref+2, pos+2, vel+2, invdt,
2544                   settle_mOrmT, settle_mHrmT, settle_ra,
2545                   settle_rb, settle_rc, settle_rra);
2546           ref[0] = ref[2];
2547           pos[0] = pos[2];
2548           vel[0] = vel[2];
2549           ref[2] = lp_ref;
2550           pos[2] = lp_pos;
2551           vel[2] = lp_vel;
2552           // determine for LP updated pos and vel
2553           swm4_omrepos(ref, pos, vel, invdt);
2554         }
2555         else {
2556             settle1(ref, pos, vel, invdt,
2557                     settle_mOrmT, settle_mHrmT, settle_ra,
2558                     settle_rb, settle_rc, settle_rra);            
2559           if (simParams->watmodel == WAT_TIP4) {
2560             tip4_omrepos(ref, pos, vel, invdt);
2561           }
2562         }
2563
2564         // which slab the hydrogen group will belong to
2565         // for pprofile calculations.
2566         int ppoffset, partition;
2567         if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
2568           atom[ig+i].position = pos[i];
2569         } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
2570           atom[ig+i].velocity = vel[i];
2571         } else for ( i = 0; i < wathgsize; ++i ) {
2572           Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
2573           Tensor vir = outer(df, ref[i]);
2574           wc += vir;
2575           f[Results::normal][ig+i] += df;
2576           atom[ig+i].velocity = vel[i];
2577           if (ppreduction) {
2578             // put all the atoms from a water in the same slab.  Atom 0
2579             // should be the parent atom.
2580             if (!i) {
2581               BigReal z = pos[i].z;
2582               partition = atom[ig].partition;
2583               int slab = (int)floor((z-zmin)*idz);
2584               if (slab < 0) slab += nslabs;
2585               else if (slab >= nslabs) slab -= nslabs;
2586               ppoffset = 3*(slab + nslabs*partition);
2587             }
2588             ppreduction->item(ppoffset  ) += vir.xx;
2589             ppreduction->item(ppoffset+1) += vir.yy;
2590             ppreduction->item(ppoffset+2) += vir.zz;
2591           }
2592         }
2593         continue;
2594       }
2595       if ( !(fixed[1] && fixed[2]) ) {
2596         dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
2597       }
2598     }
2599     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
2600       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2601         if ( !(fixed[0] && fixed[i]) ) {
2602           dsq[icnt] = tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
2603         }
2604       }
2605     }
2606     if ( icnt == 0 ) continue;  // no constraints
2607     for ( i = 0; i < icnt; ++i ) {
2608       refab[i] = ref[ial[i]] - ref[ibl[i]];
2609     }
2610     for ( i = 0; i < hgs; ++i ) {
2611       netdp[i] = 0.;
2612     }
2613     int done;
2614     int consFailure;
2615     for ( iter = 0; iter < maxiter; ++iter ) {
2616 //if (iter > 0) CkPrintf("iteration %d\n", iter);
2617       done = 1;
2618       consFailure = 0;
2619       for ( i = 0; i < icnt; ++i ) {
2620         int a = ial[i];  int b = ibl[i];
2621         Vector pab = pos[a] - pos[b];
2622               BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
2623         BigReal rabsq = dsq[i];
2624         BigReal diffsq = rabsq - pabsq;
2625         if ( fabs(diffsq) > (rabsq * tol2) ) {
2626           Vector &rab = refab[i];
2627           BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
2628           if ( rpab < ( rabsq * 1.0e-6 ) ) {
2629             done = 0;
2630             consFailure = 1;
2631             continue;
2632           }
2633           BigReal rma = rmass[a];
2634           BigReal rmb = rmass[b];
2635           BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
2636           Vector dp = rab * gab;
2637           pos[a] += rma * dp;
2638           pos[b] -= rmb * dp;
2639           if ( invdt != 0. ) {
2640             dp *= invdt;
2641             netdp[a] += dp;
2642             netdp[b] -= dp;
2643           }
2644           done = 0;
2645         }
2646       }
2647       if ( done ) break;
2648     }
2649
2650     if ( consFailure ) {
2651       if ( dieOnError ) {
2652         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
2653         << (atom[ig].id + 1) << "!\n" << endi;
2654         return -1;  // triggers early exit
2655       } else {
2656         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
2657         << (atom[ig].id + 1) << "!\n" << endi;
2658       }
2659     } else if ( ! done ) {
2660       if ( dieOnError ) {
2661         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
2662         << (atom[ig].id + 1) << "!\n" << endi;
2663         return -1;  // triggers early exit
2664       } else {
2665         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
2666         << (atom[ig].id + 1) << "!\n" << endi;
2667       }
2668     }
2669
2670     // store data back to patch
2671     int ppoffset, partition;
2672     if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
2673       atom[ig+i].position = pos[i];
2674     } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
2675       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2676     } else for ( i = 0; i < hgs; ++i ) {
2677       Force df = netdp[i] * invdt;
2678       Tensor vir = outer(df, ref[i]);
2679       wc += vir;
2680       f[Results::normal][ig+i] += df;
2681       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
2682       if (ppreduction) {
2683         if (!i) {
2684           BigReal z = pos[i].z;
2685           int partition = atom[ig].partition;
2686           int slab = (int)floor((z-zmin)*idz);
2687           if (slab < 0) slab += nslabs;
2688           else if (slab >= nslabs) slab -= nslabs;
2689           ppoffset = 3*(slab + nslabs*partition);
2690         }
2691         ppreduction->item(ppoffset  ) += vir.xx;
2692         ppreduction->item(ppoffset+1) += vir.yy;
2693         ppreduction->item(ppoffset+2) += vir.zz;
2694       }
2695     }
2696   }
2697   if ( dt && virial ) *virial += wc;
2698
2699   return 0;
2700 }
2701
2702 //  RATTLE algorithm from Allen & Tildesley
2703 void HomePatch::rattle2(const BigReal timestep, Tensor *virial)
2704 {
2705   Molecule *mol = Node::Object()->molecule;
2706   SimParameters *simParams = Node::Object()->simParameters;
2707   const int fixedAtomsOn = simParams->fixedAtomsOn;
2708   const int useSettle = simParams->useSettle;
2709   const BigReal dt = timestep / TIMEFACTOR;
2710   Tensor wc;  // constraint virial
2711   BigReal tol = simParams->rigidTol;
2712   int maxiter = simParams->rigidIter;
2713   int dieOnError = simParams->rigidDie;
2714   int i, iter;
2715   BigReal dsqi[10], tmp;
2716   int ial[10], ibl[10];
2717   Vector ref[10];  // reference position
2718   Vector refab[10];  // reference vector
2719   Vector vel[10];  // new velocity
2720   BigReal rmass[10];  // 1 / mass
2721   BigReal redmass[10];  // reduced mass
2722   int fixed[10];  // is atom fixed?
2723   
2724   // Size of a hydrogen group for water
2725   int wathgsize = 3;
2726   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
2727   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
2728
2729   //  CkPrintf("In rattle2!\n");
2730   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2731     //    CkPrintf("ig=%d\n",ig);
2732     int hgs = atom[ig].hydrogenGroupSize;
2733     if ( hgs == 1 ) continue;  // only one atom in group
2734     // cache data in local arrays and integrate positions normally
2735     int anyfixed = 0;
2736     for ( i = 0; i < hgs; ++i ) {
2737       ref[i] = atom[ig+i].position;
2738       vel[i] = atom[ig+i].velocity;
2739       rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
2740       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2741       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
2742     }
2743     int icnt = 0;
2744     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
2745       if (hgs != wathgsize) {
2746         NAMD_bug("Hydrogen group error caught in rattle2().");
2747       }
2748       // Use SETTLE for water unless some of the water atoms are fixed,
2749       if (useSettle && !anyfixed) {
2750         if (simParams->watmodel == WAT_SWM4) {
2751           // SWM4 ordering:  O D LP H1 H2
2752           // do swap(O,LP) and call settle with subarray O H1 H2
2753           // swap back after we return
2754           Vector lp_ref = ref[2];
2755           // Vector lp_vel = vel[2];
2756           ref[2] = ref[0];
2757           vel[2] = vel[0];
2758           settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
2759           ref[0] = ref[2];
2760           vel[0] = vel[2];
2761           ref[2] = lp_ref;
2762           // vel[2] = vel[0];  // set LP vel to O vel
2763         }
2764         else {
2765           settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
2766           if (simParams->watmodel == WAT_TIP4) {
2767             vel[3] = vel[0];
2768           }
2769         }
2770         for (i=0; i<hgs; i++) {
2771           atom[ig+i].velocity = vel[i];
2772         }
2773         continue;
2774       }
2775       if ( !(fixed[1] && fixed[2]) ) {
2776         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
2777         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
2778       }
2779     }
2780     //    CkPrintf("Loop 2\n");
2781     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
2782       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2783         if ( !(fixed[0] && fixed[i]) ) {
2784           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
2785           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
2786           ibl[icnt] = i;  ++icnt;
2787         }
2788       }
2789     }
2790     if ( icnt == 0 ) continue;  // no constraints
2791     //    CkPrintf("Loop 3\n");
2792     for ( i = 0; i < icnt; ++i ) {
2793       refab[i] = ref[ial[i]] - ref[ibl[i]];
2794     }
2795     //    CkPrintf("Loop 4\n");
2796     int done;
2797     for ( iter = 0; iter < maxiter; ++iter ) {
2798       done = 1;
2799       for ( i = 0; i < icnt; ++i ) {
2800         int a = ial[i];  int b = ibl[i];
2801         Vector vab = vel[a] - vel[b];
2802         Vector &rab = refab[i];
2803         BigReal rabsqi = dsqi[i];
2804         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
2805         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
2806           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
2807           wc += outer(dp,rab);
2808           vel[a] += rmass[a] * dp;
2809           vel[b] -= rmass[b] * dp;
2810           done = 0;
2811         }
2812       }
2813       if ( done ) break;
2814       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
2815     }
2816     if ( ! done ) {
2817       if ( dieOnError ) {
2818         NAMD_die("Exceeded maximum number of iterations in rattle2().");
2819       } else {
2820         iout << iWARN <<
2821           "Exceeded maximum number of iterations in rattle2().\n" << endi;
2822       }
2823     }
2824     // store data back to patch
2825     for ( i = 0; i < hgs; ++i ) {
2826       atom[ig+i].velocity = vel[i];
2827     }
2828   }
2829   //  CkPrintf("Leaving rattle2!\n");
2830   // check that there isn't a constant needed here!
2831   *virial += wc / ( 0.5 * dt );
2832
2833 }
2834
2835
2836 //  Adjust gradients for minimizer
2837 void HomePatch::minimize_rattle2(const BigReal timestep, Tensor *virial, bool forces)
2838 {
2839   Molecule *mol = Node::Object()->molecule;
2840   SimParameters *simParams = Node::Object()->simParameters;
2841   Force *f1 = f[Results::normal].begin();
2842   const int fixedAtomsOn = simParams->fixedAtomsOn;
2843   const int useSettle = simParams->useSettle;
2844   const BigReal dt = timestep / TIMEFACTOR;
2845   Tensor wc;  // constraint virial
2846   BigReal tol = simParams->rigidTol;
2847   int maxiter = simParams->rigidIter;
2848   int dieOnError = simParams->rigidDie;
2849   int i, iter;
2850   BigReal dsqi[10], tmp;
2851   int ial[10], ibl[10];
2852   Vector ref[10];  // reference position
2853   Vector refab[10];  // reference vector
2854   Vector vel[10];  // new velocity
2855   BigReal rmass[10];  // 1 / mass
2856   BigReal redmass[10];  // reduced mass
2857   int fixed[10];  // is atom fixed?
2858   
2859   // Size of a hydrogen group for water
2860   int wathgsize = 3;
2861   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
2862   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
2863
2864   //  CkPrintf("In rattle2!\n");
2865   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
2866     //    CkPrintf("ig=%d\n",ig);
2867     int hgs = atom[ig].hydrogenGroupSize;
2868     if ( hgs == 1 ) continue;  // only one atom in group
2869     // cache data in local arrays and integrate positions normally
2870     int anyfixed = 0;
2871     for ( i = 0; i < hgs; ++i ) {
2872       ref[i] = atom[ig+i].position;
2873       vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
2874       rmass[i] = 1.0;
2875       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
2876       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
2877     }
2878     int icnt = 0;
2879     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
2880       if (hgs != wathgsize) {
2881         NAMD_bug("Hydrogen group error caught in rattle2().");
2882       }
2883       // Use SETTLE for water unless some of the water atoms are fixed,
2884       if (useSettle && !anyfixed) {
2885         if (simParams->watmodel == WAT_SWM4) {
2886           // SWM4 ordering:  O D LP H1 H2
2887           // do swap(O,LP) and call settle with subarray O H1 H2
2888           // swap back after we return
2889           Vector lp_ref = ref[2];
2890           // Vector lp_vel = vel[2];
2891           ref[2] = ref[0];
2892           vel[2] = vel[0];
2893           settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
2894           ref[0] = ref[2];
2895           vel[0] = vel[2];
2896           ref[2] = lp_ref;
2897           // vel[2] = vel[0];  // set LP vel to O vel
2898         }
2899         else {
2900           settle2(1.0, 1.0, ref, vel, dt, virial);
2901           if (simParams->watmodel == WAT_TIP4) {
2902             vel[3] = vel[0];
2903           }
2904         }
2905         for (i=0; i<hgs; i++) {
2906           ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
2907         }
2908         continue;
2909       }
2910       if ( !(fixed[1] && fixed[2]) ) {
2911         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
2912         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
2913       }
2914     }
2915     //    CkPrintf("Loop 2\n");
2916     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
2917       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
2918         if ( !(fixed[0] && fixed[i]) ) {
2919           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
2920           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
2921           ibl[icnt] = i;  ++icnt;
2922         }
2923       }
2924     }
2925     if ( icnt == 0 ) continue;  // no constraints
2926     //    CkPrintf("Loop 3\n");
2927     for ( i = 0; i < icnt; ++i ) {
2928       refab[i] = ref[ial[i]] - ref[ibl[i]];
2929     }
2930     //    CkPrintf("Loop 4\n");
2931     int done;
2932     for ( iter = 0; iter < maxiter; ++iter ) {
2933       done = 1;
2934       for ( i = 0; i < icnt; ++i ) {
2935         int a = ial[i];  int b = ibl[i];
2936         Vector vab = vel[a] - vel[b];
2937         Vector &rab = refab[i];
2938         BigReal rabsqi = dsqi[i];
2939         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
2940         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
2941           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
2942           wc += outer(dp,rab);
2943           vel[a] += rmass[a] * dp;
2944           vel[b] -= rmass[b] * dp;
2945           done = 0;
2946         }
2947       }
2948       if ( done ) break;
2949       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
2950     }
2951     if ( ! done ) {
2952       if ( dieOnError ) {
2953         NAMD_die("Exceeded maximum number of iterations in rattle2().");
2954       } else {
2955         iout << iWARN <<
2956           "Exceeded maximum number of iterations in rattle2().\n" << endi;
2957       }
2958     }
2959     // store data back to patch
2960     for ( i = 0; i < hgs; ++i ) {
2961       ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
2962     }
2963   }
2964   //  CkPrintf("Leaving rattle2!\n");
2965   // check that there isn't a constant needed here!
2966   *virial += wc / ( 0.5 * dt );
2967
2968 }
2969
2970
2971 // BEGIN LA
2972 void HomePatch::loweAndersenVelocities()
2973 {
2974     DebugM(2, "loweAndersenVelocities\n");
2975     Molecule *mol = Node::Object()->molecule;
2976     SimParameters *simParams = Node::Object()->simParameters;
2977     v.resize(numAtoms);
2978     for (int i = 0; i < numAtoms; ++i) {
2979         //v[i] = p[i];
2980         // co-opt CompAtom structure to pass velocity and mass information
2981         v[i].position = atom[i].velocity;
2982         v[i].charge = atom[i].mass;
2983     }
2984     DebugM(2, "loweAndersenVelocities\n");
2985 }
2986
2987 void HomePatch::loweAndersenFinish()
2988 {
2989     DebugM(2, "loweAndersenFinish\n");
2990     v.resize(0);
2991 }
2992 // END LA
2993
2994 //LCPO
2995 void HomePatch::setLcpoType() {
2996   Molecule *mol = Node::Object()->molecule;
2997   const int *lcpoParamType = mol->getLcpoParamType();
2998
2999   lcpoType.resize(numAtoms);
3000   for (int i = 0; i < numAtoms; i++) {
3001     lcpoType[i] = lcpoParamType[pExt[i].id];
3002   }
3003 }
3004
3005 //set intrinsic radii of atom when doMigration
3006 void HomePatch::setGBISIntrinsicRadii() {
3007   intRad.resize(numAtoms*2);
3008   intRad.setall(0);
3009   Molecule *mol = Node::Object()->molecule;
3010   SimParameters *simParams = Node::Object()->simParameters;
3011   Real offset = simParams->coulomb_radius_offset;
3012   for (int i = 0; i < numAtoms; i++) {
3013     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
3014     Real screen = MassToScreen(atom[i].mass);//same
3015     intRad[2*i+0] = rad - offset;//r0
3016     intRad[2*i+1] = screen*(rad - offset);//s0
3017   }
3018 }
3019
3020 //compute born radius after phase 1, before phase 2
3021 void HomePatch::gbisComputeAfterP1() {
3022
3023   SimParameters *simParams = Node::Object()->simParameters;
3024   BigReal alphaMax = simParams->alpha_max;
3025   BigReal delta = simParams->gbis_delta;
3026   BigReal beta = simParams->gbis_beta;
3027   BigReal gamma = simParams->gbis_gamma;
3028   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
3029
3030   BigReal rhoi;
3031   BigReal rhoi0;
3032   //calculate bornRad from psiSum
3033   for (int i = 0; i < numAtoms; i++) {
3034     rhoi0 = intRad[2*i];
3035     rhoi = rhoi0+coulomb_radius_offset;
3036     psiFin[i] += psiSum[i];
3037     psiFin[i] *= rhoi0;
3038     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
3039     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
3040 #ifdef PRINT_COMP
3041     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
3042 #endif
3043   }
3044
3045   gbisP2Ready();
3046 }
3047
3048 //compute dHdrPrefix after phase 2, before phase 3
3049 void HomePatch::gbisComputeAfterP2() {
3050
3051   SimParameters *simParams = Node::Object()->simParameters;
3052   BigReal delta = simParams->gbis_delta;
3053   BigReal beta = simParams->gbis_beta;
3054   BigReal gamma = simParams->gbis_gamma;
3055   BigReal epsilon_s = simParams->solvent_dielectric;
3056   BigReal epsilon_p = simParams->dielectric;
3057   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
3058   BigReal epsilon_p_i = 1/simParams->dielectric;
3059   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
3060   BigReal kappa = simParams->kappa;
3061   BigReal fij, expkappa, Dij, dEdai, dedasum;
3062   BigReal rhoi, rhoi0, psii, nbetapsi;
3063   BigReal gammapsi2, tanhi, daidr;
3064   for (int i = 0; i < numAtoms; i++) {
3065     //add diagonal dEda term
3066     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
3067     fij = bornRad[i];//inf
3068     expkappa = exp(-kappa*fij);//0
3069     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
3070     //calculate dHij prefix
3071     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
3072                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
3073     dHdrPrefix[i] += dEdai;
3074     dedasum = dHdrPrefix[i];
3075
3076     rhoi0 = intRad[2*i];
3077     rhoi = rhoi0+coulomb_radius_offset;
3078     psii = psiFin[i];
3079     nbetapsi = -beta*psii;
3080     gammapsi2 = gamma*psii*psii;
3081     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
3082     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
3083            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
3084     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
3085 #ifdef PRINT_COMP
3086     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
3087 #endif
3088   }
3089   gbisP3Ready();
3090 }
3091
3092 //send born radius to proxies to begin phase 2
3093 void HomePatch::gbisP2Ready() {
3094   if (proxy.size() > 0) {
3095     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3096     for (int i = 0; i < proxy.size(); i++) {
3097       int node = proxy[i];
3098       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
3099       msg->patch = patchID;
3100       msg->origPe = CkMyPe();
3101       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
3102       msg->destPe = node;
3103       int seq = flags.sequence;
3104       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
3105       SET_PRIORITY(msg,seq,priority);
3106       cp[node].recvData(msg);
3107     }
3108   }
3109   Patch::gbisP2Ready();
3110 }
3111
3112 //send dHdrPrefix to proxies to begin phase 3
3113 void HomePatch::gbisP3Ready() {
3114   if (proxy.size() > 0) {
3115     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
3116     //only nonzero message should be sent for doFullElec
3117     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
3118     for (int i = 0; i < proxy.size(); i++) {
3119       int node = proxy[i];
3120       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
3121       msg->patch = patchID;
3122       msg->dHdrPrefixLen = msgAtoms;
3123       msg->origPe = CkMyPe();
3124       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
3125       msg->destPe = node;
3126       int seq = flags.sequence;
3127       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
3128       SET_PRIORITY(msg,seq,priority);
3129       cp[node].recvData(msg);
3130     }
3131   }
3132   Patch::gbisP3Ready();
3133 }
3134
3135 //receive proxy results from phase 1
3136 void HomePatch::receiveResult(ProxyGBISP1ResultMsg *msg) {
3137   ++numGBISP1Arrived;
3138     for ( int i = 0; i < msg->psiSumLen; ++i ) {
3139       psiFin[i] += msg->psiSum[i];
3140     }
3141   delete msg;
3142
3143   if (flags.doNonbonded) {
3144     //awaken if phase 1 done
3145     if (phase1BoxClosedCalled == true &&
3146         numGBISP1Arrived==proxy.size() ) {
3147         sequencer->awaken();
3148     }
3149   } else {
3150     //awaken if all phases done on noWork step
3151     if (boxesOpen == 0 &&
3152         numGBISP1Arrived == proxy.size() &&
3153         numGBISP2Arrived == proxy.size() &&
3154         numGBISP3Arrived == proxy.size()) {
3155       sequencer->awaken();
3156     }
3157   }
3158 }
3159
3160 //receive proxy results from phase 2
3161 void HomePatch::receiveResult(ProxyGBISP2ResultMsg *msg) {
3162   ++numGBISP2Arrived;
3163   //accumulate dEda
3164   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
3165     dHdrPrefix[i] += msg->dEdaSum[i];
3166   }
3167   delete msg;
3168
3169   if (flags.doNonbonded) {
3170     //awaken if phase 2 done
3171     if (phase2BoxClosedCalled == true &&
3172         numGBISP2Arrived==proxy.size() ) {
3173       sequencer->awaken();
3174     }
3175   } else {
3176     //awaken if all phases done on noWork step
3177     if (boxesOpen == 0 &&
3178         numGBISP1Arrived == proxy.size() &&
3179         numGBISP2Arrived == proxy.size() &&
3180         numGBISP3Arrived == proxy.size()) {
3181       sequencer->awaken();
3182     }
3183   }
3184 }
3185
3186 //  MOLLY algorithm part 1
3187 void HomePatch::mollyAverage()
3188 {
3189   Molecule *mol = Node::Object()->molecule;
3190   SimParameters *simParams = Node::Object()->simParameters;
3191   BigReal tol = simParams->mollyTol;
3192   int maxiter = simParams->mollyIter;
3193   int i, iter;
3194   HGArrayBigReal dsq;
3195   BigReal tmp;
3196   HGArrayInt ial, ibl;
3197   HGArrayVector ref;  // reference position
3198   HGArrayVector refab;  // reference vector
3199   HGArrayBigReal rmass;  // 1 / mass
3200   BigReal *lambda;  // Lagrange multipliers
3201   CompAtom *avg;  // averaged position
3202   int numLambdas = 0;
3203   HGArrayInt fixed;  // is atom fixed?
3204
3205   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
3206   p_avg.resize(numAtoms);
3207   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
3208
3209   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3210     int hgs = atom[ig].hydrogenGroupSize;
3211     if ( hgs == 1 ) continue;  // only one atom in group
3212         for ( i = 0; i < hgs; ++i ) {
3213           ref[i] = atom[ig+i].position;
3214           rmass[i] = 1. / atom[ig+i].mass;
3215           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
3216           if ( fixed[i] ) rmass[i] = 0.;
3217         }
3218         avg = &(p_avg[ig]);
3219         int icnt = 0;
3220
3221   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
3222           if ( hgs != 3 ) {
3223             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
3224           }
3225           if ( !(fixed[1] && fixed[2]) ) {
3226             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
3227           }
3228         }
3229         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
3230     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
3231             if ( !(fixed[0] && fixed[i]) ) {
3232               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
3233             }
3234           }
3235         }
3236         if ( icnt == 0 ) continue;  // no constraints
3237         numLambdas += icnt;
3238         molly_lambda.resize(numLambdas);
3239         lambda = &(molly_lambda[numLambdas - icnt]);
3240         for ( i = 0; i < icnt; ++i ) {
3241           refab[i] = ref[ial[i]] - ref[ibl[i]];
3242         }
3243         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
3244         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
3245         if ( iter == maxiter ) {
3246           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
3247         }
3248   }
3249
3250   // for ( i=0; i<numAtoms; ++i ) {
3251   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
3252   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
3253   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
3254   //    }
3255   // }
3256
3257 }
3258
3259
3260 //  MOLLY algorithm part 2
3261 void HomePatch::mollyMollify(Tensor *virial)
3262 {
3263   Molecule *mol = Node::Object()->molecule;
3264   SimParameters *simParams = Node::Object()->simParameters;
3265   Tensor wc;  // constraint virial
3266   int i;
3267   HGArrayInt ial, ibl;
3268   HGArrayVector ref;  // reference position
3269   CompAtom *avg;  // averaged position
3270   HGArrayVector refab;  // reference vector
3271   HGArrayVector force;  // new force
3272   HGArrayBigReal rmass;  // 1 / mass
3273   BigReal *lambda;  // Lagrange multipliers
3274   int numLambdas = 0;
3275   HGArrayInt fixed;  // is atom fixed?
3276
3277   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
3278     int hgs = atom[ig].hydrogenGroupSize;
3279     if (hgs == 1 ) continue;  // only one atom in group
3280         for ( i = 0; i < hgs; ++i ) {
3281           ref[i] = atom[ig+i].position;
3282           force[i] = f[Results::slow][ig+i];
3283           rmass[i] = 1. / atom[ig+i].mass;
3284           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
3285           if ( fixed[i] ) rmass[i] = 0.;
3286         }
3287         int icnt = 0;
3288         // c-ji I'm only going to mollify water for now
3289   if ( atom[ig].rigidBondLength > 0 ) {  // for water
3290           if ( hgs != 3 ) {
3291             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
3292           }
3293           if ( !(fixed[1] && fixed[2]) ) {
3294             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
3295           }
3296         }
3297         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
3298     if ( atom[ig+i].rigidBondLength > 0 ) {
3299             if ( !(fixed[0] && fixed[i]) ) {
3300               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
3301             }
3302           }
3303         }
3304
3305         if ( icnt == 0 ) continue;  // no constraints
3306         lambda = &(molly_lambda[numLambdas]);
3307         numLambdas += icnt;
3308         for ( i = 0; i < icnt; ++i ) {
3309           refab[i] = ref[ial[i]] - ref[ibl[i]];
3310         }
3311         avg = &(p_avg[ig]);
3312         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
3313         // store data back to patch
3314         for ( i = 0; i < hgs; ++i ) {
3315           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
3316           f[Results::slow][ig+i] = force[i];
3317         }
3318   }
3319   // check that there isn't a constant needed here!
3320   *virial += wc;
3321   p_avg.resize(0);
3322 }
3323
3324 void HomePatch::checkpoint(void) {
3325   checkpoint_atom.copy(atom);
3326   checkpoint_lattice = lattice;
3327
3328   // DMK - Atom Separation (water vs. non-water)
3329   #if NAMD_SeparateWaters != 0
3330     checkpoint_numWaterAtoms = numWaterAtoms;
3331   #endif
3332 }
3333
3334 void HomePatch::revert(void) {
3335   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3336
3337   atom.copy(checkpoint_atom);
3338   numAtoms = atom.size();
3339   lattice = checkpoint_lattice;
3340
3341   doAtomUpdate = true;
3342   rattleListValid = false;
3343
3344   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3345
3346   // DMK - Atom Separation (water vs. non-water)
3347   #if NAMD_SeparateWaters != 0
3348     numWaterAtoms = checkpoint_numWaterAtoms;
3349   #endif
3350 }
3351
3352 void HomePatch::exchangeCheckpoint(int scriptTask, int &bpc) {  // initiating replica
3353   SimParameters *simParams = Node::Object()->simParameters;
3354   checkpoint_task = scriptTask;
3355   const int remote = simParams->scriptIntArg1;
3356   const char *key = simParams->scriptStringArg1;
3357   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
3358 }
3359
3360 void HomePatch::recvCheckpointReq(int task, const char *key, int replica, int pe) {  // responding replica
3361   if ( task == SCRIPT_CHECKPOINT_FREE ) {
3362     if ( ! checkpoints.count(key) ) {
3363       NAMD_die("Unable to free checkpoint, requested key was never stored.");
3364     }
3365     delete checkpoints[key];
3366     checkpoints.erase(key);
3367     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
3368     return;
3369   }
3370   CheckpointAtomsMsg *msg;
3371   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
3372     if ( ! checkpoints.count(key) ) {
3373       NAMD_die("Unable to load checkpoint, requested key was never stored.");
3374     }
3375     checkpoint_t &cp = *checkpoints[key];
3376     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
3377     msg->lattice = cp.lattice;
3378     msg->berendsenPressure_count = cp.berendsenPressure_count;
3379     msg->numAtoms = cp.numAtoms;
3380     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
3381   } else {
3382     msg = new (0,1,0) CheckpointAtomsMsg;
3383   }
3384   msg->pid = patchID;
3385   msg->replica = CmiMyPartition();
3386   msg->pe = CkMyPe();
3387   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
3388 }
3389
3390 void HomePatch::recvCheckpointLoad(CheckpointAtomsMsg *msg) {  // initiating replica
3391   SimParameters *simParams = Node::Object()->simParameters;
3392   const int remote = simParams->scriptIntArg1;
3393   const char *key = simParams->scriptStringArg1;
3394   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
3395     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
3396   }
3397   if ( msg->replica != remote ) {
3398     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
3399   }
3400   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
3401     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
3402     strcpy(newmsg->key,key);
3403     newmsg->lattice = lattice;
3404     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
3405     newmsg->pid = patchID;
3406     newmsg->pe = CkMyPe();
3407     newmsg->replica = CmiMyPartition();
3408     newmsg->numAtoms = numAtoms;
3409     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
3410     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
3411   }
3412   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
3413     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3414     lattice = msg->lattice;
3415     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
3416     numAtoms = msg->numAtoms;
3417     atom.resize(numAtoms);
3418     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
3419     doAtomUpdate = true;
3420     rattleListValid = false;
3421     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3422   }
3423   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
3424     recvCheckpointAck();
3425   }
3426   delete msg;
3427 }
3428
3429 void HomePatch::recvCheckpointStore(CheckpointAtomsMsg *msg) {  // responding replica
3430   if ( ! checkpoints.count(msg->key) ) {
3431     checkpoints[msg->key] = new checkpoint_t;
3432   }
3433   checkpoint_t &cp = *checkpoints[msg->key];
3434   cp.lattice = msg->lattice;
3435   cp.berendsenPressure_count = msg->berendsenPressure_count;
3436   cp.numAtoms = msg->numAtoms;
3437   cp.atoms.resize(cp.numAtoms);
3438   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
3439   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
3440   delete msg;
3441 }
3442
3443 void HomePatch::recvCheckpointAck() {  // initiating replica
3444   CkpvAccess(_qd)->process();
3445 }
3446
3447
3448 void HomePatch::exchangeAtoms(int scriptTask) {
3449   SimParameters *simParams = Node::Object()->simParameters;
3450   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
3451   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
3452     exchange_dst = (int) simParams->scriptArg1;
3453     // create and save outgoing message
3454     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
3455     exchange_msg->lattice = lattice;
3456     exchange_msg->pid = patchID;
3457     exchange_msg->numAtoms = numAtoms;
3458     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
3459     if ( exchange_req >= 0 ) {
3460       recvExchangeReq(exchange_req);
3461     }
3462   }
3463   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
3464     exchange_src = (int) simParams->scriptArg2;
3465     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
3466   }
3467 }
3468
3469 void HomePatch::recvExchangeReq(int req) {
3470   exchange_req = req;
3471   if ( exchange_msg ) {
3472     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
3473     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
3474     exchange_msg = 0;
3475     exchange_req = -1;
3476     CkpvAccess(_qd)->process();
3477   }
3478 }
3479
3480 void HomePatch::recvExchangeMsg(ExchangeAtomsMsg *msg) {
3481   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
3482   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
3483   lattice = msg->lattice;
3484   numAtoms = msg->numAtoms;
3485   atom.resize(numAtoms);
3486   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
3487   delete msg;
3488   CkpvAccess(_qd)->process();
3489   doAtomUpdate = true;
3490   rattleListValid = false;
3491   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
3492 }
3493
3494 void HomePatch::submitLoadStats(int timestep)
3495 {
3496   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
3497 }
3498
3499
3500 void HomePatch::doPairlistCheck()
3501 {
3502   SimParameters *simParams = Node::Object()->simParameters;
3503
3504   if ( numAtoms == 0 || ! flags.usePairlists ) {
3505     flags.pairlistTolerance = 0.;
3506     flags.maxAtomMovement = 99999.;
3507     return;
3508   }
3509
3510   int i; int n = numAtoms;
3511   CompAtom *p_i = p.begin();
3512
3513   if ( flags.savePairlists ) {
3514     flags.pairlistTolerance = doPairlistCheck_newTolerance;
3515     flags.maxAtomMovement = 0.;
3516     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
3517     doPairlistCheck_lattice = lattice;
3518     doPairlistCheck_positions.resize(numAtoms);
3519     CompAtom *psave_i = doPairlistCheck_positions.begin();
3520     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
3521     return;
3522   }
3523
3524   Lattice &lattice_old = doPairlistCheck_lattice;
3525   Position center_cur = lattice.unscale(center);
3526   Position center_old = lattice_old.unscale(center);
3527   Vector center_delta = center_cur - center_old;
3528   
3529   // find max deviation to corner (any neighbor shares a corner)
3530   BigReal max_cd = 0.;
3531   for ( i=0; i<2; ++i ) {
3532     for ( int j=0; j<2; ++j ) {
3533       for ( int k=0; k<2; ++k ) {
3534         ScaledPosition corner(  i ? min.x : max.x ,
3535                                 j ? min.y : max.y ,
3536                                 k ? min.z : max.z );
3537         Vector corner_delta =
3538                 lattice.unscale(corner) - lattice_old.unscale(corner);
3539         corner_delta -= center_delta;
3540         BigReal cd = corner_delta.length2();
3541         if ( cd > max_cd ) max_cd = cd;
3542       }
3543     }
3544   }
3545   max_cd = sqrt(max_cd);
3546
3547   // find max deviation of atoms relative to center
3548   BigReal max_pd = 0.;
3549   CompAtom *p_old_i = doPairlistCheck_positions.begin();
3550   for ( i=0; i<n; ++i ) {
3551     Vector p_delta = p_i[i].position - p_old_i[i].position;
3552     p_delta -= center_delta;
3553     BigReal pd = p_delta.length2();
3554     if ( pd > max_pd ) max_pd = pd;
3555   }
3556   max_pd = sqrt(max_pd);
3557
3558   BigReal max_tol = max_pd + max_cd;
3559
3560   flags.maxAtomMovement = max_tol;
3561
3562   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
3563
3564   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
3565                                 doPairlistCheck_newTolerance ) ) {
3566     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
3567   }
3568
3569   if ( max_tol > doPairlistCheck_newTolerance ) {
3570     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
3571   }
3572
3573 }
3574
3575 void HomePatch::doGroupSizeCheck()
3576 {
3577   if ( ! flags.doNonbonded ) return;
3578
3579   SimParameters *simParams = Node::Object()->simParameters;
3580   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
3581   BigReal maxrad2 = 0.;
3582
3583   FullAtomList::iterator p_i = atom.begin();
3584   FullAtomList::iterator p_e = atom.end();
3585
3586   while ( p_i != p_e ) {
3587     const int hgs = p_i->hydrogenGroupSize;
3588     if ( ! hgs ) break;  // avoid infinite loop on bug
3589     int ngs = hgs;
3590     if ( ngs > 5 ) ngs = 5;  // limit to at most 5 atoms per group
3591     BigReal x = p_i->position.x;
3592     BigReal y = p_i->position.y;
3593     BigReal z = p_i->position.z;
3594     int i;
3595     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
3596       p_i[i].nonbondedGroupSize = 0;
3597       BigReal dx = p_i[i].position.x - x;
3598       BigReal dy = p_i[i].position.y - y;
3599       BigReal dz = p_i[i].position.z - z;
3600       BigReal r2 = dx * dx + dy * dy + dz * dz;
3601       if ( r2 > hgcut ) break;
3602       else if ( r2 > maxrad2 ) maxrad2 = r2;
3603     }
3604     p_i->nonbondedGroupSize = i;
3605     for ( ; i < hgs; ++i ) {
3606       p_i[i].nonbondedGroupSize = 1;
3607     }
3608     p_i += hgs;
3609   }
3610
3611   if ( p_i != p_e ) {
3612     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
3613   }
3614
3615   flags.maxGroupRadius = sqrt(maxrad2);
3616
3617 }
3618
3619 void HomePatch::doMarginCheck()
3620 {
3621   SimParameters *simParams = Node::Object()->simParameters;
3622
3623   BigReal sysdima = lattice.a_r().unit() * lattice.a();
3624   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
3625   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
3626
3627   BigReal minSize = simParams->patchDimension - simParams->margin;
3628
3629   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
3630        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
3631        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
3632
3633     NAMD_die("Periodic cell has become too small for original patch grid!\n"
3634       "Possible solutions are to restart from a recent checkpoint,\n"
3635       "increase margin, or disable useFlexibleCell for liquid simulation.");
3636   }
3637
3638   BigReal cutoff = simParams->cutoff;
3639
3640   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
3641   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
3642   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
3643
3644   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
3645     NAMD_die("Periodic cell has become too small for original patch grid!\n"
3646       "There are probably many margin violations already reported.\n"
3647       "Possible solutions are to restart from a recent checkpoint,\n"
3648       "increase margin, or disable useFlexibleCell for liquid simulation.");
3649   }
3650
3651   BigReal minx = min.x - margina;
3652   BigReal miny = min.y - marginb;
3653   BigReal minz = min.z - marginc;
3654   BigReal maxx = max.x + margina;
3655   BigReal maxy = max.y + marginb;
3656   BigReal maxz = max.z + marginc;
3657
3658   int xdev, ydev, zdev;
3659   int problemCount = 0;
3660
3661   FullAtomList::iterator p_i = atom.begin();
3662   FullAtomList::iterator p_e = atom.end();
3663   for ( ; p_i != p_e; ++p_i ) {
3664
3665     ScaledPosition s = lattice.scale(p_i->position);
3666
3667     // check if atom is within bounds
3668     if (s.x < minx) xdev = 0;
3669     else if (maxx <= s.x) xdev = 2; 
3670     else xdev = 1;
3671
3672     if (s.y < miny) ydev = 0;
3673     else if (maxy <= s.y) ydev = 2; 
3674     else ydev = 1;
3675
3676     if (s.z < minz) zdev = 0;
3677     else if (maxz <= s.z) zdev = 2; 
3678     else zdev = 1;
3679
3680     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
3681         ++problemCount;
3682     }
3683
3684   }
3685
3686   marginViolations = problemCount;
3687   // if ( problemCount ) {
3688   //     iout << iERROR <<
3689   //       "Found " << problemCount << " margin violations!\n" << endi;
3690   // } 
3691
3692 }
3693
3694
3695 void
3696 HomePatch::doAtomMigration()
3697 {
3698   int i;
3699
3700   for (i=0; i<numNeighbors; i++) {
3701     realInfo[i].mList.resize(0);
3702   }
3703
3704   // Purge the AtomMap
3705   atomMapper-&