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

Generated on Tue Nov 21 01:17:13 2017 for NAMD by  doxygen 1.4.7