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 /*
01278  * Redistribute forces from lone pair site onto its host atoms.
01279  *
01280  * Atoms are "labeled" i, j, k, l, where atom i is the lone pair.
01281  * Positions of atoms are ri, rj, rk, rl.
01282  * Forces of atoms are fi, fj, fk, fl; updated forces are returned.
01283  * Accumulate updates to the virial.
01284  *
01285  * The forces on the atoms are updated so that:
01286  *   - the force fi on the lone pair site is 0
01287  *   - the net force fi+fj+fk+fl is conserved
01288  *   - the net torque cross(ri,fi)+cross(rj,fj)+cross(rk,fk)+cross(rl,fl)
01289  *     is conserved
01290  *
01291  * If "midpt" is true (nonzero), then use the midpoint of rk and rl
01292  * (e.g. rk and rl are the hydrogen atoms for water).  Otherwise use
01293  * planes defined by ri,rj,rk and rj,rk,rl.
01294  *
01295  * Having "midpt" set true corresponds in CHARMM to having a negative
01296  * distance parameter in Lphost.
01297  *
01298  * Use FIX_FOR_WATER below to fix the case that occurs when the lone pair
01299  * site for water lies on the perpendicular bisector of rk and rl, making
01300  * the cross(r,s) zero.
01301  */
01302 #define FIX_FOR_WATER
01303 void HomePatch::redistrib_lp_force(
01304     Vector& fi, Vector& fj, Vector& fk, Vector& fl,
01305     const Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
01306     Tensor *virial, int midpt) {
01307 #ifdef DEBUG_REDISTRIB_FORCE
01308   Vector foldnet, toldnet;  // old net force, old net torque
01309   foldnet = fi + fj + fk + fl;
01310   toldnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
01311 #endif
01312   Vector fja(0), fka(0), fla(0);
01313 
01314   Vector r = ri - rj;
01315   BigReal invr2 = r.rlength();
01316   invr2 *= invr2;
01317   BigReal fdot = (fi*r) * invr2;
01318   Vector fr = r * fdot;
01319 
01320   fja += fr;
01321 
01322   Vector s, t;
01323   if (midpt) {
01324     s = rj - 0.5*(rk + rl);
01325     t = 0.5*(rk - rl);
01326   }
01327   else {
01328     s = rj - rk;
01329     t = rk - rl;
01330   }
01331   BigReal invs2 = s.rlength();
01332   invs2 *= invs2;
01333 
01334   Vector p = cross(r,s);
01335 #if !defined(FIX_FOR_WATER)
01336   BigReal invp = p.rlength();
01337 #else
01338   BigReal p2 = p.length2();  // fix division by zero above
01339 #endif
01340 
01341   Vector q = cross(s,t);
01342   BigReal invq = q.rlength();
01343 
01344 #if !defined(FIX_FOR_WATER)
01345   BigReal fpdot = (fi*p) * invp;
01346   Vector fp = p * fpdot;
01347   Vector ft = fi - fr - fp;
01348 #else
01349   BigReal fpdot;
01350   Vector fp, ft;
01351   if (p2 < 1e-6) {  // vector is near zero, assume no fp contribution to force
01352     fpdot = 0;
01353     fp = 0;
01354     ft = fi - fr;
01355   }
01356   else {
01357     fpdot = (fi*p) / sqrt(p2);
01358     fp = p * fpdot;
01359     ft = fi - fr - fp;
01360   }
01361 #endif
01362 
01363   fja += ft;
01364   Vector v = cross(r,ft);  // torque
01365   ft = cross(s,v) * invs2;
01366   fja -= ft;
01367 
01368   if (midpt) {
01369     fka += 0.5 * ft;
01370     fla += 0.5 * ft;
01371   }
01372   else {
01373     fka += ft;
01374   }
01375 
01376   BigReal srdot = (s*r) * invs2;
01377   Vector rr = r - s*srdot;
01378   BigReal rrdot = rr.length();
01379   BigReal stdot = (s*t) * invs2;
01380   Vector tt = t - s*stdot;
01381   BigReal invtt = tt.rlength();
01382   BigReal fact = rrdot*fpdot*invtt*invq;
01383   Vector fq = q * fact;
01384 
01385   fla += fq;
01386   fja += fp*(1+srdot) + fq*stdot;
01387 
01388   ft = fq*(1+stdot) + fp*srdot;
01389 
01390   if (midpt) {
01391     fka += -0.5*ft;
01392     fla += -0.5*ft;
01393   }
01394   else {
01395     fka -= ft;
01396   }
01397 
01398   if (virial) {
01399     Tensor va = outer(fja,rj);
01400     va += outer(fka,rk);
01401     va += outer(fla,rl);
01402     va -= outer(fi,ri);
01403     *virial += va;
01404   }
01405 
01406   fi = 0;  // lone pair has zero force
01407   fj += fja;
01408   fk += fka;
01409   fl += fla;
01410 
01411 #ifdef DEBUG_REDISTRIB_FORCE
01412 #define TOL_REDISTRIB  1e-4
01413   Vector fnewnet, tnewnet;  // new net force, new net torque
01414   fnewnet = fi + fj + fk + fl;
01415   tnewnet = cross(ri,fi) + cross(rj,fj) + cross(rk,fk) + cross(rl,fl);
01416   Vector fdiff = fnewnet - foldnet;
01417   Vector tdiff = tnewnet - toldnet;
01418   if (fdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
01419     printf("Error:  force redistribution for water exceeded tolerance:  "
01420         "fdiff=(%f, %f, %f)\n", fdiff.x, fdiff.y, fdiff.z);
01421   }
01422   if (tdiff.length2() > TOL_REDISTRIB*TOL_REDISTRIB) {
01423     printf("Error:  torque redistribution for water exceeded tolerance:  "
01424         "tdiff=(%f, %f, %f)\n", tdiff.x, tdiff.y, tdiff.z);
01425   }
01426 #endif
01427 }
01428 
01429 
01430 /* Redistribute forces from the massless lonepair charge particle onto
01431  * the other atoms of the water.
01432  *
01433  * This is done using the same algorithm as charmm uses for TIP4P lonepairs.
01434  *
01435  * Pass by reference the forces (O H1 H2 LP) to be modified,
01436  * pass by constant reference the corresponding positions,
01437  * and a pointer to virial.
01438  */
01439 void HomePatch::redistrib_lp_water_force(
01440     Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
01441     const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
01442     const Vector& p_lp, Tensor *virial) {
01443 
01444 #ifdef DEBUG_REDISTRIB_FORCE 
01445   // Debug information to check against results at end
01446 
01447   // total force and torque relative to origin
01448   Vector totforce, tottorque;
01449   totforce = f_ox + f_h1 + f_h2 + f_lp;
01450   tottorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
01451   //printf("Torque without LP is %f/%f/%f\n",
01452   //    tottorque.x, tottorque.y, tottorque.z);
01453   tottorque += cross(f_lp, p_lp);
01454   //printf("Torque with LP is %f/%f/%f\n",
01455   //    tottorque.x, tottorque.y, tottorque.z);
01456 #endif
01457 
01458   // accumulate force adjustments
01459   Vector fad_ox(0), fad_h(0);
01460 
01461   // Calculate the radial component of the force and add it to the oxygen
01462   Vector r_ox_lp = p_lp - p_ox;
01463   BigReal invlen2_r_ox_lp = r_ox_lp.rlength();
01464   invlen2_r_ox_lp *= invlen2_r_ox_lp;
01465   BigReal rad_factor = (f_lp * r_ox_lp) * invlen2_r_ox_lp;
01466   Vector f_rad = r_ox_lp * rad_factor;
01467 
01468   fad_ox += f_rad;
01469 
01470   // Calculate the angular component
01471   Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
01472   Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5; // half of r_h2_h1
01473 
01474   // deviation from collinearity of charge site
01475   //Vector r_oop = cross(r_ox_lp, r_hcom_ox);
01476   //
01477   // vector out of o-h-h plane
01478   //Vector r_perp = cross(r_hcom_ox, r_h2_h1_2);
01479 
01480   // Here we assume that Ox/Lp/Hcom are linear
01481   // If you want to correct for deviations, this is the place
01482 
01483 //  printf("Deviation from linearity for ox %i: %f/%f/%f\n", oxind, r_oop.x, r_oop.y, r_oop.z);
01484 
01485   Vector f_ang = f_lp - f_rad; // leave the angular component
01486 
01487   // now split this component onto the other atoms
01488   BigReal len_r_ox_lp = r_ox_lp.length();
01489   BigReal invlen_r_hcom_ox = r_hcom_ox.rlength();
01490   BigReal oxcomp = (r_hcom_ox.length() - len_r_ox_lp) * invlen_r_hcom_ox;
01491   BigReal hydcomp = 0.5 * len_r_ox_lp * invlen_r_hcom_ox;
01492 
01493   fad_ox += (f_ang * oxcomp);
01494   fad_h += (f_ang * hydcomp);  // adjustment for both hydrogens
01495 
01496   // Add virial contributions
01497   if (virial) {
01498     Tensor vir = outer(fad_ox, p_ox);
01499     vir += outer(fad_h, p_h1);
01500     vir += outer(fad_h, p_h2);
01501     vir -= outer(f_lp, p_lp);
01502     *virial += vir;
01503   }
01504 
01505   //Vector zerovec(0.0, 0.0, 0.0);
01506   f_lp = 0;
01507   f_ox += fad_ox;
01508   f_h1 += fad_h;
01509   f_h2 += fad_h;
01510 
01511 #ifdef DEBUG_REDISTRIB_FORCE 
01512   // Check that the total force and torque come out right
01513   Vector newforce, newtorque;
01514   newforce = f_ox + f_h1 + f_h2;
01515   newtorque = cross(f_ox, p_ox) + cross(f_h1, p_h1) + cross(f_h2, p_h2);
01516   Vector fdiff = newforce - totforce;
01517   Vector tdiff = newtorque - tottorque;
01518   BigReal error = fdiff.length();
01519   if (error > 0.0001) {
01520      printf("Error:  Force redistribution for water "
01521          "exceeded force tolerance:  error=%f\n", error);
01522   }
01523 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
01524   printf("Error in net force:  %f\n", error);
01525 #endif
01526 
01527   error = tdiff.length();
01528   if (error > 0.0001) {
01529      printf("Error:  Force redistribution for water "
01530          "exceeded torque tolerance:  error=%f\n", error);
01531   }
01532 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
01533   printf("Error in net torque:  %f\n", error);
01534 #endif
01535 #endif /* DEBUG */
01536 }
01537 
01538 void HomePatch::reposition_lonepair(
01539     Vector& ri, const Vector& rj, const Vector& rk, const Vector& rl,
01540     Real distance, Real angle, Real dihedral)
01541 {
01542   if ( (rj-rk).length2() > 100. || (rj-rl).length2() > 100. ) {
01543     iout << iWARN << "Large distance between lonepair reference atoms: ("
01544       << rj << ") (" << rk << ") (" << rl << ")\n" << endi;
01545   }
01546   BigReal r, t, p, cst, snt, csp, snp, invlen;
01547   Vector v, w, a, b, c;
01548 
01549   if (distance >= 0) {
01550     v = rk;
01551     r = distance;
01552   }
01553   else {
01554     v = 0.5*(rk + rl);
01555     r = -distance;
01556   }
01557 
01558   t = angle;
01559   p = dihedral;
01560   cst = cos(t);
01561   snt = sin(t);
01562   csp = cos(p);
01563   snp = sin(p);
01564   a = v - rj;
01565   b = rl - v;
01566   invlen = a.rlength();
01567   a *= invlen;
01568   c = cross(b, a);
01569   invlen = c.rlength();
01570   c *= invlen;
01571   b = cross(a, c);
01572   w.x = r*cst;
01573   w.y = r*snt*csp;
01574   w.z = r*snt*snp;
01575   ri.x = rj.x + w.x*a.x + w.y*b.x + w.z*c.x;
01576   ri.y = rj.y + w.x*a.y + w.y*b.y + w.z*c.y;
01577   ri.z = rj.z + w.x*a.z + w.y*b.z + w.z*c.z;
01578 }
01579 
01580 void HomePatch::reposition_all_lonepairs(void) {
01581   // ASSERT: simParams->lonepairs == TRUE
01582   for (int i=0;  i < numAtoms;  i++) {
01583     if (atom[i].mass < 0.01) {
01584       // found a lone pair
01585       AtomID aid = atom[i].id;  // global atom ID of lp
01586       Lphost *lph = Node::Object()->molecule->get_lphost(aid);  // its lphost
01587       if (lph == NULL) {
01588         char errmsg[512];
01589         sprintf(errmsg, "reposition lone pairs: "
01590             "no Lphost exists for LP %d\n", aid);
01591         NAMD_die(errmsg);
01592       }
01593       LocalID j = AtomMap::Object()->localID(lph->atom2);
01594       LocalID k = AtomMap::Object()->localID(lph->atom3);
01595       LocalID l = AtomMap::Object()->localID(lph->atom4);
01596       if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
01597         char errmsg[512];
01598         sprintf(errmsg, "reposition lone pairs: "
01599             "LP %d has some Lphost atom off patch\n", aid);
01600         NAMD_die(errmsg);
01601       }
01602       // reposition this lone pair
01603       reposition_lonepair(atom[i].position, atom[j.index].position,
01604           atom[k.index].position, atom[l.index].position,
01605           lph->distance, lph->angle, lph->dihedral);
01606     }
01607   }
01608 }
01609 
01610 void HomePatch::swm4_omrepos(Vector *ref, Vector *pos, Vector *vel,
01611     BigReal invdt) {
01612   // Reposition lonepair (Om) particle of Drude SWM4 water.
01613   // Same comments apply as to tip4_omrepos(), but the ordering of atoms
01614   // is different: O, D, LP, H1, H2.
01615   pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
01616   // Now, adjust velocity of particle to get it to appropriate place
01617   // during next integration "drift-step"
01618   if (invdt != 0) {
01619     vel[2] = (pos[2] - ref[2]) * invdt;
01620   }
01621   // No virial correction needed since lonepair is massless
01622 }
01623 
01624 void HomePatch::tip4_omrepos(Vector* ref, Vector* pos, Vector* vel, BigReal invdt) {
01625   /* Reposition the om particle of a tip4p water
01626    * A little geometry shows that the appropriate position is given by
01627    * R_O + (1 / 2 r_ohc) * ( 0.5 (R_H1 + R_H2) - R_O ) 
01628    * Here r_om is the distance from the oxygen to Om site, and r_ohc
01629    * is the altitude from the oxygen to the hydrogen center of mass
01630    * Those quantities are precalculated upon initialization of HomePatch
01631    *
01632    * Ordering of TIP4P atoms: O, H1, H2, LP.
01633    */
01634 
01635   //printf("rom/rohc are %f %f and invdt is %f\n", r_om, r_ohc, invdt);
01636   //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);
01637   pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc); 
01638   //printf("New position for lp is %f %f %f\n", pos[3].x, pos[3].y, pos[3].z);
01639 
01640   // Now, adjust the velocity of the particle to get it to the appropriate place
01641   if (invdt != 0) {
01642     vel[3] = (pos[3] - ref[3]) * invdt;
01643   }
01644 
01645   // No virial correction needed, since this is a massless particle
01646   return;
01647 }
01648 
01649 void HomePatch::redistrib_lonepair_forces(const int ftag, Tensor *virial) {
01650   // ASSERT: simParams->lonepairs == TRUE
01651   ForceList *f_mod = f;
01652   for (int i = 0;  i < numAtoms;  i++) {
01653     if (atom[i].mass < 0.01) {
01654       // found a lone pair
01655       AtomID aid = atom[i].id;  // global atom ID of lp
01656       Lphost *lph = Node::Object()->molecule->get_lphost(aid);  // its lphost
01657       if (lph == NULL) {
01658         char errmsg[512];
01659         sprintf(errmsg, "redistrib lone pair forces: "
01660             "no Lphost exists for LP %d\n", aid);
01661         NAMD_die(errmsg);
01662       }
01663       LocalID j = AtomMap::Object()->localID(lph->atom2);
01664       LocalID k = AtomMap::Object()->localID(lph->atom3);
01665       LocalID l = AtomMap::Object()->localID(lph->atom4);
01666       if (j.pid != patchID || k.pid != patchID || l.pid != patchID) {
01667         char errmsg[512];
01668         sprintf(errmsg, "redistrib lone pair forces: "
01669             "LP %d has some Lphost atom off patch\n", aid);
01670         NAMD_die(errmsg);
01671       }
01672       // redistribute forces from this lone pair
01673       int midpt = (lph->distance < 0);
01674       redistrib_lp_force(f_mod[ftag][i], f_mod[ftag][j.index],
01675           f_mod[ftag][k.index], f_mod[ftag][l.index],
01676           atom[i].position, atom[j.index].position,
01677           atom[k.index].position, atom[l.index].position, virial, midpt);
01678     }
01679   }
01680 }
01681 
01682 void HomePatch::redistrib_swm4_forces(const int ftag, Tensor *virial) {
01683   // Loop over the patch's atoms and apply the appropriate corrections
01684   // to get all forces off of lone pairs
01685   ForceList *f_mod = f;
01686   for (int i = 0;  i < numAtoms;  i++) {
01687     if (atom[i].mass < 0.01) {
01688       // found lonepair
01689       redistrib_lp_water_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
01690           f_mod[ftag][i+2], f_mod[ftag][i],
01691           atom[i-2].position, atom[i+1].position,
01692           atom[i+2].position, atom[i].position, virial);
01693     }
01694   }
01695 }
01696 
01697 void HomePatch::redistrib_tip4p_forces(const int ftag, Tensor* virial) {
01698   // Loop over the patch's atoms and apply the appropriate corrections
01699   // to get all forces off of lone pairs
01700   // Atom ordering:  O H1 H2 LP
01701 
01702   ForceList *f_mod =f;
01703   for (int i=0; i<numAtoms; i++) {
01704     if (atom[i].mass < 0.01) {
01705       // found lonepair
01706       redistrib_lp_water_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
01707           f_mod[ftag][i-1], f_mod[ftag][i],
01708           atom[i-3].position, atom[i-2].position,
01709           atom[i-1].position, atom[i].position, virial);
01710     }
01711   }
01712 }
01713 
01714 
01715 void HomePatch::addForceToMomentum(
01716     FullAtom       * __restrict atom_arr,
01717     const Force    * __restrict force_arr,
01718     const BigReal    dt,
01719     int              num_atoms
01720     ) {
01721   SimParameters *simParams = Node::Object()->simParameters;
01722 
01723   if ( simParams->fixedAtomsOn ) {
01724     for ( int i = 0; i < num_atoms; ++i ) {
01725       if ( atom_arr[i].atomFixed ) {
01726         atom_arr[i].velocity = 0;
01727       } else {
01728         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01729         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01730         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01731         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01732       }
01733     }
01734   } else {
01735     for ( int i = 0; i < num_atoms; ++i ) {
01736       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01737       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01738       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01739       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01740     }
01741   }
01742 }
01743 
01744 void HomePatch::addForceToMomentum3(
01745     FullAtom       * __restrict atom_arr,
01746     const Force    * __restrict force_arr1,
01747     const Force    * __restrict force_arr2,
01748     const Force    * __restrict force_arr3,
01749     const BigReal    dt1,
01750     const BigReal    dt2,
01751     const BigReal    dt3,
01752     int              num_atoms
01753     ) {
01754   SimParameters *simParams = Node::Object()->simParameters;
01755 
01756   if ( simParams->fixedAtomsOn ) {
01757     for ( int i = 0; i < num_atoms; ++i ) {
01758       if ( atom_arr[i].atomFixed ) {
01759         atom_arr[i].velocity = 0;
01760       } else {
01761         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01762         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01763             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01764         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01765             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01766         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01767             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01768       }
01769     }
01770   } else {
01771     for ( int i = 0; i < num_atoms; ++i ) {
01772       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01773       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01774           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01775       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01776           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01777       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01778           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01779     }
01780   }
01781 }
01782 
01783 void HomePatch::addVelocityToPosition(
01784     FullAtom       * __restrict atom_arr,
01785     const BigReal    dt,
01786     int              num_atoms
01787     ) {
01788   SimParameters *simParams = Node::Object()->simParameters;
01789   if ( simParams->fixedAtomsOn ) {
01790     for ( int i = 0; i < num_atoms; ++i ) {
01791       if ( ! atom_arr[i].atomFixed ) {
01792         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01793         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01794         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01795       }
01796     }
01797   } else {
01798     for ( int i = 0; i < num_atoms; ++i ) {
01799       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01800       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01801       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01802     }
01803   }
01804 }
01805 
01806 int HomePatch::hardWallDrude(const BigReal timestep, Tensor *virial,
01807     SubmitReduction *ppreduction)
01808 {
01809   Molecule *mol = Node::Object()->molecule;
01810   SimParameters *simParams = Node::Object()->simParameters;
01811   const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
01812   const int fixedAtomsOn = simParams->fixedAtomsOn;
01813   const BigReal dt = timestep / TIMEFACTOR;
01814   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
01815   int i, ia, ib, j;
01816   int dieOnError = simParams->rigidDie;
01817   Tensor wc;  // constraint virial
01818   BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
01819   int nslabs;
01820 
01821   // start data for hard wall boundary between drude and its host atom
01822   // static int Count=0;
01823   int Idx;
01824   double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
01825   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;
01826   double dot_v_r_1, dot_v_r_2;
01827   double vb_cm, dr_a, dr_b;
01828   // end data for hard wall boundary between drude and its host atom
01829 
01830   // start calculation of hard wall boundary between drude and its host atom
01831   if (simParams->drudeHardWallOn) {
01832     if (ppreduction) {
01833       nslabs = simParams->pressureProfileSlabs;
01834       idz = nslabs/lattice.c().z;
01835       zmin = lattice.origin().z - 0.5*lattice.c().z;
01836     }
01837 
01838     r_wall = simParams->drudeBondLen;
01839     r_wall_SQ = r_wall*r_wall;
01840     // Count++;
01841     for (i=1; i<numAtoms; i++)  {
01842       if ( (atom[i].mass > 0.01) && ((atom[i].mass < 1.0)) ) { // drude particle
01843         ia = i-1;
01844         ib = i;
01845 
01846         v_ab = atom[ib].position - atom[ia].position;
01847         rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
01848 
01849         if (rab_SQ > r_wall_SQ) {  // to impose the hard wall constraint
01850           rab = sqrt(rab_SQ);
01851           if ( (rab > (2.0*r_wall)) && dieOnError ) {  // unexpected situation
01852             iout << iERROR << "HardWallDrude> "
01853               << "The drude is too far away from atom "
01854               << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
01855             return -1;  // triggers early exit
01856           }
01857 
01858           v_ab.x /= rab;
01859           v_ab.y /= rab;
01860           v_ab.z /= rab;
01861 
01862           if ( fixedAtomsOn && atom[ia].atomFixed ) {  // the heavy atom is fixed
01863             if (atom[ib].atomFixed) {  // the drude is fixed too
01864               continue;
01865             }
01866             else {  // only the heavy atom is fixed
01867               dot_v_r_2 = atom[ib].velocity.x*v_ab.x
01868                 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
01869               vb_2 = dot_v_r_2 * v_ab;
01870               vp_2 = atom[ib].velocity - vb_2;
01871 
01872               dr = rab - r_wall;
01873               if(dot_v_r_2 == 0.0) {
01874                 delta_T = maxtime;
01875               }
01876               else {
01877                 delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
01878                 if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
01879               }
01880 
01881               dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
01882 
01883               vb_2 = dot_v_r_2 * v_ab;
01884 
01885               new_vel_a = atom[ia].velocity;
01886               new_vel_b = vp_2 + vb_2;
01887 
01888               dr_b = -dr + delta_T*dot_v_r_2;  // L = L_0 + dT *v_new, v was flipped
01889 
01890               new_pos_a = atom[ia].position;
01891               new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
01892             }
01893           }
01894           else {
01895             mass_a = atom[ia].mass;
01896             mass_b = atom[ib].mass;
01897             mass_sum = mass_a+mass_b;
01898 
01899             dot_v_r_1 = atom[ia].velocity.x*v_ab.x
01900               + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
01901             vb_1 = dot_v_r_1 * v_ab;
01902             vp_1 = atom[ia].velocity - vb_1;
01903 
01904             dot_v_r_2 = atom[ib].velocity.x*v_ab.x
01905               + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
01906             vb_2 = dot_v_r_2 * v_ab;
01907             vp_2 = atom[ib].velocity - vb_2;
01908 
01909             vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
01910 
01911             dot_v_r_1 -= vb_cm;
01912             dot_v_r_2 -= vb_cm;
01913 
01914             dr = rab - r_wall;
01915 
01916             if(dot_v_r_2 == dot_v_r_1) {
01917                 delta_T = maxtime;
01918             }
01919             else {
01920               delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);  // the time since the collision occurs
01921               if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
01922             }
01923             
01924             // the relative velocity between ia and ib. Drawn according to T_Drude
01925             v_Bond = sqrt(kbt/mass_b);
01926 
01927             // reflect the velocity along bond vector and scale down
01928             dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
01929             dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
01930 
01931             dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
01932             dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
01933 
01934             new_pos_a = atom[ia].position + dr_a*v_ab;  // correct the position
01935             new_pos_b = atom[ib].position + dr_b*v_ab;
01936             // atom[ia].position += (dr_a*v_ab);  // correct the position
01937             // atom[ib].position += (dr_b*v_ab);
01938             
01939             dot_v_r_1 += vb_cm;
01940             dot_v_r_2 += vb_cm;
01941 
01942             vb_1 = dot_v_r_1 * v_ab;
01943             vb_2 = dot_v_r_2 * v_ab;
01944 
01945             new_vel_a = vp_1 + vb_1;
01946             new_vel_b = vp_2 + vb_2;
01947           }
01948 
01949           int ppoffset, partition;
01950           if ( invdt == 0 ) {
01951             atom[ia].position = new_pos_a;
01952             atom[ib].position = new_pos_b;
01953           }
01954           else if ( virial == 0 ) {
01955             atom[ia].velocity = new_vel_a;
01956             atom[ib].velocity = new_vel_b;
01957           }
01958           else {
01959             for ( j = 0; j < 2; j++ ) {
01960               if (j==0) {  // atom ia, heavy atom
01961                 Idx = ia;
01962                 new_pos = &new_pos_a;
01963                 new_vel = &new_vel_a;
01964               }
01965               else if (j==1) {  // atom ib, drude
01966                 Idx = ib;
01967                 new_pos = &new_pos_b;
01968                 new_vel = &new_vel_b;
01969               }
01970               Force df = (*new_vel - atom[Idx].velocity) *
01971                 ( atom[Idx].mass * invdt );
01972               Tensor vir = outer(df, atom[Idx].position);
01973               wc += vir;
01974               atom[Idx].velocity = *new_vel;
01975               atom[Idx].position = *new_pos;
01976 
01977               if (ppreduction) {
01978                 if (!j) {
01979                   BigReal z = new_pos->z;
01980                   int partition = atom[Idx].partition;
01981                   int slab = (int)floor((z-zmin)*idz);
01982                   if (slab < 0) slab += nslabs;
01983                   else if (slab >= nslabs) slab -= nslabs;
01984                   ppoffset = 3*(slab + nslabs*partition);
01985                 }
01986                 ppreduction->item(ppoffset  ) += vir.xx;
01987                 ppreduction->item(ppoffset+1) += vir.yy;
01988                 ppreduction->item(ppoffset+2) += vir.zz;
01989               }
01990 
01991             }
01992           }
01993         }                               
01994       }
01995     }
01996 
01997     // if ( (Count>10000) && (Count%10==0) ) {
01998     //   v_ab = atom[1].position - atom[0].position;
01999     //   rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
02000     //   iout << "DBG_R: " << Count << "  " << sqrt(rab_SQ) << "\n" << endi;
02001     // }
02002 
02003   }
02004 
02005   // end calculation of hard wall boundary between drude and its host atom
02006 
02007   if ( dt && virial ) *virial += wc;
02008 
02009   return 0;
02010 }
02011 
02012 void HomePatch::buildRattleList() {
02013 
02014   SimParameters *simParams = Node::Object()->simParameters;
02015   const int fixedAtomsOn = simParams->fixedAtomsOn;
02016   const int useSettle = simParams->useSettle;
02017 
02018   // Re-size to containt numAtoms elements
02019   velNew.resize(numAtoms);
02020   posNew.resize(numAtoms);
02021 
02022   // Size of a hydrogen group for water
02023   int wathgsize = 3;
02024   int watmodel = simParams->watmodel;
02025   if (watmodel == WAT_TIP4) wathgsize = 4;
02026   else if (watmodel == WAT_SWM4) wathgsize = 5;
02027 
02028   // Initialize the settle algorithm with water parameters
02029   // settle1() assumes all waters are identical,
02030   // and will generate bad results if they are not.
02031   // XXX this will move to Molecule::build_atom_status when that 
02032   // version is debugged
02033   if ( ! settle_initialized ) {
02034     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02035       // find a water
02036       if (atom[ig].rigidBondLength > 0) {
02037         int oatm;
02038         if (simParams->watmodel == WAT_SWM4) {
02039           oatm = ig+3;  // skip over Drude and Lonepair
02040           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02041           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02042         }
02043         else {
02044           oatm = ig+1;
02045           // Avoid using the Om site to set this by mistake
02046           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02047             oatm += 1;
02048           }
02049         }
02050 
02051         // initialize settle water parameters
02052         settle1init(atom[ig].mass, atom[oatm].mass, 
02053                     atom[ig].rigidBondLength,
02054                     atom[oatm].rigidBondLength,
02055                     settle_mOrmT, settle_mHrmT, settle_ra,
02056                     settle_rb, settle_rc, settle_rra);
02057         settle_initialized = 1;
02058         break; // done with init
02059       }
02060     }
02061   }
02062   
02063   Vector ref[10];
02064   BigReal rmass[10];
02065   BigReal dsq[10];
02066   int fixed[10];
02067   int ial[10];
02068   int ibl[10];
02069 
02070   int numSettle = 0;
02071   int numRattle = 0;
02072   int posRattleParam = 0;
02073 
02074   settleList.clear();
02075   rattleList.clear();
02076   noconstList.clear();
02077   rattleParam.clear();
02078 
02079   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02080     int hgs = atom[ig].hydrogenGroupSize;
02081     if ( hgs == 1 ) {
02082       // only one atom in group
02083       noconstList.push_back(ig);
02084       continue;
02085     }
02086     int anyfixed = 0;
02087     for (int i = 0; i < hgs; ++i ) {
02088       ref[i] = atom[ig+i].position;
02089       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02090       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02091       if ( fixed[i] ) {
02092         anyfixed = 1;
02093         rmass[i] = 0.;
02094       }
02095     }
02096     int icnt = 0;
02097     BigReal tmp = atom[ig].rigidBondLength;
02098     if (tmp > 0.0) {  // for water
02099       if (hgs != wathgsize) {
02100         char errmsg[256];
02101         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02102                          "but the specified water model requires %d atoms.\n",
02103                          atom[ig].id+1, hgs, wathgsize);
02104         NAMD_die(errmsg);
02105       }
02106       // Use SETTLE for water unless some of the water atoms are fixed,
02107       if (useSettle && !anyfixed) {
02108         // Store to Settle -list
02109         settleList.push_back(ig);
02110         continue;
02111       }
02112       if ( !(fixed[1] && fixed[2]) ) {
02113         dsq[icnt] = tmp * tmp;
02114         ial[icnt] = 1;
02115         ibl[icnt] = 2;
02116         ++icnt;
02117       }
02118     } // if (tmp > 0.0)
02119     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02120       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02121         if ( !(fixed[0] && fixed[i]) ) {
02122           dsq[icnt] = tmp * tmp;
02123           ial[icnt] = 0;
02124           ibl[icnt] = i;
02125           ++icnt;
02126         }
02127       }
02128     }
02129     if ( icnt == 0 ) {
02130       // no constraints
02131       noconstList.push_back(ig);
02132       continue;  
02133     }
02134     // Store to Rattle -list
02135     RattleList rattleListElem;
02136     rattleListElem.ig  = ig;
02137     rattleListElem.icnt = icnt;
02138     rattleList.push_back(rattleListElem);
02139     for (int i = 0; i < icnt; ++i ) {
02140       int a = ial[i];
02141       int b = ibl[i];
02142       RattleParam rattleParamElem;
02143       rattleParamElem.ia = a;
02144       rattleParamElem.ib = b;
02145       rattleParamElem.dsq = dsq[i];
02146       rattleParamElem.rma = rmass[a];
02147       rattleParamElem.rmb = rmass[b];
02148       rattleParam.push_back(rattleParamElem);
02149     }
02150   }
02151 
02152 }
02153 
02154 void HomePatch::addRattleForce(const BigReal invdt, Tensor& wc) {
02155   for (int ig = 0; ig < numAtoms; ++ig ) {
02156     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
02157     Tensor vir = outer(df, atom[ig].position);
02158     wc += vir;
02159     f[Results::normal][ig] += df;
02160     atom[ig].velocity = velNew[ig];
02161   }
02162 }
02163 
02164 int HomePatch::rattle1(const BigReal timestep, Tensor *virial, 
02165   SubmitReduction *ppreduction) {
02166 
02167   SimParameters *simParams = Node::Object()->simParameters;
02168   if (simParams->watmodel != WAT_TIP3 || ppreduction) {
02169     // Call old rattle1 -method instead
02170     return rattle1old(timestep, virial, ppreduction);
02171   }
02172 
02173   if (!rattleListValid) {
02174     buildRattleList();
02175     rattleListValid = true;
02176   }
02177 
02178   const int fixedAtomsOn = simParams->fixedAtomsOn;
02179   const int useSettle = simParams->useSettle;
02180   const BigReal dt = timestep / TIMEFACTOR;
02181   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02182   const BigReal tol2 = 2.0 * simParams->rigidTol;
02183   int maxiter = simParams->rigidIter;
02184   int dieOnError = simParams->rigidDie;
02185 
02186   Vector ref[10];  // reference position
02187   Vector pos[10];  // new position
02188   Vector vel[10];  // new velocity
02189 
02190   // Manual un-roll
02191   int n = (settleList.size()/2)*2;
02192   for (int j=0;j < n;j+=2) {
02193     int ig;
02194     ig = settleList[j];
02195     for (int i = 0; i < 3; ++i ) {
02196       ref[i] = atom[ig+i].position;
02197       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02198     }
02199     ig = settleList[j+1];
02200     for (int i = 0; i < 3; ++i ) {
02201       ref[i+3] = atom[ig+i].position;
02202       pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
02203     }
02204     settle1_SIMD<2>(ref, pos,
02205       settle_mOrmT, settle_mHrmT, settle_ra,
02206       settle_rb, settle_rc, settle_rra);
02207 
02208     ig = settleList[j];
02209     for (int i = 0; i < 3; ++i ) {
02210       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02211       posNew[ig+i] = pos[i];
02212     }
02213     ig = settleList[j+1];
02214     for (int i = 0; i < 3; ++i ) {
02215       velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
02216       posNew[ig+i] = pos[i+3];
02217     }
02218 
02219   }
02220 
02221   if (settleList.size() % 2) {
02222     int ig = settleList[settleList.size()-1];
02223     for (int i = 0; i < 3; ++i ) {
02224       ref[i] = atom[ig+i].position;
02225       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02226     }
02227     settle1_SIMD<1>(ref, pos,
02228             settle_mOrmT, settle_mHrmT, settle_ra,
02229             settle_rb, settle_rc, settle_rra);        
02230     for (int i = 0; i < 3; ++i ) {
02231       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02232       posNew[ig+i] = pos[i];
02233     }
02234   }
02235 
02236   int posParam = 0;
02237   for (int j=0;j < rattleList.size();++j) {
02238 
02239     BigReal refx[10];
02240     BigReal refy[10];
02241     BigReal refz[10];
02242 
02243     BigReal posx[10];
02244     BigReal posy[10];
02245     BigReal posz[10];
02246 
02247     int ig = rattleList[j].ig;
02248     int icnt = rattleList[j].icnt;
02249     int hgs = atom[ig].hydrogenGroupSize;
02250     for (int i = 0; i < hgs; ++i ) {
02251       ref[i] = atom[ig+i].position;
02252       pos[i] = atom[ig+i].position;
02253       if (!(fixedAtomsOn && atom[ig+i].atomFixed))
02254         pos[i] += atom[ig+i].velocity * dt;
02255       refx[i] = ref[i].x;
02256       refy[i] = ref[i].y;
02257       refz[i] = ref[i].z;
02258       posx[i] = pos[i].x;
02259       posy[i] = pos[i].y;
02260       posz[i] = pos[i].z;
02261     }
02262 
02263 
02264     bool done;
02265     bool consFailure;
02266     if (icnt == 1) {
02267       rattlePair<1>(&rattleParam[posParam],
02268         refx, refy, refz,
02269         posx, posy, posz);
02270       done = true;
02271       consFailure = false;
02272     } else {
02273       rattleN(icnt, &rattleParam[posParam],
02274         refx, refy, refz,
02275         posx, posy, posz,
02276         tol2, maxiter,
02277         done, consFailure);
02278     }
02279 
02280     // Advance position in rattleParam
02281     posParam += icnt;
02282 
02283     for (int i = 0; i < hgs; ++i ) {
02284       pos[i].x = posx[i];
02285       pos[i].y = posy[i];
02286       pos[i].z = posz[i];
02287     }
02288 
02289     for (int i = 0; i < hgs; ++i ) {
02290       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02291       posNew[ig+i] = pos[i];
02292     }
02293 
02294     if ( consFailure ) {
02295       if ( dieOnError ) {
02296         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02297         << (atom[ig].id + 1) << "!\n" << endi;
02298         return -1;  // triggers early exit
02299       } else {
02300         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02301         << (atom[ig].id + 1) << "!\n" << endi;
02302       }
02303     } else if ( ! done ) {
02304       if ( dieOnError ) {
02305         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02306         << (atom[ig].id + 1) << "!\n" << endi;
02307         return -1;  // triggers early exit
02308       } else {
02309         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02310         << (atom[ig].id + 1) << "!\n" << endi;
02311       }
02312     }
02313 
02314   }
02315 
02316   // Finally, we have to go through atoms that are not involved in rattle just so that we have
02317   // their positions and velocities up-to-date in posNew and velNew
02318   for (int j=0;j < noconstList.size();++j) {
02319     int ig = noconstList[j];
02320     int hgs = atom[ig].hydrogenGroupSize;
02321     for (int i = 0; i < hgs; ++i ) {
02322       velNew[ig+i] = atom[ig+i].velocity;
02323       posNew[ig+i] = atom[ig+i].position;
02324     }
02325   }
02326 
02327   if ( invdt == 0 ) {
02328     for (int ig = 0; ig < numAtoms; ++ig )
02329       atom[ig].position = posNew[ig];
02330   } else if ( virial == 0 ) {
02331     for (int ig = 0; ig < numAtoms; ++ig )
02332       atom[ig].velocity = velNew[ig];
02333   } else {
02334     Tensor wc;  // constraint virial
02335     addRattleForce(invdt, wc);
02336     *virial += wc;
02337   }
02338 
02339   return 0;
02340 }
02341 
02342 //  RATTLE algorithm from Allen & Tildesley
02343 int HomePatch::rattle1old(const BigReal timestep, Tensor *virial, 
02344     SubmitReduction *ppreduction)
02345 {
02346   Molecule *mol = Node::Object()->molecule;
02347   SimParameters *simParams = Node::Object()->simParameters;
02348   const int fixedAtomsOn = simParams->fixedAtomsOn;
02349   const int useSettle = simParams->useSettle;
02350   const BigReal dt = timestep / TIMEFACTOR;
02351   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02352   BigReal tol2 = 2.0 * simParams->rigidTol;
02353   int maxiter = simParams->rigidIter;
02354   int dieOnError = simParams->rigidDie;
02355   int i, iter;
02356   BigReal dsq[10], tmp;
02357   int ial[10], ibl[10];
02358   Vector ref[10];  // reference position
02359   Vector refab[10];  // reference vector
02360   Vector pos[10];  // new position
02361   Vector vel[10];  // new velocity
02362   Vector netdp[10];  // total momentum change from constraint
02363   BigReal rmass[10];  // 1 / mass
02364   int fixed[10];  // is atom fixed?
02365   Tensor wc;  // constraint virial
02366   BigReal idz, zmin;
02367   int nslabs;
02368 
02369   // Initialize the settle algorithm with water parameters
02370   // settle1() assumes all waters are identical,
02371   // and will generate bad results if they are not.
02372   // XXX this will move to Molecule::build_atom_status when that 
02373   // version is debugged
02374   if ( ! settle_initialized ) {
02375     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02376       // find a water
02377       if (atom[ig].rigidBondLength > 0) {
02378         int oatm;
02379         if (simParams->watmodel == WAT_SWM4) {
02380           oatm = ig+3;  // skip over Drude and Lonepair
02381           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02382           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02383         }
02384         else {
02385           oatm = ig+1;
02386           // Avoid using the Om site to set this by mistake
02387           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02388             oatm += 1;
02389           }
02390         }
02391 
02392         // initialize settle water parameters
02393         settle1init(atom[ig].mass, atom[oatm].mass, 
02394                     atom[ig].rigidBondLength,
02395                     atom[oatm].rigidBondLength,
02396                     settle_mOrmT, settle_mHrmT, settle_ra,
02397                     settle_rb, settle_rc, settle_rra);
02398         settle_initialized = 1;
02399         break; // done with init
02400       }
02401     }
02402   }
02403 
02404   if (ppreduction) {
02405     nslabs = simParams->pressureProfileSlabs;
02406     idz = nslabs/lattice.c().z;
02407     zmin = lattice.origin().z - 0.5*lattice.c().z;
02408   }
02409 
02410   // Size of a hydrogen group for water
02411   int wathgsize = 3;
02412   int watmodel = simParams->watmodel;
02413   if (watmodel == WAT_TIP4) wathgsize = 4;
02414   else if (watmodel == WAT_SWM4) wathgsize = 5;
02415   
02416   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02417     int hgs = atom[ig].hydrogenGroupSize;
02418     if ( hgs == 1 ) continue;  // only one atom in group
02419     // cache data in local arrays and integrate positions normally
02420     int anyfixed = 0;
02421     for ( i = 0; i < hgs; ++i ) {
02422       ref[i] = atom[ig+i].position;
02423       pos[i] = atom[ig+i].position;
02424       vel[i] = atom[ig+i].velocity;
02425       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02426       //printf("rmass of %i is %f\n", ig+i, rmass[i]);
02427       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02428       //printf("fixed status of %i is %i\n", i, fixed[i]);
02429       // undo addVelocityToPosition to get proper reference coordinates
02430       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
02431     }
02432     int icnt = 0;
02433     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02434       if (hgs != wathgsize) {
02435         char errmsg[256];
02436         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02437                          "but the specified water model requires %d atoms.\n",
02438                          atom[ig].id+1, hgs, wathgsize);
02439         NAMD_die(errmsg);
02440       }
02441       // Use SETTLE for water unless some of the water atoms are fixed,
02442       if (useSettle && !anyfixed) {
02443         if (simParams->watmodel == WAT_SWM4) {
02444           // SWM4 ordering:  O D LP H1 H2
02445           // do swap(O,LP) and call settle with subarray O H1 H2
02446           // swap back after we return
02447           Vector lp_ref = ref[2];
02448           Vector lp_pos = pos[2];
02449           Vector lp_vel = vel[2];
02450           ref[2] = ref[0];
02451           pos[2] = pos[0];
02452           vel[2] = vel[0];
02453           settle1(ref+2, pos+2, vel+2, invdt,
02454                   settle_mOrmT, settle_mHrmT, settle_ra,
02455                   settle_rb, settle_rc, settle_rra);
02456           ref[0] = ref[2];
02457           pos[0] = pos[2];
02458           vel[0] = vel[2];
02459           ref[2] = lp_ref;
02460           pos[2] = lp_pos;
02461           vel[2] = lp_vel;
02462           // determine for LP updated pos and vel
02463           swm4_omrepos(ref, pos, vel, invdt);
02464         }
02465         else {
02466             settle1(ref, pos, vel, invdt,
02467                     settle_mOrmT, settle_mHrmT, settle_ra,
02468                     settle_rb, settle_rc, settle_rra);            
02469           if (simParams->watmodel == WAT_TIP4) {
02470             tip4_omrepos(ref, pos, vel, invdt);
02471           }
02472         }
02473 
02474         // which slab the hydrogen group will belong to
02475         // for pprofile calculations.
02476         int ppoffset, partition;
02477         if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02478           atom[ig+i].position = pos[i];
02479         } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02480           atom[ig+i].velocity = vel[i];
02481         } else for ( i = 0; i < wathgsize; ++i ) {
02482           Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
02483           Tensor vir = outer(df, ref[i]);
02484           wc += vir;
02485           f[Results::normal][ig+i] += df;
02486           atom[ig+i].velocity = vel[i];
02487           if (ppreduction) {
02488             // put all the atoms from a water in the same slab.  Atom 0
02489             // should be the parent atom.
02490             if (!i) {
02491               BigReal z = pos[i].z;
02492               partition = atom[ig].partition;
02493               int slab = (int)floor((z-zmin)*idz);
02494               if (slab < 0) slab += nslabs;
02495               else if (slab >= nslabs) slab -= nslabs;
02496               ppoffset = 3*(slab + nslabs*partition);
02497             }
02498             ppreduction->item(ppoffset  ) += vir.xx;
02499             ppreduction->item(ppoffset+1) += vir.yy;
02500             ppreduction->item(ppoffset+2) += vir.zz;
02501           }
02502         }
02503         continue;
02504       }
02505       if ( !(fixed[1] && fixed[2]) ) {
02506         dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02507       }
02508     }
02509     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02510       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02511         if ( !(fixed[0] && fixed[i]) ) {
02512           dsq[icnt] = tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
02513         }
02514       }
02515     }
02516     if ( icnt == 0 ) continue;  // no constraints
02517     for ( i = 0; i < icnt; ++i ) {
02518       refab[i] = ref[ial[i]] - ref[ibl[i]];
02519     }
02520     for ( i = 0; i < hgs; ++i ) {
02521       netdp[i] = 0.;
02522     }
02523     int done;
02524     int consFailure;
02525     for ( iter = 0; iter < maxiter; ++iter ) {
02526 //if (iter > 0) CkPrintf("iteration %d\n", iter);
02527       done = 1;
02528       consFailure = 0;
02529       for ( i = 0; i < icnt; ++i ) {
02530         int a = ial[i];  int b = ibl[i];
02531         Vector pab = pos[a] - pos[b];
02532               BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
02533         BigReal rabsq = dsq[i];
02534         BigReal diffsq = rabsq - pabsq;
02535         if ( fabs(diffsq) > (rabsq * tol2) ) {
02536           Vector &rab = refab[i];
02537           BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
02538           if ( rpab < ( rabsq * 1.0e-6 ) ) {
02539             done = 0;
02540             consFailure = 1;
02541             continue;
02542           }
02543           BigReal rma = rmass[a];
02544           BigReal rmb = rmass[b];
02545           BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
02546           Vector dp = rab * gab;
02547           pos[a] += rma * dp;
02548           pos[b] -= rmb * dp;
02549           if ( invdt != 0. ) {
02550             dp *= invdt;
02551             netdp[a] += dp;
02552             netdp[b] -= dp;
02553           }
02554           done = 0;
02555         }
02556       }
02557       if ( done ) break;
02558     }
02559 
02560     if ( consFailure ) {
02561       if ( dieOnError ) {
02562         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02563         << (atom[ig].id + 1) << "!\n" << endi;
02564         return -1;  // triggers early exit
02565       } else {
02566         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02567         << (atom[ig].id + 1) << "!\n" << endi;
02568       }
02569     } else if ( ! done ) {
02570       if ( dieOnError ) {
02571         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02572         << (atom[ig].id + 1) << "!\n" << endi;
02573         return -1;  // triggers early exit
02574       } else {
02575         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02576         << (atom[ig].id + 1) << "!\n" << endi;
02577       }
02578     }
02579 
02580     // store data back to patch
02581     int ppoffset, partition;
02582     if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
02583       atom[ig+i].position = pos[i];
02584     } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
02585       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02586     } else for ( i = 0; i < hgs; ++i ) {
02587       Force df = netdp[i] * invdt;
02588       Tensor vir = outer(df, ref[i]);
02589       wc += vir;
02590       f[Results::normal][ig+i] += df;
02591       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02592       if (ppreduction) {
02593         if (!i) {
02594           BigReal z = pos[i].z;
02595           int partition = atom[ig].partition;
02596           int slab = (int)floor((z-zmin)*idz);
02597           if (slab < 0) slab += nslabs;
02598           else if (slab >= nslabs) slab -= nslabs;
02599           ppoffset = 3*(slab + nslabs*partition);
02600         }
02601         ppreduction->item(ppoffset  ) += vir.xx;
02602         ppreduction->item(ppoffset+1) += vir.yy;
02603         ppreduction->item(ppoffset+2) += vir.zz;
02604       }
02605     }
02606   }
02607   if ( dt && virial ) *virial += wc;
02608 
02609   return 0;
02610 }
02611 
02612 //  RATTLE algorithm from Allen & Tildesley
02613 void HomePatch::rattle2(const BigReal timestep, Tensor *virial)
02614 {
02615   Molecule *mol = Node::Object()->molecule;
02616   SimParameters *simParams = Node::Object()->simParameters;
02617   const int fixedAtomsOn = simParams->fixedAtomsOn;
02618   const int useSettle = simParams->useSettle;
02619   const BigReal dt = timestep / TIMEFACTOR;
02620   Tensor wc;  // constraint virial
02621   BigReal tol = simParams->rigidTol;
02622   int maxiter = simParams->rigidIter;
02623   int dieOnError = simParams->rigidDie;
02624   int i, iter;
02625   BigReal dsqi[10], tmp;
02626   int ial[10], ibl[10];
02627   Vector ref[10];  // reference position
02628   Vector refab[10];  // reference vector
02629   Vector vel[10];  // new velocity
02630   BigReal rmass[10];  // 1 / mass
02631   BigReal redmass[10];  // reduced mass
02632   int fixed[10];  // is atom fixed?
02633   
02634   // Size of a hydrogen group for water
02635   int wathgsize = 3;
02636   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02637   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02638 
02639   //  CkPrintf("In rattle2!\n");
02640   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02641     //    CkPrintf("ig=%d\n",ig);
02642     int hgs = atom[ig].hydrogenGroupSize;
02643     if ( hgs == 1 ) continue;  // only one atom in group
02644     // cache data in local arrays and integrate positions normally
02645     int anyfixed = 0;
02646     for ( i = 0; i < hgs; ++i ) {
02647       ref[i] = atom[ig+i].position;
02648       vel[i] = atom[ig+i].velocity;
02649       rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
02650       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02651       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02652     }
02653     int icnt = 0;
02654     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02655       if (hgs != wathgsize) {
02656         NAMD_bug("Hydrogen group error caught in rattle2().");
02657       }
02658       // Use SETTLE for water unless some of the water atoms are fixed,
02659       if (useSettle && !anyfixed) {
02660         if (simParams->watmodel == WAT_SWM4) {
02661           // SWM4 ordering:  O D LP H1 H2
02662           // do swap(O,LP) and call settle with subarray O H1 H2
02663           // swap back after we return
02664           Vector lp_ref = ref[2];
02665           // Vector lp_vel = vel[2];
02666           ref[2] = ref[0];
02667           vel[2] = vel[0];
02668           settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
02669           ref[0] = ref[2];
02670           vel[0] = vel[2];
02671           ref[2] = lp_ref;
02672           // vel[2] = vel[0];  // set LP vel to O vel
02673         }
02674         else {
02675           settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
02676           if (simParams->watmodel == WAT_TIP4) {
02677             vel[3] = vel[0];
02678           }
02679         }
02680         for (i=0; i<hgs; i++) {
02681           atom[ig+i].velocity = vel[i];
02682         }
02683         continue;
02684       }
02685       if ( !(fixed[1] && fixed[2]) ) {
02686         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02687         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02688       }
02689     }
02690     //    CkPrintf("Loop 2\n");
02691     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02692       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02693         if ( !(fixed[0] && fixed[i]) ) {
02694           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02695           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02696           ibl[icnt] = i;  ++icnt;
02697         }
02698       }
02699     }
02700     if ( icnt == 0 ) continue;  // no constraints
02701     //    CkPrintf("Loop 3\n");
02702     for ( i = 0; i < icnt; ++i ) {
02703       refab[i] = ref[ial[i]] - ref[ibl[i]];
02704     }
02705     //    CkPrintf("Loop 4\n");
02706     int done;
02707     for ( iter = 0; iter < maxiter; ++iter ) {
02708       done = 1;
02709       for ( i = 0; i < icnt; ++i ) {
02710         int a = ial[i];  int b = ibl[i];
02711         Vector vab = vel[a] - vel[b];
02712         Vector &rab = refab[i];
02713         BigReal rabsqi = dsqi[i];
02714         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02715         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02716           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02717           wc += outer(dp,rab);
02718           vel[a] += rmass[a] * dp;
02719           vel[b] -= rmass[b] * dp;
02720           done = 0;
02721         }
02722       }
02723       if ( done ) break;
02724       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02725     }
02726     if ( ! done ) {
02727       if ( dieOnError ) {
02728         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02729       } else {
02730         iout << iWARN <<
02731           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02732       }
02733     }
02734     // store data back to patch
02735     for ( i = 0; i < hgs; ++i ) {
02736       atom[ig+i].velocity = vel[i];
02737     }
02738   }
02739   //  CkPrintf("Leaving rattle2!\n");
02740   // check that there isn't a constant needed here!
02741   *virial += wc / ( 0.5 * dt );
02742 
02743 }
02744 
02745 
02746 //  Adjust gradients for minimizer
02747 void HomePatch::minimize_rattle2(const BigReal timestep, Tensor *virial, bool forces)
02748 {
02749   Molecule *mol = Node::Object()->molecule;
02750   SimParameters *simParams = Node::Object()->simParameters;
02751   Force *f1 = f[Results::normal].begin();
02752   const int fixedAtomsOn = simParams->fixedAtomsOn;
02753   const int useSettle = simParams->useSettle;
02754   const BigReal dt = timestep / TIMEFACTOR;
02755   Tensor wc;  // constraint virial
02756   BigReal tol = simParams->rigidTol;
02757   int maxiter = simParams->rigidIter;
02758   int dieOnError = simParams->rigidDie;
02759   int i, iter;
02760   BigReal dsqi[10], tmp;
02761   int ial[10], ibl[10];
02762   Vector ref[10];  // reference position
02763   Vector refab[10];  // reference vector
02764   Vector vel[10];  // new velocity
02765   BigReal rmass[10];  // 1 / mass
02766   BigReal redmass[10];  // reduced mass
02767   int fixed[10];  // is atom fixed?
02768   
02769   // Size of a hydrogen group for water
02770   int wathgsize = 3;
02771   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02772   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02773 
02774   //  CkPrintf("In rattle2!\n");
02775   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02776     //    CkPrintf("ig=%d\n",ig);
02777     int hgs = atom[ig].hydrogenGroupSize;
02778     if ( hgs == 1 ) continue;  // only one atom in group
02779     // cache data in local arrays and integrate positions normally
02780     int anyfixed = 0;
02781     for ( i = 0; i < hgs; ++i ) {
02782       ref[i] = atom[ig+i].position;
02783       vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
02784       rmass[i] = 1.0;
02785       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02786       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02787     }
02788     int icnt = 0;
02789     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02790       if (hgs != wathgsize) {
02791         NAMD_bug("Hydrogen group error caught in rattle2().");
02792       }
02793       // Use SETTLE for water unless some of the water atoms are fixed,
02794       if (useSettle && !anyfixed) {
02795         if (simParams->watmodel == WAT_SWM4) {
02796           // SWM4 ordering:  O D LP H1 H2
02797           // do swap(O,LP) and call settle with subarray O H1 H2
02798           // swap back after we return
02799           Vector lp_ref = ref[2];
02800           // Vector lp_vel = vel[2];
02801           ref[2] = ref[0];
02802           vel[2] = vel[0];
02803           settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
02804           ref[0] = ref[2];
02805           vel[0] = vel[2];
02806           ref[2] = lp_ref;
02807           // vel[2] = vel[0];  // set LP vel to O vel
02808         }
02809         else {
02810           settle2(1.0, 1.0, ref, vel, dt, virial);
02811           if (simParams->watmodel == WAT_TIP4) {
02812             vel[3] = vel[0];
02813           }
02814         }
02815         for (i=0; i<hgs; i++) {
02816           ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
02817         }
02818         continue;
02819       }
02820       if ( !(fixed[1] && fixed[2]) ) {
02821         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02822         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02823       }
02824     }
02825     //    CkPrintf("Loop 2\n");
02826     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02827       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02828         if ( !(fixed[0] && fixed[i]) ) {
02829           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02830           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02831           ibl[icnt] = i;  ++icnt;
02832         }
02833       }
02834     }
02835     if ( icnt == 0 ) continue;  // no constraints
02836     //    CkPrintf("Loop 3\n");
02837     for ( i = 0; i < icnt; ++i ) {
02838       refab[i] = ref[ial[i]] - ref[ibl[i]];
02839     }
02840     //    CkPrintf("Loop 4\n");
02841     int done;
02842     for ( iter = 0; iter < maxiter; ++iter ) {
02843       done = 1;
02844       for ( i = 0; i < icnt; ++i ) {
02845         int a = ial[i];  int b = ibl[i];
02846         Vector vab = vel[a] - vel[b];
02847         Vector &rab = refab[i];
02848         BigReal rabsqi = dsqi[i];
02849         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02850         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02851           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02852           wc += outer(dp,rab);
02853           vel[a] += rmass[a] * dp;
02854           vel[b] -= rmass[b] * dp;
02855           done = 0;
02856         }
02857       }
02858       if ( done ) break;
02859       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02860     }
02861     if ( ! done ) {
02862       if ( dieOnError ) {
02863         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02864       } else {
02865         iout << iWARN <<
02866           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02867       }
02868     }
02869     // store data back to patch
02870     for ( i = 0; i < hgs; ++i ) {
02871       ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
02872     }
02873   }
02874   //  CkPrintf("Leaving rattle2!\n");
02875   // check that there isn't a constant needed here!
02876   *virial += wc / ( 0.5 * dt );
02877 
02878 }
02879 
02880 
02881 // BEGIN LA
02882 void HomePatch::loweAndersenVelocities()
02883 {
02884     DebugM(2, "loweAndersenVelocities\n");
02885     Molecule *mol = Node::Object()->molecule;
02886     SimParameters *simParams = Node::Object()->simParameters;
02887     v.resize(numAtoms);
02888     for (int i = 0; i < numAtoms; ++i) {
02889         //v[i] = p[i];
02890         // co-opt CompAtom structure to pass velocity and mass information
02891         v[i].position = atom[i].velocity;
02892         v[i].charge = atom[i].mass;
02893     }
02894     DebugM(2, "loweAndersenVelocities\n");
02895 }
02896 
02897 void HomePatch::loweAndersenFinish()
02898 {
02899     DebugM(2, "loweAndersenFinish\n");
02900     v.resize(0);
02901 }
02902 // END LA
02903 
02904 //LCPO
02905 void HomePatch::setLcpoType() {
02906   Molecule *mol = Node::Object()->molecule;
02907   const int *lcpoParamType = mol->getLcpoParamType();
02908 
02909   lcpoType.resize(numAtoms);
02910   for (int i = 0; i < numAtoms; i++) {
02911     lcpoType[i] = lcpoParamType[pExt[i].id];
02912   }
02913 }
02914 
02915 //set intrinsic radii of atom when doMigration
02916 void HomePatch::setGBISIntrinsicRadii() {
02917   intRad.resize(numAtoms*2);
02918   intRad.setall(0);
02919   Molecule *mol = Node::Object()->molecule;
02920   SimParameters *simParams = Node::Object()->simParameters;
02921   Real offset = simParams->coulomb_radius_offset;
02922   for (int i = 0; i < numAtoms; i++) {
02923     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
02924     Real screen = MassToScreen(atom[i].mass);//same
02925     intRad[2*i+0] = rad - offset;//r0
02926     intRad[2*i+1] = screen*(rad - offset);//s0
02927   }
02928 }
02929 
02930 //compute born radius after phase 1, before phase 2
02931 void HomePatch::gbisComputeAfterP1() {
02932 
02933   SimParameters *simParams = Node::Object()->simParameters;
02934   BigReal alphaMax = simParams->alpha_max;
02935   BigReal delta = simParams->gbis_delta;
02936   BigReal beta = simParams->gbis_beta;
02937   BigReal gamma = simParams->gbis_gamma;
02938   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
02939 
02940   BigReal rhoi;
02941   BigReal rhoi0;
02942   //calculate bornRad from psiSum
02943   for (int i = 0; i < numAtoms; i++) {
02944     rhoi0 = intRad[2*i];
02945     rhoi = rhoi0+coulomb_radius_offset;
02946     psiFin[i] += psiSum[i];
02947     psiFin[i] *= rhoi0;
02948     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
02949     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
02950 #ifdef PRINT_COMP
02951     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
02952 #endif
02953   }
02954 
02955   gbisP2Ready();
02956 }
02957 
02958 //compute dHdrPrefix after phase 2, before phase 3
02959 void HomePatch::gbisComputeAfterP2() {
02960 
02961   SimParameters *simParams = Node::Object()->simParameters;
02962   BigReal delta = simParams->gbis_delta;
02963   BigReal beta = simParams->gbis_beta;
02964   BigReal gamma = simParams->gbis_gamma;
02965   BigReal epsilon_s = simParams->solvent_dielectric;
02966   BigReal epsilon_p = simParams->dielectric;
02967   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
02968   BigReal epsilon_p_i = 1/simParams->dielectric;
02969   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
02970   BigReal kappa = simParams->kappa;
02971   BigReal fij, expkappa, Dij, dEdai, dedasum;
02972   BigReal rhoi, rhoi0, psii, nbetapsi;
02973   BigReal gammapsi2, tanhi, daidr;
02974   for (int i = 0; i < numAtoms; i++) {
02975     //add diagonal dEda term
02976     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
02977     fij = bornRad[i];//inf
02978     expkappa = exp(-kappa*fij);//0
02979     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
02980     //calculate dHij prefix
02981     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
02982                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
02983     dHdrPrefix[i] += dEdai;
02984     dedasum = dHdrPrefix[i];
02985 
02986     rhoi0 = intRad[2*i];
02987     rhoi = rhoi0+coulomb_radius_offset;
02988     psii = psiFin[i];
02989     nbetapsi = -beta*psii;
02990     gammapsi2 = gamma*psii*psii;
02991     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
02992     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
02993            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
02994     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
02995 #ifdef PRINT_COMP
02996     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
02997 #endif
02998   }
02999   gbisP3Ready();
03000 }
03001 
03002 //send born radius to proxies to begin phase 2
03003 void HomePatch::gbisP2Ready() {
03004   if (proxy.size() > 0) {
03005     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03006     for (int i = 0; i < proxy.size(); i++) {
03007       int node = proxy[i];
03008       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
03009       msg->patch = patchID;
03010       msg->origPe = CkMyPe();
03011       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
03012       msg->destPe = node;
03013       int seq = flags.sequence;
03014       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03015       SET_PRIORITY(msg,seq,priority);
03016       cp[node].recvData(msg);
03017     }
03018   }
03019   Patch::gbisP2Ready();
03020 }
03021 
03022 //send dHdrPrefix to proxies to begin phase 3
03023 void HomePatch::gbisP3Ready() {
03024   if (proxy.size() > 0) {
03025     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03026     //only nonzero message should be sent for doFullElec
03027     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
03028     for (int i = 0; i < proxy.size(); i++) {
03029       int node = proxy[i];
03030       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
03031       msg->patch = patchID;
03032       msg->dHdrPrefixLen = msgAtoms;
03033       msg->origPe = CkMyPe();
03034       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
03035       msg->destPe = node;
03036       int seq = flags.sequence;
03037       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03038       SET_PRIORITY(msg,seq,priority);
03039       cp[node].recvData(msg);
03040     }
03041   }
03042   Patch::gbisP3Ready();
03043 }
03044 
03045 //receive proxy results from phase 1
03046 void HomePatch::receiveResult(ProxyGBISP1ResultMsg *msg) {
03047   ++numGBISP1Arrived;
03048     for ( int i = 0; i < msg->psiSumLen; ++i ) {
03049       psiFin[i] += msg->psiSum[i];
03050     }
03051   delete msg;
03052 
03053   if (flags.doNonbonded) {
03054     //awaken if phase 1 done
03055     if (phase1BoxClosedCalled == true &&
03056         numGBISP1Arrived==proxy.size() ) {
03057         sequencer->awaken();
03058     }
03059   } else {
03060     //awaken if all phases done on noWork step
03061     if (boxesOpen == 0 &&
03062         numGBISP1Arrived == proxy.size() &&
03063         numGBISP2Arrived == proxy.size() &&
03064         numGBISP3Arrived == proxy.size()) {
03065       sequencer->awaken();
03066     }
03067   }
03068 }
03069 
03070 //receive proxy results from phase 2
03071 void HomePatch::receiveResult(ProxyGBISP2ResultMsg *msg) {
03072   ++numGBISP2Arrived;
03073   //accumulate dEda
03074   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
03075     dHdrPrefix[i] += msg->dEdaSum[i];
03076   }
03077   delete msg;
03078 
03079   if (flags.doNonbonded) {
03080     //awaken if phase 2 done
03081     if (phase2BoxClosedCalled == true &&
03082         numGBISP2Arrived==proxy.size() ) {
03083       sequencer->awaken();
03084     }
03085   } else {
03086     //awaken if all phases done on noWork step
03087     if (boxesOpen == 0 &&
03088         numGBISP1Arrived == proxy.size() &&
03089         numGBISP2Arrived == proxy.size() &&
03090         numGBISP3Arrived == proxy.size()) {
03091       sequencer->awaken();
03092     }
03093   }
03094 }
03095 
03096 //  MOLLY algorithm part 1
03097 void HomePatch::mollyAverage()
03098 {
03099   Molecule *mol = Node::Object()->molecule;
03100   SimParameters *simParams = Node::Object()->simParameters;
03101   BigReal tol = simParams->mollyTol;
03102   int maxiter = simParams->mollyIter;
03103   int i, iter;
03104   HGArrayBigReal dsq;
03105   BigReal tmp;
03106   HGArrayInt ial, ibl;
03107   HGArrayVector ref;  // reference position
03108   HGArrayVector refab;  // reference vector
03109   HGArrayBigReal rmass;  // 1 / mass
03110   BigReal *lambda;  // Lagrange multipliers
03111   CompAtom *avg;  // averaged position
03112   int numLambdas = 0;
03113   HGArrayInt fixed;  // is atom fixed?
03114 
03115   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
03116   p_avg.resize(numAtoms);
03117   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
03118 
03119   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03120     int hgs = atom[ig].hydrogenGroupSize;
03121     if ( hgs == 1 ) continue;  // only one atom in group
03122         for ( i = 0; i < hgs; ++i ) {
03123           ref[i] = atom[ig+i].position;
03124           rmass[i] = 1. / atom[ig+i].mass;
03125           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03126           if ( fixed[i] ) rmass[i] = 0.;
03127         }
03128         avg = &(p_avg[ig]);
03129         int icnt = 0;
03130 
03131   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
03132           if ( hgs != 3 ) {
03133             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
03134           }
03135           if ( !(fixed[1] && fixed[2]) ) {
03136             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03137           }
03138         }
03139         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03140     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
03141             if ( !(fixed[0] && fixed[i]) ) {
03142               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03143             }
03144           }
03145         }
03146         if ( icnt == 0 ) continue;  // no constraints
03147         numLambdas += icnt;
03148         molly_lambda.resize(numLambdas);
03149         lambda = &(molly_lambda[numLambdas - icnt]);
03150         for ( i = 0; i < icnt; ++i ) {
03151           refab[i] = ref[ial[i]] - ref[ibl[i]];
03152         }
03153         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
03154         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
03155         if ( iter == maxiter ) {
03156           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
03157         }
03158   }
03159 
03160   // for ( i=0; i<numAtoms; ++i ) {
03161   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
03162   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
03163   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
03164   //    }
03165   // }
03166 
03167 }
03168 
03169 
03170 //  MOLLY algorithm part 2
03171 void HomePatch::mollyMollify(Tensor *virial)
03172 {
03173   Molecule *mol = Node::Object()->molecule;
03174   SimParameters *simParams = Node::Object()->simParameters;
03175   Tensor wc;  // constraint virial
03176   int i;
03177   HGArrayInt ial, ibl;
03178   HGArrayVector ref;  // reference position
03179   CompAtom *avg;  // averaged position
03180   HGArrayVector refab;  // reference vector
03181   HGArrayVector force;  // new force
03182   HGArrayBigReal rmass;  // 1 / mass
03183   BigReal *lambda;  // Lagrange multipliers
03184   int numLambdas = 0;
03185   HGArrayInt fixed;  // is atom fixed?
03186 
03187   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03188     int hgs = atom[ig].hydrogenGroupSize;
03189     if (hgs == 1 ) continue;  // only one atom in group
03190         for ( i = 0; i < hgs; ++i ) {
03191           ref[i] = atom[ig+i].position;
03192           force[i] = f[Results::slow][ig+i];
03193           rmass[i] = 1. / atom[ig+i].mass;
03194           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03195           if ( fixed[i] ) rmass[i] = 0.;
03196         }
03197         int icnt = 0;
03198         // c-ji I'm only going to mollify water for now
03199   if ( atom[ig].rigidBondLength > 0 ) {  // for water
03200           if ( hgs != 3 ) {
03201             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
03202           }
03203           if ( !(fixed[1] && fixed[2]) ) {
03204             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03205           }
03206         }
03207         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03208     if ( atom[ig+i].rigidBondLength > 0 ) {
03209             if ( !(fixed[0] && fixed[i]) ) {
03210               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03211             }
03212           }
03213         }
03214 
03215         if ( icnt == 0 ) continue;  // no constraints
03216         lambda = &(molly_lambda[numLambdas]);
03217         numLambdas += icnt;
03218         for ( i = 0; i < icnt; ++i ) {
03219           refab[i] = ref[ial[i]] - ref[ibl[i]];
03220         }
03221         avg = &(p_avg[ig]);
03222         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
03223         // store data back to patch
03224         for ( i = 0; i < hgs; ++i ) {
03225           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
03226           f[Results::slow][ig+i] = force[i];
03227         }
03228   }
03229   // check that there isn't a constant needed here!
03230   *virial += wc;
03231   p_avg.resize(0);
03232 }
03233 
03234 void HomePatch::checkpoint(void) {
03235   checkpoint_atom.copy(atom);
03236   checkpoint_lattice = lattice;
03237 
03238   // DMK - Atom Separation (water vs. non-water)
03239   #if NAMD_SeparateWaters != 0
03240     checkpoint_numWaterAtoms = numWaterAtoms;
03241   #endif
03242 }
03243 
03244 void HomePatch::revert(void) {
03245   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03246 
03247   atom.copy(checkpoint_atom);
03248   numAtoms = atom.size();
03249   lattice = checkpoint_lattice;
03250 
03251   doAtomUpdate = true;
03252   rattleListValid = false;
03253 
03254   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03255 
03256   // DMK - Atom Separation (water vs. non-water)
03257   #if NAMD_SeparateWaters != 0
03258     numWaterAtoms = checkpoint_numWaterAtoms;
03259   #endif
03260 }
03261 
03262 void HomePatch::exchangeCheckpoint(int scriptTask, int &bpc) {  // initiating replica
03263   SimParameters *simParams = Node::Object()->simParameters;
03264   checkpoint_task = scriptTask;
03265   const int remote = simParams->scriptIntArg1;
03266   const char *key = simParams->scriptStringArg1;
03267   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
03268 }
03269 
03270 void HomePatch::recvCheckpointReq(int task, const char *key, int replica, int pe) {  // responding replica
03271   if ( task == SCRIPT_CHECKPOINT_FREE ) {
03272     if ( ! checkpoints.count(key) ) {
03273       NAMD_die("Unable to free checkpoint, requested key was never stored.");
03274     }
03275     delete checkpoints[key];
03276     checkpoints.erase(key);
03277     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
03278     return;
03279   }
03280   CheckpointAtomsMsg *msg;
03281   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
03282     if ( ! checkpoints.count(key) ) {
03283       NAMD_die("Unable to load checkpoint, requested key was never stored.");
03284     }
03285     checkpoint_t &cp = *checkpoints[key];
03286     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
03287     msg->lattice = cp.lattice;
03288     msg->berendsenPressure_count = cp.berendsenPressure_count;
03289     msg->numAtoms = cp.numAtoms;
03290     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
03291   } else {
03292     msg = new (0,1,0) CheckpointAtomsMsg;
03293   }
03294   msg->pid = patchID;
03295   msg->replica = CmiMyPartition();
03296   msg->pe = CkMyPe();
03297   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
03298 }
03299 
03300 void HomePatch::recvCheckpointLoad(CheckpointAtomsMsg *msg) {  // initiating replica
03301   SimParameters *simParams = Node::Object()->simParameters;
03302   const int remote = simParams->scriptIntArg1;
03303   const char *key = simParams->scriptStringArg1;
03304   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
03305     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
03306   }
03307   if ( msg->replica != remote ) {
03308     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
03309   }
03310   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03311     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
03312     strcpy(newmsg->key,key);
03313     newmsg->lattice = lattice;
03314     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
03315     newmsg->pid = patchID;
03316     newmsg->pe = CkMyPe();
03317     newmsg->replica = CmiMyPartition();
03318     newmsg->numAtoms = numAtoms;
03319     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03320     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
03321   }
03322   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03323     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03324     lattice = msg->lattice;
03325     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
03326     numAtoms = msg->numAtoms;
03327     atom.resize(numAtoms);
03328     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03329     doAtomUpdate = true;
03330     rattleListValid = false;
03331     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03332   }
03333   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
03334     recvCheckpointAck();
03335   }
03336   delete msg;
03337 }
03338 
03339 void HomePatch::recvCheckpointStore(CheckpointAtomsMsg *msg) {  // responding replica
03340   if ( ! checkpoints.count(msg->key) ) {
03341     checkpoints[msg->key] = new checkpoint_t;
03342   }
03343   checkpoint_t &cp = *checkpoints[msg->key];
03344   cp.lattice = msg->lattice;
03345   cp.berendsenPressure_count = msg->berendsenPressure_count;
03346   cp.numAtoms = msg->numAtoms;
03347   cp.atoms.resize(cp.numAtoms);
03348   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
03349   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
03350   delete msg;
03351 }
03352 
03353 void HomePatch::recvCheckpointAck() {  // initiating replica
03354   CkpvAccess(_qd)->process();
03355 }
03356 
03357 
03358 void HomePatch::exchangeAtoms(int scriptTask) {
03359   SimParameters *simParams = Node::Object()->simParameters;
03360   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
03361   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
03362     exchange_dst = (int) simParams->scriptArg1;
03363     // create and save outgoing message
03364     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
03365     exchange_msg->lattice = lattice;
03366     exchange_msg->pid = patchID;
03367     exchange_msg->numAtoms = numAtoms;
03368     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03369     if ( exchange_req >= 0 ) {
03370       recvExchangeReq(exchange_req);
03371     }
03372   }
03373   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
03374     exchange_src = (int) simParams->scriptArg2;
03375     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
03376   }
03377 }
03378 
03379 void HomePatch::recvExchangeReq(int req) {
03380   exchange_req = req;
03381   if ( exchange_msg ) {
03382     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
03383     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
03384     exchange_msg = 0;
03385     exchange_req = -1;
03386     CkpvAccess(_qd)->process();
03387   }
03388 }
03389 
03390 void HomePatch::recvExchangeMsg(ExchangeAtomsMsg *msg) {
03391   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
03392   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03393   lattice = msg->lattice;
03394   numAtoms = msg->numAtoms;
03395   atom.resize(numAtoms);
03396   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03397   delete msg;
03398   CkpvAccess(_qd)->process();
03399   doAtomUpdate = true;
03400   rattleListValid = false;
03401   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03402 }
03403 
03404 void HomePatch::submitLoadStats(int timestep)
03405 {
03406   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
03407 }
03408 
03409 
03410 void HomePatch::doPairlistCheck()
03411 {
03412   SimParameters *simParams = Node::Object()->simParameters;
03413 
03414   if ( numAtoms == 0 || ! flags.usePairlists ) {
03415     flags.pairlistTolerance = 0.;
03416     flags.maxAtomMovement = 99999.;
03417     return;
03418   }
03419 
03420   int i; int n = numAtoms;
03421   CompAtom *p_i = p.begin();
03422 
03423   if ( flags.savePairlists ) {
03424     flags.pairlistTolerance = doPairlistCheck_newTolerance;
03425     flags.maxAtomMovement = 0.;
03426     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
03427     doPairlistCheck_lattice = lattice;
03428     doPairlistCheck_positions.resize(numAtoms);
03429     CompAtom *psave_i = doPairlistCheck_positions.begin();
03430     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
03431     return;
03432   }
03433 
03434   Lattice &lattice_old = doPairlistCheck_lattice;
03435   Position center_cur = lattice.unscale(center);
03436   Position center_old = lattice_old.unscale(center);
03437   Vector center_delta = center_cur - center_old;
03438   
03439   // find max deviation to corner (any neighbor shares a corner)
03440   BigReal max_cd = 0.;
03441   for ( i=0; i<2; ++i ) {
03442     for ( int j=0; j<2; ++j ) {
03443       for ( int k=0; k<2; ++k ) {
03444         ScaledPosition corner(  i ? min.x : max.x ,
03445                                 j ? min.y : max.y ,
03446                                 k ? min.z : max.z );
03447         Vector corner_delta =
03448                 lattice.unscale(corner) - lattice_old.unscale(corner);
03449         corner_delta -= center_delta;
03450         BigReal cd = corner_delta.length2();
03451         if ( cd > max_cd ) max_cd = cd;
03452       }
03453     }
03454   }
03455   max_cd = sqrt(max_cd);
03456 
03457   // find max deviation of atoms relative to center
03458   BigReal max_pd = 0.;
03459   CompAtom *p_old_i = doPairlistCheck_positions.begin();
03460   for ( i=0; i<n; ++i ) {
03461     Vector p_delta = p_i[i].position - p_old_i[i].position;
03462     p_delta -= center_delta;
03463     BigReal pd = p_delta.length2();
03464     if ( pd > max_pd ) max_pd = pd;
03465   }
03466   max_pd = sqrt(max_pd);
03467 
03468   BigReal max_tol = max_pd + max_cd;
03469 
03470   flags.maxAtomMovement = max_tol;
03471 
03472   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
03473 
03474   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
03475                                 doPairlistCheck_newTolerance ) ) {
03476     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
03477   }
03478 
03479   if ( max_tol > doPairlistCheck_newTolerance ) {
03480     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
03481   }
03482 
03483 }
03484 
03485 void HomePatch::doGroupSizeCheck()
03486 {
03487   if ( ! flags.doNonbonded ) return;
03488 
03489   SimParameters *simParams = Node::Object()->simParameters;
03490   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
03491   BigReal maxrad2 = 0.;
03492 
03493   FullAtomList::iterator p_i = atom.begin();
03494   FullAtomList::iterator p_e = atom.end();
03495 
03496   while ( p_i != p_e ) {
03497     const int hgs = p_i->hydrogenGroupSize;
03498     if ( ! hgs ) break;  // avoid infinite loop on bug
03499     int ngs = hgs;
03500     if ( ngs > 5 ) ngs = 5;  // limit to at most 5 atoms per group
03501     BigReal x = p_i->position.x;
03502     BigReal y = p_i->position.y;
03503     BigReal z = p_i->position.z;
03504     int i;
03505     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
03506       p_i[i].nonbondedGroupSize = 0;
03507       BigReal dx = p_i[i].position.x - x;
03508       BigReal dy = p_i[i].position.y - y;
03509       BigReal dz = p_i[i].position.z - z;
03510       BigReal r2 = dx * dx + dy * dy + dz * dz;
03511       if ( r2 > hgcut ) break;
03512       else if ( r2 > maxrad2 ) maxrad2 = r2;
03513     }
03514     p_i->nonbondedGroupSize = i;
03515     for ( ; i < hgs; ++i ) {
03516       p_i[i].nonbondedGroupSize = 1;
03517     }
03518     p_i += hgs;
03519   }
03520 
03521   if ( p_i != p_e ) {
03522     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
03523   }
03524 
03525   flags.maxGroupRadius = sqrt(maxrad2);
03526 
03527 }
03528 
03529 void HomePatch::doMarginCheck()
03530 {
03531   SimParameters *simParams = Node::Object()->simParameters;
03532 
03533   BigReal sysdima = lattice.a_r().unit() * lattice.a();
03534   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
03535   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
03536 
03537   BigReal minSize = simParams->patchDimension - simParams->margin;
03538 
03539   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
03540        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
03541        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
03542 
03543     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03544       "Possible solutions are to restart from a recent checkpoint,\n"
03545       "increase margin, or disable useFlexibleCell for liquid simulation.");
03546   }
03547 
03548   BigReal cutoff = simParams->cutoff;
03549 
03550   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
03551   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
03552   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
03553 
03554   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
03555     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03556       "There are probably many margin violations already reported.\n"
03557       "Possible solutions are to restart from a recent checkpoint,\n"
03558       "increase margin, or disable useFlexibleCell for liquid simulation.");
03559   }
03560 
03561   BigReal minx = min.x - margina;
03562   BigReal miny = min.y - marginb;
03563   BigReal minz = min.z - marginc;
03564   BigReal maxx = max.x + margina;
03565   BigReal maxy = max.y + marginb;
03566   BigReal maxz = max.z + marginc;
03567 
03568   int xdev, ydev, zdev;
03569   int problemCount = 0;
03570 
03571   FullAtomList::iterator p_i = atom.begin();
03572   FullAtomList::iterator p_e = atom.end();
03573   for ( ; p_i != p_e; ++p_i ) {
03574 
03575     ScaledPosition s = lattice.scale(p_i->position);
03576 
03577     // check if atom is within bounds
03578     if (s.x < minx) xdev = 0;
03579     else if (maxx <= s.x) xdev = 2; 
03580     else xdev = 1;
03581 
03582     if (s.y < miny) ydev = 0;
03583     else if (maxy <= s.y) ydev = 2; 
03584     else ydev = 1;
03585 
03586     if (s.z < minz) zdev = 0;
03587     else if (maxz <= s.z) zdev = 2; 
03588     else zdev = 1;
03589 
03590     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
03591         ++problemCount;
03592     }
03593 
03594   }
03595 
03596   marginViolations = problemCount;
03597   // if ( problemCount ) {
03598   //     iout << iERROR <<
03599   //       "Found " << problemCount << " margin violations!\n" << endi;
03600   // } 
03601 
03602 }
03603 
03604 
03605 void
03606 HomePatch::doAtomMigration()
03607 {
03608   int i;
03609 
03610   for (i=0; i<numNeighbors; i++) {
03611     realInfo[i].mList.resize(0);
03612   }
03613 
03614   // Purge the AtomMap
03615   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03616 
03617   // Determine atoms that need to migrate
03618 
03619   BigReal minx = min.x;
03620   BigReal miny = min.y;
03621   BigReal minz = min.z;
03622   BigReal maxx = max.x;
03623   BigReal maxy = max.y;
03624   BigReal maxz = max.z;
03625 
03626   int xdev, ydev, zdev;
03627   int delnum = 0;
03628 
03629   FullAtomList::iterator atom_i = atom.begin();
03630   FullAtomList::iterator atom_e = atom.end();
03631 
03632   // DMK - Atom Separation (water vs. non-water)
03633   #if NAMD_SeparateWaters != 0
03634     FullAtomList::iterator atom_first = atom_i;
03635     int numLostWaterAtoms = 0;
03636   #endif
03637 
03638   while ( atom_i != atom_e ) {
03639     if ( atom_i->migrationGroupSize ) {
03640       Position pos = atom_i->position;
03641       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
03642         int mgs = atom_i->migrationGroupSize;
03643         int c = 1;
03644         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
03645                                 j+=(atom_i+j)->hydrogenGroupSize ) {
03646           pos += (atom_i+j)->position;
03647           ++c;
03648         }
03649         pos *= 1./c;
03650         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
03651       }
03652 
03653       ScaledPosition s = lattice.scale(pos);
03654 
03655       // check if atom is within bounds
03656       if (s.x < minx) xdev = 0;
03657       else if (maxx <= s.x) xdev = 2;
03658       else xdev = 1;
03659 
03660       if (s.y < miny) ydev = 0;
03661       else if (maxy <= s.y) ydev = 2;
03662       else ydev = 1;
03663 
03664       if (s.z < minz) zdev = 0;
03665       else if (maxz <= s.z) zdev = 2;
03666       else zdev = 1;
03667 
03668     }
03669 
03670     if (mInfo[xdev][ydev][zdev]) { // process atom for migration
03671                                     // Don't migrate if destination is myself
03672 
03673       // See if we have a migration list already
03674       MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
03675       DebugM(3,"Migrating atom " << atom_i->id << " from patch "
03676                 << patchID << " with position " << atom_i->position << "\n");
03677       mCur.add(*atom_i);
03678 
03679       ++delnum;
03680 
03681 
03682       // DMK - Atom Separation (water vs. non-water)
03683       #if NAMD_SeparateWaters != 0
03684         // Check to see if this atom is part of a water molecule.  If
03685         //   so, numWaterAtoms needs to be adjusted to refect the lost of
03686         //   this atom.
03687         // NOTE: The atom separation code assumes that if the oxygen
03688         //   atom of the hydrogen group making up the water molecule is
03689         //   migrated to another HomePatch, the hydrogens will also
03690         //   move!!!
03691         int atomIndex = atom_i - atom_first;
03692         if (atomIndex < numWaterAtoms)
03693           numLostWaterAtoms++;
03694       #endif
03695 
03696 
03697     } else {
03698 
03699       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
03700 
03701     }
03702 
03703     ++atom_i;
03704   }
03705 
03706   // DMK - Atom Separation (water vs. non-water)
03707   #if NAMD_SeparateWaters != 0
03708     numWaterAtoms -= numLostWaterAtoms;
03709   #endif
03710 
03711 
03712   int delpos = numAtoms - delnum;
03713   DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
03714   atom.del(delpos,delnum);
03715 
03716   numAtoms = atom.size();
03717 
03718   PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
03719 
03720   // signal depositMigration() that we are inMigration mode
03721   inMigration = true;
03722 
03723   // Drain the migration message buffer
03724   for (i=0; i<numMlBuf; i++) {
03725      DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
03726      depositMigration(msgbuf[i]);
03727   }
03728   numMlBuf = 0;
03729      
03730   if (!allMigrationIn) {
03731     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
03732     migrationSuspended = true;
03733     sequencer->suspend();
03734     migrationSuspended = false;
03735   }
03736   allMigrationIn = false;
03737 
03738   inMigration = false;
03739   marginViolations = 0;
03740 }
03741 
03742 void 
03743 HomePatch::depositMigration(MigrateAtomsMsg *msg)
03744 {
03745 
03746   if (!inMigration) { // We have to buffer changes due to migration
03747                       // until our patch is in migration mode
03748     msgbuf[numMlBuf++] = msg;
03749     return;
03750   } 
03751 
03752 
03753   // DMK - Atom Separation (water vs. non-water)
03754   #if NAMD_SeparateWaters != 0
03755 
03756 
03757     // Merge the incoming list of atoms with the current list of
03758     //   atoms.  Note that mergeSeparatedAtomList() will apply any
03759     //   required transformations to the incoming atoms as it is
03760     //   separating them.
03761     mergeAtomList(msg->migrationList);
03762 
03763 
03764   #else
03765 
03766     // Merge the incoming list of atoms with the current list of
03767     // atoms.  Apply transformations to the incoming atoms as they are
03768     // added to this patch's list.
03769     {
03770       MigrationList &migrationList = msg->migrationList;
03771       MigrationList::iterator mi;
03772       Transform mother_transform;
03773       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
03774         DebugM(1,"Migrating atom " << mi->id << " to patch "
03775                   << patchID << " with position " << mi->position << "\n"); 
03776         if ( mi->migrationGroupSize ) {
03777           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
03778             Position pos = mi->position;
03779             int mgs = mi->migrationGroupSize;
03780             int c = 1;
03781             for ( int j=mi->hydrogenGroupSize; j<mgs;
03782                                 j+=(mi+j)->hydrogenGroupSize ) {
03783               pos += (mi+j)->position;
03784               ++c;
03785             }
03786             pos *= 1./c;
03787             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
03788             mother_transform = mi->transform;
03789             pos = lattice.nearest(pos,center,&mother_transform);
03790             mi->position = lattice.reverse_transform(mi->position,mi->transform);
03791             mi->position = lattice.apply_transform(mi->position,mother_transform);
03792             mi->transform = mother_transform;
03793           } else {
03794             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
03795             mother_transform = mi->transform;
03796           }
03797         } else {
03798           mi->position = lattice.reverse_transform(mi->position,mi->transform);
03799           mi->position = lattice.apply_transform(mi->position,mother_transform);
03800           mi->transform = mother_transform;
03801         }
03802         atom.add(*mi);
03803       }
03804     }
03805 
03806 
03807   #endif // if (NAMD_SeparateWaters != 0)
03808 
03809 
03810   numAtoms = atom.size();
03811   delete msg;
03812 
03813   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
03814   if (!--patchMigrationCounter) {
03815     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
03816     allMigrationIn = true;
03817     patchMigrationCounter = numNeighbors;
03818     if (migrationSuspended) {
03819       DebugM(3,"patch "<<patchID<<" is being awakened\n");
03820       migrationSuspended = false;
03821       sequencer->awaken();
03822       return;
03823     }
03824     else {
03825        DebugM(3,"patch "<<patchID<<" was not suspended\n");
03826     }
03827   }
03828 }
03829 
03830 
03831 
03832 // DMK - Atom Separation (water vs. non-water)
03833 #if NAMD_SeparateWaters != 0
03834 
03835 // This function will separate waters from non-waters in the current
03836 //   atom list (regardless of whether or not the atom list is has been
03837 //   sorted yet or not).
03838 void HomePatch::separateAtoms() {
03839 
03840   // Basic Idea:  Iterate through all the atoms in the current list
03841   //   of atoms.  Pack the waters in the current atoms list and move
03842   //   the non-waters to the scratch list.  Once the atoms have all
03843   //   been separated, copy the non-waters to the end of the waters.
03844   // NOTE:  This code does not assume that the atoms list has been
03845   //   separated in any manner.
03846 
03847   // NOTE: Sanity check - Doesn't look like the default constructor actually
03848   //   adds any atoms but does set numAtoms. ???
03849   if (atom.size() < 0) return;  // Nothing to do.
03850 
03851   // Resize the scratch FullAtomList (tempAtom)
03852   tempAtom.resize(numAtoms);  // NOTE: Worst case: all non-water
03853 
03854   // Define size of a water hydrogen group
03855   int wathgsize = 3;
03856   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
03857   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
03858 
03859   // Iterate through all the atoms
03860   int i = 0;
03861   int waterIndex = 0;
03862   int nonWaterIndex = 0;
03863   while (i < numAtoms) {
03864 
03865     FullAtom &atom_i = atom[i];
03866     Mass mass = atom_i.mass;
03867     int hgs = atom_i.hydrogenGroupSize; 
03868     // Check to see if this hydrogen group is a water molecule
03869     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
03870 
03871       // Move this hydrogen group up in the current atom list
03872       if (waterIndex != i) {
03873         atom[waterIndex    ] = atom[i    ];  // Oxygen
03874         atom[waterIndex + 1] = atom[i + 1];  // Hydrogen
03875         atom[waterIndex + 2] = atom[i + 2];  // Hydrogen
03876         if (wathgsize > 3) atom[waterIndex + 3] = atom[i + 3];  // lonepair
03877         if (wathgsize > 4) atom[waterIndex + 4] = atom[i + 4];  // drude
03878           // actual Drude water is arranged:  O D LP H H
03879       }
03880 
03881       waterIndex += wathgsize;
03882       i += wathgsize;
03883 
03884     } else {
03885 
03886       // Move this hydrogen group into non-water (scratch) atom list
03887       for (int j = 0; j < hgs; j++)
03888         tempAtom[nonWaterIndex + j] = atom[i + j];
03889 
03890       nonWaterIndex += hgs;
03891       i += hgs;
03892     }
03893 
03894   } // end iterating through atoms
03895 
03896   // Iterate through the non-water (scratch) atom list, adding the
03897   //   atoms to the end of the atom list.
03898   // NOTE: This could be done with a straight memcpy if the internal
03899   //   data structures of ResizeArray could be accessed directly.
03900   //   Or, perhaps add a member to ResizeArray that can add a consecutive
03901   //   list of elements starting at a particular index (would be cleaner).
03902   for (i = 0; i < nonWaterIndex; i++)
03903     atom[waterIndex + i] = tempAtom[i];
03904 
03905   // Set numWaterAtoms
03906   numWaterAtoms = waterIndex;
03907 }
03908 
03909 
03910 // This function will merge the given list of atoms (not assumed to
03911 //   be separated) with the current list of atoms (already assumed
03912 //   to be separated).
03913 // NOTE: This function applies the transformations to the incoming
03914 //   atoms as it is separating them.
03915 void HomePatch::mergeAtomList(FullAtomList &al) {
03916 
03917   // Sanity check
03918   if (al.size() <= 0) return;  // Nothing to do
03919 
03920   const int orig_atomSize = atom.size();
03921   const int orig_alSize = al.size();
03922 
03923   // Resize the atom list (will eventually hold contents of both lists)
03924   atom.resize(orig_atomSize + orig_alSize); // NOTE: Will have contents of both
03925 
03926 
03927   #if 0  // version where non-waters are moved to scratch first
03928 
03929   
03930   // Basic Idea:  The current list is separated already so copy the
03931   //   non-water atoms out of it into the scratch atom array.  Then
03932   //   separate the incoming/given list (al), adding the waters to the
03933   //   end of the waters in atom list and non-waters to the end of the
03934   //   scratch list.  At this point, all waters are in atom list and all
03935   //   non-waters are in the scratch list so just copy the scratch list
03936   //   to the end of the atom list.
03937   // NOTE: If al is already separated and the number of waters in it
03938   //   is know, could simply move the non-waters in atoms back by that
03939   //   amount and directly copy the waters in al into the created gap
03940   //   and the non-waters in al to the end.  Leave this as an
03941   //   optimization for later since I'm not sure if this would actually
03942   //   do better as the combining code (for combining migration
03943   //   messages) would also have to merge the contents of the atom lists
03944   //   they carry.  Generally speaking, there is probably a faster way
03945   //   to do this, but this will get it working.
03946 
03947   // Copy all the non-waters in the current atom list into the
03948   //   scratch atom list.
03949   const int orig_atom_numNonWaters = orig_atomSize - numWaterAtoms;
03950   tempAtom.resize(orig_atom_numNonWaters + al.size()); // NOTE: Worst case
03951   for (int i = 0; i < orig_atom_numNonWaters; i++)
03952     tempAtom[i] = atom[numWaterAtoms + i];
03953 
03954   // Separate the contents of the given atom list (applying the
03955   // transforms as needed)
03956   int atom_waterIndex = numWaterAtoms;
03957   int atom_nonWaterIndex = orig_atom_numNonWaters;
03958   int i = 0;
03959   while (i < orig_alSize) {
03960 
03961     FullAtom &atom_i = al[i];
03962     int hgs = atom_i.hydrogenGroupSize;
03963     if ( hgs != atom_i.migrationGroupSize ) {
03964       NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
03965     }
03966     Mass mass = atom_i.mass;
03967 
03968     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
03969 
03970       // Apply the transforms
03971 
03972       // Oxygen (@ +0)
03973       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
03974       Transform mother_transform = al[i].transform;
03975 
03976       // Hydrogen (@ +1)
03977       al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
03978       al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
03979       al[i+1].transform = mother_transform;
03980 
03981       // Hydrogen (@ +2)
03982       al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
03983       al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
03984       al[i+2].transform = mother_transform;
03985 
03986       // Add to the end of the waters in the current list of atoms
03987       atom[atom_waterIndex    ] = al[i    ];
03988       atom[atom_waterIndex + 1] = al[i + 1];
03989       atom[atom_waterIndex + 2] = al[i + 2];
03990 
03991       atom_waterIndex += 3;
03992       i += 3;
03993 
03994     } else {
03995 
03996       // Apply the transforms
03997 
03998       // Non-Hydrogen (@ +0)
03999       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
04000       Transform mother_transform = al[i].transform;
04001 
04002       // Hydrogens (@ +1 -> +(hgs-1))
04003       for (int j = 1; j < hgs; j++) {
04004         al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
04005         al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
04006         al[i+j].transform = mother_transform;
04007       }
04008 
04009       // Add to the end of the non-waters (scratch) atom list
04010       for (int j = 0; j < hgs; j++)
04011         tempAtom[atom_nonWaterIndex + j] = al[i + j];
04012 
04013       atom_nonWaterIndex += hgs;
04014       i += hgs;
04015     }
04016 
04017   } // end while iterating through given atom list
04018 
04019   // Copy all the non-waters to the end of the current atom list
04020   for (int i = 0; i < atom_nonWaterIndex; i++)
04021     atom[atom_waterIndex + i] = tempAtom[i];
04022 
04023   // Set numWaterAtoms and numAtoms
04024   numWaterAtoms = atom_waterIndex;
04025   numAtoms = atom.size();
04026 
04027 
04028   #else
04029 
04030 
04031   // Basic Idea:  Count the number of water atoms in the incoming atom
04032   //   list then move the non-waters back in the current atom list to
04033   //   make room for the incoming waters.  Once there is room in the
04034   //   current list, separate the incoming list as the atoms are being
04035   //   added to the current list.
04036   // NOTE:  Since the incoming atom list is likely to be small,
04037   //   iterating over its hydrogen groups twice should not be too bad.
04038   // NOTE:  This code assumes the current list is already separated,
04039   //   the incoming list may not be separated, and the transforms are
04040   //   applied to the incoming atoms as the separation occurs.
04041 
04042   // size of a water hydrogen group
04043   int wathgsize = 3;
04044   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
04045   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
04046 
04047   // Count the number of waters in the given atom list
04048   int al_numWaterAtoms = 0;
04049   int i = 0;
04050   while (i < orig_alSize) {
04051 
04052     FullAtom &atom_i = al[i];
04053     int hgs = atom_i.hydrogenGroupSize;
04054     Mass mass = atom_i.mass;
04055 
04056     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
04057       al_numWaterAtoms += wathgsize;
04058     }
04059 
04060     i += hgs;
04061   }
04062 
04063   // Move all of the non-waters in the current atom list back (to a
04064   //   higher index) by the number of waters in the given list.
04065   if (al_numWaterAtoms > 0) {
04066     for (i = orig_atomSize - 1; i >= numWaterAtoms; i--) {
04067       atom[i + al_numWaterAtoms] = atom[i];
04068     }
04069   }
04070 
04071   // Separte the atoms in the given atom list.  Perform the
04072   //   transformations on them and then add them to the appropriate
04073   //   location in the current atom list.
04074   int atom_waterIndex = numWaterAtoms;
04075   int atom_nonWaterIndex = orig_atomSize + al_numWaterAtoms;
04076   i = 0;
04077   while (i < orig_alSize) {
04078 
04079     FullAtom &atom_i = al[i];
04080     int hgs = atom_i.hydrogenGroupSize;
04081     if ( hgs != atom_i.migrationGroupSize ) {
04082       NAMD_bug("HomePatch::mergeAtomList() not updated for migration groups!");
04083     }
04084     Mass mass = atom_i.mass;
04085 
04086     if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
04087 
04088       // Apply the transforms
04089 
04090       // Oxygen (@ +0)
04091       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
04092       Transform mother_transform = al[i].transform;
04093 
04094       // Hydrogen (@ +1)
04095       al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
04096       al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
04097       al[i+1].transform = mother_transform;
04098 
04099       // Hydrogen (@ +2)
04100       al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
04101       al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
04102       al[i+2].transform = mother_transform;
04103 
04104       // Add to the end of the waters in the current list of atoms
04105       atom[atom_waterIndex    ] = al[i    ];
04106       atom[atom_waterIndex + 1] = al[i + 1];
04107       atom[atom_waterIndex + 2] = al[i + 2];
04108 
04109       if (wathgsize > 3) atom[atom_waterIndex + 3] = al[i + 3];
04110 
04111       atom_waterIndex += wathgsize;
04112       i += wathgsize;
04113 
04114     } else {
04115 
04116       // Apply the transforms
04117 
04118       // Non-Hydrogen (@ +0)
04119       al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
04120       Transform mother_transform = al[i].transform;
04121 
04122       // Hydrogens (@ +1 -> +(hgs-1))
04123       for (int j = 1; j < hgs; j++) {
04124         al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
04125         al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
04126         al[i+j].transform = mother_transform;
04127       }
04128 
04129       // Add to the end of the non-waters (scratch) atom list
04130       for (int j = 0; j < hgs; j++)
04131         atom[atom_nonWaterIndex + j] = al[i + j];
04132 
04133       atom_nonWaterIndex += hgs;
04134       i += hgs;
04135     }
04136 
04137   } // end while iterating through given atom list
04138 
04139   // Set numWaterAtoms and numAtoms
04140   numWaterAtoms = atom_waterIndex;
04141   numAtoms = atom_nonWaterIndex;
04142 
04143   #endif
04144 }
04145 
04146 #endif
04147 
04148 
04149 
04150 inline void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx,
04151                                               HGArrayBigReal &b)
04152 {
04153         int i,ii=-1,ip,j;
04154         double sum;
04155 
04156         for (i=0;i<n;i++) {
04157                 ip=indx[i];
04158                 sum=b[ip];
04159                 b[ip]=b[i];
04160                 if (ii >= 0)
04161                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
04162                 else if (sum) ii=i;
04163                 b[i]=sum;
04164         }
04165         for (i=n-1;i>=0;i--) {
04166                 sum=b[i];
04167                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
04168                 b[i]=sum/a[i][i];
04169         }
04170 }
04171 
04172 
04173 inline void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
04174 {
04175 
04176         int i,imax,j,k;
04177         double big,dum,sum,temp;
04178         HGArrayBigReal vv;
04179         *d=1.0;
04180         for (i=0;i<n;i++) {
04181                 big=0.0;
04182                 for (j=0;j<n;j++)
04183                         if ((temp=fabs(a[i][j])) > big) big=temp;
04184                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
04185                 vv[i]=1.0/big;
04186         }
04187         for (j=0;j<n;j++) {
04188                 for (i=0;i<j;i++) {
04189                         sum=a[i][j];
04190                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
04191                         a[i][j]=sum;
04192                 }
04193                 big=0.0;
04194                 for (i=j;i<n;i++) {
04195                         sum=a[i][j];
04196                         for (k=0;k<j;k++)
04197                                 sum -= a[i][k]*a[k][j];
04198                         a[i][j]=sum;
04199                         if ( (dum=vv[i]*fabs(sum)) >= big) {
04200                                 big=dum;
04201                                 imax=i;
04202                         }
04203                 }
04204                 if (j != imax) {
04205                         for (k=0;k<n;k++) {
04206                                 dum=a[imax][k];
04207                                 a[imax][k]=a[j][k];
04208                                 a[j][k]=dum;
04209                         }
04210                         *d = -(*d);
04211                         vv[imax]=vv[j];
04212                 }
04213                 indx[j]=imax;
04214                 if (a[j][j] == 0.0) a[j][j]=TINY;
04215                 if (j != n-1) {
04216                         dum=1.0/(a[j][j]);
04217                         for (i=j+1;i<n;i++) a[i][j] *= dum;
04218                 }
04219         }
04220 }
04221 
04222 
04223 inline void G_q(const HGArrayVector &refab, HGMatrixVector &gqij,
04224      const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl) {
04225   int i; 
04226   // step through the rows of the matrix
04227   for(i=0;i<m;i++) {
04228     gqij[i][ial[i]]=2.0*refab[i];
04229     gqij[i][ibl[i]]=-gqij[i][ial[i]];
04230   }
04231 };
04232 
04233 
04234 // c-ji code for MOLLY 7-31-99
04235 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) {
04236   //  input:  n = length of hydrogen group to be averaged (shaked)
04237   //          q[n] = vector array of original positions
04238   //          m = number of constraints
04239   //          imass[n] = inverse mass for each atom
04240   //          length2[m] = square of reference bond length for each constraint
04241   //          ial[m] = atom a in each constraint 
04242   //          ibl[m] = atom b in each constraint 
04243   //          refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
04244   //          tolf = function error tolerance for Newton's iteration
04245   //          ntrial = max number of Newton's iterations
04246   //  output: lambda[m] = double array of lagrange multipliers (used by mollify)
04247   //          qtilde[n] = vector array of averaged (shaked) positions
04248 
04249   int k,k1,i,j;
04250   BigReal errx,errf,d,tolx;
04251 
04252   HGArrayInt indx;
04253   HGArrayBigReal p;
04254   HGArrayBigReal fvec;
04255   HGMatrixBigReal fjac;
04256   HGArrayVector avgab;
04257   HGMatrixVector grhs;
04258   HGMatrixVector auxrhs;
04259   HGMatrixVector glhs;
04260 
04261   //  iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
04262   tolx=tolf; 
04263   
04264   // initialize lambda, globalGrhs
04265 
04266   for (i=0;i<m;i++) {
04267     lambda[i]=0.0;
04268   }
04269 
04270   // define grhs, auxrhs for all iterations
04271   // grhs= g_x(q)
04272   //
04273   G_q(refab,grhs,n,m,ial,ibl);
04274   for (k=1;k<=ntrial;k++) {
04275     //    usrfun(qtilde,q0,lambda,fvec,fjac,n,water); 
04276     HGArrayBigReal gij;
04277     // this used to be main loop of usrfun
04278     // compute qtilde given q0, lambda, IMASSes
04279     {
04280       BigReal multiplier;
04281       HGArrayVector tmp;
04282       for (i=0;i<m;i++) {
04283         multiplier = lambda[i];
04284         // auxrhs = M^{-1}grhs^{T}
04285         for (j=0;j<n;j++) {
04286           auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
04287         }
04288       }
04289       for (j=0;j<n;j++) {
04290         //      tmp[j]=0.0;      
04291         for (i=0;i<m;i++) {
04292           tmp[j]+=auxrhs[i][j];
04293         }
04294       }
04295  
04296       for (j=0;j<n;j++) {
04297         qtilde[j].position=q[j]+tmp[j];
04298       }
04299       //      delete [] tmp;
04300     }
04301   
04302     for ( i = 0; i < m; i++ ) {
04303       avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04304     }
04305 
04306     //  iout<<iINFO << "Calling Jac" << std::endl<<endi;
04307     //  Jac(qtilde, q0, fjac,n,water);
04308     {
04309       //  Vector glhs[3*n+3];
04310 
04311       HGMatrixVector grhs2;
04312 
04313       G_q(avgab,glhs,n,m,ial,ibl);
04314 #ifdef DEBUG0
04315       iout<<iINFO << "G_q:" << std::endl<<endi;
04316       for (i=0;i<m;i++) {
04317         iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
04318       }
04319 #endif
04320       //      G_q(refab,grhs2,m,ial,ibl);
04321       // update with the masses
04322       for (j=0; j<n; j++) { // number of atoms
04323         for (i=0; i<m;i++) { // number of constraints
04324           grhs2[i][j] = grhs[i][j]*imass[j];
04325         }
04326       }
04327 
04328       // G_q(qtilde) * M^-1 G_q'(q0) =
04329       // G_q(qtilde) * grhs'
04330       for (i=0;i<m;i++) { // number of constraints
04331         for (j=0;j<m;j++) { // number of constraints
04332           fjac[i][j] = 0; 
04333           for (k1=0;k1<n;k1++) {
04334             fjac[i][j] += glhs[i][k1]*grhs2[j][k1]; 
04335           }
04336         }
04337       }
04338 #ifdef DEBUG0  
04339       iout<<iINFO << "glhs" <<endi;
04340       for(i=0;i<9;i++) {
04341         iout<<iINFO << glhs[i] << ","<<endi;
04342       }
04343       iout<<iINFO << std::endl<<endi;
04344       for(i=0;i<9;i++) {
04345         iout<<iINFO << grhs2[i] << ","<<endi;
04346       }
04347       iout<<iINFO << std::endl<<endi;
04348 #endif
04349       //      delete[] grhs2;
04350     }
04351     // end of Jac calculation
04352 #ifdef DEBUG0
04353     iout<<iINFO << "Jac" << std::endl<<endi;
04354     for (i=0;i<m;i++) 
04355       for (j=0;j<m;j++)
04356         iout<<iINFO << fjac[i][j] << " "<<endi;
04357     iout<< std::endl<<endi;
04358 #endif
04359     // calculate constraints in gij for n constraints this being a water
04360     //  G(qtilde, gij, n, water);
04361     for (i=0;i<m;i++) {
04362       gij[i]=avgab[i]*avgab[i]-length2[i];
04363     }
04364 #ifdef DEBUG0
04365     iout<<iINFO << "G" << std::endl<<endi;
04366     iout<<iINFO << "( "<<endi;
04367     for(i=0;i<m-1;i++) {
04368       iout<<iINFO << gij[i] << ", "<<endi;
04369     }
04370     iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
04371 #endif
04372     // fill the return vector
04373     for(i=0;i<m;i++) {
04374       fvec[i] = gij[i];
04375     }
04376     // free up the constraints
04377     //    delete[] gij;
04378     // continue Newton's iteration    
04379     errf=0.0;
04380     for (i=0;i<m;i++) errf += fabs(fvec[i]);
04381 #ifdef DEBUG0
04382     iout<<iINFO << "errf: " << errf << std::endl<<endi;
04383 #endif
04384     if (errf <= tolf) {
04385       break;
04386     }
04387     for (i=0;i<m;i++) p[i] = -fvec[i];
04388     //    iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
04389     ludcmp(fjac,m,indx,&d);
04390     lubksb(fjac,m,indx,p);
04391 
04392     errx=0.0;
04393     for (i=0;i<m;i++) {
04394       errx += fabs(p[i]);
04395     }
04396     for (i=0;i<m;i++)  
04397       lambda[i] += p[i];
04398 
04399 #ifdef DEBUG0
04400     iout<<iINFO << "lambda:" << lambda[0] 
04401          << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04402     iout<<iINFO << "errx: " << errx << std::endl<<endi;
04403 #endif
04404     if (errx <= tolx) break;
04405 #ifdef DEBUG0
04406     iout<<iINFO << "Qtilde:" << std::endl<<endi;
04407     iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi; 
04408 #endif
04409   }
04410 #ifdef DEBUG
04411   iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04412 #endif
04413 
04414   return k; // 
04415 }
04416 
04417 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) {
04418   int i,j,k;
04419   BigReal d;
04420   HGMatrixBigReal fjac;
04421   Vector zero(0.0,0.0,0.0);
04422   
04423   HGArrayVector tmpforce;
04424   HGArrayVector tmpforce2;
04425   HGArrayVector y;
04426   HGMatrixVector grhs;
04427   HGMatrixVector glhs;
04428   HGArrayBigReal aux;
04429   HGArrayInt indx;
04430 
04431   for(i=0;i<n;i++) {
04432     tmpforce[i]=imass[i]*force[i];
04433   }
04434 
04435   HGMatrixVector grhs2;
04436   HGArrayVector avgab;
04437 
04438   for ( i = 0; i < m; i++ ) {
04439         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04440   }
04441 
04442   G_q(avgab,glhs,n,m,ial,ibl);
04443   G_q(refab,grhs,n,m,ial,ibl);
04444   // update with the masses
04445   for (j=0; j<n; j++) { // number of atoms
04446         for (i=0; i<m;i++) { // number of constraints
04447           grhs2[i][j] = grhs[i][j]*imass[j];
04448         }
04449   }
04450 
04451   // G_q(qtilde) * M^-1 G_q'(q0) =
04452   // G_q(qtilde) * grhs'
04453   for (i=0;i<m;i++) { // number of constraints
04454         for (j=0;j<m;j++) { // number of constraints
04455           fjac[j][i] = 0; 
04456           for (k=0;k<n;k++) {
04457             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
04458           }
04459         }
04460   }
04461 
04462   // aux=gqij*tmpforce
04463   //  globalGrhs::computeGlobalGrhs(q0,n,water);
04464   //  G_q(refab,grhs,m,ial,ibl);
04465   for(i=0;i<m;i++) {
04466     aux[i]=0.0;
04467     for(j=0;j<n;j++) {
04468       aux[i]+=grhs[i][j]*tmpforce[j];
04469     }
04470   }
04471 
04472   ludcmp(fjac,m,indx,&d);
04473   lubksb(fjac,m,indx,aux);
04474 
04475   for(j=0;j<n;j++) {
04476     y[j] = zero;
04477     for(i=0;i<m;i++) {
04478       y[j] += aux[i]*glhs[i][j];
04479     }
04480   }
04481   for(i=0;i<n;i++) {
04482     y[i]=force[i]-y[i];
04483   }
04484     
04485   // gqq12*y
04486   for(i=0;i<n;i++) {
04487     tmpforce2[i]=imass[i]*y[i];
04488   }
04489 
04490   // here we assume that tmpforce is initialized to zero.
04491   for (i=0;i<n;i++) {
04492     tmpforce[i]=zero;
04493   }
04494   
04495   for (j=0;j<m;j++) {
04496     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
04497     tmpforce[ial[j]] += tmpf;
04498     tmpforce[ibl[j]] -= tmpf;
04499   }
04500   // c-ji the other bug for 2 constraint water was this line (2-4-99)
04501   //  for(i=0;i<m;i++) {
04502   for(i=0;i<n;i++) {
04503     force[i]=tmpforce[i]+y[i];
04504   }
04505 
04506 }

Generated on Wed Jun 20 01:17:13 2018 for NAMD by  doxygen 1.4.7