HomePatch.C

Go to the documentation of this file.
00001 
00008 /*
00009    HomePatch owns the actual atoms of a Patch of space
00010    Proxy(s) get messages via ProxyMgr from HomePatch(es)
00011    to update lists of atoms and their coordinates
00012    HomePatch(es) also have a Sequencer bound to them
00013 
00014    superclass:  Patch           
00015 */
00016 
00017 #include "time.h"
00018 #include <math.h>
00019 #include "charm++.h"
00020 #include "qd.h"
00021 
00022 #include "SimParameters.h"
00023 #include "HomePatch.h"
00024 #include "AtomMap.h"
00025 #include "Node.h"
00026 #include "PatchMap.inl"
00027 #include "main.h"
00028 #include "ProxyMgr.decl.h"
00029 #include "ProxyMgr.h"
00030 #include "Migration.h"
00031 #include "Molecule.h"
00032 #include "PatchMgr.h"
00033 #include "Sequencer.h"
00034 #include "Broadcasts.h"
00035 #include "LdbCoordinator.h"
00036 #include "ReductionMgr.h"
00037 #include "Sync.h"
00038 #include "Random.h"
00039 #include "Priorities.h"
00040 #include "ComputeNonbondedUtil.h"
00041 #include "ComputeGBIS.inl"
00042 #include "Priorities.h"
00043 #include "SortAtoms.h"
00044 
00045 #include "ComputeQM.h"
00046 #include "ComputeQMMgr.decl.h"
00047 
00048 //#define PRINT_COMP
00049 #define TINY 1.0e-20;
00050 #define MAXHGS 10
00051 #define MIN_DEBUG_LEVEL 2
00052 //#define DEBUGM
00053 #include "Debug.h"
00054 
00055 #include <vector>
00056 #include <algorithm>
00057 using namespace std;
00058 
00059 typedef int HGArrayInt[MAXHGS];
00060 typedef BigReal HGArrayBigReal[MAXHGS];
00061 typedef zVector HGArrayVector[MAXHGS];
00062 typedef BigReal HGMatrixBigReal[MAXHGS][MAXHGS];
00063 typedef zVector HGMatrixVector[MAXHGS][MAXHGS];
00064 
00065 #include "ComputeNonbondedMICKernel.h"
00066 
00067 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);
00068 
00069 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);
00070 
00071 #define MASS_EPSILON  (1.0e-35)  //a very small floating point number
00072 
00073 
00074 // DMK - Atom Separation (water vs. non-water)
00075 #if NAMD_SeparateWaters != 0
00076 
00077 // Macro to test if a hydrogen group represents a water molecule.
00078 // NOTE: This test is the same test in Molecule.C for setting the
00079 //   OxygenAtom flag in status.
00080 // hgtype should be the number of atoms in a water hydrogen group
00081 // It must now be set based on simulation parameters because we might
00082 // be using tip4p
00083 
00084 // DJH: This will give false positive for full Drude model,
00085 //      e.g. O D H is not water but has hgs==3
00086 #define IS_HYDROGEN_GROUP_WATER(hgs, mass)                 \
00087   ((hgs >= 3) && ((mass >= 14.0) && (mass <= 18.0)))
00088 
00089 #endif
00090 
00091 
00092 HomePatch::HomePatch(PatchID pd, FullAtomList &al) : Patch(pd)
00093 // DMK - Atom Separation (water vs. non-water)
00094 #if NAMD_SeparateWaters != 0
00095   ,tempAtom()
00096 #endif
00097 { 
00098   atom.swap(al);
00099   settle_initialized = 0;
00100 
00101   doAtomUpdate = true;
00102   rattleListValid = false;
00103 
00104   exchange_msg = 0;
00105   exchange_req = -1;
00106 
00107   numGBISP1Arrived = 0;
00108   numGBISP2Arrived = 0;
00109   numGBISP3Arrived = 0;
00110   phase1BoxClosedCalled = false;
00111   phase2BoxClosedCalled = false;
00112   phase3BoxClosedCalled = false;
00113 
00114   min.x = PatchMap::Object()->min_a(patchID);
00115   min.y = PatchMap::Object()->min_b(patchID);
00116   min.z = PatchMap::Object()->min_c(patchID);
00117   max.x = PatchMap::Object()->max_a(patchID);
00118   max.y = PatchMap::Object()->max_b(patchID);
00119   max.z = PatchMap::Object()->max_c(patchID);
00120   center = 0.5*(min+max);
00121 
00122   int aAway = PatchMap::Object()->numaway_a();
00123   if ( PatchMap::Object()->periodic_a() ||
00124        PatchMap::Object()->gridsize_a() > aAway + 1 ) {
00125     aAwayDist = (max.x - min.x) * aAway;
00126   } else {
00127     aAwayDist = Node::Object()->simParameters->patchDimension;
00128   }
00129   int bAway = PatchMap::Object()->numaway_b();
00130   if ( PatchMap::Object()->periodic_b() ||
00131        PatchMap::Object()->gridsize_b() > bAway + 1 ) {
00132     bAwayDist = (max.y - min.y) * bAway;
00133   } else {
00134     bAwayDist = Node::Object()->simParameters->patchDimension;
00135   }
00136   int cAway = PatchMap::Object()->numaway_c();
00137   if ( PatchMap::Object()->periodic_c() ||
00138        PatchMap::Object()->gridsize_c() > cAway + 1 ) {
00139     cAwayDist = (max.z - min.z) * cAway;
00140   } else {
00141     cAwayDist = Node::Object()->simParameters->patchDimension;
00142   }
00143 
00144   migrationSuspended = false;
00145   allMigrationIn = false;
00146   marginViolations = 0;
00147   patchMapRead = 0; // We delay read of PatchMap data
00148                     // to make sure it is really valid
00149   inMigration = false;
00150   numMlBuf = 0;
00151   flags.sequence = -1;
00152   flags.maxForceUsed = -1;
00153 
00154   numAtoms = atom.size();
00155   replacementForces = 0;
00156 
00157   SimParameters *simParams = Node::Object()->simParameters;
00158   doPairlistCheck_newTolerance = 
00159         0.5 * ( simParams->pairlistDist - simParams->cutoff );
00160 
00161 
00162   numFixedAtoms = 0;
00163   if ( simParams->fixedAtomsOn ) {
00164     for ( int i = 0; i < numAtoms; ++i ) {
00165       numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
00166     }
00167   }
00168 
00169 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00170   ptnTree.resize(0);
00171   /*children = NULL;
00172   numChild = 0;*/
00173 #else
00174   child =  new int[proxySpanDim];
00175   nChild = 0;   // number of proxy spanning tree children
00176 #endif
00177 
00178 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00179   nphs = 0;
00180   localphs = NULL;
00181   isProxyChanged = 0;
00182 #endif
00183 
00184 
00185   // DMK - Atom Separation (water vs. non-water)
00186   #if NAMD_SeparateWaters != 0
00187 
00188     // Create the scratch memory for separating atoms
00189     tempAtom.resize(numAtoms);
00190     numWaterAtoms = -1;
00191 
00192     // Separate the current list of atoms
00193     separateAtoms();
00194 
00195   #endif
00196     
00197   // Handle unusual water models here
00198   if (simParams->watmodel == WAT_TIP4) init_tip4();
00199   else if (simParams->watmodel == WAT_SWM4) init_swm4();
00200 
00201   isNewProxyAdded = 0;
00202 }
00203 
00204 void HomePatch::write_tip4_props() {
00205   printf("Writing r_om and r_ohc: %f | %f\n", r_om, r_ohc);
00206 }
00207 
00208 void HomePatch::init_tip4() {
00209   // initialize the distances needed for the tip4p water model
00210   Molecule *mol = Node::Object()->molecule;
00211   r_om = mol->r_om;
00212   r_ohc = mol->r_ohc;
00213 }
00214 
00215 
00216 void ::HomePatch::init_swm4() {
00217   // initialize the distances needed for the SWM4 water model
00218   Molecule *mol = Node::Object()->molecule;
00219   r_om = mol->r_om;
00220   r_ohc = mol->r_ohc;
00221 }
00222 
00223 
00224 void HomePatch::reinitAtoms(FullAtomList &al) {
00225   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
00226 
00227   atom.swap(al);
00228   numAtoms = atom.size();
00229 
00230   doAtomUpdate = true;
00231   rattleListValid = false;
00232 
00233   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
00234 
00235   // DMK - Atom Separation (water vs. non-water)
00236   #if NAMD_SeparateWaters != 0
00237 
00238     // Reset the numWaterAtoms value
00239     numWaterAtoms = -1;
00240 
00241     // Separate the atoms
00242     separateAtoms();
00243 
00244   #endif
00245 }
00246 
00247 // Bind a Sequencer to this HomePatch
00248 void HomePatch::useSequencer(Sequencer *sequencerPtr)
00249 { sequencer=sequencerPtr; }
00250 
00251 // start simulation over this Patch of atoms
00252 void HomePatch::runSequencer(void)
00253 { sequencer->run(); }
00254 
00255 void HomePatch::readPatchMap() {
00256   // iout << "Patch " << patchID << " has " << proxy.size() << " proxies.\n" << endi;
00257   PatchMap *p = PatchMap::Object();
00258   PatchID nnPatchID[PatchMap::MaxOneAway];
00259 
00260   patchMigrationCounter = numNeighbors 
00261     = PatchMap::Object()->oneAwayNeighbors(patchID, nnPatchID);
00262   DebugM( 1, "NumNeighbors for pid " <<patchID<<" is "<< numNeighbors << "\n");
00263   int n;
00264   for (n=0; n<numNeighbors; n++) {
00265     realInfo[n].destNodeID = p->node(realInfo[n].destPatchID = nnPatchID[n]);
00266      DebugM( 1, " nnPatchID=" <<nnPatchID[n]<<" nnNodeID="<< realInfo[n].destNodeID<< "\n");
00267     realInfo[n].mList.resize(0);
00268   }
00269 
00270   // Make mapping from the 3x3x3 cube of pointers to real migration info
00271   for (int i=0; i<3; i++)
00272     for (int j=0; j<3; j++)
00273       for (int k=0; k<3; k++)
00274       {
00275         int pid =  p->pid(p->index_a(patchID)+i-1, 
00276             p->index_b(patchID)+j-1, p->index_c(patchID)+k-1);
00277         if (pid < 0) {
00278            DebugM(5, "ERROR, for patchID " << patchID <<" I got neigh pid = " << pid << "\n");
00279         }
00280         if (pid == patchID && ! (
00281                 ( (i-1) && p->periodic_a() ) ||
00282                 ( (j-1) && p->periodic_b() ) ||
00283                 ( (k-1) && p->periodic_c() ) )) {
00284           mInfo[i][j][k] = NULL;
00285         }
00286         else {
00287           // Does not work as expected for periodic with only two patches.
00288           // Also need to check which image we want, but OK for now.  -JCP
00289           for (n = 0; n<numNeighbors; n++) {
00290             if (pid == realInfo[n].destPatchID) {
00291               mInfo[i][j][k] = &realInfo[n];
00292               break;
00293             }
00294           }
00295           if (n == numNeighbors) { // disaster! 
00296             DebugM(4,"BAD News, I could not find PID " << pid << "\n");
00297           }
00298         }
00299       }
00300 
00301   DebugM(1,"Patch("<<patchID<<") # of neighbors = " << numNeighbors << "\n");
00302 }
00303 
00304 HomePatch::~HomePatch()
00305 {
00306     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
00307 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00308     ptnTree.resize(0);
00309     #ifdef USE_NODEPATCHMGR
00310     delete [] nodeChildren;    
00311     #endif
00312 #endif
00313   delete [] child;
00314 }
00315 
00316 
00317 void HomePatch::boxClosed(int box) {
00318   // begin gbis
00319   if (box == 5) {// end of phase 1
00320     phase1BoxClosedCalled = true;
00321     if (!psiSumBox.isOpen() && numGBISP1Arrived == proxy.size()) {
00322       if (flags.doGBIS && flags.doNonbonded) {
00323         sequencer->awaken();
00324       }
00325     } else {
00326       //need to wait until proxies arrive before awakening
00327     }
00328   } else if (box == 6) {// intRad
00329     //do nothing
00330   } else if (box == 7) {// bornRad
00331     //do nothing
00332   } else if (box == 8) {// end of phase 2
00333     phase2BoxClosedCalled = true;
00334     //if no proxies, AfterP1 can't be called from receive
00335     //so it will be called from here
00336     if (!dEdaSumBox.isOpen() && numGBISP2Arrived == proxy.size()) {
00337       if (flags.doGBIS && flags.doNonbonded) {
00338         sequencer->awaken();
00339       }
00340     } else {
00341       //need to wait until proxies arrive before awakening
00342     }
00343   } else if (box == 9) {
00344     //do nothing
00345   } else if (box == 10) {
00346     //lcpoType Box closed: do nothing
00347   } else {
00348     //do nothing
00349   }
00350   // end gbis
00351 
00352   if ( ! --boxesOpen ) {
00353     if ( replacementForces ) {
00354       for ( int i = 0; i < numAtoms; ++i ) {
00355         if ( replacementForces[i].replace ) {
00356           for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
00357           f[Results::normal][i] = replacementForces[i].force;
00358         }
00359       }
00360       replacementForces = 0;
00361     }
00362     DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
00363         << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
00364     // only awaken suspended threads.  Then say it is suspended.
00365 
00366     phase3BoxClosedCalled = true;
00367     if (flags.doGBIS) {
00368       if (flags.doNonbonded) {
00369         sequencer->awaken();
00370       } else {
00371         if (numGBISP1Arrived == proxy.size() &&
00372           numGBISP2Arrived == proxy.size() &&
00373           numGBISP3Arrived == proxy.size()) {
00374           sequencer->awaken();//all boxes closed and all proxies arrived
00375         }
00376       }
00377     } else {//non-gbis awaken
00378       sequencer->awaken();
00379     }
00380   } else {
00381     DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
00382   }
00383 }
00384 
00385 void HomePatch::registerProxy(RegisterProxyMsg *msg) {
00386   DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
00387   proxy.add(msg->node);
00388   forceBox.clientAdd();
00389 
00390   isNewProxyAdded = 1;
00391 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00392   isProxyChanged = 1;
00393 #endif
00394 
00395   Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
00396   delete msg;
00397 }
00398 
00399 void HomePatch::unregisterProxy(UnregisterProxyMsg *msg) {
00400 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00401   isProxyChanged = 1;
00402 #endif
00403   int n = msg->node;
00404   NodeID *pe = proxy.begin();
00405   for ( ; *pe != n; ++pe );
00406   forceBox.clientRemove();
00407   proxy.del(pe - proxy.begin());
00408   delete msg;
00409 }
00410 
00411 #if USE_TOPOMAP && USE_SPANNING_TREE
00412 
00413 int HomePatch::findSubroots(int dim, int* subroots, int psize, int* pidscopy){
00414   int nChild = 0;
00415   int cones[6][proxySpanDim*proxySpanDim+proxySpanDim];
00416   int conesizes[6] = {0,0,0,0,0,0};
00417   int conecounters[6] = {0,0,0,0,0,0};
00418   int childcounter = 0;
00419   nChild = (psize>proxySpanDim)?proxySpanDim:psize;
00420   TopoManager tmgr;
00421   for(int i=0;i<psize;i++){
00422     int cone = tmgr.getConeNumberForRank(pidscopy[i]);
00423     cones[cone][conesizes[cone]++] = pidscopy[i];
00424   }
00425 
00426   while(childcounter<nChild){
00427     for(int i=0;i<6;i++){
00428       if(conecounters[i]<conesizes[i]){
00429         subroots[childcounter++] = cones[i][conecounters[i]++];
00430       }
00431     }
00432   }
00433   for(int i=nChild;i<proxySpanDim;i++)
00434     subroots[i] = -1;
00435   return nChild;
00436 }
00437 #endif // USE_TOPOMAP 
00438 
00439 static int compDistance(const void *a, const void *b)
00440 {
00441   int d1 = abs(*(int *)a - CkMyPe());
00442   int d2 = abs(*(int *)b - CkMyPe());
00443   if (d1 < d2) 
00444     return -1;
00445   else if (d1 == d2) 
00446     return 0;
00447   else 
00448     return 1;
00449 }
00450 
00451 void HomePatch::sendProxies()
00452 {
00453 #if USE_NODEPATCHMGR
00454         CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
00455         NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
00456         npm->sendProxyList(patchID, proxy.begin(), proxy.size());
00457 #else
00458         ProxyMgr::Object()->sendProxies(patchID, proxy.begin(), proxy.size());
00459 #endif
00460 }
00461 
00462 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00463 void HomePatch::buildNodeAwareSpanningTree(void){
00464 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00465         DebugFileTrace *dft = DebugFileTrace::Object();
00466         dft->openTrace();
00467         dft->writeTrace("HomePatch[%d] has %d proxy on proc[%d] node[%d]\n", patchID, proxy.size(), CkMyPe(), CkMyNode());
00468         dft->writeTrace("Proxies are: ");
00469         for(int i=0; i<proxy.size(); i++) dft->writeTrace("%d(%d), ", proxy[i], CkNodeOf(proxy[i]));
00470         dft->writeTrace("\n");
00471         dft->closeTrace();
00472 #endif
00473  
00474     //build the naive spanning tree for this home patch
00475     if(! proxy.size()) {
00476         //this case will not happen in practice.
00477         //In debugging state where spanning tree is enforced, then this could happen
00478         //Chao Mei        
00479        return;
00480     }
00481     ProxyMgr::buildSinglePatchNodeAwareSpanningTree(patchID, proxy, ptnTree);
00482     //optimize on the naive spanning tree
00483 
00484     //setup the children
00485     setupChildrenFromProxySpanningTree();
00486     //send down to children
00487     sendNodeAwareSpanningTree();
00488 }
00489 
00490 void HomePatch::setupChildrenFromProxySpanningTree(){
00491 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00492     isProxyChanged = 1;
00493 #endif
00494     if(ptnTree.size()==0) {
00495         nChild = 0;
00496         delete [] child;
00497         child = NULL;
00498         #ifdef USE_NODEPATCHMGR
00499         numNodeChild = 0;
00500         delete [] nodeChildren;
00501         nodeChildren = NULL;        
00502         #endif
00503         return;
00504     }
00505     proxyTreeNode *rootnode = &ptnTree.item(0);
00506     CmiAssert(rootnode->peIDs[0] == CkMyPe());
00507     //set up children
00508     //1. add external children (the first proc inside the proxy tree node)    
00509     //2. add internal children (with threshold that would enable spanning    
00510     int internalChild = rootnode->numPes-1;
00511     int externalChild = ptnTree.size()-1;
00512     externalChild = (proxySpanDim>externalChild)?externalChild:proxySpanDim;
00513     int internalSlots = proxySpanDim-externalChild;
00514     if(internalChild>0){
00515       if(internalSlots==0) {
00516          //at least having one internal child
00517         internalChild = 1;
00518       }else{
00519         internalChild = (internalSlots>internalChild)?internalChild:internalSlots;
00520       }
00521     }
00522     
00523     nChild = externalChild+internalChild;
00524     CmiAssert(nChild>0);
00525 
00526     //exclude the root node        
00527     delete [] child;
00528     child = new int[nChild];    
00529 
00530     for(int i=0; i<externalChild; i++) {
00531         child[i] = ptnTree.item(i+1).peIDs[0];
00532     }
00533     for(int i=externalChild, j=1; i<nChild; i++, j++) {
00534         child[i] = rootnode->peIDs[j];
00535     }
00536 
00537 #ifdef USE_NODEPATCHMGR
00538     //only register the cores that have proxy patches. The HomePach's core
00539     //doesn't need to be registered.
00540     CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
00541     NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
00542     if(rootnode->numPes==1){
00543         npm->registerPatch(patchID, 0, NULL);        
00544     }
00545     else{
00546         npm->registerPatch(patchID, rootnode->numPes-1, &rootnode->peIDs[1]);        
00547     }
00548 
00549     //set up childrens in terms of node ids
00550     numNodeChild = externalChild;
00551     if(internalChild) numNodeChild++;
00552     delete [] nodeChildren;
00553     nodeChildren = new int[numNodeChild];    
00554     for(int i=0; i<externalChild; i++) {
00555         nodeChildren[i] = ptnTree.item(i+1).nodeID;        
00556     }
00557     //the last entry always stores this node id if there are proxies
00558     //on other cores of the same node for this patch.
00559     if(internalChild)
00560         nodeChildren[numNodeChild-1] = rootnode->nodeID;
00561 #endif
00562     
00563 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00564     DebugFileTrace *dft = DebugFileTrace::Object();
00565     dft->openTrace();
00566     dft->writeTrace("HomePatch[%d] has %d children: ", patchID, nChild);
00567     for(int i=0; i<nChild; i++)
00568         dft->writeTrace("%d ", child[i]);
00569     dft->writeTrace("\n");
00570     dft->closeTrace();
00571 #endif   
00572 }
00573 #endif
00574 
00575 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00576 //This is not an entry method, but takes an argument of message type
00577 void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){
00578     //set up the whole tree ptnTree
00579     int treesize = msg->numNodesWithProxies;    
00580     ptnTree.resize(treesize);    
00581     int *pAllPes = msg->allPes;
00582     for(int i=0; i<treesize; i++) {
00583         proxyTreeNode *oneNode = &ptnTree.item(i);
00584         delete [] oneNode->peIDs;
00585         oneNode->numPes = msg->numPesOfNode[i];
00586         oneNode->nodeID = CkNodeOf(*pAllPes);
00587         oneNode->peIDs = new int[oneNode->numPes];
00588         for(int j=0; j<oneNode->numPes; j++) {
00589             oneNode->peIDs[j] = *pAllPes;
00590             pAllPes++;
00591         }
00592     }
00593     //setup children
00594     setupChildrenFromProxySpanningTree();
00595     //send down to children
00596     sendNodeAwareSpanningTree();
00597 }
00598 
00599 void HomePatch::sendNodeAwareSpanningTree(){
00600     if(ptnTree.size()==0) return;    
00601     ProxyNodeAwareSpanningTreeMsg *msg = 
00602         ProxyNodeAwareSpanningTreeMsg::getANewMsg(patchID, CkMyPe(), ptnTree.begin(), ptnTree.size());
00603    
00604     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00605     msg->printOut("HP::sendST");
00606     #endif
00607 
00608     CmiAssert(CkMyPe() == msg->allPes[0]);
00609     ProxyMgr::Object()->sendNodeAwareSpanningTree(msg);
00610 
00611 }
00612 #else
00613 void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){}
00614 void HomePatch::sendNodeAwareSpanningTree(){}
00615 #endif
00616 
00617 #ifndef NODEAWARE_PROXY_SPANNINGTREE
00618 // recv a spanning tree from load balancer
00619 void HomePatch::recvSpanningTree(int *t, int n)
00620 {
00621   int i;
00622   nChild=0;
00623   tree.resize(n);
00624   for (i=0; i<n; i++) {
00625     tree[i] = t[i];
00626   }
00627 
00628   for (i=1; i<=proxySpanDim; i++) {
00629     if (tree.size() <= i) break;
00630     child[i-1] = tree[i];
00631     nChild++;
00632   }
00633 
00634 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00635   isProxyChanged = 1;
00636 #endif
00637 
00638   // send down to children
00639   sendSpanningTree();
00640 }
00641 
00642 void HomePatch::sendSpanningTree()
00643 {
00644   if (tree.size() == 0) return;
00645   ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
00646   msg->patch = patchID;
00647   msg->node = CkMyPe();
00648   msg->tree.copy(tree);  // copy data for thread safety
00649   ProxyMgr::Object()->sendSpanningTree(msg);  
00650 }
00651 #else
00652 void HomePatch::recvSpanningTree(int *t, int n){}
00653 void HomePatch::sendSpanningTree(){}
00654 #endif
00655 
00656 #ifndef NODEAWARE_PROXY_SPANNINGTREE
00657 void HomePatch::buildSpanningTree(void)
00658 {
00659   nChild = 0;
00660   int psize = proxy.size();
00661   if (psize == 0) return;
00662   NodeIDList oldtree;  oldtree.swap(tree);
00663   int oldsize = oldtree.size();
00664   tree.resize(psize + 1);
00665   tree.setall(-1);
00666   tree[0] = CkMyPe();
00667   int s=1, e=psize+1;
00668   NodeIDList::iterator pli;
00669   int patchNodesLast =
00670     ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00671   int nNonPatch = 0;
00672 #if 1
00673     // try to put it to the same old tree
00674   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
00675   {
00676     int oldindex = oldtree.find(*pli);
00677     if (oldindex != -1 && oldindex < psize) {
00678       tree[oldindex] = *pli;
00679     }
00680   }
00681   s=1; e=psize;
00682   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
00683   {
00684     if (tree.find(*pli) != -1) continue;    // already assigned
00685     if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
00686       while (tree[e] != -1) { e--; if (e==-1) e = psize; }
00687       tree[e] = *pli;
00688     } else {
00689       while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
00690       tree[s] = *pli;
00691       nNonPatch++;
00692     }
00693   }
00694 #if 1
00695   if (oldsize==0 && nNonPatch) {
00696     // first time, sort by distance
00697     qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
00698   }
00699 #endif
00700 
00701   //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
00702 
00703 #if USE_TOPOMAP && USE_SPANNING_TREE
00704   
00705   //Right now only works for spanning trees with two levels
00706   int *treecopy = new int [psize];
00707   int subroots[proxySpanDim];
00708   int subsizes[proxySpanDim];
00709   int subtrees[proxySpanDim][proxySpanDim];
00710   int idxes[proxySpanDim];
00711   int i = 0;
00712 
00713   for(i=0;i<proxySpanDim;i++){
00714     subsizes[i] = 0;
00715     idxes[i] = i;
00716   }
00717   
00718   for(i=0;i<psize;i++){
00719     treecopy[i] = tree[i+1];
00720   }
00721   
00722   TopoManager tmgr;
00723   tmgr.sortRanksByHops(treecopy,nNonPatch);
00724   tmgr.sortRanksByHops(treecopy+nNonPatch,
00725                                                 psize-nNonPatch);  
00726   
00727   /* build tree and subtrees */
00728   nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
00729   delete [] treecopy;
00730   
00731   for(int i=1;i<psize+1;i++){
00732     int isSubroot=0;
00733     for(int j=0;j<nChild;j++)
00734       if(tree[i]==subroots[j]){
00735         isSubroot=1;
00736         break;
00737       }
00738     if(isSubroot) continue;
00739     
00740     int bAdded = 0;
00741     tmgr.sortIndexByHops(tree[i], subroots,
00742                                                   idxes, proxySpanDim);
00743     for(int j=0;j<proxySpanDim;j++){
00744       if(subsizes[idxes[j]]<proxySpanDim){
00745         subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
00746         bAdded = 1; 
00747         break;
00748       }
00749     }
00750     if( psize > proxySpanDim && ! bAdded ) {
00751       NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
00752     }
00753   }
00754 
00755 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
00756   
00757   for (int i=1; i<=proxySpanDim; i++) {
00758     if (tree.size() <= i) break;
00759     child[i-1] = tree[i];
00760     nChild++;
00761   }
00762 #endif
00763 #endif
00764   
00765 #if 0
00766   // for debugging
00767   CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
00768   for (int i=0; i<psize+1; i++) {
00769     CkPrintf("%d ", tree[i]);
00770   }
00771   CkPrintf("\n");
00772 #endif
00773   // send to children nodes
00774   sendSpanningTree();
00775 }
00776 #endif
00777 
00778 
00779 void HomePatch::receiveResults(ProxyResultVarsizeMsg *msg) {
00780 
00781     numGBISP3Arrived++;
00782     DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00783     int n = msg->node;
00784     Results *r = forceBox.clientOpen();
00785 
00786     char *iszeroPtr = msg->isZero;
00787     Force *msgFPtr = msg->forceArr;
00788 
00789     for ( int k = 0; k < Results::maxNumForces; ++k )
00790     {
00791       Force *rfPtr = r->f[k];
00792       for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
00793           if((*iszeroPtr)!=1) {
00794               *rfPtr += *msgFPtr;
00795               msgFPtr++;
00796           }
00797       }      
00798     }
00799     forceBox.clientClose();
00800     delete msg;
00801 }
00802 
00803 void HomePatch::receiveResults(ProxyResultMsg *msg) {
00804   numGBISP3Arrived++;
00805   DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00806   int n = msg->node;
00807   Results *r = forceBox.clientOpen();
00808   for ( int k = 0; k < Results::maxNumForces; ++k )
00809   {
00810     Force *f = r->f[k];
00811     register ForceList::iterator f_i, f_e;
00812     f_i = msg->forceList[k]->begin();
00813     f_e = msg->forceList[k]->end();
00814     for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
00815   }
00816   forceBox.clientClose();
00817   delete msg;
00818 }
00819 
00820 void HomePatch::receiveResults(ProxyCombinedResultRawMsg* msg)
00821 {
00822     numGBISP3Arrived++;
00823   DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodeSize<<")\n");
00824     Results *r = forceBox.clientOpen(msg->nodeSize);
00825       register char* isNonZero = msg->isForceNonZero;
00826       register Force* f_i = msg->forceArr;
00827       for ( int k = 0; k < Results::maxNumForces; ++k )
00828       {
00829         Force *f = r->f[k];
00830                 int nf = msg->flLen[k];
00831 #ifdef ARCH_POWERPC
00832 #pragma disjoint (*f_i, *f)
00833 #endif
00834         for (int count = 0; count < nf; count++) {
00835           if(*isNonZero){
00836                         f[count].x += f_i->x;
00837                         f[count].y += f_i->y;
00838                         f[count].z += f_i->z;
00839                         f_i++;
00840           }
00841           isNonZero++;
00842         }
00843       }
00844     forceBox.clientClose(msg->nodeSize);
00845 
00846     delete msg;
00847 }
00848 
00849 void HomePatch::qmSwapAtoms()
00850 {
00851     // This is used for LSS in QM/MM simulations.
00852     // Changes atom labels so that we effectively exchange solvent
00853     // molecules between classical and quantum modes.
00854     
00855     SimParameters *simParams = Node::Object()->simParameters;
00856     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
00857     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
00858     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
00859     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
00860     
00861     ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
00862             CkpvAccess(BOCclass_group).computeQMMgr) ;
00863     
00864     FullAtom *a_i = atom.begin();
00865     
00866     for (int i=0; i<numAtoms; ++i ) { 
00867         
00868         LSSSubsDat *subP = lssSubs(mgrP).find( LSSSubsDat(a_i[i].id) ) ;
00869         
00870         if ( subP != NULL ) {
00871             a_i[i].id = subP->newID ;
00872             a_i[i].vdwType = subP->newVdWType ;
00873             
00874             // If we are swappign a classical atom with a QM one, the charge
00875             // will need extra handling.
00876             if (qmAtomGroup[subP->newID] > 0 && simParams->PMEOn) {
00877                 // We make sure that the last atom charge calculated for the 
00878                 // QM atom being transfered here is available for PME
00879                 // in the next step.
00880                 
00881                 // Loops over all QM atoms (in all QM groups) comparing their 
00882                 // global indices (the sequential atom ID from NAMD).
00883                 for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
00884                     
00885                     if (qmAtmIndx[qmIter] == subP->newID) {
00886                         qmAtmChrg[qmIter] = subP->newCharge;
00887                         break;
00888                     }
00889                     
00890                 }
00891                 
00892                 // For QM atoms, the charge in the full atom structure is zero.
00893                 // Electrostatic interactions between QM atoms and their 
00894                 // environment is handled in ComputeQM.
00895                 a_i[i].charge = 0;
00896             }
00897             else {
00898                 // If we are swappign a QM atom with a Classical one, only the charge
00899                 // in the full atomstructure needs updating, since it used to be zero.
00900                 a_i[i].charge = subP->newCharge ;
00901             }
00902         }
00903     }
00904     
00905     return;
00906 }
00907 
00908 void HomePatch::positionsReady(int doMigration)
00909 {    
00910   SimParameters *simParams = Node::Object()->simParameters;
00911 
00912   flags.sequence++;
00913 
00914   if (!patchMapRead) { readPatchMap(); }
00915       
00916   if (numNeighbors && ! simParams->staticAtomAssignment) {
00917     if (doMigration) {
00918       rattleListValid = false;
00919       doAtomMigration();
00920     } else {
00921       doMarginCheck();
00922     }
00923   }
00924   
00925   if (doMigration && simParams->qmLSSOn)
00926       qmSwapAtoms();
00927 
00928 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
00929   if ( doMigration ) {
00930     int n = numAtoms;
00931     FullAtom *a_i = atom.begin();
00932     #if defined(NAMD_CUDA) || (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
00933     int *ao = new int[n];
00934     int nfree;
00935     if ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces ) {
00936       int k = 0;
00937       int k2 = n;
00938       for ( int j=0; j<n; ++j ) {
00939         // put fixed atoms at end
00940         if ( a_i[j].atomFixed ) ao[--k2] = j;
00941         else ao[k++] = j;
00942       }
00943       nfree = k;
00944     } else {
00945       nfree = n;
00946       for ( int j=0; j<n; ++j ) {
00947         ao[j] = j;
00948       }
00949     }
00950 
00951     sortAtomsForCUDA(ao,a_i,nfree,n);
00952   
00953     for ( int i=0; i<n; ++i ) { 
00954       a_i[i].sortOrder = ao[i];
00955     }
00956     delete [] ao;
00957     #else
00958       for (int i = 0; i < n; ++i) {
00959         a_i[i].sortOrder = i;
00960       }
00961     #endif
00962   }
00963 
00964   { 
00965     const double charge_scaling = sqrt(COULOMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1);
00966     const Vector ucenter = lattice.unscale(center);
00967     const BigReal ucenter_x = ucenter.x;
00968     const BigReal ucenter_y = ucenter.y;
00969     const BigReal ucenter_z = ucenter.z;
00970     const int n = numAtoms;
00971     #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
00972       int n_16 = n;
00973       n_16 = (n + 15) & (~15);
00974       cudaAtomList.resize(n_16);
00975       CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
00976       mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
00977       mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
00978       mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
00979       mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
00980     #else
00981     cudaAtomList.resize(n);
00982     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
00983     #endif
00984     const FullAtom *a = atom.begin();
00985     for ( int k=0; k<n; ++k ) {
00986       #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
00987         int j = a[k].sortOrder;
00988         atom_x[k] = a[j].position.x - ucenter_x;
00989         atom_y[k] = a[j].position.y - ucenter_y;
00990         atom_z[k] = a[j].position.z - ucenter_z;
00991         atom_q[k] = charge_scaling * a[j].charge;
00992       #else
00993       int j = a[k].sortOrder;
00994       ac[k].x = a[j].position.x - ucenter_x;
00995       ac[k].y = a[j].position.y - ucenter_y;
00996       ac[k].z = a[j].position.z - ucenter_z;
00997       ac[k].q = charge_scaling * a[j].charge;
00998       #endif
00999     }
01000   }
01001 #else
01002   doMigration = doMigration && numNeighbors;
01003 #endif
01004   doMigration = doMigration || ! patchMapRead;
01005 
01006   doMigration = doMigration || doAtomUpdate;
01007   doAtomUpdate = false;
01008 
01009   // Workaround for oversize groups
01010   doGroupSizeCheck();
01011 
01012   // Copy information needed by computes and proxys to Patch::p.
01013   p.resize(numAtoms);
01014   CompAtom *p_i = p.begin();
01015   pExt.resize(numAtoms);
01016   CompAtomExt *pExt_i = pExt.begin();
01017   FullAtom *a_i = atom.begin();
01018   int i; int n = numAtoms;
01019   for ( i=0; i<n; ++i ) { 
01020     p_i[i] = a_i[i]; 
01021     pExt_i[i] = a_i[i];
01022   }
01023 
01024   // Measure atom movement to test pairlist validity
01025   doPairlistCheck();
01026 
01027   if (flags.doMolly) mollyAverage();
01028   // BEGIN LA
01029   if (flags.doLoweAndersen) loweAndersenVelocities();
01030   // END LA
01031 
01032     if (flags.doGBIS) {
01033       //reset for next time step
01034       numGBISP1Arrived = 0;
01035       phase1BoxClosedCalled = false;
01036       numGBISP2Arrived = 0;
01037       phase2BoxClosedCalled = false;
01038       numGBISP3Arrived = 0;
01039       phase3BoxClosedCalled = false;
01040       if (doMigration || isNewProxyAdded)
01041         setGBISIntrinsicRadii();
01042     }
01043 
01044   if (flags.doLCPO) {
01045     if (doMigration || isNewProxyAdded) {
01046       setLcpoType();
01047     }
01048   }
01049 
01050   // Must Add Proxy Changes when migration completed!
01051   NodeIDList::iterator pli;
01052   int *pids = NULL;
01053   int pidsPreAllocated = 1;
01054   int npid;
01055   if (proxySendSpanning == 0) {
01056     npid = proxy.size();
01057     pids = new int[npid];
01058     pidsPreAllocated = 0;
01059     int *pidi = pids;
01060     int *pide = pids + proxy.size();
01061     int patchNodesLast =
01062       ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
01063     for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
01064     {
01065       if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
01066         *(--pide) = *pli;
01067       } else {
01068         *(pidi++) = *pli;
01069       }
01070     }
01071   }
01072   else {
01073 #ifdef NODEAWARE_PROXY_SPANNINGTREE
01074     #ifdef USE_NODEPATCHMGR
01075     npid = numNodeChild;
01076     pids = nodeChildren;
01077     #else
01078     npid = nChild;
01079     pids = child;
01080     #endif
01081 #else
01082     npid = nChild;
01083     pidsPreAllocated = 0;
01084     pids = new int[proxySpanDim];
01085     for (int i=0; i<nChild; i++) pids[i] = child[i];
01086 #endif
01087   }
01088   if (npid) { //have proxies
01089     int seq = flags.sequence;
01090     int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
01091     //begin to prepare proxy msg and send it
01092     int pdMsgPLLen = p.size();
01093     int pdMsgAvgPLLen = 0;
01094     if(flags.doMolly) {
01095         pdMsgAvgPLLen = p_avg.size();
01096     }
01097     // BEGIN LA
01098     int pdMsgVLLen = 0;
01099     if (flags.doLoweAndersen) {
01100         pdMsgVLLen = v.size();
01101     }
01102     // END LA
01103 
01104     int intRadLen = 0;
01105     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
01106             intRadLen = numAtoms * 2;
01107     }
01108 
01109     //LCPO
01110     int lcpoTypeLen = 0;
01111     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
01112             lcpoTypeLen = numAtoms;
01113     }
01114 
01115     int pdMsgPLExtLen = 0;
01116     if(doMigration || isNewProxyAdded) {
01117         pdMsgPLExtLen = pExt.size();
01118     }
01119 
01120     int cudaAtomLen = 0;
01121 #ifdef NAMD_CUDA
01122     cudaAtomLen = numAtoms;
01123 #endif
01124 
01125     #ifdef NAMD_MIC
01126       #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
01127         cudaAtomLen = numAtoms;
01128       #else
01129         cudaAtomLen = (numAtoms + 15) & (~15);
01130       #endif
01131     #endif
01132 
01133     ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
01134       lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
01135 
01136     SET_PRIORITY(nmsg,seq,priority);
01137     nmsg->patch = patchID;
01138     nmsg->flags = flags;
01139     nmsg->plLen = pdMsgPLLen;                
01140     //copying data to the newly created msg
01141     memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
01142     nmsg->avgPlLen = pdMsgAvgPLLen;        
01143     if(flags.doMolly) {
01144         memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
01145     }
01146     // BEGIN LA
01147     nmsg->vlLen = pdMsgVLLen;
01148     if (flags.doLoweAndersen) {
01149         memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
01150     }
01151     // END LA
01152 
01153     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
01154       for (int i = 0; i < numAtoms * 2; i++) {
01155         nmsg->intRadList[i] = intRad[i];
01156       }
01157     }
01158 
01159     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
01160       for (int i = 0; i < numAtoms; i++) {
01161         nmsg->lcpoTypeList[i] = lcpoType[i];
01162       }
01163     }
01164 
01165     nmsg->plExtLen = pdMsgPLExtLen;
01166     if(doMigration || isNewProxyAdded){     
01167         memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
01168     }
01169 
01170 // DMK
01171 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
01172     memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
01173 #endif
01174     
01175 #if NAMD_SeparateWaters != 0
01176     //DMK - Atom Separation (water vs. non-water)
01177     nmsg->numWaterAtoms = numWaterAtoms;
01178 #endif
01179 
01180 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
01181     nmsg->isFromImmMsgCall = 0;
01182 #endif
01183     
01184     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
01185     DebugFileTrace *dft = DebugFileTrace::Object();
01186     dft->openTrace();
01187     dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
01188     for(int i=0; i<npid; i++) {
01189         dft->writeTrace("%d ", pids[i]);
01190     }
01191     dft->writeTrace("\n");
01192     dft->closeTrace();
01193     #endif
01194 
01195 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
01196     if (isProxyChanged || localphs == NULL)
01197     {
01198 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
01199      //CmiAssert(isProxyChanged);
01200      if (nphs) {
01201        for (int i=0; i<nphs; i++) {
01202          CmiDestoryPersistent(localphs[i]);
01203        }
01204        delete [] localphs;
01205      }
01206      localphs = new PersistentHandle[npid];
01207      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;
01208      for (int i=0; i<npid; i++) {
01209 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
01210        if (proxySendSpanning)
01211            localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
01212        else
01213 #endif
01214        localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
01215      }
01216      nphs = npid;
01217     }
01218     CmiAssert(nphs == npid && localphs != NULL);
01219     CmiUsePersistentHandle(localphs, nphs);
01220 #endif
01221     if(doMigration || isNewProxyAdded) {
01222         ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
01223     }else{
01224         ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
01225     }
01226 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
01227     CmiUsePersistentHandle(NULL, 0);
01228 #endif
01229     isNewProxyAdded = 0;
01230   }
01231   isProxyChanged = 0;
01232   if(!pidsPreAllocated) delete [] pids;
01233   DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
01234 
01235 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
01236   positionPtrBegin = p.begin();
01237   positionPtrEnd = p.end();
01238 #endif
01239 
01240   if(flags.doMolly) {
01241       avgPositionPtrBegin = p_avg.begin();
01242       avgPositionPtrEnd = p_avg.end();
01243   }
01244   // BEGIN LA
01245   if (flags.doLoweAndersen) {
01246       velocityPtrBegin = v.begin();
01247       velocityPtrEnd = v.end();
01248   }
01249   // END LA
01250 
01251   Patch::positionsReady(doMigration);
01252 
01253   patchMapRead = 1;
01254 
01255   // gzheng
01256   Sync::Object()->PatchReady();
01257 }
01258 
01259 void HomePatch::replaceForces(ExtForce *f)
01260 {
01261   replacementForces = f;
01262 }
01263 
01264 void HomePatch::saveForce(const int ftag)
01265 {
01266   f_saved[ftag].resize(numAtoms);
01267   for ( int i = 0; i < numAtoms; ++i )
01268   {
01269     f_saved[ftag][i] = f[ftag][i];
01270   }
01271 }
01272 
01273 
01274 #undef DEBUG_REDISTRIB_FORCE 
01275 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
01276 //#define DEBUG_REDISTRIB_FORCE
01277 
01292 void HomePatch::redistrib_colinear_lp_force(
01293     Vector& fi, Vector& fj, Vector& fk,
01294     const Vector& ri, const Vector& rj, const Vector& rk,
01295     Real distance_f, Real scale_f) {
01296   BigReal distance = distance_f;
01297   BigReal scale = scale_f;
01298   Vector r_jk = rj - rk;
01299   // TODO: Add a check for small distances?
01300   BigReal r_jk_rlength = r_jk.rlength();
01301   distance *= r_jk_rlength;
01302   BigReal fdot = distance*(fi*r_jk)*r_jk_rlength*r_jk_rlength;
01303   fj += (1. + scale + distance)*fi - r_jk*fdot;
01304   fk -= (scale + distance)*fi + r_jk*fdot;
01305 }
01306 
01332 #define FIX_FOR_WATER
01333 void HomePatch::redistrib_relative_lp_force(
01334     Vector& fi, Vector& fj, Vector& fk, Vector& fl,
01335     const Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
01336     Tensor *virial, int midpt) {
01337 #ifdef DEBUG_REDISTRIB_FORCE
01338   Vector foldnet, toldnet;  // old net force, old net torque
01339   foldnet = fi + fj + fk + fl;
01340   toldnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
01341 #endif
01342   Vector fja(0), fka(0), fla(0);
01343 
01344   Vector r = ri - rj;
01345   BigReal invr2 = r.rlength();
01346   invr2 *= invr2;
01347   BigReal fdot = (fi*r) * invr2;
01348   Vector fr = r * fdot;
01349 
01350   fja += fr;
01351 
01352   Vector s, t;
01353   if (midpt) {
01354     s = rj - 0.5*(rk + rl);
01355     t = 0.5*(rk - rl);
01356   }
01357   else {
01358     s = rj - rk;
01359     t = rk - rl;
01360   }
01361   BigReal invs2 = s.rlength();
01362   invs2 *= invs2;
01363 
01364   Vector p = cross(r,s);
01365 #if !defined(FIX_FOR_WATER)
01366   BigReal invp = p.rlength();
01367 #else
01368   BigReal p2 = p.length2();  // fix division by zero above
01369 #endif
01370 
01371   Vector q = cross(s,t);
01372   BigReal invq = q.rlength();
01373 
01374 #if !defined(FIX_FOR_WATER)
01375   BigReal fpdot = (fi*p) * invp;
01376   Vector fp = p * fpdot;
01377   Vector ft = fi - fr - fp;
01378 #else
01379   BigReal fpdot;
01380   Vector fp, ft;
01381   if (p2 < 1e-6) {  // vector is near zero, assume no fp contribution to force
01382     fpdot = 0;
01383     fp = 0;
01384     ft = fi - fr;
01385   }
01386   else {
01387     fpdot = (fi*p) / sqrt(p2);
01388     fp = p * fpdot;
01389     ft = fi - fr - fp;
01390   }
01391 #endif
01392 
01393   fja += ft;
01394   Vector v = cross(r,ft);  // torque
01395   ft = cross(s,v) * invs2;
01396   fja -= ft;
01397 
01398   if (midpt) {
01399     fka += 0.5 * ft;
01400     fla += 0.5 * ft;
01401   }
01402   else {
01403     fka += ft;
01404   }
01405 
01406   BigReal srdot = (s*r) * invs2;
01407   Vector rr = r - s*srdot;
01408   BigReal rrdot = rr.length();
01409   BigReal stdot = (s*t) * invs2;
01410   Vector tt = t - s*stdot;
01411   BigReal invtt = tt.rlength();
01412   BigReal fact = rrdot*fpdot*invtt*invq;
01413   Vector fq = q * fact;
01414 
01415   fla += fq;
01416   fja += fp*(1+srdot) + fq*stdot;
01417 
01418   ft = fq*(1+stdot) + fp*srdot;
01419 
01420   if (midpt) {
01421     fka += -0.5*ft;
01422     fla += -0.5*ft;
01423   }
01424   else {
01425     fka -= ft;
01426   }
01427 
01428   if (virial) {
01429     Tensor va = outer(fja,rj);
01430     va += outer(fka,rk);
01431     va += outer(fla,rl);
01432     va -= outer(fi,ri);
01433     *virial += va;
01434   }
01435 
01436   fi = 0;  // lone pair has zero force
01437   fj += fja;
01438   fk += fka;
01439   fl += fla;
01440 
01441 #ifdef DEBUG_REDISTRIB_FORCE
01442 #define TOL_REDISTRIB  1e-4
01443   Vector fnewnet, tnewnet;  // new net force, new net torque
01444   fnewnet = fi + fj + fk + fl;
01445   tnewnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
01446   Vector fdiff = fnewnet - foldnet;
01447   Vector tdiff = tnewnet - toldnet;
01448   if (fdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
01449     printf("Error:  force redistribution for water exceeded tolerance:  "
01450         "fdiff=(%f, %f, %f)\n", fdiff.x, fdiff.y, fdiff.z);
01451   }
01452   if (tdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
01453     printf("Error:  torque redistribution for water exceeded tolerance:  "
01454         "tdiff=(%f, %f, %f)\n", tdiff.x, tdiff.y, tdiff.z);
01455   }
01456 #endif
01457 }
01458 
01459 
01460 /* Redistribute forces from the massless lonepair charge particle onto
01461  * the other atoms of the water.
01462  *
01463  * This is done using the same algorithm as charmm uses for TIP4P lonepairs.
01464  *
01465  * Pass by reference the forces (O H1 H2 LP) to be modified,
01466  * pass by constant reference the corresponding positions,
01467  * and a pointer to virial.
01468  */
01469 void HomePatch::redistrib_lp_water_force(
01470     Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
01471     const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
01472     const Vector& p_lp, Tensor *virial) {
01473 
01474 #ifdef DEBUG_REDISTRIB_FORCE 
01475   // Debug information to check against results at end
01476 
01477   // total force and torque relative to origin
01478   Vector totforce, tottorque;
01479   totforce = f_ox + f_h1 + f_h2 + f_lp;
01480   tottorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
01481   //printf("Torque without LP is %f/%f/%f\n",
01482   //    tottorque.x, tottorque.y, tottorque.z);
01483   tottorque += cross(f_lp, p_lp);
01484   //printf("Torque with LP is %f/%f/%f\n",
01485   //    tottorque.x, tottorque.y, tottorque.z);
01486 #endif
01487 
01488   // accumulate force adjustments
01489   Vector fad_ox(0), fad_h(0);
01490 
01491   // Calculate the radial component of the force and add it to the oxygen
01492   Vector r_ox_lp = p_lp - p_ox;
01493   BigReal invlen2_r_ox_lp = r_ox_lp.rlength();
01494   invlen2_r_ox_lp *= invlen2_r_ox_lp;
01495   BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
01496   Vector f_rad = r_ox_lp * rad_factor;
01497 
01498   fad_ox += f_rad;
01499 
01500   // Calculate the angular component
01501   Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
01502   Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5; // half of r_h2_h1
01503 
01504   // deviation from collinearity of charge site
01505   //Vector r_oop = cross(r_ox_lp, r_hcom_ox);
01506   //
01507   // vector out of o-h-h plane
01508   //Vector r_perp = cross(r_hcom_ox, r_h2_h1_2);
01509 
01510   // Here we assume that Ox/Lp/Hcom are linear
01511   // If you want to correct for deviations, this is the place
01512 
01513 //  printf("Deviation from linearity for ox %i: %f/%f/%f\n", oxind, r_oop.x, r_oop.y, r_oop.z);
01514 
01515   Vector f_ang = f_lp - f_rad; // leave the angular component
01516 
01517   // now split this component onto the other atoms
01518   BigReal len_r_ox_lp = r_ox_lp.length();
01519   BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
01520   BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
01521   BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
01522 
01523   fad_ox += (f_ang * oxcomp);
01524   fad_h += (f_ang * hydcomp);  // adjustment for both hydrogens
01525 
01526   // Add virial contributions
01527   if (virial) {
01528     Tensor vir = outer(fad_ox, p_ox);
01529     vir += outer(fad_h, p_h1);
01530     vir += outer(fad_h, p_h2);
01531     vir -= outer(f_lp, p_lp);
01532     *virial += vir;
01533   }
01534 
01535   //Vector zerovec(0.0, 0.0, 0.0);
01536   f_lp = 0;
01537   f_ox += fad_ox;
01538   f_h1 += fad_h;
01539   f_h2 += fad_h;
01540 
01541 #ifdef DEBUG_REDISTRIB_FORCE 
01542   // Check that the total force and torque come out right
01543   Vector newforce, newtorque;
01544   newforce = f_ox + f_h1 + f_h2;
01545   newtorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
01546   Vector fdiff = newforce - totforce;
01547   Vector tdiff = newtorque - tottorque;
01548   BigReal error = fdiff.length();
01549   if (error > 0.0001) {
01550      printf("Error:  Force redistribution for water "
01551          "exceeded force tolerance:  error=%f\n", error);
01552   }
01553 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
01554   printf("Error in net force:  %f\n", error);
01555 #endif
01556 
01557   error = tdiff.length();
01558   if (error > 0.0001) {
01559      printf("Error:  Force redistribution for water "
01560          "exceeded torque tolerance:  error=%f\n", error);
01561   }
01562 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
01563   printf("Error in net torque:  %f\n", error);
01564 #endif
01565 #endif /* DEBUG */
01566 }
01567 
01586 void HomePatch::reposition_colinear_lonepair(
01587     Vector& ri, const Vector& rj, const Vector& rk,
01588     Real distance_f, Real scale_f)
01589 {
01590   BigReal distance = distance_f;
01591   BigReal scale = scale_f;
01592   Vector r_jk = rj - rk;
01593   BigReal r2 = r_jk.length2();
01594   if (r2 < 1e-10 || 100. < r2) { // same low tolerance as used in CHARMM
01595     iout << iWARN << "Large/small distance between lonepair reference atoms: ("
01596          << rj << ") (" << rk << ")\n" << endi;
01597   }
01598   ri = rj + (scale + distance*r_jk.rlength())*r_jk;
01599 }
01600 
01615 void HomePatch::reposition_relative_lonepair(
01616     Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
01617     Real distance, Real angle, Real dihedral)
01618 {
01619   if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
01620     iout << iWARN << "Large distance between lonepair reference atoms: ("
01621       << rj << ") (" << rk << ") (" << rl << ")\n" << endi;
01622   }
01623   BigReal r, t, p, cst, snt, csp, snp, invlen;
01624   Vector v, w, a, b, c;
01625 
01626   if (distance >= 0) {
01627     v = rk;
01628     r = distance;
01629   }
01630   else {
01631     v = 0.5*(rk + rl);
01632     r = -distance;
01633   }
01634 
01635   t = angle;
01636   p = dihedral;
01637   cst = cos(t);
01638   snt = sin(t);
01639   csp = cos(p);
01640   snp = sin(p);
01641   a = v - rj;
01642   b = rl - v;
01643   invlen = a.rlength();
01644   a *= invlen;
01645   c = cross(b, a);
01646   invlen = c.rlength();
01647   c *= invlen;
01648   b = cross(a, c);
01649   w.x = r*cst;
01650   w.y = r*snt*csp;
01651   w.z = r*snt*snp;
01652   ri.x = rj.x + w.x*a.x + w.y*b.x + w.z*c.x;
01653   ri.y = rj.y + w.x*a.y + w.y*b.y + w.z*c.y;
01654   ri.z = rj.z + w.x*a.z + w.y*b.z + w.z*c.z;
01655 }
01656 
01657 void HomePatch::reposition_all_lonepairs(void) {
01658   // ASSERT: simParams->lonepairs == TRUE
01659   for (int i=0;  i < numAtoms;  i++) {
01660     if (atom[i].mass < 0.01) {
01661       // found a lone pair
01662       AtomID aid = atom[i].id;  // global atom ID of lp
01663       Lphost *lph = Node::Object()->molecule->get_lphost(aid);  // its lphost
01664       if (lph == NULL) {
01665         char errmsg[512];
01666         sprintf(errmsg, "reposition lone pairs: "
01667             "no Lphost exists for LP %d\n", aid);
01668         NAMD_die(errmsg);
01669       }
01670       LocalID j = AtomMap::Object()->localID(lph->atom2);
01671       LocalID k = AtomMap::Object()->localID(lph->atom3);
01672       LocalID l = AtomMap::Object()->localID(lph->atom4);
01673       if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
01674         char errmsg[512];
01675         sprintf(errmsg, "reposition lone pairs: "
01676             "LP %d has some Lphost atom off patch\n", aid);
01677         NAMD_die(errmsg);
01678       }
01679       // reposition this lone pair
01680       if (lph->numhosts == 2) {
01681         reposition_colinear_lonepair(atom[i].position, atom[j.index].position,
01682             atom[k.index].position, lph->distance, lph->angle);
01683       }
01684       else if (lph->numhosts == 3) {
01685         reposition_relative_lonepair(atom[i].position, atom[j.index].position,
01686             atom[k.index].position, atom[l.index].position,
01687             lph->distance, lph->angle, lph->dihedral);
01688       }
01689     }
01690   }
01691 }
01692 
01693 void HomePatch::swm4_omrepos(Vector *ref, Vector *pos, Vector *vel,
01694     BigReal invdt) {
01695   // Reposition lonepair (Om) particle of Drude SWM4 water.
01696   // Same comments apply as to tip4_omrepos(), but the ordering of atoms
01697   // is different: O, D, LP, H1, H2.
01698   pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
01699   // Now, adjust velocity of particle to get it to appropriate place
01700   // during next integration "drift-step"
01701   if (invdt != 0) {
01702     vel[2] = (pos[2] - ref[2]) * invdt;
01703   }
01704   // No virial correction needed since lonepair is massless
01705 }
01706 
01707 void HomePatch::tip4_omrepos(Vector* ref, Vector* pos, Vector* vel, BigReal invdt) {
01708   /* Reposition the om particle of a tip4p water
01709    * A little geometry shows that the appropriate position is given by
01710    * R_O + (1 / 2 r_ohc) * ( 0.5 (R_H1 + R_H2) - R_O ) 
01711    * Here r_om is the distance from the oxygen to Om site, and r_ohc
01712    * is the altitude from the oxygen to the hydrogen center of mass
01713    * Those quantities are precalculated upon initialization of HomePatch
01714    *
01715    * Ordering of TIP4P atoms: O, H1, H2, LP.
01716    */
01717 
01718   //printf("rom/rohc are %f %f and invdt is %f\n", r_om, r_ohc, invdt);
01719   //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);
01720   pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc); 
01721   //printf("New position for lp is %f %f %f\n", pos[3].x, pos[3].y, pos[3].z);
01722 
01723   // Now, adjust the velocity of the particle to get it to the appropriate place
01724   if (invdt != 0) {
01725     vel[3] = (pos[3] - ref[3]) * invdt;
01726   }
01727 
01728   // No virial correction needed, since this is a massless particle
01729   return;
01730 }
01731 
01732 void HomePatch::redistrib_lonepair_forces(const int ftag, Tensor *virial) {
01733   // ASSERT: simParams->lonepairs == TRUE
01734   ForceList *f_mod = f;
01735   for (int i = 0;  i < numAtoms;  i++) {
01736     if (atom[i].mass < 0.01) {
01737       // found a lone pair
01738       AtomID aid = atom[i].id;  // global atom ID of lp
01739       Lphost *lph = Node::Object()->molecule->get_lphost(aid);  // its lphost
01740       if (lph == NULL) {
01741         char errmsg[512];
01742         sprintf(errmsg, "redistrib lone pair forces: "
01743             "no Lphost exists for LP %d\n", aid);
01744         NAMD_die(errmsg);
01745       }
01746       LocalID j = AtomMap::Object()->localID(lph->atom2);
01747       LocalID k = AtomMap::Object()->localID(lph->atom3);
01748       LocalID l = AtomMap::Object()->localID(lph->atom4);
01749       if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
01750         char errmsg[512];
01751         sprintf(errmsg, "redistrib lone pair forces: "
01752             "LP %d has some Lphost atom off patch\n", aid);
01753         NAMD_die(errmsg);
01754       }
01755       // redistribute forces from this lone pair
01756       if (lph->numhosts == 2) {
01757         redistrib_colinear_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
01758             f_mod[ftag][k.index], atom[i].position, atom[j.index].position,
01759             atom[k.index].position, lph->distance, lph->angle);
01760       }
01761       else if (lph->numhosts == 3) {
01762         int midpt = (lph->distance < 0);
01763         redistrib_relative_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
01764             f_mod[ftag][k.index], f_mod[ftag][l.index],
01765             atom[i].position, atom[j.index].position,
01766             atom[k.index].position, atom[l.index].position, virial, midpt);
01767       }
01768     }
01769   }
01770 }
01771 
01772 void HomePatch::redistrib_swm4_forces(const int ftag, Tensor *virial) {
01773   // Loop over the patch's atoms and apply the appropriate corrections
01774   // to get all forces off of lone pairs
01775   ForceList *f_mod = f;
01776   for (int i = 0;  i < numAtoms;  i++) {
01777     if (atom[i].mass < 0.01) {
01778       // found lonepair
01779       redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
01780           f_mod[ftag][i+2], f_mod[ftag][i],
01781           atom[i-2].position, atom[i+1].position,
01782           atom[i+2].position, atom[i].position, virial);
01783     }
01784   }
01785 }
01786 
01787 void HomePatch::redistrib_tip4p_forces(const int ftag, Tensor* virial) {
01788   // Loop over the patch's atoms and apply the appropriate corrections
01789   // to get all forces off of lone pairs
01790   // Atom ordering:  O H1 H2 LP
01791 
01792   ForceList *f_mod =f;
01793   for (int i=0; i<numAtoms; i++) {
01794     if (atom[i].mass < 0.01) {
01795       // found lonepair
01796       redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
01797           f_mod[ftag][i-1], f_mod[ftag][i],
01798           atom[i-3].position, atom[i-2].position,
01799           atom[i-1].position, atom[i].position, virial);
01800     }
01801   }
01802 }
01803 
01804 
01805 void HomePatch::addForceToMomentum(
01806     FullAtom       * __restrict atom_arr,
01807     const Force    * __restrict force_arr,
01808     const BigReal    dt,
01809     int              num_atoms
01810     ) {
01811   SimParameters *simParams = Node::Object()->simParameters;
01812 
01813   if ( simParams->fixedAtomsOn ) {
01814     for ( int i = 0; i < num_atoms; ++i ) {
01815       if ( atom_arr[i].atomFixed ) {
01816         atom_arr[i].velocity = 0;
01817       } else {
01818         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01819         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01820         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01821         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01822       }
01823     }
01824   } else {
01825     for ( int i = 0; i < num_atoms; ++i ) {
01826       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01827       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01828       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01829       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01830     }
01831   }
01832 }
01833 
01834 void HomePatch::addForceToMomentum3(
01835     FullAtom       * __restrict atom_arr,
01836     const Force    * __restrict force_arr1,
01837     const Force    * __restrict force_arr2,
01838     const Force    * __restrict force_arr3,
01839     const BigReal    dt1,
01840     const BigReal    dt2,
01841     const BigReal    dt3,
01842     int              num_atoms
01843     ) {
01844   SimParameters *simParams = Node::Object()->simParameters;
01845 
01846   if ( simParams->fixedAtomsOn ) {
01847     for ( int i = 0; i < num_atoms; ++i ) {
01848       if ( atom_arr[i].atomFixed ) {
01849         atom_arr[i].velocity = 0;
01850       } else {
01851         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01852         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01853             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01854         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01855             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01856         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01857             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01858       }
01859     }
01860   } else {
01861     for ( int i = 0; i < num_atoms; ++i ) {
01862       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01863       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01864           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01865       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01866           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01867       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01868           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01869     }
01870   }
01871 }
01872 
01873 void HomePatch::addVelocityToPosition(
01874     FullAtom       * __restrict atom_arr,
01875     const BigReal    dt,
01876     int              num_atoms
01877     ) {
01878   SimParameters *simParams = Node::Object()->simParameters;
01879   if ( simParams->fixedAtomsOn ) {
01880     for ( int i = 0; i < num_atoms; ++i ) {
01881       if ( ! atom_arr[i].atomFixed ) {
01882         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01883         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01884         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01885       }
01886     }
01887   } else {
01888     for ( int i = 0; i < num_atoms; ++i ) {
01889       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01890       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01891       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01892     }
01893   }
01894 }
01895 
01896 int HomePatch::hardWallDrude(const BigReal timestep, Tensor *virial,
01897     SubmitReduction *ppreduction)
01898 {
01899   Molecule *mol = Node::Object()->molecule;
01900   SimParameters *simParams = Node::Object()->simParameters;
01901   const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
01902   const int fixedAtomsOn = simParams->fixedAtomsOn;
01903   const BigReal dt = timestep / TIMEFACTOR;
01904   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
01905   int i, ia, ib, j;
01906   int dieOnError = simParams->rigidDie;
01907   Tensor wc;  // constraint virial
01908   BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
01909   int nslabs;
01910 
01911   // start data for hard wall boundary between drude and its host atom
01912   // static int Count=0;
01913   int Idx;
01914   double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
01915   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;
01916   double dot_v_r_1, dot_v_r_2;
01917   double vb_cm, dr_a, dr_b;
01918   // end data for hard wall boundary between drude and its host atom
01919 
01920   // start calculation of hard wall boundary between drude and its host atom
01921   if (simParams->drudeHardWallOn) {
01922     if (ppreduction) {
01923       nslabs = simParams->pressureProfileSlabs;
01924       idz = nslabs/lattice.c().z;
01925       zmin = lattice.origin().z - 0.5*lattice.c().z;
01926     }
01927 
01928     r_wall = simParams->drudeBondLen;
01929     r_wall_SQ = r_wall*r_wall;
01930     // Count++;
01931     for (i=1; i<numAtoms; i++)  {
01932       if ( (atom[i].mass > 0.01) && ((atom[i].mass < 1.0)) ) { // drude particle
01933         ia = i-1;
01934         ib = i;
01935 
01936         v_ab = atom[ib].position - atom[ia].position;
01937         rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
01938 
01939         if (rab_SQ > r_wall_SQ) {  // to impose the hard wall constraint
01940           rab = sqrt(rab_SQ);
01941           if ( (rab > (2.0*r_wall)) && dieOnError ) {  // unexpected situation
01942             iout << iERROR << "HardWallDrude> "
01943               << "The drude is too far away from atom "
01944               << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
01945             return -1;  // triggers early exit
01946           }
01947 
01948           v_ab.x /= rab;
01949           v_ab.y /= rab;
01950           v_ab.z /= rab;
01951 
01952           if ( fixedAtomsOn && atom[ia].atomFixed ) {  // the heavy atom is fixed
01953             if (atom[ib].atomFixed) {  // the drude is fixed too
01954               continue;
01955             }
01956             else {  // only the heavy atom is fixed
01957               dot_v_r_2 = atom[ib].velocity.x*v_ab.x
01958                 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
01959               vb_2 = dot_v_r_2 * v_ab;
01960               vp_2 = atom[ib].velocity - vb_2;
01961 
01962               dr = rab - r_wall;
01963               if(dot_v_r_2 == 0.0) {
01964                 delta_T = maxtime;
01965               }
01966               else {
01967                 delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
01968                 if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
01969               }
01970 
01971               dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
01972 
01973               vb_2 = dot_v_r_2 * v_ab;
01974 
01975               new_vel_a = atom[ia].velocity;
01976               new_vel_b = vp_2 + vb_2;
01977 
01978               dr_b = -dr + delta_T*dot_v_r_2;  // L = L_0 + dT *v_new, v was flipped
01979 
01980               new_pos_a = atom[ia].position;
01981               new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
01982             }
01983           }
01984           else {
01985             mass_a = atom[ia].mass;
01986             mass_b = atom[ib].mass;
01987             mass_sum = mass_a+mass_b;
01988 
01989             dot_v_r_1 = atom[ia].velocity.x*v_ab.x
01990               + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
01991             vb_1 = dot_v_r_1 * v_ab;
01992             vp_1 = atom[ia].velocity - vb_1;
01993 
01994             dot_v_r_2 = atom[ib].velocity.x*v_ab.x
01995               + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
01996             vb_2 = dot_v_r_2 * v_ab;
01997             vp_2 = atom[ib].velocity - vb_2;
01998 
01999             vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
02000 
02001             dot_v_r_1 -= vb_cm;
02002             dot_v_r_2 -= vb_cm;
02003 
02004             dr = rab - r_wall;
02005 
02006             if(dot_v_r_2 == dot_v_r_1) {
02007                 delta_T = maxtime;
02008             }
02009             else {
02010               delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);  // the time since the collision occurs
02011               if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
02012             }
02013             
02014             // the relative velocity between ia and ib. Drawn according to T_Drude
02015             v_Bond = sqrt(kbt/mass_b);
02016 
02017             // reflect the velocity along bond vector and scale down
02018             dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
02019             dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
02020 
02021             dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
02022             dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
02023 
02024             new_pos_a = atom[ia].position + dr_a*v_ab;  // correct the position
02025             new_pos_b = atom[ib].position + dr_b*v_ab;
02026             // atom[ia].position += (dr_a*v_ab);  // correct the position
02027             // atom[ib].position += (dr_b*v_ab);
02028             
02029             dot_v_r_1 += vb_cm;
02030             dot_v_r_2 += vb_cm;
02031 
02032             vb_1 = dot_v_r_1 * v_ab;
02033             vb_2 = dot_v_r_2 * v_ab;
02034 
02035             new_vel_a = vp_1 + vb_1;
02036             new_vel_b = vp_2 + vb_2;
02037           }
02038 
02039           int ppoffset, partition;
02040           if ( invdt == 0 ) {
02041             atom[ia].position = new_pos_a;
02042             atom[ib].position = new_pos_b;
02043           }
02044           else if ( virial == 0 ) {
02045             atom[ia].velocity = new_vel_a;
02046             atom[ib].velocity = new_vel_b;
02047           }
02048           else {
02049             for ( j = 0; j < 2; j++ ) {
02050               if (j==0) {  // atom ia, heavy atom
02051                 Idx = ia;
02052                 new_pos = &new_pos_a;
02053                 new_vel = &new_vel_a;
02054               }
02055               else if (j==1) {  // atom ib, drude
02056                 Idx = ib;
02057                 new_pos = &new_pos_b;
02058                 new_vel = &new_vel_b;
02059               }
02060               Force df = (*new_vel - atom[Idx].velocity) *
02061                 ( atom[Idx].mass * invdt );
02062               Tensor vir = outer(df, atom[Idx].position);
02063               wc += vir;
02064               atom[Idx].velocity = *new_vel;
02065               atom[Idx].position = *new_pos;
02066 
02067               if (ppreduction) {
02068                 if (!j) {
02069                   BigReal z = new_pos->z;
02070                   int partition = atom[Idx].partition;
02071                   int slab = (int)floor((z-zmin)*idz);
02072                   if (slab < 0) slab += nslabs;
02073                   else if (slab >= nslabs) slab -= nslabs;
02074                   ppoffset = 3*(slab + nslabs*partition);
02075                 }
02076                 ppreduction->item(ppoffset  ) += vir.xx;
02077                 ppreduction->item(ppoffset+1) += vir.yy;
02078                 ppreduction->item(ppoffset+2) += vir.zz;
02079               }
02080 
02081             }
02082           }
02083         }                               
02084       }
02085     }
02086 
02087     // if ( (Count>10000) && (Count%10==0) ) {
02088     //   v_ab = atom[1].position - atom[0].position;
02089     //   rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
02090     //   iout << "DBG_R: " << Count << "  " << sqrt(rab_SQ) << "\n" << endi;
02091     // }
02092 
02093   }
02094 
02095   // end calculation of hard wall boundary between drude and its host atom
02096 
02097   if ( dt && virial ) *virial += wc;
02098 
02099   return 0;
02100 }
02101 
02102 void HomePatch::buildRattleList() {
02103 
02104   SimParameters *simParams = Node::Object()->simParameters;
02105   const int fixedAtomsOn = simParams->fixedAtomsOn;
02106   const int useSettle = simParams->useSettle;
02107 
02108   // Re-size to containt numAtoms elements
02109   velNew.resize(numAtoms);
02110   posNew.resize(numAtoms);
02111 
02112   // Size of a hydrogen group for water
02113   int wathgsize = 3;
02114   int watmodel = simParams->watmodel;
02115   if (watmodel == WAT_TIP4) wathgsize = 4;
02116   else if (watmodel == WAT_SWM4) wathgsize = 5;
02117 
02118   // Initialize the settle algorithm with water parameters
02119   // settle1() assumes all waters are identical,
02120   // and will generate bad results if they are not.
02121   // XXX this will move to Molecule::build_atom_status when that 
02122   // version is debugged
02123   if ( ! settle_initialized ) {
02124     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02125       // find a water
02126       if (atom[ig].rigidBondLength > 0) {
02127         int oatm;
02128         if (simParams->watmodel == WAT_SWM4) {
02129           oatm = ig+3;  // skip over Drude and Lonepair
02130           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02131           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02132         }
02133         else {
02134           oatm = ig+1;
02135           // Avoid using the Om site to set this by mistake
02136           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02137             oatm += 1;
02138           }
02139         }
02140 
02141         // initialize settle water parameters
02142         settle1init(atom[ig].mass, atom[oatm].mass, 
02143                     atom[ig].rigidBondLength,
02144                     atom[oatm].rigidBondLength,
02145                     settle_mOrmT, settle_mHrmT, settle_ra,
02146                     settle_rb, settle_rc, settle_rra);
02147         settle_initialized = 1;
02148         break; // done with init
02149       }
02150     }
02151   }
02152   
02153   Vector ref[10];
02154   BigReal rmass[10];
02155   BigReal dsq[10];
02156   int fixed[10];
02157   int ial[10];
02158   int ibl[10];
02159 
02160   int numSettle = 0;
02161   int numRattle = 0;
02162   int posRattleParam = 0;
02163 
02164   settleList.clear();
02165   rattleList.clear();
02166   noconstList.clear();
02167   rattleParam.clear();
02168 
02169   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02170     int hgs = atom[ig].hydrogenGroupSize;
02171     if ( hgs == 1 ) {
02172       // only one atom in group
02173       noconstList.push_back(ig);
02174       continue;
02175     }
02176     int anyfixed = 0;
02177     for (int i = 0; i < hgs; ++i ) {
02178       ref[i] = atom[ig+i].position;
02179       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02180       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02181       if ( fixed[i] ) {
02182         anyfixed = 1;
02183         rmass[i] = 0.;
02184       }
02185     }
02186     int icnt = 0;
02187     BigReal tmp = atom[ig].rigidBondLength;
02188     if (tmp > 0.0) {  // for water
02189       if (hgs != wathgsize) {
02190         char errmsg[256];
02191         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02192                          "but the specified water model requires %d atoms.\n",
02193                          atom[ig].id+1, hgs, wathgsize);
02194         NAMD_die(errmsg);
02195       }
02196       // Use SETTLE for water unless some of the water atoms are fixed,
02197       if (useSettle && !anyfixed) {
02198         // Store to Settle -list
02199         settleList.push_back(ig);
02200         continue;
02201       }
02202       if ( !(fixed[1] && fixed[2]) ) {
02203         dsq[icnt] = tmp * tmp;
02204         ial[icnt] = 1;
02205         ibl[icnt] = 2;
02206         ++icnt;
02207       }
02208     } // if (tmp > 0.0)
02209     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02210       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02211         if ( !(fixed[0] && fixed[i]) ) {
02212           dsq[icnt] = tmp * tmp;
02213           ial[icnt] = 0;
02214           ibl[icnt] = i;
02215           ++icnt;
02216         }
02217       }
02218     }
02219     if ( icnt == 0 ) {
02220       // no constraints
02221       noconstList.push_back(ig);
02222       continue;  
02223     }
02224     // Store to Rattle -list
02225     RattleList rattleListElem;
02226     rattleListElem.ig  = ig;
02227     rattleListElem.icnt = icnt;
02228     rattleList.push_back(rattleListElem);
02229     for (int i = 0; i < icnt; ++i ) {
02230       int a = ial[i];
02231       int b = ibl[i];
02232       RattleParam rattleParamElem;
02233       rattleParamElem.ia = a;
02234       rattleParamElem.ib = b;
02235       rattleParamElem.dsq = dsq[i];
02236       rattleParamElem.rma = rmass[a];
02237       rattleParamElem.rmb = rmass[b];
02238       rattleParam.push_back(rattleParamElem);
02239     }
02240   }
02241 
02242 }
02243 
02244 void HomePatch::addRattleForce(const BigReal invdt, Tensor& wc) {
02245   for (int ig = 0; ig < numAtoms; ++ig ) {
02246     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
02247     Tensor vir = outer(df, atom[ig].position);
02248     wc += vir;
02249     f[Results::normal][ig] += df;
02250     atom[ig].velocity = velNew[ig];
02251   }
02252 }
02253 
02254 int HomePatch::rattle1(const BigReal timestep, Tensor *virial, 
02255   SubmitReduction *ppreduction) {
02256 
02257   SimParameters *simParams = Node::Object()->simParameters;
02258   if (simParams->watmodel != WAT_TIP3 || ppreduction) {
02259     // Call old rattle1 -method instead
02260     return rattle1old(timestep, virial, ppreduction);
02261   }
02262 
02263   if (!rattleListValid) {
02264     buildRattleList();
02265     rattleListValid = true;
02266   }
02267 
02268   const int fixedAtomsOn = simParams->fixedAtomsOn;
02269   const int useSettle = simParams->useSettle;
02270   const BigReal dt = timestep / TIMEFACTOR;
02271   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02272   const BigReal tol2 = 2.0 * simParams->rigidTol;
02273   int maxiter = simParams->rigidIter;
02274   int dieOnError = simParams->rigidDie;
02275 
02276   Vector ref[10];  // reference position
02277   Vector pos[10];  // new position
02278   Vector vel[10];  // new velocity
02279 
02280   // Manual un-roll
02281   int n = (settleList.size()/2)*2;
02282   for (int j=0;j < n;j+=2) {
02283     int ig;
02284     ig = settleList[j];
02285     for (int i = 0; i < 3; ++i ) {
02286       ref[i] = atom[ig+i].position;
02287       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02288     }
02289     ig = settleList[j+1];
02290     for (int i = 0; i < 3; ++i ) {
02291       ref[i+3] = atom[ig+i].position;
02292       pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
02293     }
02294     settle1_SIMD<2>(ref, pos,
02295       settle_mOrmT, settle_mHrmT, settle_ra,
02296       settle_rb, settle_rc, settle_rra);
02297 
02298     ig = settleList[j];
02299     for (int i = 0; i < 3; ++i ) {
02300       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02301       posNew[ig+i] = pos[i];
02302     }
02303     ig = settleList[j+1];
02304     for (int i = 0; i < 3; ++i ) {
02305       velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
02306       posNew[ig+i] = pos[i+3];
02307     }
02308 
02309   }
02310 
02311   if (settleList.size() % 2) {
02312     int ig = settleList[settleList.size()-1];
02313     for (int i = 0; i < 3; ++i ) {
02314       ref[i] = atom[ig+i].position;
02315       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02316     }
02317     settle1_SIMD<1>(ref, pos,
02318             settle_mOrmT, settle_mHrmT, settle_ra,
02319             settle_rb, settle_rc, settle_rra);        
02320     for (int i = 0; i < 3; ++i ) {
02321       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02322       posNew[ig+i] = pos[i];
02323     }
02324   }
02325 
02326   int posParam = 0;
02327   for (int j=0;j < rattleList.size();++j) {
02328 
02329     BigReal refx[10];
02330     BigReal refy[10];
02331     BigReal refz[10];
02332 
02333     BigReal posx[10];
02334     BigReal posy[10];
02335     BigReal posz[10];
02336 
02337     int ig = rattleList[j].ig;
02338     int icnt = rattleList[j].icnt;
02339     int hgs = atom[ig].hydrogenGroupSize;
02340     for (int i = 0; i < hgs; ++i ) {
02341       ref[i] = atom[ig+i].position;
02342       pos[i] = atom[ig+i].position;
02343       if (!(fixedAtomsOn && atom[ig+i].atomFixed))
02344         pos[i] += atom[ig+i].velocity * dt;
02345       refx[i] = ref[i].x;
02346       refy[i] = ref[i].y;
02347       refz[i] = ref[i].z;
02348       posx[i] = pos[i].x;
02349       posy[i] = pos[i].y;
02350       posz[i] = pos[i].z;
02351     }
02352 
02353 
02354     bool done;
02355     bool consFailure;
02356     if (icnt == 1) {
02357       rattlePair<1>(&rattleParam[posParam],
02358         refx, refy, refz,
02359         posx, posy, posz,
02360         consFailure);
02361       done = true;
02362     } else {
02363       rattleN(icnt, &rattleParam[posParam],
02364         refx, refy, refz,
02365         posx, posy, posz,
02366         tol2, maxiter,
02367         done, consFailure);
02368     }
02369 
02370     // Advance position in rattleParam
02371     posParam += icnt;
02372 
02373     for (int i = 0; i < hgs; ++i ) {
02374       pos[i].x = posx[i];
02375       pos[i].y = posy[i];
02376       pos[i].z = posz[i];
02377     }
02378 
02379     for (int i = 0; i < hgs; ++i ) {
02380       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02381       posNew[ig+i] = pos[i];
02382     }
02383 
02384     if ( consFailure ) {
02385       if ( dieOnError ) {
02386         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02387         << (atom[ig].id + 1) << "!\n" << endi;
02388         return -1;  // triggers early exit
02389       } else {
02390         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02391         << (atom[ig].id + 1) << "!\n" << endi;
02392       }
02393     } else if ( ! done ) {
02394       if ( dieOnError ) {
02395         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02396         << (atom[ig].id + 1) << "!\n" << endi;
02397         return -1;  // triggers early exit
02398       } else {
02399         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02400         << (atom[ig].id + 1) << "!\n" << endi;
02401       }
02402     }
02403 
02404   }
02405 
02406   // Finally, we have to go through atoms that are not involved in rattle just so that we have
02407   // their positions and velocities up-to-date in posNew and velNew
02408   for (int j=0;j < noconstList.size();++j) {
02409     int ig = noconstList[j];
02410     int hgs = atom[ig].hydrogenGroupSize;
02411     for (int i = 0; i < hgs; ++i ) {
02412       velNew[ig+i] = atom[ig+i].velocity;
02413       posNew[ig+i] = atom[ig+i].position;
02414     }
02415   }
02416 
02417   if ( invdt == 0 ) {
02418     for (int ig = 0; ig < numAtoms; ++ig )
02419       atom[ig].position = posNew[ig];
02420   } else if ( virial == 0 ) {
02421     for (int ig = 0; ig < numAtoms; ++ig )
02422       atom[ig].velocity = velNew[ig];
02423   } else {
02424     Tensor wc;  // constraint virial
02425     addRattleForce(invdt, wc);
02426     *virial += wc;
02427   }
02428 
02429   return 0;
02430 }
02431 
02432 //  RATTLE algorithm from Allen & Tildesley
02433 int HomePatch::rattle1old(const BigReal timestep, Tensor *virial, 
02434     SubmitReduction *ppreduction)
02435 {
02436   Molecule *mol = Node::Object()->molecule;
02437   SimParameters *simParams = Node::Object()->simParameters;
02438   const int fixedAtomsOn = simParams->fixedAtomsOn;
02439   const int useSettle = simParams->useSettle;
02440   const BigReal dt = timestep / TIMEFACTOR;
02441   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02442   BigReal tol2 = 2.0 * simParams->rigidTol;
02443   int maxiter = simParams->rigidIter;
02444   int dieOnError = simParams->rigidDie;
02445   int i, iter;
02446   BigReal dsq[10], tmp;
02447   int ial[10], ibl[10];
02448   Vector ref[10];  // reference position
02449   Vector refab[10];  // reference vector
02450   Vector pos[10];  // new position
02451   Vector vel[10];  // new velocity
02452   Vector netdp[10];  // total momentum change from constraint
02453   BigReal rmass[10];  // 1 / mass
02454   int fixed[10];  // is atom fixed?
02455   Tensor wc;  // constraint virial
02456   BigReal idz, zmin;
02457   int nslabs;
02458 
02459   // Initialize the settle algorithm with water parameters
02460   // settle1() assumes all waters are identical,
02461   // and will generate bad results if they are not.
02462   // XXX this will move to Molecule::build_atom_status when that 
02463   // version is debugged
02464   if ( ! settle_initialized ) {
02465     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02466       // find a water
02467       if (atom[ig].rigidBondLength > 0) {
02468         int oatm;
02469         if (simParams->watmodel == WAT_SWM4) {
02470           oatm = ig+3;  // skip over Drude and Lonepair
02471           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02472           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02473         }
02474         else {
02475           oatm = ig+1;
02476           // Avoid using the Om site to set this by mistake
02477           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02478             oatm += 1;
02479           }
02480         }
02481 
02482         // initialize settle water parameters
02483         settle1init(atom[ig].mass, atom[oatm].mass, 
02484                     atom[ig].rigidBondLength,
02485                     atom[oatm].rigidBondLength,
02486                     settle_mOrmT, settle_mHrmT, settle_ra,
02487                     settle_rb, settle_rc, settle_rra);
02488         settle_initialized = 1;
02489         break; // done with init
02490       }
02491     }
02492   }
02493 
02494   if (ppreduction) {
02495     nslabs = simParams->pressureProfileSlabs;
02496     idz = nslabs/lattice.c().z;
02497     zmin = lattice.origin().z - 0.5*lattice.c().z;
02498   }
02499 
02500   // Size of a hydrogen group for water
02501   int wathgsize = 3;
02502   int watmodel = simParams->watmodel;
02503   if (watmodel == WAT_TIP4) wathgsize = 4;
02504   else if (watmodel == WAT_SWM4) wathgsize = 5;
02505   
02506   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02507     int hgs = atom[ig].hydrogenGroupSize;
02508     if ( hgs == 1 ) continue;  // only one atom in group
02509     // cache data in local arrays and integrate positions normally
02510     int anyfixed = 0;
02511     for ( i = 0; i < hgs; ++i ) {
02512       ref[i] = atom[ig+i].position;
02513       pos[i] = atom[ig+i].position;
02514       vel[i] = atom[ig+i].velocity;
02515       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02516       //printf("rmass of %i is %f\n", ig+i, rmass[i]);
02517       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02518       //printf("fixed status of %i is %i\n", i, fixed[i]);
02519       // undo addVelocityToPosition to get proper reference coordinates
02520       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
02521     }
02522     int icnt = 0;
02523     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02524       if (hgs != wathgsize) {
02525         char errmsg[256];
02526         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02527                          "but the specified water model requires %d atoms.\n",
02528                          atom[ig].id+1, hgs, wathgsize);
02529         NAMD_die(errmsg);
02530       }
02531       // Use SETTLE for water unless some of the water atoms are fixed,
02532       if (useSettle && !anyfixed) {
02533         if (simParams->watmodel == WAT_SWM4) {
02534           // SWM4 ordering:  O D LP H1 H2
02535           // do swap(O,LP) and call settle with subarray O H1 H2
02536           // swap back after we return
02537           Vector lp_ref = ref[2];
02538           Vector lp_pos = pos[2];
02539           Vector lp_vel = vel[2];
02540           ref[2] = ref[0];
02541           pos[2] = pos[0];
02542           vel[2] = vel[0];
02543           settle1(ref+2, pos+2, vel+2, invdt,
02544                   settle_mOrmT, settle_mHrmT, settle_ra,
02545                   settle_rb, settle_rc, settle_rra);
02546           ref[0] = ref[2];
02547           pos[0] = pos[2];
02548           vel[0] = vel[2];
02549           ref[2] = lp_ref;
02550           pos[2] = lp_pos;
02551           vel[2] = lp_vel;
02552           // determine for LP updated pos and vel
02553           swm4_omrepos(ref, pos, vel, invdt);
02554         }
02555         else {
02556             settle1(ref, pos, vel, invdt,
02557                     settle_mOrmT, settle_mHrmT, settle_ra,
02558                     settle_rb, settle_rc, settle_rra);            
02559           if (simParams->watmodel == WAT_TIP4) {
02560             tip4_omrepos(ref, pos, vel, invdt);
02561           }
02562         }
02563 
02564         // which slab the hydrogen group will belong to
02565         // for pprofile calculations.
02566         int ppoffset, partition;
02567         if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02568           atom[ig+i].position = pos[i];
02569         } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02570           atom[ig+i].velocity = vel[i];
02571         } else for ( i = 0; i < wathgsize; ++i ) {
02572           Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
02573           Tensor vir = outer(df, ref[i]);
02574           wc += vir;
02575           f[Results::normal][ig+i] += df;
02576           atom[ig+i].velocity = vel[i];
02577           if (ppreduction) {
02578             // put all the atoms from a water in the same slab.  Atom 0
02579             // should be the parent atom.
02580             if (!i) {
02581               BigReal z = pos[i].z;
02582               partition = atom[ig].partition;
02583               int slab = (int)floor((z-zmin)*idz);
02584               if (slab < 0) slab += nslabs;
02585               else if (slab >= nslabs) slab -= nslabs;
02586               ppoffset = 3*(slab + nslabs*partition);
02587             }
02588             ppreduction->item(ppoffset  ) += vir.xx;
02589             ppreduction->item(ppoffset+1) += vir.yy;
02590             ppreduction->item(ppoffset+2) += vir.zz;
02591           }
02592         }
02593         continue;
02594       }
02595       if ( !(fixed[1] && fixed[2]) ) {
02596         dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02597       }
02598     }
02599     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02600       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02601         if ( !(fixed[0] && fixed[i]) ) {
02602           dsq[icnt] = tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
02603         }
02604       }
02605     }
02606     if ( icnt == 0 ) continue;  // no constraints
02607     for ( i = 0; i < icnt; ++i ) {
02608       refab[i] = ref[ial[i]] - ref[ibl[i]];
02609     }
02610     for ( i = 0; i < hgs; ++i ) {
02611       netdp[i] = 0.;
02612     }
02613     int done;
02614     int consFailure;
02615     for ( iter = 0; iter < maxiter; ++iter ) {
02616 //if (iter > 0) CkPrintf("iteration %d\n", iter);
02617       done = 1;
02618       consFailure = 0;
02619       for ( i = 0; i < icnt; ++i ) {
02620         int a = ial[i];  int b = ibl[i];
02621         Vector pab = pos[a] - pos[b];
02622               BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
02623         BigReal rabsq = dsq[i];
02624         BigReal diffsq = rabsq - pabsq;
02625         if ( fabs(diffsq) > (rabsq * tol2) ) {
02626           Vector &rab = refab[i];
02627           BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
02628           if ( rpab < ( rabsq * 1.0e-6 ) ) {
02629             done = 0;
02630             consFailure = 1;
02631             continue;
02632           }
02633           BigReal rma = rmass[a];
02634           BigReal rmb = rmass[b];
02635           BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
02636           Vector dp = rab * gab;
02637           pos[a] += rma * dp;
02638           pos[b] -= rmb * dp;
02639           if ( invdt != 0. ) {
02640             dp *= invdt;
02641             netdp[a] += dp;
02642             netdp[b] -= dp;
02643           }
02644           done = 0;
02645         }
02646       }
02647       if ( done ) break;
02648     }
02649 
02650     if ( consFailure ) {
02651       if ( dieOnError ) {
02652         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02653         << (atom[ig].id + 1) << "!\n" << endi;
02654         return -1;  // triggers early exit
02655       } else {
02656         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02657         << (atom[ig].id + 1) << "!\n" << endi;
02658       }
02659     } else if ( ! done ) {
02660       if ( dieOnError ) {
02661         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02662         << (atom[ig].id + 1) << "!\n" << endi;
02663         return -1;  // triggers early exit
02664       } else {
02665         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02666         << (atom[ig].id + 1) << "!\n" << endi;
02667       }
02668     }
02669 
02670     // store data back to patch
02671     int ppoffset, partition;
02672     if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
02673       atom[ig+i].position = pos[i];
02674     } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
02675       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02676     } else for ( i = 0; i < hgs; ++i ) {
02677       Force df = netdp[i] * invdt;
02678       Tensor vir = outer(df, ref[i]);
02679       wc += vir;
02680       f[Results::normal][ig+i] += df;
02681       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02682       if (ppreduction) {
02683         if (!i) {
02684           BigReal z = pos[i].z;
02685           int partition = atom[ig].partition;
02686           int slab = (int)floor((z-zmin)*idz);
02687           if (slab < 0) slab += nslabs;
02688           else if (slab >= nslabs) slab -= nslabs;
02689           ppoffset = 3*(slab + nslabs*partition);
02690         }
02691         ppreduction->item(ppoffset  ) += vir.xx;
02692         ppreduction->item(ppoffset+1) += vir.yy;
02693         ppreduction->item(ppoffset+2) += vir.zz;
02694       }
02695     }
02696   }
02697   if ( dt && virial ) *virial += wc;
02698 
02699   return 0;
02700 }
02701 
02702 //  RATTLE algorithm from Allen & Tildesley
02703 void HomePatch::rattle2(const BigReal timestep, Tensor *virial)
02704 {
02705   Molecule *mol = Node::Object()->molecule;
02706   SimParameters *simParams = Node::Object()->simParameters;
02707   const int fixedAtomsOn = simParams->fixedAtomsOn;
02708   const int useSettle = simParams->useSettle;
02709   const BigReal dt = timestep / TIMEFACTOR;
02710   Tensor wc;  // constraint virial
02711   BigReal tol = simParams->rigidTol;
02712   int maxiter = simParams->rigidIter;
02713   int dieOnError = simParams->rigidDie;
02714   int i, iter;
02715   BigReal dsqi[10], tmp;
02716   int ial[10], ibl[10];
02717   Vector ref[10];  // reference position
02718   Vector refab[10];  // reference vector
02719   Vector vel[10];  // new velocity
02720   BigReal rmass[10];  // 1 / mass
02721   BigReal redmass[10];  // reduced mass
02722   int fixed[10];  // is atom fixed?
02723   
02724   // Size of a hydrogen group for water
02725   int wathgsize = 3;
02726   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02727   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02728 
02729   //  CkPrintf("In rattle2!\n");
02730   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02731     //    CkPrintf("ig=%d\n",ig);
02732     int hgs = atom[ig].hydrogenGroupSize;
02733     if ( hgs == 1 ) continue;  // only one atom in group
02734     // cache data in local arrays and integrate positions normally
02735     int anyfixed = 0;
02736     for ( i = 0; i < hgs; ++i ) {
02737       ref[i] = atom[ig+i].position;
02738       vel[i] = atom[ig+i].velocity;
02739       rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
02740       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02741       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02742     }
02743     int icnt = 0;
02744     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02745       if (hgs != wathgsize) {
02746         NAMD_bug("Hydrogen group error caught in rattle2().");
02747       }
02748       // Use SETTLE for water unless some of the water atoms are fixed,
02749       if (useSettle && !anyfixed) {
02750         if (simParams->watmodel == WAT_SWM4) {
02751           // SWM4 ordering:  O D LP H1 H2
02752           // do swap(O,LP) and call settle with subarray O H1 H2
02753           // swap back after we return
02754           Vector lp_ref = ref[2];
02755           // Vector lp_vel = vel[2];
02756           ref[2] = ref[0];
02757           vel[2] = vel[0];
02758           settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
02759           ref[0] = ref[2];
02760           vel[0] = vel[2];
02761           ref[2] = lp_ref;
02762           // vel[2] = vel[0];  // set LP vel to O vel
02763         }
02764         else {
02765           settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
02766           if (simParams->watmodel == WAT_TIP4) {
02767             vel[3] = vel[0];
02768           }
02769         }
02770         for (i=0; i<hgs; i++) {
02771           atom[ig+i].velocity = vel[i];
02772         }
02773         continue;
02774       }
02775       if ( !(fixed[1] && fixed[2]) ) {
02776         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02777         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02778       }
02779     }
02780     //    CkPrintf("Loop 2\n");
02781     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02782       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02783         if ( !(fixed[0] && fixed[i]) ) {
02784           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02785           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02786           ibl[icnt] = i;  ++icnt;
02787         }
02788       }
02789     }
02790     if ( icnt == 0 ) continue;  // no constraints
02791     //    CkPrintf("Loop 3\n");
02792     for ( i = 0; i < icnt; ++i ) {
02793       refab[i] = ref[ial[i]] - ref[ibl[i]];
02794     }
02795     //    CkPrintf("Loop 4\n");
02796     int done;
02797     for ( iter = 0; iter < maxiter; ++iter ) {
02798       done = 1;
02799       for ( i = 0; i < icnt; ++i ) {
02800         int a = ial[i];  int b = ibl[i];
02801         Vector vab = vel[a] - vel[b];
02802         Vector &rab = refab[i];
02803         BigReal rabsqi = dsqi[i];
02804         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02805         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02806           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02807           wc += outer(dp,rab);
02808           vel[a] += rmass[a] * dp;
02809           vel[b] -= rmass[b] * dp;
02810           done = 0;
02811         }
02812       }
02813       if ( done ) break;
02814       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02815     }
02816     if ( ! done ) {
02817       if ( dieOnError ) {
02818         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02819       } else {
02820         iout << iWARN <<
02821           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02822       }
02823     }
02824     // store data back to patch
02825     for ( i = 0; i < hgs; ++i ) {
02826       atom[ig+i].velocity = vel[i];
02827     }
02828   }
02829   //  CkPrintf("Leaving rattle2!\n");
02830   // check that there isn't a constant needed here!
02831   *virial += wc / ( 0.5 * dt );
02832 
02833 }
02834 
02835 
02836 //  Adjust gradients for minimizer
02837 void HomePatch::minimize_rattle2(const BigReal timestep, Tensor *virial, bool forces)
02838 {
02839   Molecule *mol = Node::Object()->molecule;
02840   SimParameters *simParams = Node::Object()->simParameters;
02841   Force *f1 = f[Results::normal].begin();
02842   const int fixedAtomsOn = simParams->fixedAtomsOn;
02843   const int useSettle = simParams->useSettle;
02844   const BigReal dt = timestep / TIMEFACTOR;
02845   Tensor wc;  // constraint virial
02846   BigReal tol = simParams->rigidTol;
02847   int maxiter = simParams->rigidIter;
02848   int dieOnError = simParams->rigidDie;
02849   int i, iter;
02850   BigReal dsqi[10], tmp;
02851   int ial[10], ibl[10];
02852   Vector ref[10];  // reference position
02853   Vector refab[10];  // reference vector
02854   Vector vel[10];  // new velocity
02855   BigReal rmass[10];  // 1 / mass
02856   BigReal redmass[10];  // reduced mass
02857   int fixed[10];  // is atom fixed?
02858   
02859   // Size of a hydrogen group for water
02860   int wathgsize = 3;
02861   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02862   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02863 
02864   //  CkPrintf("In rattle2!\n");
02865   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02866     //    CkPrintf("ig=%d\n",ig);
02867     int hgs = atom[ig].hydrogenGroupSize;
02868     if ( hgs == 1 ) continue;  // only one atom in group
02869     // cache data in local arrays and integrate positions normally
02870     int anyfixed = 0;
02871     for ( i = 0; i < hgs; ++i ) {
02872       ref[i] = atom[ig+i].position;
02873       vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
02874       rmass[i] = 1.0;
02875       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02876       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02877     }
02878     int icnt = 0;
02879     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02880       if (hgs != wathgsize) {
02881         NAMD_bug("Hydrogen group error caught in rattle2().");
02882       }
02883       // Use SETTLE for water unless some of the water atoms are fixed,
02884       if (useSettle && !anyfixed) {
02885         if (simParams->watmodel == WAT_SWM4) {
02886           // SWM4 ordering:  O D LP H1 H2
02887           // do swap(O,LP) and call settle with subarray O H1 H2
02888           // swap back after we return
02889           Vector lp_ref = ref[2];
02890           // Vector lp_vel = vel[2];
02891           ref[2] = ref[0];
02892           vel[2] = vel[0];
02893           settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
02894           ref[0] = ref[2];
02895           vel[0] = vel[2];
02896           ref[2] = lp_ref;
02897           // vel[2] = vel[0];  // set LP vel to O vel
02898         }
02899         else {
02900           settle2(1.0, 1.0, ref, vel, dt, virial);
02901           if (simParams->watmodel == WAT_TIP4) {
02902             vel[3] = vel[0];
02903           }
02904         }
02905         for (i=0; i<hgs; i++) {
02906           ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
02907         }
02908         continue;
02909       }
02910       if ( !(fixed[1] && fixed[2]) ) {
02911         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02912         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02913       }
02914     }
02915     //    CkPrintf("Loop 2\n");
02916     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02917       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02918         if ( !(fixed[0] && fixed[i]) ) {
02919           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02920           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02921           ibl[icnt] = i;  ++icnt;
02922         }
02923       }
02924     }
02925     if ( icnt == 0 ) continue;  // no constraints
02926     //    CkPrintf("Loop 3\n");
02927     for ( i = 0; i < icnt; ++i ) {
02928       refab[i] = ref[ial[i]] - ref[ibl[i]];
02929     }
02930     //    CkPrintf("Loop 4\n");
02931     int done;
02932     for ( iter = 0; iter < maxiter; ++iter ) {
02933       done = 1;
02934       for ( i = 0; i < icnt; ++i ) {
02935         int a = ial[i];  int b = ibl[i];
02936         Vector vab = vel[a] - vel[b];
02937         Vector &rab = refab[i];
02938         BigReal rabsqi = dsqi[i];
02939         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02940         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02941           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02942           wc += outer(dp,rab);
02943           vel[a] += rmass[a] * dp;
02944           vel[b] -= rmass[b] * dp;
02945           done = 0;
02946         }
02947       }
02948       if ( done ) break;
02949       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02950     }
02951     if ( ! done ) {
02952       if ( dieOnError ) {
02953         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02954       } else {
02955         iout << iWARN <<
02956           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02957       }
02958     }
02959     // store data back to patch
02960     for ( i = 0; i < hgs; ++i ) {
02961       ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
02962     }
02963   }
02964   //  CkPrintf("Leaving rattle2!\n");
02965   // check that there isn't a constant needed here!
02966   *virial += wc / ( 0.5 * dt );
02967 
02968 }
02969 
02970 
02971 // BEGIN LA
02972 void HomePatch::loweAndersenVelocities()
02973 {
02974     DebugM(2, "loweAndersenVelocities\n");
02975     Molecule *mol = Node::Object()->molecule;
02976     SimParameters *simParams = Node::Object()->simParameters;
02977     v.resize(numAtoms);
02978     for (int i = 0; i < numAtoms; ++i) {
02979         //v[i] = p[i];
02980         // co-opt CompAtom structure to pass velocity and mass information
02981         v[i].position = atom[i].velocity;
02982         v[i].charge = atom[i].mass;
02983     }
02984     DebugM(2, "loweAndersenVelocities\n");
02985 }
02986 
02987 void HomePatch::loweAndersenFinish()
02988 {
02989     DebugM(2, "loweAndersenFinish\n");
02990     v.resize(0);
02991 }
02992 // END LA
02993 
02994 //LCPO
02995 void HomePatch::setLcpoType() {
02996   Molecule *mol = Node::Object()->molecule;
02997   const int *lcpoParamType = mol->getLcpoParamType();
02998 
02999   lcpoType.resize(numAtoms);
03000   for (int i = 0; i < numAtoms; i++) {
03001     lcpoType[i] = lcpoParamType[pExt[i].id];
03002   }
03003 }
03004 
03005 //set intrinsic radii of atom when doMigration
03006 void HomePatch::setGBISIntrinsicRadii() {
03007   intRad.resize(numAtoms*2);
03008   intRad.setall(0);
03009   Molecule *mol = Node::Object()->molecule;
03010   SimParameters *simParams = Node::Object()->simParameters;
03011   Real offset = simParams->coulomb_radius_offset;
03012   for (int i = 0; i < numAtoms; i++) {
03013     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
03014     Real screen = MassToScreen(atom[i].mass);//same
03015     intRad[2*i+0] = rad - offset;//r0
03016     intRad[2*i+1] = screen*(rad - offset);//s0
03017   }
03018 }
03019 
03020 //compute born radius after phase 1, before phase 2
03021 void HomePatch::gbisComputeAfterP1() {
03022 
03023   SimParameters *simParams = Node::Object()->simParameters;
03024   BigReal alphaMax = simParams->alpha_max;
03025   BigReal delta = simParams->gbis_delta;
03026   BigReal beta = simParams->gbis_beta;
03027   BigReal gamma = simParams->gbis_gamma;
03028   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03029 
03030   BigReal rhoi;
03031   BigReal rhoi0;
03032   //calculate bornRad from psiSum
03033   for (int i = 0; i < numAtoms; i++) {
03034     rhoi0 = intRad[2*i];
03035     rhoi = rhoi0+coulomb_radius_offset;
03036     psiFin[i] += psiSum[i];
03037     psiFin[i] *= rhoi0;
03038     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
03039     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
03040 #ifdef PRINT_COMP
03041     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
03042 #endif
03043   }
03044 
03045   gbisP2Ready();
03046 }
03047 
03048 //compute dHdrPrefix after phase 2, before phase 3
03049 void HomePatch::gbisComputeAfterP2() {
03050 
03051   SimParameters *simParams = Node::Object()->simParameters;
03052   BigReal delta = simParams->gbis_delta;
03053   BigReal beta = simParams->gbis_beta;
03054   BigReal gamma = simParams->gbis_gamma;
03055   BigReal epsilon_s = simParams->solvent_dielectric;
03056   BigReal epsilon_p = simParams->dielectric;
03057   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
03058   BigReal epsilon_p_i = 1/simParams->dielectric;
03059   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03060   BigReal kappa = simParams->kappa;
03061   BigReal fij, expkappa, Dij, dEdai, dedasum;
03062   BigReal rhoi, rhoi0, psii, nbetapsi;
03063   BigReal gammapsi2, tanhi, daidr;
03064   for (int i = 0; i < numAtoms; i++) {
03065     //add diagonal dEda term
03066     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
03067     fij = bornRad[i];//inf
03068     expkappa = exp(-kappa*fij);//0
03069     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
03070     //calculate dHij prefix
03071     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
03072                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
03073     dHdrPrefix[i] += dEdai;
03074     dedasum = dHdrPrefix[i];
03075 
03076     rhoi0 = intRad[2*i];
03077     rhoi = rhoi0+coulomb_radius_offset;
03078     psii = psiFin[i];
03079     nbetapsi = -beta*psii;
03080     gammapsi2 = gamma*psii*psii;
03081     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
03082     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
03083            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
03084     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
03085 #ifdef PRINT_COMP
03086     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
03087 #endif
03088   }
03089   gbisP3Ready();
03090 }
03091 
03092 //send born radius to proxies to begin phase 2
03093 void HomePatch::gbisP2Ready() {
03094   if (proxy.size() > 0) {
03095     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03096     for (int i = 0; i < proxy.size(); i++) {
03097       int node = proxy[i];
03098       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
03099       msg->patch = patchID;
03100       msg->origPe = CkMyPe();
03101       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
03102       msg->destPe = node;
03103       int seq = flags.sequence;
03104       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03105       SET_PRIORITY(msg,seq,priority);
03106       cp[node].recvData(msg);
03107     }
03108   }
03109   Patch::gbisP2Ready();
03110 }
03111 
03112 //send dHdrPrefix to proxies to begin phase 3
03113 void HomePatch::gbisP3Ready() {
03114   if (proxy.size() > 0) {
03115     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03116     //only nonzero message should be sent for doFullElec
03117     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
03118     for (int i = 0; i < proxy.size(); i++) {
03119       int node = proxy[i];
03120       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
03121       msg->patch = patchID;
03122       msg->dHdrPrefixLen = msgAtoms;
03123       msg->origPe = CkMyPe();
03124       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
03125       msg->destPe = node;
03126       int seq = flags.sequence;
03127       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03128       SET_PRIORITY(msg,seq,priority);
03129       cp[node].recvData(msg);
03130     }
03131   }
03132   Patch::gbisP3Ready();
03133 }
03134 
03135 //receive proxy results from phase 1
03136 void HomePatch::receiveResult(ProxyGBISP1ResultMsg *msg) {
03137   ++numGBISP1Arrived;
03138     for ( int i = 0; i < msg->psiSumLen; ++i ) {
03139       psiFin[i] += msg->psiSum[i];
03140     }
03141   delete msg;
03142 
03143   if (flags.doNonbonded) {
03144     //awaken if phase 1 done
03145     if (phase1BoxClosedCalled == true &&
03146         numGBISP1Arrived==proxy.size() ) {
03147         sequencer->awaken();
03148     }
03149   } else {
03150     //awaken if all phases done on noWork step
03151     if (boxesOpen == 0 &&
03152         numGBISP1Arrived == proxy.size() &&
03153         numGBISP2Arrived == proxy.size() &&
03154         numGBISP3Arrived == proxy.size()) {
03155       sequencer->awaken();
03156     }
03157   }
03158 }
03159 
03160 //receive proxy results from phase 2
03161 void HomePatch::receiveResult(ProxyGBISP2ResultMsg *msg) {
03162   ++numGBISP2Arrived;
03163   //accumulate dEda
03164   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
03165     dHdrPrefix[i] += msg->dEdaSum[i];
03166   }
03167   delete msg;
03168 
03169   if (flags.doNonbonded) {
03170     //awaken if phase 2 done
03171     if (phase2BoxClosedCalled == true &&
03172         numGBISP2Arrived==proxy.size() ) {
03173       sequencer->awaken();
03174     }
03175   } else {
03176     //awaken if all phases done on noWork step
03177     if (boxesOpen == 0 &&
03178         numGBISP1Arrived == proxy.size() &&
03179         numGBISP2Arrived == proxy.size() &&
03180         numGBISP3Arrived == proxy.size()) {
03181       sequencer->awaken();
03182     }
03183   }
03184 }
03185 
03186 //  MOLLY algorithm part 1
03187 void HomePatch::mollyAverage()
03188 {
03189   Molecule *mol = Node::Object()->molecule;
03190   SimParameters *simParams = Node::Object()->simParameters;
03191   BigReal tol = simParams->mollyTol;
03192   int maxiter = simParams->mollyIter;
03193   int i, iter;
03194   HGArrayBigReal dsq;
03195   BigReal tmp;
03196   HGArrayInt ial, ibl;
03197   HGArrayVector ref;  // reference position
03198   HGArrayVector refab;  // reference vector
03199   HGArrayBigReal rmass;  // 1 / mass
03200   BigReal *lambda;  // Lagrange multipliers
03201   CompAtom *avg;  // averaged position
03202   int numLambdas = 0;
03203   HGArrayInt fixed;  // is atom fixed?
03204 
03205   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
03206   p_avg.resize(numAtoms);
03207   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
03208 
03209   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03210     int hgs = atom[ig].hydrogenGroupSize;
03211     if ( hgs == 1 ) continue;  // only one atom in group
03212         for ( i = 0; i < hgs; ++i ) {
03213           ref[i] = atom[ig+i].position;
03214           rmass[i] = 1. / atom[ig+i].mass;
03215           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03216           if ( fixed[i] ) rmass[i] = 0.;
03217         }
03218         avg = &(p_avg[ig]);
03219         int icnt = 0;
03220 
03221   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
03222           if ( hgs != 3 ) {
03223             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
03224           }
03225           if ( !(fixed[1] && fixed[2]) ) {
03226             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03227           }
03228         }
03229         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03230     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
03231             if ( !(fixed[0] && fixed[i]) ) {
03232               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03233             }
03234           }
03235         }
03236         if ( icnt == 0 ) continue;  // no constraints
03237         numLambdas += icnt;
03238         molly_lambda.resize(numLambdas);
03239         lambda = &(molly_lambda[numLambdas - icnt]);
03240         for ( i = 0; i < icnt; ++i ) {
03241           refab[i] = ref[ial[i]] - ref[ibl[i]];
03242         }
03243         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
03244         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
03245         if ( iter == maxiter ) {
03246           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
03247         }
03248   }
03249 
03250   // for ( i=0; i<numAtoms; ++i ) {
03251   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
03252   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
03253   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
03254   //    }
03255   // }
03256 
03257 }
03258 
03259 
03260 //  MOLLY algorithm part 2
03261 void HomePatch::mollyMollify(Tensor *virial)
03262 {
03263   Molecule *mol = Node::Object()->molecule;
03264   SimParameters *simParams = Node::Object()->simParameters;
03265   Tensor wc;  // constraint virial
03266   int i;
03267   HGArrayInt ial, ibl;
03268   HGArrayVector ref;  // reference position
03269   CompAtom *avg;  // averaged position
03270   HGArrayVector refab;  // reference vector
03271   HGArrayVector force;  // new force
03272   HGArrayBigReal rmass;  // 1 / mass
03273   BigReal *lambda;  // Lagrange multipliers
03274   int numLambdas = 0;
03275   HGArrayInt fixed;  // is atom fixed?
03276 
03277   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03278     int hgs = atom[ig].hydrogenGroupSize;
03279     if (hgs == 1 ) continue;  // only one atom in group
03280         for ( i = 0; i < hgs; ++i ) {
03281           ref[i] = atom[ig+i].position;
03282           force[i] = f[Results::slow][ig+i];
03283           rmass[i] = 1. / atom[ig+i].mass;
03284           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03285           if ( fixed[i] ) rmass[i] = 0.;
03286         }
03287         int icnt = 0;
03288         // c-ji I'm only going to mollify water for now
03289   if ( atom[ig].rigidBondLength > 0 ) {  // for water
03290           if ( hgs != 3 ) {
03291             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
03292           }
03293           if ( !(fixed[1] && fixed[2]) ) {
03294             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03295           }
03296         }
03297         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03298     if ( atom[ig+i].rigidBondLength > 0 ) {
03299             if ( !(fixed[0] && fixed[i]) ) {
03300               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03301             }
03302           }
03303         }
03304 
03305         if ( icnt == 0 ) continue;  // no constraints
03306         lambda = &(molly_lambda[numLambdas]);
03307         numLambdas += icnt;
03308         for ( i = 0; i < icnt; ++i ) {
03309           refab[i] = ref[ial[i]] - ref[ibl[i]];
03310         }
03311         avg = &(p_avg[ig]);
03312         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
03313         // store data back to patch
03314         for ( i = 0; i < hgs; ++i ) {
03315           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
03316           f[Results::slow][ig+i] = force[i];
03317         }
03318   }
03319   // check that there isn't a constant needed here!
03320   *virial += wc;
03321   p_avg.resize(0);
03322 }
03323 
03324 void HomePatch::checkpoint(void) {
03325   checkpoint_atom.copy(atom);
03326   checkpoint_lattice = lattice;
03327 
03328   // DMK - Atom Separation (water vs. non-water)
03329   #if NAMD_SeparateWaters != 0
03330     checkpoint_numWaterAtoms = numWaterAtoms;
03331   #endif
03332 }
03333 
03334 void HomePatch::revert(void) {
03335   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03336 
03337   atom.copy(checkpoint_atom);
03338   numAtoms = atom.size();
03339   lattice = checkpoint_lattice;
03340 
03341   doAtomUpdate = true;
03342   rattleListValid = false;
03343 
03344   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03345 
03346   // DMK - Atom Separation (water vs. non-water)
03347   #if NAMD_SeparateWaters != 0
03348     numWaterAtoms = checkpoint_numWaterAtoms;
03349   #endif
03350 }
03351 
03352 void HomePatch::exchangeCheckpoint(int scriptTask, int &bpc) {  // initiating replica
03353   SimParameters *simParams = Node::Object()->simParameters;
03354   checkpoint_task = scriptTask;
03355   const int remote = simParams->scriptIntArg1;
03356   const char *key = simParams->scriptStringArg1;
03357   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
03358 }
03359 
03360 void HomePatch::recvCheckpointReq(int task, const char *key, int replica, int pe) {  // responding replica
03361   if ( task == SCRIPT_CHECKPOINT_FREE ) {
03362     if ( ! checkpoints.count(key) ) {
03363       NAMD_die("Unable to free checkpoint, requested key was never stored.");
03364     }
03365     delete checkpoints[key];
03366     checkpoints.erase(key);
03367     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
03368     return;
03369   }
03370   CheckpointAtomsMsg *msg;
03371   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
03372     if ( ! checkpoints.count(key) ) {
03373       NAMD_die("Unable to load checkpoint, requested key was never stored.");
03374     }
03375     checkpoint_t &cp = *checkpoints[key];
03376     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
03377     msg->lattice = cp.lattice;
03378     msg->berendsenPressure_count = cp.berendsenPressure_count;
03379     msg->numAtoms = cp.numAtoms;
03380     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
03381   } else {
03382     msg = new (0,1,0) CheckpointAtomsMsg;
03383   }
03384   msg->pid = patchID;
03385   msg->replica = CmiMyPartition();
03386   msg->pe = CkMyPe();
03387   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
03388 }
03389 
03390 void HomePatch::recvCheckpointLoad(CheckpointAtomsMsg *msg) {  // initiating replica
03391   SimParameters *simParams = Node::Object()->simParameters;
03392   const int remote = simParams->scriptIntArg1;
03393   const char *key = simParams->scriptStringArg1;
03394   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
03395     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
03396   }
03397   if ( msg->replica != remote ) {
03398     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
03399   }
03400   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03401     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
03402     strcpy(newmsg->key,key);
03403     newmsg->lattice = lattice;
03404     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
03405     newmsg->pid = patchID;
03406     newmsg->pe = CkMyPe();
03407     newmsg->replica = CmiMyPartition();
03408     newmsg->numAtoms = numAtoms;
03409     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03410     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
03411   }
03412   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03413     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03414     lattice = msg->lattice;
03415     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
03416     numAtoms = msg->numAtoms;
03417     atom.resize(numAtoms);
03418     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03419     doAtomUpdate = true;
03420     rattleListValid = false;
03421     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03422   }
03423   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
03424     recvCheckpointAck();
03425   }
03426   delete msg;
03427 }
03428 
03429 void HomePatch::recvCheckpointStore(CheckpointAtomsMsg *msg) {  // responding replica
03430   if ( ! checkpoints.count(msg->key) ) {
03431     checkpoints[msg->key] = new checkpoint_t;
03432   }
03433   checkpoint_t &cp = *checkpoints[msg->key];
03434   cp.lattice = msg->lattice;
03435   cp.berendsenPressure_count = msg->berendsenPressure_count;
03436   cp.numAtoms = msg->numAtoms;
03437   cp.atoms.resize(cp.numAtoms);
03438   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
03439   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
03440   delete msg;
03441 }
03442 
03443 void HomePatch::recvCheckpointAck() {  // initiating replica
03444   CkpvAccess(_qd)->process();
03445 }
03446 
03447 
03448 void HomePatch::exchangeAtoms(int scriptTask) {
03449   SimParameters *simParams = Node::Object()->simParameters;
03450   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
03451   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
03452     exchange_dst = (int) simParams->scriptArg1;
03453     // create and save outgoing message
03454     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
03455     exchange_msg->lattice = lattice;
03456     exchange_msg->pid = patchID;
03457     exchange_msg->numAtoms = numAtoms;
03458     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03459     if ( exchange_req >= 0 ) {
03460       recvExchangeReq(exchange_req);
03461     }
03462   }
03463   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
03464     exchange_src = (int) simParams->scriptArg2;
03465     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
03466   }
03467 }
03468 
03469 void HomePatch::recvExchangeReq(int req) {
03470   exchange_req = req;
03471   if ( exchange_msg ) {
03472     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
03473     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
03474     exchange_msg = 0;
03475     exchange_req = -1;
03476     CkpvAccess(_qd)->process();
03477   }
03478 }
03479 
03480 void HomePatch::recvExchangeMsg(ExchangeAtomsMsg *msg) {
03481   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
03482   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03483   lattice = msg->lattice;
03484   numAtoms = msg->numAtoms;
03485   atom.resize(numAtoms);
03486   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03487   delete msg;
03488   CkpvAccess(_qd)->process();
03489   doAtomUpdate = true;
03490   rattleListValid = false;
03491   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03492 }
03493 
03494 void HomePatch::submitLoadStats(int timestep)
03495 {
03496   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
03497 }
03498 
03499 
03500 void HomePatch::doPairlistCheck()
03501 {
03502   SimParameters *simParams = Node::Object()->simParameters;
03503 
03504   if ( numAtoms == 0 || ! flags.usePairlists ) {
03505     flags.pairlistTolerance = 0.;
03506     flags.maxAtomMovement = 99999.;
03507     return;
03508   }
03509 
03510   int i; int n = numAtoms;
03511   CompAtom *p_i = p.begin();
03512 
03513   if ( flags.savePairlists ) {
03514     flags.pairlistTolerance = doPairlistCheck_newTolerance;
03515     flags.maxAtomMovement = 0.;
03516     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
03517     doPairlistCheck_lattice = lattice;
03518     doPairlistCheck_positions.resize(numAtoms);
03519     CompAtom *psave_i = doPairlistCheck_positions.begin();
03520     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
03521     return;
03522   }
03523 
03524   Lattice &lattice_old = doPairlistCheck_lattice;
03525   Position center_cur = lattice.unscale(center);
03526   Position center_old = lattice_old.unscale(center);
03527   Vector center_delta = center_cur - center_old;
03528   
03529   // find max deviation to corner (any neighbor shares a corner)
03530   BigReal max_cd = 0.;
03531   for ( i=0; i<2; ++i ) {
03532     for ( int j=0; j<2; ++j ) {
03533       for ( int k=0; k<2; ++k ) {
03534         ScaledPosition corner(  i ? min.x : max.x ,
03535                                 j ? min.y : max.y ,
03536                                 k ? min.z : max.z );
03537         Vector corner_delta =
03538                 lattice.unscale(corner) - lattice_old.unscale(corner);
03539         corner_delta -= center_delta;
03540         BigReal cd = corner_delta.length2();
03541         if ( cd > max_cd ) max_cd = cd;
03542       }
03543     }
03544   }
03545   max_cd = sqrt(max_cd);
03546 
03547   // find max deviation of atoms relative to center
03548   BigReal max_pd = 0.;
03549   CompAtom *p_old_i = doPairlistCheck_positions.begin();
03550   for ( i=0; i<n; ++i ) {
03551     Vector p_delta = p_i[i].position - p_old_i[i].position;
03552     p_delta -= center_delta;
03553     BigReal pd = p_delta.length2();
03554     if ( pd > max_pd ) max_pd = pd;
03555   }
03556   max_pd = sqrt(max_pd);
03557 
03558   BigReal max_tol = max_pd + max_cd;
03559 
03560   flags.maxAtomMovement = max_tol;
03561 
03562   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
03563 
03564   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
03565                                 doPairlistCheck_newTolerance ) ) {
03566     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
03567   }
03568 
03569   if ( max_tol > doPairlistCheck_newTolerance ) {
03570     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
03571   }
03572 
03573 }
03574 
03575 void HomePatch::doGroupSizeCheck()
03576 {
03577   if ( ! flags.doNonbonded ) return;
03578 
03579   SimParameters *simParams = Node::Object()->simParameters;
03580   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
03581   BigReal maxrad2 = 0.;
03582 
03583   FullAtomList::iterator p_i = atom.begin();
03584   FullAtomList::iterator p_e = atom.end();
03585 
03586   while ( p_i != p_e ) {
03587     const int hgs = p_i->hydrogenGroupSize;
03588     if ( ! hgs ) break;  // avoid infinite loop on bug
03589     int ngs = hgs;
03590     if ( ngs > 5 ) ngs = 5;  // limit to at most 5 atoms per group
03591     BigReal x = p_i->position.x;
03592     BigReal y = p_i->position.y;
03593     BigReal z = p_i->position.z;
03594     int i;
03595     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
03596       p_i[i].nonbondedGroupSize = 0;
03597       BigReal dx = p_i[i].position.x - x;
03598       BigReal dy = p_i[i].position.y - y;
03599       BigReal dz = p_i[i].position.z - z;
03600       BigReal r2 = dx * dx + dy * dy + dz * dz;
03601       if ( r2 > hgcut ) break;
03602       else if ( r2 > maxrad2 ) maxrad2 = r2;
03603     }
03604     p_i->nonbondedGroupSize = i;
03605     for ( ; i < hgs; ++i ) {
03606       p_i[i].nonbondedGroupSize = 1;
03607     }
03608     p_i += hgs;
03609   }
03610 
03611   if ( p_i != p_e ) {
03612     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
03613   }
03614 
03615   flags.maxGroupRadius = sqrt(maxrad2);
03616 
03617 }
03618 
03619 void HomePatch::doMarginCheck()
03620 {
03621   SimParameters *simParams = Node::Object()->simParameters;
03622 
03623   BigReal sysdima = lattice.a_r().unit() * lattice.a();
03624   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
03625   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
03626 
03627   BigReal minSize = simParams->patchDimension - simParams->margin;
03628 
03629   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
03630        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
03631        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
03632 
03633     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03634       "Possible solutions are to restart from a recent checkpoint,\n"
03635       "increase margin, or disable useFlexibleCell for liquid simulation.");
03636   }
03637 
03638   BigReal cutoff = simParams->cutoff;
03639 
03640   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
03641   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
03642   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
03643 
03644   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
03645     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03646       "There are probably many margin violations already reported.\n"
03647       "Possible solutions are to restart from a recent checkpoint,\n"
03648       "increase margin, or disable useFlexibleCell for liquid simulation.");
03649   }
03650 
03651   BigReal minx = min.x - margina;
03652   BigReal miny = min.y - marginb;
03653   BigReal minz = min.z - marginc;
03654   BigReal maxx = max.x + margina;
03655   BigReal maxy = max.y + marginb;
03656   BigReal maxz = max.z + marginc;
03657 
03658   int xdev, ydev, zdev;
03659   int problemCount = 0;
03660 
03661   FullAtomList::iterator p_i = atom.begin();
03662   FullAtomList::iterator p_e = atom.end();
03663   for ( ; p_i != p_e; ++p_i ) {
03664 
03665     ScaledPosition s = lattice.scale(p_i->position);
03666 
03667     // check if atom is within bounds
03668     if (s.x < minx) xdev = 0;
03669     else if (maxx <= s.x) xdev = 2; 
03670     else xdev = 1;
03671 
03672     if (s.y < miny) ydev = 0;
03673     else if (maxy <= s.y) ydev = 2; 
03674     else ydev = 1;
03675 
03676     if (s.z < minz) zdev = 0;
03677     else if (maxz <= s.z) zdev = 2; 
03678     else zdev = 1;
03679 
03680     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
03681         ++problemCount;
03682     }
03683 
03684   }
03685 
03686   marginViolations = problemCount;
03687   // if ( problemCount ) {
03688   //     iout << iERROR <<
03689   //       "Found " << problemCount << " margin violations!\n" << endi;
03690   // } 
03691 
03692 }
03693 
03694 
03695 void
03696 HomePatch::doAtomMigration()
03697 {
03698   int i;
03699 
03700   for (i=0; i<numNeighbors; i++) {
03701     realInfo[i].mList.resize(0);
03702   }
03703 
03704   // Purge the AtomMap
03705   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03706 
03707   // Determine atoms that need to migrate
03708 
03709   BigReal minx = min.x;
03710   BigReal miny = min.y;
03711   BigReal minz = min.z;
03712   BigReal maxx = max.x;
03713   BigReal maxy = max.y;
03714   BigReal maxz = max.z;
03715 
03716   int xdev, ydev, zdev;
03717   int delnum = 0;
03718 
03719   FullAtomList::iterator atom_i = atom.begin();
03720   FullAtomList::iterator atom_e = atom.end();
03721 
03722   // DMK - Atom Separation (water vs. non-water)
03723   #if NAMD_SeparateWaters != 0
03724     FullAtomList::iterator atom_first = atom_i;
03725     int numLostWaterAtoms = 0;
03726   #endif
03727 
03728   while ( atom_i != atom_e ) {
03729     if ( atom_i->migrationGroupSize ) {
03730       Position pos = atom_i->position;
03731       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
03732         int mgs = atom_i->migrationGroupSize;
03733         int c = 1;
03734         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
03735                                 j+=(atom_i+j)->hydrogenGroupSize ) {
03736           pos += (atom_i+j)->position;
03737           ++c;
03738         }
03739         pos *= 1./c;
03740         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
03741       }
03742 
03743       ScaledPosition s = lattice.scale(pos);
03744 
03745       // check if atom is within bounds
03746       if (s.x < minx) xdev = 0;
03747       else if (maxx <= s.x) xdev = 2;
03748       else xdev = 1;
03749 
03750       if (s.y < miny) ydev = 0;
03751       else if (maxy <= s.y) ydev = 2;
03752       else ydev = 1;
03753 
03754       if (s.z < minz) zdev = 0;
03755       else if (maxz <= s.z) zdev = 2;
03756       else zdev = 1;
03757 
03758     }
03759 
03760     if (mInfo[xdev][ydev][zdev]) { // process atom for migration
03761                                     // Don't migrate if destination is myself
03762 
03763       // See if we have a migration list already
03764       MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
03765       DebugM(3,"Migrating atom " << atom_i->id << " from patch "
03766                 << patchID << " with position " << atom_i->position << "\n");
03767       mCur.add(*atom_i);
03768 
03769       ++delnum;
03770 
03771 
03772       // DMK - Atom Separation (water vs. non-water)
03773       #if NAMD_SeparateWaters != 0
03774         // Check to see if this atom is part of a water molecule.  If
03775         //   so, numWaterAtoms needs to be adjusted to refect the lost of
03776         //   this atom.
03777         // NOTE: The atom separation code assumes that if the oxygen
03778         //   atom of the hydrogen group making up the water molecule is
03779         //   migrated to another HomePatch, the hydrogens will also
03780         //   move!!!
03781         int atomIndex = atom_i - atom_first;
03782         if (atomIndex < numWaterAtoms)
03783           numLostWaterAtoms++;
03784       #endif
03785 
03786 
03787     } else {
03788 
03789       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
03790 
03791     }
03792 
03793     ++atom_i;
03794   }
03795 
03796   // DMK - Atom Separation (water vs. non-water)
03797   #if NAMD_SeparateWaters != 0
03798     numWaterAtoms -= numLostWaterAtoms;
03799   #endif
03800 
03801 
03802   int delpos = numAtoms - delnum;
03803   DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
03804   atom.del(delpos,delnum);
03805 
03806   numAtoms = atom.size();
03807 
03808   PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
03809 
03810   // signal depositMigration() that we are inMigration mode
03811   inMigration = true;
03812 
03813   // Drain the migration message buffer
03814   for (i=0; i<numMlBuf; i++) {
03815      DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
03816      depositMigration(msgbuf[i]);
03817   }
03818   numMlBuf = 0;
03819      
03820   if (!allMigrationIn) {
03821     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
03822     migrationSuspended = true;
03823     sequencer->suspend();
03824     migrationSuspended = false;
03825   }
03826   allMigrationIn = false;
03827 
03828   inMigration = false;
03829   marginViolations = 0;
03830 }
03831 
03832 void 
03833 HomePatch::depositMigration(MigrateAtomsMsg *msg)
03834 {
03835 
03836   if (!inMigration) { // We have to buffer changes due to migration
03837                       // until our patch is in migration mode
03838     msgbuf[numMlBuf++] = msg;
03839     return;
03840   } 
03841 
03842 
03843   // DMK - Atom Separation (water vs. non-water)
03844   #if NAMD_SeparateWaters != 0
03845 
03846 
03847     // Merge the incoming list of atoms with the current list of
03848     //   atoms.  Note that mergeSeparatedAtomList() will apply any
03849     //   required transformations to the incoming atoms as it is
03850     //   separating them.
03851     mergeAtomList(msg->migrationList);
03852 
03853 
03854   #else
03855 
03856     // Merge the incoming list of atoms with the current list of
03857     // atoms.  Apply transformations to the incoming atoms as they are
03858     // added to this patch's list.
03859     {
03860       MigrationList &migrationList = msg->migrationList;
03861       MigrationList::iterator mi;
03862       Transform mother_transform;
03863       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
03864         DebugM(1,"Migrating atom " << mi->id << " to patch "
03865                   << patchID << " with position " << mi->position << "\n"); 
03866         if ( mi->migrationGroupSize ) {
03867           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
03868             Position pos = mi->position;
03869             int mgs = mi->migrationGroupSize;
03870             int c = 1;
03871             for ( int j=mi->hydrogenGroupSize; j<mgs;
03872                                 j+=(mi+j)->hydrogenGroupSize ) {
03873               pos += (mi+j)->position;
03874               ++c;
03875             }
03876             pos *= 1./c;
03877             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
03878             mother_transform = mi->transform;
03879             pos = lattice.nearest(pos,center,&mother_transform);
03880             mi->position = lattice.reverse_transform(mi->position,mi->transform);
03881             mi->position = lattice.apply_transform(mi->position,mother_transform);
03882             mi->transform = mother_transform;
03883           } else {
03884             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
03885             mother_transform = mi->transform;
03886           }
03887         } else {
03888           mi->position = lattice.reverse_transform(mi->position,mi->transform);
03889           mi->position = lattice.apply_transform(mi->position,mother_transform);
03890           mi->transform = mother_transform;
03891         }
03892         atom.add(*mi);
03893       }
03894     }
03895 
03896 
03897   #endif // if (NAMD_SeparateWaters != 0)
03898 
03899 
03900   numAtoms = atom.size();
03901   delete msg;
03902 
03903   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
03904   if (!--patchMigrationCounter) {
03905     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
03906     allMigrationIn = true;
03907     patchMigrationCounter = numNeighbors;
03908     if (migrationSuspended) {
03909       DebugM(3,"patch "<<patchID<<" is being awakened\n");
03910       migrationSuspended = false;
03911       sequencer->awaken();
03912       return;
03913     }
03914     else {
03915        DebugM(3,"patch "<<patchID<<" was not suspended\n");
03916     }
03917   }
03918 }
03919 
03920 
03921 
03922 // DMK - Atom Separation (water vs. non-water)
03923 #if NAMD_SeparateWaters != 0
03924 
03925 // This function will separate waters from non-waters in the current
03926 //   atom list (regardless of whether or not the atom list is has been
03927 //   sorted yet or not).
03928 void HomePatch::separateAtoms() {
03929 
03930   // Basic Idea:  Iterate through all the atoms in the current list
03931   //   of atoms.  Pack the waters in the current atoms list and move
03932   //   the non-waters to the scratch list.  Once the atoms have all
03933   //   been separated, copy the non-waters to the end of the waters.
03934   // NOTE:  This code does not assume that the atoms list has been
03935   //   separated in any manner.
03936 
03937   // NOTE: Sanity check - Doesn't look like the default constructor actually
03938   //   adds any atoms but does set numAtoms. ???
03939   if (atom.size() < 0) return;  // Nothing to do.
03940 
03941   // Resize the scratch FullAtomList (tempAtom)
03942   tempAtom.resize(numAtoms);  // NOTE: Worst case: all non-water
03943 
03944   // Define size of a water hydrogen group
03945   int wathgsize = 3;
03946   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
03947   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
03948 
03949   // Iterate through all the atoms
03950   int i = 0;
03951   int waterIndex = 0;
03952   int nonWaterIndex = 0;
03953   while (i < numAtoms) {
03954 
03955     FullAtom &atom_i = atom[i];
03956     Mass mass = atom_i.mass;
03957     int hgs = atom_i.hydrogenGroupSize; 
03958     // Check to see if this hydrogen group is a water molecule
03959     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
03960 
03961       // Move this hydrogen group up in the current atom list
03962       if (waterIndex != i) {
03963         atom[waterIndex    ] = atom[i    ];  // Oxygen
03964         atom[waterIndex + 1] = atom[i + 1];  // Hydrogen
03965         atom[waterIndex + 2] = atom[i + 2];  // Hydrogen
03966         if (wathgsize > 3) atom[waterIndex + 3] = atom[i + 3];  // lonepair
03967         if (wathgsize > 4) atom[waterIndex + 4] = atom[i + 4];  // drude
03968           // actual Drude water is arranged:  O D LP H H
03969       }
03970 
03971       waterIndex += wathgsize;
03972       i += wathgsize;
03973 
03974     } else {
03975 
03976       // Move this hydrogen group into non-water (scratch) atom list
03977       for (int j = 0; j < hgs; j++)
03978         tempAtom[nonWaterIndex + j] = atom[i + j];
03979 
03980       nonWaterIndex += hgs;
03981       i += hgs;
03982     }
03983 
03984   } // end iterating through atoms
03985 
03986   // Iterate through the non-water (scratch) atom list, adding the
03987   //   atoms to the end of the atom list.
03988   // NOTE: This could be done with a straight memcpy if the internal
03989   //   data structures of ResizeArray could be accessed directly.
03990   //   Or, perhaps add a member to ResizeArray that can add a consecutive
03991   //   list of elements starting at a particular index (would be cleaner).
03992   for (i = 0; i < nonWaterIndex; i++)
03993     atom[waterIndex + i] = tempAtom[i];
03994 
03995   // Set numWaterAtoms
03996   numWaterAtoms = waterIndex;
03997 }
03998 
03999 
04000 // This function will merge the given list of atoms (not assumed to
04001 //   be separated) with the current list of atoms (already assumed
04002 //   to be separated).
04003 // NOTE: This function applies the transformations to the incoming
04004 //   atoms as it is separating them.
04005 void HomePatch::mergeAtomList(FullAtomList &al) {
04006 
04007   // Sanity check
04008   if (al.size() <= 0) return;  // Nothing to do
04009 
04010   const int orig_atomSize = atom.size();
04011   const int orig_alSize = al.size();
04012 
04013   // Resize the atom list (will eventually hold contents of both lists)
04014   atom.resize(orig_atomSize + orig_alSize); // NOTE: Will have contents of both
04015 
04016 
04017   #if 0  // version where non-waters are moved to scratch first
04018 
04019   
04020   // Basic Idea:  The current list is separated already so copy the
04021   //   non-water atoms out of it into the scratch atom array.  Then
04022   //   separate the incoming/given list (al), adding the waters to the
04023   //   end of the waters in atom list and non-waters to the end of the
04024   //   scratch list.  At this point, all waters are in atom list and all
04025   //   non-waters are in the scratch list so just copy the scratch list
04026   //   to the end of the atom list.
04027   // NOTE: If al is already separated and the number of waters in it
04028   //   is know, could simply move the non-waters in atoms back by that
04029   //   amount and directly copy the waters in al into the created gap
04030   //   and the non-waters in al to the end.  Leave this as an
04031   //   optimization for later since I'm not sure if this would actually
04032   //   do better as the combining code (for combining migration
04033   //   messages) would also have to merge the contents of the atom lists
04034   //   they carry.  Generally speaking, there is probably a faster way
04035   //   to do this, but this will get it working.
04036 
04037   // Copy all the non-waters in the current atom list into the
04038   //   scratch atom list.
04039   const int orig_atom_numNonWaters = orig_atomSize - numWaterAtoms;
04040   tempAtom.resize(orig_atom_numNonWaters + al.size()); // NOTE: Worst case
04041   for (int i = 0; i < orig_atom_numNonWaters; i++)
04042     tempAtom[i] = atom[numWaterAtoms + i];
04043 
04044   // Separate the contents of the given atom list (applying the
04045   // transforms as needed)
04046   int atom_waterIndex = numWaterAtoms;
04047   int atom_nonWaterIndex = orig_atom_numNonWaters;
04048   int i = 0;
04049   while (i < orig_alSize) {
04050 
04051     FullAtom &atom_i = al[i];
04052     int hgs = atom_i.hydrogenGroupSize;
04053     if ( hgs != atom_i.migrationGroupSize ) {
04054       NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
04055     }
04056     Mass mass = atom_i.mass;
04057 
04058     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
04059 
04060       // Apply the transforms
04061 
04062       // Oxygen (@ +0)
04063       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
04064       Transform mother_transform = al[i].transform;
04065 
04066       // Hydrogen (@ +1)
04067       al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
04068       al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
04069       al[i+1].transform = mother_transform;
04070 
04071       // Hydrogen (@ +2)
04072       al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
04073       al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
04074       al[i+2].transform = mother_transform;
04075 
04076       // Add to the end of the waters in the current list of atoms
04077       atom[atom_waterIndex    ] = al[i    ];
04078       atom[atom_waterIndex + 1] = al[i + 1];
04079       atom[atom_waterIndex + 2] = al[i + 2];
04080 
04081       atom_waterIndex += 3;
04082       i += 3;
04083 
04084     } else {
04085 
04086       // Apply the transforms
04087 
04088       // Non-Hydrogen (@ +0)
04089       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
04090       Transform mother_transform = al[i].transform;
04091 
04092       // Hydrogens (@ +1 -> +(hgs-1))
04093       for (int j = 1; j < hgs; j++) {
04094         al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
04095         al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
04096         al[i+j].transform = mother_transform;
04097       }
04098 
04099       // Add to the end of the non-waters (scratch) atom list
04100       for (int j = 0; j < hgs; j++)
04101         tempAtom[atom_nonWaterIndex + j] = al[i + j];
04102 
04103       atom_nonWaterIndex += hgs;
04104       i += hgs;
04105     }
04106 
04107   } // end while iterating through given atom list
04108 
04109   // Copy all the non-waters to the end of the current atom list
04110   for (int i = 0; i < atom_nonWaterIndex; i++)
04111     atom[atom_waterIndex + i] = tempAtom[i];
04112 
04113   // Set numWaterAtoms and numAtoms
04114   numWaterAtoms = atom_waterIndex;
04115   numAtoms = atom.size();
04116 
04117 
04118   #else
04119 
04120 
04121   // Basic Idea:  Count the number of water atoms in the incoming atom
04122   //   list then move the non-waters back in the current atom list to
04123   //   make room for the incoming waters.  Once there is room in the
04124   //   current list, separate the incoming list as the atoms are being
04125   //   added to the current list.
04126   // NOTE:  Since the incoming atom list is likely to be small,
04127   //   iterating over its hydrogen groups twice should not be too bad.
04128   // NOTE:  This code assumes the current list is already separated,
04129   //   the incoming list may not be separated, and the transforms are
04130   //   applied to the incoming atoms as the separation occurs.
04131 
04132   // size of a water hydrogen group
04133   int wathgsize = 3;
04134   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
04135   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
04136 
04137   // Count the number of waters in the given atom list
04138   int al_numWaterAtoms = 0;
04139   int i = 0;
04140   while (i < orig_alSize) {
04141 
04142     FullAtom &atom_i = al[i];
04143     int hgs = atom_i.hydrogenGroupSize;
04144     Mass mass = atom_i.mass;
04145 
04146     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
04147       al_numWaterAtoms += wathgsize;
04148     }
04149 
04150     i += hgs;
04151   }
04152 
04153   // Move all of the non-waters in the current atom list back (to a
04154   //   higher index) by the number of waters in the given list.
04155   if (al_numWaterAtoms > 0) {
04156     for (i = orig_atomSize - 1; i >= numWaterAtoms; i--) {
04157       atom[i + al_numWaterAtoms] = atom[i];
04158     }
04159   }
04160 
04161   // Separte the atoms in the given atom list.  Perform the
04162   //   transformations on them and then add them to the appropriate
04163   //   location in the current atom list.
04164   int atom_waterIndex = numWaterAtoms;
04165   int atom_nonWaterIndex = orig_atomSize + al_numWaterAtoms;
04166   i = 0;
04167   while (i < orig_alSize) {
04168 
04169     FullAtom &atom_i = al[i];
04170     int hgs = atom_i.hydrogenGroupSize;
04171     if ( hgs != atom_i.migrationGroupSize ) {
04172       NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
04173     }
04174     Mass mass = atom_i.mass;
04175 
04176     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
04177 
04178       // Apply the transforms
04179 
04180       // Oxygen (@ +0)
04181       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
04182       Transform mother_transform = al[i].transform;
04183 
04184       // Hydrogen (@ +1)
04185       al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
04186       al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
04187       al[i+1].transform = mother_transform;
04188 
04189       // Hydrogen (@ +2)
04190       al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
04191       al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
04192       al[i+2].transform = mother_transform;
04193 
04194       // Add to the end of the waters in the current list of atoms
04195       atom[atom_waterIndex    ] = al[i    ];
04196       atom[atom_waterIndex + 1] = al[i + 1];
04197       atom[atom_waterIndex + 2] = al[i + 2];
04198 
04199       if (wathgsize > 3) atom[atom_waterIndex + 3] = al[i + 3];
04200 
04201       atom_waterIndex += wathgsize;
04202       i += wathgsize;
04203 
04204     } else {
04205 
04206       // Apply the transforms
04207 
04208       // Non-Hydrogen (@ +0)
04209       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
04210       Transform mother_transform = al[i].transform;
04211 
04212       // Hydrogens (@ +1 -> +(hgs-1))
04213       for (int j = 1; j < hgs; j++) {
04214         al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
04215         al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
04216         al[i+j].transform = mother_transform;
04217       }
04218 
04219       // Add to the end of the non-waters (scratch) atom list
04220       for (int j = 0; j < hgs; j++)
04221         atom[atom_nonWaterIndex + j] = al[i + j];
04222 
04223       atom_nonWaterIndex += hgs;
04224       i += hgs;
04225     }
04226 
04227   } // end while iterating through given atom list
04228 
04229   // Set numWaterAtoms and numAtoms
04230   numWaterAtoms = atom_waterIndex;
04231   numAtoms = atom_nonWaterIndex;
04232 
04233   #endif
04234 }
04235 
04236 #endif
04237 
04238 
04239 
04240 inline void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx,
04241                                               HGArrayBigReal &b)
04242 {
04243         int i,ii=-1,ip,j;
04244         double sum;
04245 
04246         for (i=0;i<n;i++) {
04247                 ip=indx[i];
04248                 sum=b[ip];
04249                 b[ip]=b[i];
04250                 if (ii >= 0)
04251                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
04252                 else if (sum) ii=i;
04253                 b[i]=sum;
04254         }
04255         for (i=n-1;i>=0;i--) {
04256                 sum=b[i];
04257                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
04258                 b[i]=sum/a[i][i];
04259         }
04260 }
04261 
04262 
04263 inline void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
04264 {
04265 
04266         int i,imax,j,k;
04267         double big,dum,sum,temp;
04268         HGArrayBigReal vv;
04269         *d=1.0;
04270         for (i=0;i<n;i++) {
04271                 big=0.0;
04272                 for (j=0;j<n;j++)
04273                         if ((temp=fabs(a[i][j])) > big) big=temp;
04274                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
04275                 vv[i]=1.0/big;
04276         }
04277         for (j=0;j<n;j++) {
04278                 for (i=0;i<j;i++) {
04279                         sum=a[i][j];
04280                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
04281                         a[i][j]=sum;
04282                 }
04283                 big=0.0;
04284                 for (i=j;i<n;i++) {
04285                         sum=a[i][j];
04286                         for (k=0;k<j;k++)
04287                                 sum -= a[i][k]*a[k][j];
04288                         a[i][j]=sum;
04289                         if ( (dum=vv[i]*fabs(sum)) >= big) {
04290                                 big=dum;
04291                                 imax=i;
04292                         }
04293                 }
04294                 if (j != imax) {
04295                         for (k=0;k<n;k++) {
04296                                 dum=a[imax][k];
04297                                 a[imax][k]=a[j][k];
04298                                 a[j][k]=dum;
04299                         }
04300                         *d = -(*d);
04301                         vv[imax]=vv[j];
04302                 }
04303                 indx[j]=imax;
04304                 if (a[j][j] == 0.0) a[j][j]=TINY;
04305                 if (j != n-1) {
04306                         dum=1.0/(a[j][j]);
04307                         for (i=j+1;i<n;i++) a[i][j] *= dum;
04308                 }
04309         }
04310 }
04311 
04312 
04313 inline void G_q(const HGArrayVector &refab, HGMatrixVector &gqij,
04314      const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl) {
04315   int i; 
04316   // step through the rows of the matrix
04317   for(i=0;i<m;i++) {
04318     gqij[i][ial[i]]=2.0*refab[i];
04319     gqij[i][ibl[i]]=-gqij[i][ial[i]];
04320   }
04321 };
04322 
04323 
04324 // c-ji code for MOLLY 7-31-99
04325 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) {
04326   //  input:  n = length of hydrogen group to be averaged (shaked)
04327   //          q[n] = vector array of original positions
04328   //          m = number of constraints
04329   //          imass[n] = inverse mass for each atom
04330   //          length2[m] = square of reference bond length for each constraint
04331   //          ial[m] = atom a in each constraint 
04332   //          ibl[m] = atom b in each constraint 
04333   //          refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
04334   //          tolf = function error tolerance for Newton's iteration
04335   //          ntrial = max number of Newton's iterations
04336   //  output: lambda[m] = double array of lagrange multipliers (used by mollify)
04337   //          qtilde[n] = vector array of averaged (shaked) positions
04338 
04339   int k,k1,i,j;
04340   BigReal errx,errf,d,tolx;
04341 
04342   HGArrayInt indx;
04343   HGArrayBigReal p;
04344   HGArrayBigReal fvec;
04345   HGMatrixBigReal fjac;
04346   HGArrayVector avgab;
04347   HGMatrixVector grhs;
04348   HGMatrixVector auxrhs;
04349   HGMatrixVector glhs;
04350 
04351   //  iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
04352   tolx=tolf; 
04353   
04354   // initialize lambda, globalGrhs
04355 
04356   for (i=0;i<m;i++) {
04357     lambda[i]=0.0;
04358   }
04359 
04360   // define grhs, auxrhs for all iterations
04361   // grhs= g_x(q)
04362   //
04363   G_q(refab,grhs,n,m,ial,ibl);
04364   for (k=1;k<=ntrial;k++) {
04365     //    usrfun(qtilde,q0,lambda,fvec,fjac,n,water); 
04366     HGArrayBigReal gij;
04367     // this used to be main loop of usrfun
04368     // compute qtilde given q0, lambda, IMASSes
04369     {
04370       BigReal multiplier;
04371       HGArrayVector tmp;
04372       for (i=0;i<m;i++) {
04373         multiplier = lambda[i];
04374         // auxrhs = M^{-1}grhs^{T}
04375         for (j=0;j<n;j++) {
04376           auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
04377         }
04378       }
04379       for (j=0;j<n;j++) {
04380         //      tmp[j]=0.0;      
04381         for (i=0;i<m;i++) {
04382           tmp[j]+=auxrhs[i][j];
04383         }
04384       }
04385  
04386       for (j=0;j<n;j++) {
04387         qtilde[j].position=q[j]+tmp[j];
04388       }
04389       //      delete [] tmp;
04390     }
04391   
04392     for ( i = 0; i < m; i++ ) {
04393       avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04394     }
04395 
04396     //  iout<<iINFO << "Calling Jac" << std::endl<<endi;
04397     //  Jac(qtilde, q0, fjac,n,water);
04398     {
04399       //  Vector glhs[3*n+3];
04400 
04401       HGMatrixVector grhs2;
04402 
04403       G_q(avgab,glhs,n,m,ial,ibl);
04404 #ifdef DEBUG0
04405       iout<<iINFO << "G_q:" << std::endl<<endi;
04406       for (i=0;i<m;i++) {
04407         iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
04408       }
04409 #endif
04410       //      G_q(refab,grhs2,m,ial,ibl);
04411       // update with the masses
04412       for (j=0; j<n; j++) { // number of atoms
04413         for (i=0; i<m;i++) { // number of constraints
04414           grhs2[i][j] = grhs[i][j]*imass[j];
04415         }
04416       }
04417 
04418       // G_q(qtilde) * M^-1 G_q'(q0) =
04419       // G_q(qtilde) * grhs'
04420       for (i=0;i<m;i++) { // number of constraints
04421         for (j=0;j<m;j++) { // number of constraints
04422           fjac[i][j] = 0; 
04423           for (k1=0;k1<n;k1++) {
04424             fjac[i][j] += glhs[i][k1]*grhs2[j][k1]; 
04425           }
04426         }
04427       }
04428 #ifdef DEBUG0  
04429       iout<<iINFO << "glhs" <<endi;
04430       for(i=0;i<9;i++) {
04431         iout<<iINFO << glhs[i] << ","<<endi;
04432       }
04433       iout<<iINFO << std::endl<<endi;
04434       for(i=0;i<9;i++) {
04435         iout<<iINFO << grhs2[i] << ","<<endi;
04436       }
04437       iout<<iINFO << std::endl<<endi;
04438 #endif
04439       //      delete[] grhs2;
04440     }
04441     // end of Jac calculation
04442 #ifdef DEBUG0
04443     iout<<iINFO << "Jac" << std::endl<<endi;
04444     for (i=0;i<m;i++) 
04445       for (j=0;j<m;j++)
04446         iout<<iINFO << fjac[i][j] << " "<<endi;
04447     iout<< std::endl<<endi;
04448 #endif
04449     // calculate constraints in gij for n constraints this being a water
04450     //  G(qtilde, gij, n, water);
04451     for (i=0;i<m;i++) {
04452       gij[i]=avgab[i]*avgab[i]-length2[i];
04453     }
04454 #ifdef DEBUG0
04455     iout<<iINFO << "G" << std::endl<<endi;
04456     iout<<iINFO << "( "<<endi;
04457     for(i=0;i<m-1;i++) {
04458       iout<<iINFO << gij[i] << ", "<<endi;
04459     }
04460     iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
04461 #endif
04462     // fill the return vector
04463     for(i=0;i<m;i++) {
04464       fvec[i] = gij[i];
04465     }
04466     // free up the constraints
04467     //    delete[] gij;
04468     // continue Newton's iteration    
04469     errf=0.0;
04470     for (i=0;i<m;i++) errf += fabs(fvec[i]);
04471 #ifdef DEBUG0
04472     iout<<iINFO << "errf: " << errf << std::endl<<endi;
04473 #endif
04474     if (errf <= tolf) {
04475       break;
04476     }
04477     for (i=0;i<m;i++) p[i] = -fvec[i];
04478     //    iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
04479     ludcmp(fjac,m,indx,&d);
04480     lubksb(fjac,m,indx,p);
04481 
04482     errx=0.0;
04483     for (i=0;i<m;i++) {
04484       errx += fabs(p[i]);
04485     }
04486     for (i=0;i<m;i++)  
04487       lambda[i] += p[i];
04488 
04489 #ifdef DEBUG0
04490     iout<<iINFO << "lambda:" << lambda[0] 
04491          << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04492     iout<<iINFO << "errx: " << errx << std::endl<<endi;
04493 #endif
04494     if (errx <= tolx) break;
04495 #ifdef DEBUG0
04496     iout<<iINFO << "Qtilde:" << std::endl<<endi;
04497     iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi; 
04498 #endif
04499   }
04500 #ifdef DEBUG
04501   iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04502 #endif
04503 
04504   return k; // 
04505 }
04506 
04507 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) {
04508   int i,j,k;
04509   BigReal d;
04510   HGMatrixBigReal fjac;
04511   Vector zero(0.0,0.0,0.0);
04512   
04513   HGArrayVector tmpforce;
04514   HGArrayVector tmpforce2;
04515   HGArrayVector y;
04516   HGMatrixVector grhs;
04517   HGMatrixVector glhs;
04518   HGArrayBigReal aux;
04519   HGArrayInt indx;
04520 
04521   for(i=0;i<n;i++) {
04522     tmpforce[i]=imass[i]*force[i];
04523   }
04524 
04525   HGMatrixVector grhs2;
04526   HGArrayVector avgab;
04527 
04528   for ( i = 0; i < m; i++ ) {
04529         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04530   }
04531 
04532   G_q(avgab,glhs,n,m,ial,ibl);
04533   G_q(refab,grhs,n,m,ial,ibl);
04534   // update with the masses
04535   for (j=0; j<n; j++) { // number of atoms
04536         for (i=0; i<m;i++) { // number of constraints
04537           grhs2[i][j] = grhs[i][j]*imass[j];
04538         }
04539   }
04540 
04541   // G_q(qtilde) * M^-1 G_q'(q0) =
04542   // G_q(qtilde) * grhs'
04543   for (i=0;i<m;i++) { // number of constraints
04544         for (j=0;j<m;j++) { // number of constraints
04545           fjac[j][i] = 0; 
04546           for (k=0;k<n;k++) {
04547             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
04548           }
04549         }
04550   }
04551 
04552   // aux=gqij*tmpforce
04553   //  globalGrhs::computeGlobalGrhs(q0,n,water);
04554   //  G_q(refab,grhs,m,ial,ibl);
04555   for(i=0;i<m;i++) {
04556     aux[i]=0.0;
04557     for(j=0;j<n;j++) {
04558       aux[i]+=grhs[i][j]*tmpforce[j];
04559     }
04560   }
04561 
04562   ludcmp(fjac,m,indx,&d);
04563   lubksb(fjac,m,indx,aux);
04564 
04565   for(j=0;j<n;j++) {
04566     y[j] = zero;
04567     for(i=0;i<m;i++) {
04568       y[j] += aux[i]*glhs[i][j];
04569     }
04570   }
04571   for(i=0;i<n;i++) {
04572     y[i]=force[i]-y[i];
04573   }
04574     
04575   // gqq12*y
04576   for(i=0;i<n;i++) {
04577     tmpforce2[i]=imass[i]*y[i];
04578   }
04579 
04580   // here we assume that tmpforce is initialized to zero.
04581   for (i=0;i<n;i++) {
04582     tmpforce[i]=zero;
04583   }
04584   
04585   for (j=0;j<m;j++) {
04586     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
04587     tmpforce[ial[j]] += tmpf;
04588     tmpforce[ibl[j]] -= tmpf;
04589   }
04590   // c-ji the other bug for 2 constraint water was this line (2-4-99)
04591   //  for(i=0;i<m;i++) {
04592   for(i=0;i<n;i++) {
04593     force[i]=tmpforce[i]+y[i];
04594   }
04595 
04596 }

Generated on Sun Oct 21 01:17:13 2018 for NAMD by  doxygen 1.4.7