Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

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

Generated on Tue May 21 04:07:16 2013 for NAMD by  doxygen 1.3.9.1