HomePatch Class Reference

#include <HomePatch.h>

Inheritance diagram for HomePatch:

Patch List of all members.

Public Member Functions

 ~HomePatch ()
void registerProxy (RegisterProxyMsg *)
void unregisterProxy (UnregisterProxyMsg *)
void receiveResults (ProxyResultVarsizeMsg *msg)
void receiveResults (ProxyResultMsg *msg)
void receiveResult (ProxyGBISP1ResultMsg *msg)
void receiveResult (ProxyGBISP2ResultMsg *msg)
void receiveResults (ProxyCombinedResultRawMsg *msg)
void depositMigration (MigrateAtomsMsg *)
void useSequencer (Sequencer *sequencerPtr)
void runSequencer (void)
void positionsReady (int doMigration=0)
void saveForce (const int ftag=Results::normal)
void addForceToMomentum (const BigReal, const int ftag=Results::normal, const int useSaved=0)
void addForceToMomentum3 (const BigReal timestep1, const int ftag1, const int useSaved1, const BigReal timestep2, const int ftag2, const int useSaved2, const BigReal timestep3, const int ftag3, const int useSaved3)
void addVelocityToPosition (const BigReal)
int hardWallDrude (const BigReal, Tensor *virial, SubmitReduction *)
void addRattleForce (const BigReal invdt, Tensor &wc)
void buildRattleList ()
int rattle1old (const BigReal, Tensor *virial, SubmitReduction *)
int rattle1 (const BigReal, Tensor *virial, SubmitReduction *)
void rattle2 (const BigReal, Tensor *virial)
void minimize_rattle2 (const BigReal, Tensor *virial, bool forces=false)
void mollyAverage ()
void mollyMollify (Tensor *virial)
void loweAndersenVelocities ()
void loweAndersenFinish ()
void setGBISIntrinsicRadii ()
void gbisComputeAfterP1 ()
void gbisComputeAfterP2 ()
void gbisP2Ready ()
void gbisP3Ready ()
void setLcpoType ()
void checkpoint (void)
void revert (void)
void exchangeCheckpoint (int scriptTask, int &bpc)
void recvCheckpointReq (int task, const char *key, int replica, int pe)
void recvCheckpointLoad (CheckpointAtomsMsg *msg)
void recvCheckpointStore (CheckpointAtomsMsg *msg)
void recvCheckpointAck ()
void exchangeAtoms (int scriptTask)
void recvExchangeReq (int req)
void recvExchangeMsg (ExchangeAtomsMsg *msg)
void replaceForces (ExtForce *f)
void qmSwapAtoms ()
void submitLoadStats (int timestep)
FullAtomListgetAtomList ()
void buildSpanningTree (void)
void sendNodeAwareSpanningTree ()
void recvNodeAwareSpanningTree (ProxyNodeAwareSpanningTreeMsg *msg)
void sendSpanningTree ()
void recvSpanningTree (int *t, int n)
void sendProxies ()

Public Attributes

int marginViolations
std::vector< int > settleList
std::vector< RattleListrattleList
std::vector< RattleParamrattleParam
std::vector< int > noconstList
bool rattleListValid
std::vector< VectorvelNew
std::vector< VectorposNew
int checkpoint_task
std::map< std::string, checkpoint_t * > checkpoints
int exchange_dst
int exchange_src
int exchange_req
ExchangeAtomsMsgexchange_msg
LDObjHandle ldObjHandle

Protected Member Functions

virtual void boxClosed (int)
void doPairlistCheck ()
void doGroupSizeCheck ()
void doMarginCheck ()
void doAtomMigration ()

Protected Attributes

int inMigration
int numMlBuf
MigrateAtomsMsgmsgbuf [PatchMap::MaxOneAway]

Friends

class PatchMgr
class Sequencer
class ComputeGlobal

Classes

struct  checkpoint_t
struct  RattleList

Detailed Description

Definition at line 46 of file HomePatch.h.


Constructor & Destructor Documentation

HomePatch::~HomePatch (  ) 

Definition at line 411 of file HomePatch.C.

References Patch::atomMapper, ResizeArray< Elem >::begin(), ResizeArray< Elem >::end(), and AtomMapper::unregisterIDsFullAtom().

00412 {
00413     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
00414 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00415     ptnTree.resize(0);
00416     #ifdef USE_NODEPATCHMGR
00417     delete [] nodeChildren;    
00418     #endif
00419 #endif
00420   delete [] child;
00421 }


Member Function Documentation

void HomePatch::addForceToMomentum ( const   BigReal,
const int  ftag = Results::normal,
const int  useSaved = 0 
)

Definition at line 1822 of file HomePatch.C.

References ResizeArray< Elem >::begin(), ResizeArray< Elem >::const_begin(), Patch::f, namd_reciprocal, Node::Object(), Node::simParameters, simParams, TIMEFACTOR, FullAtom::velocity, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::addForceToMomentum().

01824 {
01825   SimParameters *simParams = Node::Object()->simParameters;
01826   const BigReal dt = timestep / TIMEFACTOR;
01827   ForceList *f_use = (useSaved ? f_saved : f);
01828 
01829   if ( simParams->fixedAtomsOn ) {
01830     for ( int i = 0; i < numAtoms; ++i ) {
01831       if ( atom[i].atomFixed ) {
01832         atom[i].velocity = 0;
01833       } else {
01834         BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.); 
01835         atom[i].velocity += f_use[ftag][i] * recip_val;
01836       }
01837     }
01838   } else {
01839     FullAtom *atom_arr  = atom.begin();
01840     const Force    *force_arr = f_use[ftag].const_begin();
01841 #ifdef ARCH_POWERPC
01842 #pragma disjoint (*force_arr, *atom_arr)
01843 #endif
01844     for ( int i = 0; i < numAtoms; ++i ) {
01845       if (atom[i].mass == 0.) continue;
01846       BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.); 
01847       //printf("Taking reciprocal of mass %f\n", atom[i].mass);
01848       atom_arr[i].velocity.x += force_arr[i].x * recip_val;
01849       atom_arr[i].velocity.y += force_arr[i].y * recip_val;
01850       atom_arr[i].velocity.z += force_arr[i].z * recip_val;
01851     }
01852   }
01853 }

void HomePatch::addForceToMomentum3 ( const BigReal  timestep1,
const int  ftag1,
const int  useSaved1,
const BigReal  timestep2,
const int  ftag2,
const int  useSaved2,
const BigReal  timestep3,
const int  ftag3,
const int  useSaved3 
)

Definition at line 1855 of file HomePatch.C.

References ResizeArray< Elem >::begin(), ResizeArray< Elem >::const_begin(), Patch::f, namd_reciprocal, Node::Object(), Node::simParameters, simParams, TIMEFACTOR, FullAtom::velocity, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::addForceToMomentum3().

01858 {
01859   SimParameters *simParams = Node::Object()->simParameters;
01860   const BigReal dt1 = timestep1 / TIMEFACTOR;
01861   const BigReal dt2 = timestep2 / TIMEFACTOR;
01862   const BigReal dt3 = timestep3 / TIMEFACTOR;
01863   ForceList *f_use1 = (useSaved1 ? f_saved : f);
01864   ForceList *f_use2 = (useSaved2 ? f_saved : f);
01865   ForceList *f_use3 = (useSaved3 ? f_saved : f);
01866 
01867   if ( simParams->fixedAtomsOn ) {
01868     for ( int i = 0; i < numAtoms; ++i ) {
01869       if ( atom[i].atomFixed ) {
01870         atom[i].velocity = 0;
01871       } else {
01872         BigReal recip_val = ( atom[i].mass > 0. ? namd_reciprocal( atom[i].mass ) : 0.); 
01873         atom[i].velocity += (f_use1[ftag1][i]*dt1 + f_use2[ftag2][i]*dt2 + f_use3[ftag3][i]*dt3)*recip_val;
01874       }
01875     }
01876   } else {
01877     FullAtom *atom_arr  = atom.begin();
01878     const Force *force_arr1 = f_use1[ftag1].const_begin();
01879     const Force *force_arr2 = f_use2[ftag2].const_begin();
01880     const Force *force_arr3 = f_use3[ftag3].const_begin();
01881 #ifdef ARCH_POWERPC
01882 #pragma disjoint (*force_arr1, *atom_arr)
01883 #pragma disjoint (*force_arr2, *atom_arr)
01884 #pragma disjoint (*force_arr3, *atom_arr)
01885 #endif
01886     for ( int i = 0; i < numAtoms; ++i ) {
01887       if (atom[i].mass == 0.) continue;
01888       BigReal recip_val = ( atom[i].mass > 0. ? namd_reciprocal( atom[i].mass ) : 0.); 
01889       //printf("Taking reciprocal of mass %f\n", atom[i].mass);
01890       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * recip_val;
01891       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * recip_val;
01892       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * recip_val;
01893     }
01894   }
01895 }

void HomePatch::addRattleForce ( const BigReal  invdt,
Tensor wc 
)

Definition at line 2263 of file HomePatch.C.

References Patch::f, Results::normal, outer(), and velNew.

Referenced by rattle1().

02263                                                               {
02264   for (int ig = 0; ig < numAtoms; ++ig ) {
02265     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
02266     Tensor vir = outer(df, atom[ig].position);
02267     wc += vir;
02268     f[Results::normal][ig] += df;
02269     atom[ig].velocity = velNew[ig];
02270   }
02271 }

void HomePatch::addVelocityToPosition ( const   BigReal  ) 

Definition at line 1897 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Node::Object(), CompAtom::position, Node::simParameters, simParams, TIMEFACTOR, FullAtom::velocity, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::addVelocityToPosition().

01898 {
01899   SimParameters *simParams = Node::Object()->simParameters;
01900   const BigReal dt = timestep / TIMEFACTOR;
01901   if ( simParams->fixedAtomsOn ) {
01902     for ( int i = 0; i < numAtoms; ++i ) {
01903       if ( ! atom[i].atomFixed ) atom[i].position += atom[i].velocity * dt;
01904     }
01905   } else {
01906     FullAtom *atom_arr  = atom.begin();
01907     for ( int i = 0; i < numAtoms; ++i ) {
01908       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01909       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01910       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01911     }
01912   }
01913 }

void HomePatch::boxClosed ( int   )  [protected, virtual]

Implements Patch.

Definition at line 424 of file HomePatch.C.

References Sequencer::awaken(), Patch::boxesOpen, DebugM, Patch::dEdaSumBox, Flags::doGBIS, Flags::doNonbonded, Patch::f, Patch::flags, j, Results::maxNumForces, Results::normal, Patch::patchID, Patch::psiSumBox, ResizeArray< Elem >::size(), and Sequencer::thread.

00424                                  {
00425   // begin gbis
00426   if (box == 5) {// end of phase 1
00427     phase1BoxClosedCalled = true;
00428     if (!psiSumBox.isOpen() && numGBISP1Arrived == proxy.size()) {
00429       if (flags.doGBIS && flags.doNonbonded) {
00430         sequencer->awaken();
00431       }
00432     } else {
00433       //need to wait until proxies arrive before awakening
00434     }
00435   } else if (box == 6) {// intRad
00436     //do nothing
00437   } else if (box == 7) {// bornRad
00438     //do nothing
00439   } else if (box == 8) {// end of phase 2
00440     phase2BoxClosedCalled = true;
00441     //if no proxies, AfterP1 can't be called from receive
00442     //so it will be called from here
00443     if (!dEdaSumBox.isOpen() && numGBISP2Arrived == proxy.size()) {
00444       if (flags.doGBIS && flags.doNonbonded) {
00445         sequencer->awaken();
00446       }
00447     } else {
00448       //need to wait until proxies arrive before awakening
00449     }
00450   } else if (box == 9) {
00451     //do nothing
00452   } else if (box == 10) {
00453     //lcpoType Box closed: do nothing
00454   } else {
00455     //do nothing
00456   }
00457   // end gbis
00458 
00459   if ( ! --boxesOpen ) {
00460     if ( replacementForces ) {
00461       for ( int i = 0; i < numAtoms; ++i ) {
00462         if ( replacementForces[i].replace ) {
00463           for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
00464           f[Results::normal][i] = replacementForces[i].force;
00465         }
00466       }
00467       replacementForces = 0;
00468     }
00469     DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
00470         << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
00471     // only awaken suspended threads.  Then say it is suspended.
00472 
00473     phase3BoxClosedCalled = true;
00474     if (flags.doGBIS) {
00475       if (flags.doNonbonded) {
00476         sequencer->awaken();
00477       } else {
00478         if (numGBISP1Arrived == proxy.size() &&
00479           numGBISP2Arrived == proxy.size() &&
00480           numGBISP3Arrived == proxy.size()) {
00481           sequencer->awaken();//all boxes closed and all proxies arrived
00482         }
00483       }
00484     } else {//non-gbis awaken
00485       sequencer->awaken();
00486     }
00487   } else {
00488     DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
00489   }
00490 }

void HomePatch::buildRattleList (  ) 

Definition at line 2121 of file HomePatch.C.

References Node::Object(), posNew, settle1init(), Node::simParameters, simParams, velNew, WAT_SWM4, and WAT_TIP4.

Referenced by rattle1().

02121                                 {
02122 
02123   SimParameters *simParams = Node::Object()->simParameters;
02124   const int fixedAtomsOn = simParams->fixedAtomsOn;
02125   const int useSettle = simParams->useSettle;
02126 
02127   // Re-size to containt numAtoms elements
02128   velNew.resize(numAtoms);
02129   posNew.resize(numAtoms);
02130 
02131   // Size of a hydrogen group for water
02132   int wathgsize = 3;
02133   int watmodel = simParams->watmodel;
02134   if (watmodel == WAT_TIP4) wathgsize = 4;
02135   else if (watmodel == WAT_SWM4) wathgsize = 5;
02136 
02137   // Initialize the settle algorithm with water parameters
02138   // settle1() assumes all waters are identical,
02139   // and will generate bad results if they are not.
02140   // XXX this will move to Molecule::build_atom_status when that 
02141   // version is debugged
02142   if ( ! settle_initialized ) {
02143     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02144       // find a water
02145       if (atom[ig].rigidBondLength > 0) {
02146         int oatm;
02147         if (simParams->watmodel == WAT_SWM4) {
02148           oatm = ig+3;  // skip over Drude and Lonepair
02149           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02150           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02151         }
02152         else {
02153           oatm = ig+1;
02154           // Avoid using the Om site to set this by mistake
02155           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02156             oatm += 1;
02157           }
02158         }
02159 
02160         // initialize settle water parameters
02161         settle1init(atom[ig].mass, atom[oatm].mass, 
02162                     atom[ig].rigidBondLength,
02163                     atom[oatm].rigidBondLength,
02164                     settle_mOrmT, settle_mHrmT, settle_ra,
02165                     settle_rb, settle_rc, settle_rra);
02166         settle_initialized = 1;
02167         break; // done with init
02168       }
02169     }
02170   }
02171   
02172   Vector ref[10];
02173   BigReal rmass[10];
02174   BigReal dsq[10];
02175   int fixed[10];
02176   int ial[10];
02177   int ibl[10];
02178 
02179   int numSettle = 0;
02180   int numRattle = 0;
02181   int posRattleParam = 0;
02182 
02183   settleList.clear();
02184   rattleList.clear();
02185   noconstList.clear();
02186   rattleParam.clear();
02187 
02188   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02189     int hgs = atom[ig].hydrogenGroupSize;
02190     if ( hgs == 1 ) {
02191       // only one atom in group
02192       noconstList.push_back(ig);
02193       continue;
02194     }
02195     int anyfixed = 0;
02196     for (int i = 0; i < hgs; ++i ) {
02197       ref[i] = atom[ig+i].position;
02198       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02199       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02200       if ( fixed[i] ) {
02201         anyfixed = 1;
02202         rmass[i] = 0.;
02203       }
02204     }
02205     int icnt = 0;
02206     BigReal tmp = atom[ig].rigidBondLength;
02207     if (tmp > 0.0) {  // for water
02208       if (hgs != wathgsize) {
02209         char errmsg[256];
02210         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02211                          "but the specified water model requires %d atoms.\n",
02212                          atom[ig].id+1, hgs, wathgsize);
02213         NAMD_die(errmsg);
02214       }
02215       // Use SETTLE for water unless some of the water atoms are fixed,
02216       if (useSettle && !anyfixed) {
02217         // Store to Settle -list
02218         settleList.push_back(ig);
02219         continue;
02220       }
02221       if ( !(fixed[1] && fixed[2]) ) {
02222         dsq[icnt] = tmp * tmp;
02223         ial[icnt] = 1;
02224         ibl[icnt] = 2;
02225         ++icnt;
02226       }
02227     } // if (tmp > 0.0)
02228     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02229       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02230         if ( !(fixed[0] && fixed[i]) ) {
02231           dsq[icnt] = tmp * tmp;
02232           ial[icnt] = 0;
02233           ibl[icnt] = i;
02234           ++icnt;
02235         }
02236       }
02237     }
02238     if ( icnt == 0 ) {
02239       // no constraints
02240       noconstList.push_back(ig);
02241       continue;  
02242     }
02243     // Store to Rattle -list
02244     RattleList rattleListElem;
02245     rattleListElem.ig  = ig;
02246     rattleListElem.icnt = icnt;
02247     rattleList.push_back(rattleListElem);
02248     for (int i = 0; i < icnt; ++i ) {
02249       int a = ial[i];
02250       int b = ibl[i];
02251       RattleParam rattleParamElem;
02252       rattleParamElem.ia = a;
02253       rattleParamElem.ib = b;
02254       rattleParamElem.dsq = dsq[i];
02255       rattleParamElem.rma = rmass[a];
02256       rattleParamElem.rmb = rmass[b];
02257       rattleParam.push_back(rattleParamElem);
02258     }
02259   }
02260 
02261 }

void HomePatch::buildSpanningTree ( void   ) 

Definition at line 764 of file HomePatch.C.

References ResizeArray< Elem >::begin(), compDistance(), ResizeArray< Elem >::end(), ResizeArray< Elem >::find(), j, NAMD_bug(), PatchMap::numNodesWithPatches(), PatchMap::Object(), Patch::patchID, proxySpanDim, ResizeArray< Elem >::resize(), sendSpanningTree(), ResizeArray< Elem >::setall(), ResizeArray< Elem >::size(), and ResizeArray< Elem >::swap().

00765 {
00766   nChild = 0;
00767   int psize = proxy.size();
00768   if (psize == 0) return;
00769   NodeIDList oldtree;  oldtree.swap(tree);
00770   int oldsize = oldtree.size();
00771   tree.resize(psize + 1);
00772   tree.setall(-1);
00773   tree[0] = CkMyPe();
00774   int s=1, e=psize+1;
00775   NodeIDList::iterator pli;
00776   int patchNodesLast =
00777     ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00778   int nNonPatch = 0;
00779 #if 1
00780     // try to put it to the same old tree
00781   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
00782   {
00783     int oldindex = oldtree.find(*pli);
00784     if (oldindex != -1 && oldindex < psize) {
00785       tree[oldindex] = *pli;
00786     }
00787   }
00788   s=1; e=psize;
00789   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
00790   {
00791     if (tree.find(*pli) != -1) continue;    // already assigned
00792     if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
00793       while (tree[e] != -1) { e--; if (e==-1) e = psize; }
00794       tree[e] = *pli;
00795     } else {
00796       while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
00797       tree[s] = *pli;
00798       nNonPatch++;
00799     }
00800   }
00801 #if 1
00802   if (oldsize==0 && nNonPatch) {
00803     // first time, sort by distance
00804     qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
00805   }
00806 #endif
00807 
00808   //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
00809 
00810 #if USE_TOPOMAP && USE_SPANNING_TREE
00811   
00812   //Right now only works for spanning trees with two levels
00813   int *treecopy = new int [psize];
00814   int subroots[proxySpanDim];
00815   int subsizes[proxySpanDim];
00816   int subtrees[proxySpanDim][proxySpanDim];
00817   int idxes[proxySpanDim];
00818   int i = 0;
00819 
00820   for(i=0;i<proxySpanDim;i++){
00821     subsizes[i] = 0;
00822     idxes[i] = i;
00823   }
00824   
00825   for(i=0;i<psize;i++){
00826     treecopy[i] = tree[i+1];
00827   }
00828   
00829   TopoManager tmgr;
00830   tmgr.sortRanksByHops(treecopy,nNonPatch);
00831   tmgr.sortRanksByHops(treecopy+nNonPatch,
00832                                                 psize-nNonPatch);  
00833   
00834   /* build tree and subtrees */
00835   nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
00836   delete [] treecopy;
00837   
00838   for(int i=1;i<psize+1;i++){
00839     int isSubroot=0;
00840     for(int j=0;j<nChild;j++)
00841       if(tree[i]==subroots[j]){
00842         isSubroot=1;
00843         break;
00844       }
00845     if(isSubroot) continue;
00846     
00847     int bAdded = 0;
00848     tmgr.sortIndexByHops(tree[i], subroots,
00849                                                   idxes, proxySpanDim);
00850     for(int j=0;j<proxySpanDim;j++){
00851       if(subsizes[idxes[j]]<proxySpanDim){
00852         subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
00853         bAdded = 1; 
00854         break;
00855       }
00856     }
00857     if( psize > proxySpanDim && ! bAdded ) {
00858       NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
00859     }
00860   }
00861 
00862 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
00863   
00864   for (int i=1; i<=proxySpanDim; i++) {
00865     if (tree.size() <= i) break;
00866     child[i-1] = tree[i];
00867     nChild++;
00868   }
00869 #endif
00870 #endif
00871   
00872 #if 0
00873   // for debugging
00874   CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
00875   for (int i=0; i<psize+1; i++) {
00876     CkPrintf("%d ", tree[i]);
00877   }
00878   CkPrintf("\n");
00879 #endif
00880   // send to children nodes
00881   sendSpanningTree();
00882 }

void HomePatch::checkpoint ( void   ) 

Definition at line 3343 of file HomePatch.C.

References ResizeArray< Elem >::copy(), and Patch::lattice.

Referenced by Sequencer::algorithm().

03343                                {
03344   checkpoint_atom.copy(atom);
03345   checkpoint_lattice = lattice;
03346 
03347   // DMK - Atom Separation (water vs. non-water)
03348   #if NAMD_SeparateWaters != 0
03349     checkpoint_numWaterAtoms = numWaterAtoms;
03350   #endif
03351 }

void HomePatch::depositMigration ( MigrateAtomsMsg  ) 

Definition at line 3852 of file HomePatch.C.

References ResizeArray< Elem >::add(), Lattice::apply_transform(), Sequencer::awaken(), ResizeArray< Elem >::begin(), DebugM, ResizeArray< Elem >::end(), inMigration, j, Patch::lattice, MigrateAtomsMsg::migrationList, msgbuf, Lattice::nearest(), numMlBuf, Patch::patchID, Lattice::reverse_transform(), and ResizeArray< Elem >::size().

Referenced by MigrateAtomsCombinedMsg::distribute(), and doAtomMigration().

03853 {
03854 
03855   if (!inMigration) { // We have to buffer changes due to migration
03856                       // until our patch is in migration mode
03857     msgbuf[numMlBuf++] = msg;
03858     return;
03859   } 
03860 
03861 
03862   // DMK - Atom Separation (water vs. non-water)
03863   #if NAMD_SeparateWaters != 0
03864 
03865 
03866     // Merge the incoming list of atoms with the current list of
03867     //   atoms.  Note that mergeSeparatedAtomList() will apply any
03868     //   required transformations to the incoming atoms as it is
03869     //   separating them.
03870     mergeAtomList(msg->migrationList);
03871 
03872 
03873   #else
03874 
03875     // Merge the incoming list of atoms with the current list of
03876     // atoms.  Apply transformations to the incoming atoms as they are
03877     // added to this patch's list.
03878     {
03879       MigrationList &migrationList = msg->migrationList;
03880       MigrationList::iterator mi;
03881       Transform mother_transform;
03882       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
03883         DebugM(1,"Migrating atom " << mi->id << " to patch "
03884                   << patchID << " with position " << mi->position << "\n"); 
03885         if ( mi->migrationGroupSize ) {
03886           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
03887             Position pos = mi->position;
03888             int mgs = mi->migrationGroupSize;
03889             int c = 1;
03890             for ( int j=mi->hydrogenGroupSize; j<mgs;
03891                                 j+=(mi+j)->hydrogenGroupSize ) {
03892               pos += (mi+j)->position;
03893               ++c;
03894             }
03895             pos *= 1./c;
03896             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
03897             mother_transform = mi->transform;
03898             pos = lattice.nearest(pos,center,&mother_transform);
03899             mi->position = lattice.reverse_transform(mi->position,mi->transform);
03900             mi->position = lattice.apply_transform(mi->position,mother_transform);
03901             mi->transform = mother_transform;
03902           } else {
03903             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
03904             mother_transform = mi->transform;
03905           }
03906         } else {
03907           mi->position = lattice.reverse_transform(mi->position,mi->transform);
03908           mi->position = lattice.apply_transform(mi->position,mother_transform);
03909           mi->transform = mother_transform;
03910         }
03911         atom.add(*mi);
03912       }
03913     }
03914 
03915 
03916   #endif // if (NAMD_SeparateWaters != 0)
03917 
03918 
03919   numAtoms = atom.size();
03920   delete msg;
03921 
03922   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
03923   if (!--patchMigrationCounter) {
03924     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
03925     allMigrationIn = true;
03926     patchMigrationCounter = numNeighbors;
03927     if (migrationSuspended) {
03928       DebugM(3,"patch "<<patchID<<" is being awakened\n");
03929       migrationSuspended = false;
03930       sequencer->awaken();
03931       return;
03932     }
03933     else {
03934        DebugM(3,"patch "<<patchID<<" was not suspended\n");
03935     }
03936   }
03937 }

void HomePatch::doAtomMigration (  )  [protected]

Definition at line 3715 of file HomePatch.C.

References ResizeArray< Elem >::add(), atomIndex, Patch::atomMapper, ResizeArray< Elem >::begin(), DebugM, ResizeArray< Elem >::del(), depositMigration(), ResizeArray< Elem >::end(), inMigration, j, Patch::lattice, marginViolations, MigrationInfo::mList, msgbuf, numMlBuf, PatchMgr::Object(), Patch::patchID, CompAtom::position, Lattice::scale(), PatchMgr::sendMigrationMsgs(), ResizeArray< Elem >::size(), Sequencer::suspend(), AtomMapper::unregisterIDsFullAtom(), Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

03716 {
03717   int i;
03718 
03719   for (i=0; i<numNeighbors; i++) {
03720     realInfo[i].mList.resize(0);
03721   }
03722 
03723   // Purge the AtomMap
03724   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03725 
03726   // Determine atoms that need to migrate
03727 
03728   BigReal minx = min.x;
03729   BigReal miny = min.y;
03730   BigReal minz = min.z;
03731   BigReal maxx = max.x;
03732   BigReal maxy = max.y;
03733   BigReal maxz = max.z;
03734 
03735   int xdev, ydev, zdev;
03736   int delnum = 0;
03737 
03738   FullAtomList::iterator atom_i = atom.begin();
03739   FullAtomList::iterator atom_e = atom.end();
03740 
03741   // DMK - Atom Separation (water vs. non-water)
03742   #if NAMD_SeparateWaters != 0
03743     FullAtomList::iterator atom_first = atom_i;
03744     int numLostWaterAtoms = 0;
03745   #endif
03746 
03747   while ( atom_i != atom_e ) {
03748     if ( atom_i->migrationGroupSize ) {
03749       Position pos = atom_i->position;
03750       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
03751         int mgs = atom_i->migrationGroupSize;
03752         int c = 1;
03753         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
03754                                 j+=(atom_i+j)->hydrogenGroupSize ) {
03755           pos += (atom_i+j)->position;
03756           ++c;
03757         }
03758         pos *= 1./c;
03759         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
03760       }
03761 
03762       ScaledPosition s = lattice.scale(pos);
03763 
03764       // check if atom is within bounds
03765       if (s.x < minx) xdev = 0;
03766       else if (maxx <= s.x) xdev = 2;
03767       else xdev = 1;
03768 
03769       if (s.y < miny) ydev = 0;
03770       else if (maxy <= s.y) ydev = 2;
03771       else ydev = 1;
03772 
03773       if (s.z < minz) zdev = 0;
03774       else if (maxz <= s.z) zdev = 2;
03775       else zdev = 1;
03776 
03777     }
03778 
03779     if (mInfo[xdev][ydev][zdev]) { // process atom for migration
03780                                     // Don't migrate if destination is myself
03781 
03782       // See if we have a migration list already
03783       MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
03784       DebugM(3,"Migrating atom " << atom_i->id << " from patch "
03785                 << patchID << " with position " << atom_i->position << "\n");
03786       mCur.add(*atom_i);
03787 
03788       ++delnum;
03789 
03790 
03791       // DMK - Atom Separation (water vs. non-water)
03792       #if NAMD_SeparateWaters != 0
03793         // Check to see if this atom is part of a water molecule.  If
03794         //   so, numWaterAtoms needs to be adjusted to refect the lost of
03795         //   this atom.
03796         // NOTE: The atom separation code assumes that if the oxygen
03797         //   atom of the hydrogen group making up the water molecule is
03798         //   migrated to another HomePatch, the hydrogens will also
03799         //   move!!!
03800         int atomIndex = atom_i - atom_first;
03801         if (atomIndex < numWaterAtoms)
03802           numLostWaterAtoms++;
03803       #endif
03804 
03805 
03806     } else {
03807 
03808       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
03809 
03810     }
03811 
03812     ++atom_i;
03813   }
03814 
03815   // DMK - Atom Separation (water vs. non-water)
03816   #if NAMD_SeparateWaters != 0
03817     numWaterAtoms -= numLostWaterAtoms;
03818   #endif
03819 
03820 
03821   int delpos = numAtoms - delnum;
03822   DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
03823   atom.del(delpos,delnum);
03824 
03825   numAtoms = atom.size();
03826 
03827   PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
03828 
03829   // signal depositMigration() that we are inMigration mode
03830   inMigration = true;
03831 
03832   // Drain the migration message buffer
03833   for (i=0; i<numMlBuf; i++) {
03834      DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
03835      depositMigration(msgbuf[i]);
03836   }
03837   numMlBuf = 0;
03838      
03839   if (!allMigrationIn) {
03840     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
03841     migrationSuspended = true;
03842     sequencer->suspend();
03843     migrationSuspended = false;
03844   }
03845   allMigrationIn = false;
03846 
03847   inMigration = false;
03848   marginViolations = 0;
03849 }

void HomePatch::doGroupSizeCheck (  )  [protected]

Definition at line 3594 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Flags::doNonbonded, ResizeArray< Elem >::end(), Patch::flags, Flags::maxGroupRadius, NAMD_bug(), Node::Object(), Node::simParameters, simParams, x, y, and z.

Referenced by positionsReady().

03595 {
03596   if ( ! flags.doNonbonded ) return;
03597 
03598   SimParameters *simParams = Node::Object()->simParameters;
03599   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
03600   BigReal maxrad2 = 0.;
03601 
03602   FullAtomList::iterator p_i = atom.begin();
03603   FullAtomList::iterator p_e = atom.end();
03604 
03605   while ( p_i != p_e ) {
03606     const int hgs = p_i->hydrogenGroupSize;
03607     if ( ! hgs ) break;  // avoid infinite loop on bug
03608     int ngs = hgs;
03609     if ( ngs > 5 ) ngs = 5;  // limit to at most 5 atoms per group
03610     BigReal x = p_i->position.x;
03611     BigReal y = p_i->position.y;
03612     BigReal z = p_i->position.z;
03613     int i;
03614     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
03615       p_i[i].nonbondedGroupSize = 0;
03616       BigReal dx = p_i[i].position.x - x;
03617       BigReal dy = p_i[i].position.y - y;
03618       BigReal dz = p_i[i].position.z - z;
03619       BigReal r2 = dx * dx + dy * dy + dz * dz;
03620       if ( r2 > hgcut ) break;
03621       else if ( r2 > maxrad2 ) maxrad2 = r2;
03622     }
03623     p_i->nonbondedGroupSize = i;
03624     for ( ; i < hgs; ++i ) {
03625       p_i[i].nonbondedGroupSize = 1;
03626     }
03627     p_i += hgs;
03628   }
03629 
03630   if ( p_i != p_e ) {
03631     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
03632   }
03633 
03634   flags.maxGroupRadius = sqrt(maxrad2);
03635 
03636 }

void HomePatch::doMarginCheck (  )  [protected]

Definition at line 3638 of file HomePatch.C.

References Lattice::a(), Lattice::a_r(), Lattice::b(), Lattice::b_r(), ResizeArray< Elem >::begin(), Lattice::c(), Lattice::c_r(), ResizeArray< Elem >::end(), Patch::lattice, marginViolations, NAMD_die(), Node::Object(), Lattice::scale(), Node::simParameters, simParams, Vector::unit(), Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

03639 {
03640   SimParameters *simParams = Node::Object()->simParameters;
03641 
03642   BigReal sysdima = lattice.a_r().unit() * lattice.a();
03643   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
03644   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
03645 
03646   BigReal minSize = simParams->patchDimension - simParams->margin;
03647 
03648   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
03649        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
03650        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
03651 
03652     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03653       "Possible solutions are to restart from a recent checkpoint,\n"
03654       "increase margin, or disable useFlexibleCell for liquid simulation.");
03655   }
03656 
03657   BigReal cutoff = simParams->cutoff;
03658 
03659   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
03660   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
03661   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
03662 
03663   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
03664     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03665       "There are probably many margin violations already reported.\n"
03666       "Possible solutions are to restart from a recent checkpoint,\n"
03667       "increase margin, or disable useFlexibleCell for liquid simulation.");
03668   }
03669 
03670   BigReal minx = min.x - margina;
03671   BigReal miny = min.y - marginb;
03672   BigReal minz = min.z - marginc;
03673   BigReal maxx = max.x + margina;
03674   BigReal maxy = max.y + marginb;
03675   BigReal maxz = max.z + marginc;
03676 
03677   int xdev, ydev, zdev;
03678   int problemCount = 0;
03679 
03680   FullAtomList::iterator p_i = atom.begin();
03681   FullAtomList::iterator p_e = atom.end();
03682   for ( ; p_i != p_e; ++p_i ) {
03683 
03684     ScaledPosition s = lattice.scale(p_i->position);
03685 
03686     // check if atom is within bounds
03687     if (s.x < minx) xdev = 0;
03688     else if (maxx <= s.x) xdev = 2; 
03689     else xdev = 1;
03690 
03691     if (s.y < miny) ydev = 0;
03692     else if (maxy <= s.y) ydev = 2; 
03693     else ydev = 1;
03694 
03695     if (s.z < minz) zdev = 0;
03696     else if (maxz <= s.z) zdev = 2; 
03697     else zdev = 1;
03698 
03699     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
03700         ++problemCount;
03701     }
03702 
03703   }
03704 
03705   marginViolations = problemCount;
03706   // if ( problemCount ) {
03707   //     iout << iERROR <<
03708   //       "Found " << problemCount << " margin violations!\n" << endi;
03709   // } 
03710 
03711 }

void HomePatch::doPairlistCheck (  )  [protected]

Definition at line 3519 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Patch::flags, j, Patch::lattice, Flags::maxAtomMovement, Node::Object(), Flags::pairlistTolerance, ResizeArray< Elem >::resize(), Flags::savePairlists, Node::simParameters, simParams, Lattice::unscale(), Flags::usePairlists, Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

03520 {
03521   SimParameters *simParams = Node::Object()->simParameters;
03522 
03523   if ( numAtoms == 0 || ! flags.usePairlists ) {
03524     flags.pairlistTolerance = 0.;
03525     flags.maxAtomMovement = 99999.;
03526     return;
03527   }
03528 
03529   int i; int n = numAtoms;
03530   CompAtom *p_i = p.begin();
03531 
03532   if ( flags.savePairlists ) {
03533     flags.pairlistTolerance = doPairlistCheck_newTolerance;
03534     flags.maxAtomMovement = 0.;
03535     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
03536     doPairlistCheck_lattice = lattice;
03537     doPairlistCheck_positions.resize(numAtoms);
03538     CompAtom *psave_i = doPairlistCheck_positions.begin();
03539     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
03540     return;
03541   }
03542 
03543   Lattice &lattice_old = doPairlistCheck_lattice;
03544   Position center_cur = lattice.unscale(center);
03545   Position center_old = lattice_old.unscale(center);
03546   Vector center_delta = center_cur - center_old;
03547   
03548   // find max deviation to corner (any neighbor shares a corner)
03549   BigReal max_cd = 0.;
03550   for ( i=0; i<2; ++i ) {
03551     for ( int j=0; j<2; ++j ) {
03552       for ( int k=0; k<2; ++k ) {
03553         ScaledPosition corner(  i ? min.x : max.x ,
03554                                 j ? min.y : max.y ,
03555                                 k ? min.z : max.z );
03556         Vector corner_delta =
03557                 lattice.unscale(corner) - lattice_old.unscale(corner);
03558         corner_delta -= center_delta;
03559         BigReal cd = corner_delta.length2();
03560         if ( cd > max_cd ) max_cd = cd;
03561       }
03562     }
03563   }
03564   max_cd = sqrt(max_cd);
03565 
03566   // find max deviation of atoms relative to center
03567   BigReal max_pd = 0.;
03568   CompAtom *p_old_i = doPairlistCheck_positions.begin();
03569   for ( i=0; i<n; ++i ) {
03570     Vector p_delta = p_i[i].position - p_old_i[i].position;
03571     p_delta -= center_delta;
03572     BigReal pd = p_delta.length2();
03573     if ( pd > max_pd ) max_pd = pd;
03574   }
03575   max_pd = sqrt(max_pd);
03576 
03577   BigReal max_tol = max_pd + max_cd;
03578 
03579   flags.maxAtomMovement = max_tol;
03580 
03581   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
03582 
03583   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
03584                                 doPairlistCheck_newTolerance ) ) {
03585     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
03586   }
03587 
03588   if ( max_tol > doPairlistCheck_newTolerance ) {
03589     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
03590   }
03591 
03592 }

void HomePatch::exchangeAtoms ( int  scriptTask  ) 

Definition at line 3467 of file HomePatch.C.

References ExchangeAtomsMsg::atoms, ResizeArray< Elem >::begin(), exchange_dst, exchange_msg, exchange_req, exchange_src, Patch::lattice, ExchangeAtomsMsg::lattice, ExchangeAtomsMsg::numAtoms, PatchMgr::Object(), Node::Object(), Patch::patchID, ExchangeAtomsMsg::pid, recvExchangeReq(), SCRIPT_ATOMRECV, SCRIPT_ATOMSEND, SCRIPT_ATOMSENDRECV, PatchMgr::sendExchangeReq(), Node::simParameters, and simParams.

Referenced by Sequencer::algorithm().

03467                                             {
03468   SimParameters *simParams = Node::Object()->simParameters;
03469   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
03470   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
03471     exchange_dst = (int) simParams->scriptArg1;
03472     // create and save outgoing message
03473     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
03474     exchange_msg->lattice = lattice;
03475     exchange_msg->pid = patchID;
03476     exchange_msg->numAtoms = numAtoms;
03477     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03478     if ( exchange_req >= 0 ) {
03479       recvExchangeReq(exchange_req);
03480     }
03481   }
03482   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
03483     exchange_src = (int) simParams->scriptArg2;
03484     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
03485   }
03486 }

void HomePatch::exchangeCheckpoint ( int  scriptTask,
int &  bpc 
)

Definition at line 3371 of file HomePatch.C.

References checkpoint_task, PatchMgr::Object(), Node::Object(), Patch::patchID, PatchMgr::sendCheckpointReq(), Node::simParameters, and simParams.

Referenced by Sequencer::algorithm().

03371                                                            {  // initiating replica
03372   SimParameters *simParams = Node::Object()->simParameters;
03373   checkpoint_task = scriptTask;
03374   const int remote = simParams->scriptIntArg1;
03375   const char *key = simParams->scriptStringArg1;
03376   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
03377 }

void HomePatch::gbisComputeAfterP1 (  ) 

Definition at line 3040 of file HomePatch.C.

References Patch::bornRad, Patch::flags, gbisP2Ready(), Patch::intRad, Node::Object(), Patch::pExt, Patch::psiFin, Patch::psiSum, Flags::sequence, Node::simParameters, and simParams.

03040                                    {
03041 
03042   SimParameters *simParams = Node::Object()->simParameters;
03043   BigReal alphaMax = simParams->alpha_max;
03044   BigReal delta = simParams->gbis_delta;
03045   BigReal beta = simParams->gbis_beta;
03046   BigReal gamma = simParams->gbis_gamma;
03047   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03048 
03049   BigReal rhoi;
03050   BigReal rhoi0;
03051   //calculate bornRad from psiSum
03052   for (int i = 0; i < numAtoms; i++) {
03053     rhoi0 = intRad[2*i];
03054     rhoi = rhoi0+coulomb_radius_offset;
03055     psiFin[i] += psiSum[i];
03056     psiFin[i] *= rhoi0;
03057     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
03058     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
03059 #ifdef PRINT_COMP
03060     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
03061 #endif
03062   }
03063 
03064   gbisP2Ready();
03065 }

void HomePatch::gbisComputeAfterP2 (  ) 

Definition at line 3068 of file HomePatch.C.

References Patch::bornRad, COULOMB, Patch::dEdaSum, Patch::dHdrPrefix, Patch::flags, gbisP3Ready(), Patch::intRad, Node::Object(), Patch::pExt, Patch::psiFin, Flags::sequence, Node::simParameters, and simParams.

03068                                    {
03069 
03070   SimParameters *simParams = Node::Object()->simParameters;
03071   BigReal delta = simParams->gbis_delta;
03072   BigReal beta = simParams->gbis_beta;
03073   BigReal gamma = simParams->gbis_gamma;
03074   BigReal epsilon_s = simParams->solvent_dielectric;
03075   BigReal epsilon_p = simParams->dielectric;
03076   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
03077   BigReal epsilon_p_i = 1/simParams->dielectric;
03078   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03079   BigReal kappa = simParams->kappa;
03080   BigReal fij, expkappa, Dij, dEdai, dedasum;
03081   BigReal rhoi, rhoi0, psii, nbetapsi;
03082   BigReal gammapsi2, tanhi, daidr;
03083   for (int i = 0; i < numAtoms; i++) {
03084     //add diagonal dEda term
03085     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
03086     fij = bornRad[i];//inf
03087     expkappa = exp(-kappa*fij);//0
03088     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
03089     //calculate dHij prefix
03090     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
03091                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
03092     dHdrPrefix[i] += dEdai;
03093     dedasum = dHdrPrefix[i];
03094 
03095     rhoi0 = intRad[2*i];
03096     rhoi = rhoi0+coulomb_radius_offset;
03097     psii = psiFin[i];
03098     nbetapsi = -beta*psii;
03099     gammapsi2 = gamma*psii*psii;
03100     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
03101     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
03102            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
03103     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
03104 #ifdef PRINT_COMP
03105     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
03106 #endif
03107   }
03108   gbisP3Ready();
03109 }

void HomePatch::gbisP2Ready (  ) 

Reimplemented from Patch.

Definition at line 3112 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Patch::bornRad, ProxyGBISP2DataMsg::bornRad, ProxyGBISP2DataMsg::destPe, Patch::flags, GB2_PROXY_DATA_PRIORITY, Patch::gbisP2Ready(), ProxyGBISP2DataMsg::origPe, ProxyGBISP2DataMsg::patch, PATCH_PRIORITY, Patch::patchID, PRIORITY_SIZE, Flags::sequence, SET_PRIORITY, and ResizeArray< Elem >::size().

Referenced by gbisComputeAfterP1().

03112                             {
03113   if (proxy.size() > 0) {
03114     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03115     for (int i = 0; i < proxy.size(); i++) {
03116       int node = proxy[i];
03117       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
03118       msg->patch = patchID;
03119       msg->origPe = CkMyPe();
03120       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
03121       msg->destPe = node;
03122       int seq = flags.sequence;
03123       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03124       SET_PRIORITY(msg,seq,priority);
03125       cp[node].recvData(msg);
03126     }
03127   }
03128   Patch::gbisP2Ready();
03129 }

void HomePatch::gbisP3Ready (  ) 

Reimplemented from Patch.

Definition at line 3132 of file HomePatch.C.

References ResizeArray< Elem >::begin(), ProxyGBISP3DataMsg::destPe, Patch::dHdrPrefix, ProxyGBISP3DataMsg::dHdrPrefix, ProxyGBISP3DataMsg::dHdrPrefixLen, Flags::doFullElectrostatics, Patch::flags, GB3_PROXY_DATA_PRIORITY, Patch::gbisP3Ready(), ProxyGBISP3DataMsg::origPe, ProxyGBISP3DataMsg::patch, PATCH_PRIORITY, Patch::patchID, PRIORITY_SIZE, Flags::sequence, SET_PRIORITY, and ResizeArray< Elem >::size().

Referenced by gbisComputeAfterP2().

03132                             {
03133   if (proxy.size() > 0) {
03134     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03135     //only nonzero message should be sent for doFullElec
03136     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
03137     for (int i = 0; i < proxy.size(); i++) {
03138       int node = proxy[i];
03139       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
03140       msg->patch = patchID;
03141       msg->dHdrPrefixLen = msgAtoms;
03142       msg->origPe = CkMyPe();
03143       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
03144       msg->destPe = node;
03145       int seq = flags.sequence;
03146       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03147       SET_PRIORITY(msg,seq,priority);
03148       cp[node].recvData(msg);
03149     }
03150   }
03151   Patch::gbisP3Ready();
03152 }

FullAtomList& HomePatch::getAtomList (  )  [inline]

Definition at line 205 of file HomePatch.h.

Referenced by ComputeHomePatch::doWork().

00205 { return (atom); }

int HomePatch::hardWallDrude ( const   BigReal,
Tensor virial,
SubmitReduction  
)

Definition at line 1915 of file HomePatch.C.

References BOLTZMANN, Lattice::c(), endi(), iERROR(), iout, SubmitReduction::item(), j, Patch::lattice, Node::molecule, Node::Object(), Lattice::origin(), outer(), partition(), Node::simParameters, simParams, TIMEFACTOR, Vector::x, Tensor::xx, Vector::y, Tensor::yy, z, Vector::z, and Tensor::zz.

Referenced by Sequencer::hardWallDrude().

01917 {
01918   Molecule *mol = Node::Object()->molecule;
01919   SimParameters *simParams = Node::Object()->simParameters;
01920   const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
01921   const int fixedAtomsOn = simParams->fixedAtomsOn;
01922   const BigReal dt = timestep / TIMEFACTOR;
01923   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
01924   int i, ia, ib, j;
01925   int dieOnError = simParams->rigidDie;
01926   Tensor wc;  // constraint virial
01927   BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
01928   int nslabs;
01929 
01930   // start data for hard wall boundary between drude and its host atom
01931   // static int Count=0;
01932   int Idx;
01933   double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
01934   Vector v_ab, vb_1, vp_1, vb_2, vp_2, new_vel_a, new_vel_b, new_pos_a, new_pos_b, *new_pos, *new_vel;
01935   double dot_v_r_1, dot_v_r_2;
01936   double vb_cm, dr_a, dr_b;
01937   // end data for hard wall boundary between drude and its host atom
01938 
01939   // start calculation of hard wall boundary between drude and its host atom
01940   if (simParams->drudeHardWallOn) {
01941     if (ppreduction) {
01942       nslabs = simParams->pressureProfileSlabs;
01943       idz = nslabs/lattice.c().z;
01944       zmin = lattice.origin().z - 0.5*lattice.c().z;
01945     }
01946 
01947     r_wall = simParams->drudeBondLen;
01948     r_wall_SQ = r_wall*r_wall;
01949     // Count++;
01950     for (i=1; i<numAtoms; i++)  {
01951       if ( (atom[i].mass > 0.01) && ((atom[i].mass < 1.0)) ) { // drude particle
01952         ia = i-1;
01953         ib = i;
01954 
01955         v_ab = atom[ib].position - atom[ia].position;
01956         rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
01957 
01958         if (rab_SQ > r_wall_SQ) {  // to impose the hard wall constraint
01959           rab = sqrt(rab_SQ);
01960           if ( (rab > (2.0*r_wall)) && dieOnError ) {  // unexpected situation
01961             iout << iERROR << "HardWallDrude> "
01962               << "The drude is too far away from atom "
01963               << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
01964             return -1;  // triggers early exit
01965           }
01966 
01967           v_ab.x /= rab;
01968           v_ab.y /= rab;
01969           v_ab.z /= rab;
01970 
01971           if ( fixedAtomsOn && atom[ia].atomFixed ) {  // the heavy atom is fixed
01972             if (atom[ib].atomFixed) {  // the drude is fixed too
01973               continue;
01974             }
01975             else {  // only the heavy atom is fixed
01976               dot_v_r_2 = atom[ib].velocity.x*v_ab.x
01977                 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
01978               vb_2 = dot_v_r_2 * v_ab;
01979               vp_2 = atom[ib].velocity - vb_2;
01980 
01981               dr = rab - r_wall;
01982               if(dot_v_r_2 == 0.0) {
01983                 delta_T = maxtime;
01984               }
01985               else {
01986                 delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
01987                 if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
01988               }
01989 
01990               dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
01991 
01992               vb_2 = dot_v_r_2 * v_ab;
01993 
01994               new_vel_a = atom[ia].velocity;
01995               new_vel_b = vp_2 + vb_2;
01996 
01997               dr_b = -dr + delta_T*dot_v_r_2;  // L = L_0 + dT *v_new, v was flipped
01998 
01999               new_pos_a = atom[ia].position;
02000               new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
02001             }
02002           }
02003           else {
02004             mass_a = atom[ia].mass;
02005             mass_b = atom[ib].mass;
02006             mass_sum = mass_a+mass_b;
02007 
02008             dot_v_r_1 = atom[ia].velocity.x*v_ab.x
02009               + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
02010             vb_1 = dot_v_r_1 * v_ab;
02011             vp_1 = atom[ia].velocity - vb_1;
02012 
02013             dot_v_r_2 = atom[ib].velocity.x*v_ab.x
02014               + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
02015             vb_2 = dot_v_r_2 * v_ab;
02016             vp_2 = atom[ib].velocity - vb_2;
02017 
02018             vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
02019 
02020             dot_v_r_1 -= vb_cm;
02021             dot_v_r_2 -= vb_cm;
02022 
02023             dr = rab - r_wall;
02024 
02025             if(dot_v_r_2 == dot_v_r_1) {
02026                 delta_T = maxtime;
02027             }
02028             else {
02029               delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);  // the time since the collision occurs
02030               if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
02031             }
02032             
02033             // the relative velocity between ia and ib. Drawn according to T_Drude
02034             v_Bond = sqrt(kbt/mass_b);
02035 
02036             // reflect the velocity along bond vector and scale down
02037             dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
02038             dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
02039 
02040             dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
02041             dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
02042 
02043             new_pos_a = atom[ia].position + dr_a*v_ab;  // correct the position
02044             new_pos_b = atom[ib].position + dr_b*v_ab;
02045             // atom[ia].position += (dr_a*v_ab);  // correct the position
02046             // atom[ib].position += (dr_b*v_ab);
02047             
02048             dot_v_r_1 += vb_cm;
02049             dot_v_r_2 += vb_cm;
02050 
02051             vb_1 = dot_v_r_1 * v_ab;
02052             vb_2 = dot_v_r_2 * v_ab;
02053 
02054             new_vel_a = vp_1 + vb_1;
02055             new_vel_b = vp_2 + vb_2;
02056           }
02057 
02058           int ppoffset, partition;
02059           if ( invdt == 0 ) {
02060             atom[ia].position = new_pos_a;
02061             atom[ib].position = new_pos_b;
02062           }
02063           else if ( virial == 0 ) {
02064             atom[ia].velocity = new_vel_a;
02065             atom[ib].velocity = new_vel_b;
02066           }
02067           else {
02068             for ( j = 0; j < 2; j++ ) {
02069               if (j==0) {  // atom ia, heavy atom
02070                 Idx = ia;
02071                 new_pos = &new_pos_a;
02072                 new_vel = &new_vel_a;
02073               }
02074               else if (j==1) {  // atom ib, drude
02075                 Idx = ib;
02076                 new_pos = &new_pos_b;
02077                 new_vel = &new_vel_b;
02078               }
02079               Force df = (*new_vel - atom[Idx].velocity) *
02080                 ( atom[Idx].mass * invdt );
02081               Tensor vir = outer(df, atom[Idx].position);
02082               wc += vir;
02083               atom[Idx].velocity = *new_vel;
02084               atom[Idx].position = *new_pos;
02085 
02086               if (ppreduction) {
02087                 if (!j) {
02088                   BigReal z = new_pos->z;
02089                   int partition = atom[Idx].partition;
02090                   int slab = (int)floor((z-zmin)*idz);
02091                   if (slab < 0) slab += nslabs;
02092                   else if (slab >= nslabs) slab -= nslabs;
02093                   ppoffset = 3*(slab + nslabs*partition);
02094                 }
02095                 ppreduction->item(ppoffset  ) += vir.xx;
02096                 ppreduction->item(ppoffset+1) += vir.yy;
02097                 ppreduction->item(ppoffset+2) += vir.zz;
02098               }
02099 
02100             }
02101           }
02102         }                               
02103       }
02104     }
02105 
02106     // if ( (Count>10000) && (Count%10==0) ) {
02107     //   v_ab = atom[1].position - atom[0].position;
02108     //   rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
02109     //   iout << "DBG_R: " << Count << "  " << sqrt(rab_SQ) << "\n" << endi;
02110     // }
02111 
02112   }
02113 
02114   // end calculation of hard wall boundary between drude and its host atom
02115 
02116   if ( dt && virial ) *virial += wc;
02117 
02118   return 0;
02119 }

void HomePatch::loweAndersenFinish (  ) 

Definition at line 3006 of file HomePatch.C.

References DebugM.

03007 {
03008     DebugM(2, "loweAndersenFinish\n");
03009     v.resize(0);
03010 }

void HomePatch::loweAndersenVelocities (  ) 

Definition at line 2991 of file HomePatch.C.

References DebugM, Node::molecule, Node::Object(), Node::simParameters, and simParams.

Referenced by positionsReady().

02992 {
02993     DebugM(2, "loweAndersenVelocities\n");
02994     Molecule *mol = Node::Object()->molecule;
02995     SimParameters *simParams = Node::Object()->simParameters;
02996     v.resize(numAtoms);
02997     for (int i = 0; i < numAtoms; ++i) {
02998         //v[i] = p[i];
02999         // co-opt CompAtom structure to pass velocity and mass information
03000         v[i].position = atom[i].velocity;
03001         v[i].charge = atom[i].mass;
03002     }
03003     DebugM(2, "loweAndersenVelocities\n");
03004 }

void HomePatch::minimize_rattle2 ( const   BigReal,
Tensor virial,
bool  forces = false 
)

Definition at line 2856 of file HomePatch.C.

References ResizeArray< Elem >::begin(), endi(), Patch::f, iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Results::normal, Node::Object(), outer(), settle2(), Node::simParameters, simParams, TIMEFACTOR, WAT_SWM4, WAT_TIP4, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::newMinimizeDirection(), and Sequencer::submitMinimizeReductions().

02857 {
02858   Molecule *mol = Node::Object()->molecule;
02859   SimParameters *simParams = Node::Object()->simParameters;
02860   Force *f1 = f[Results::normal].begin();
02861   const int fixedAtomsOn = simParams->fixedAtomsOn;
02862   const int useSettle = simParams->useSettle;
02863   const BigReal dt = timestep / TIMEFACTOR;
02864   Tensor wc;  // constraint virial
02865   BigReal tol = simParams->rigidTol;
02866   int maxiter = simParams->rigidIter;
02867   int dieOnError = simParams->rigidDie;
02868   int i, iter;
02869   BigReal dsqi[10], tmp;
02870   int ial[10], ibl[10];
02871   Vector ref[10];  // reference position
02872   Vector refab[10];  // reference vector
02873   Vector vel[10];  // new velocity
02874   BigReal rmass[10];  // 1 / mass
02875   BigReal redmass[10];  // reduced mass
02876   int fixed[10];  // is atom fixed?
02877   
02878   // Size of a hydrogen group for water
02879   int wathgsize = 3;
02880   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02881   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02882 
02883   //  CkPrintf("In rattle2!\n");
02884   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02885     //    CkPrintf("ig=%d\n",ig);
02886     int hgs = atom[ig].hydrogenGroupSize;
02887     if ( hgs == 1 ) continue;  // only one atom in group
02888     // cache data in local arrays and integrate positions normally
02889     int anyfixed = 0;
02890     for ( i = 0; i < hgs; ++i ) {
02891       ref[i] = atom[ig+i].position;
02892       vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
02893       rmass[i] = 1.0;
02894       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02895       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02896     }
02897     int icnt = 0;
02898     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02899       if (hgs != wathgsize) {
02900         NAMD_bug("Hydrogen group error caught in rattle2().");
02901       }
02902       // Use SETTLE for water unless some of the water atoms are fixed,
02903       if (useSettle && !anyfixed) {
02904         if (simParams->watmodel == WAT_SWM4) {
02905           // SWM4 ordering:  O D LP H1 H2
02906           // do swap(O,LP) and call settle with subarray O H1 H2
02907           // swap back after we return
02908           Vector lp_ref = ref[2];
02909           // Vector lp_vel = vel[2];
02910           ref[2] = ref[0];
02911           vel[2] = vel[0];
02912           settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
02913           ref[0] = ref[2];
02914           vel[0] = vel[2];
02915           ref[2] = lp_ref;
02916           // vel[2] = vel[0];  // set LP vel to O vel
02917         }
02918         else {
02919           settle2(1.0, 1.0, ref, vel, dt, virial);
02920           if (simParams->watmodel == WAT_TIP4) {
02921             vel[3] = vel[0];
02922           }
02923         }
02924         for (i=0; i<hgs; i++) {
02925           ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
02926         }
02927         continue;
02928       }
02929       if ( !(fixed[1] && fixed[2]) ) {
02930         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02931         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02932       }
02933     }
02934     //    CkPrintf("Loop 2\n");
02935     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02936       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02937         if ( !(fixed[0] && fixed[i]) ) {
02938           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02939           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02940           ibl[icnt] = i;  ++icnt;
02941         }
02942       }
02943     }
02944     if ( icnt == 0 ) continue;  // no constraints
02945     //    CkPrintf("Loop 3\n");
02946     for ( i = 0; i < icnt; ++i ) {
02947       refab[i] = ref[ial[i]] - ref[ibl[i]];
02948     }
02949     //    CkPrintf("Loop 4\n");
02950     int done;
02951     for ( iter = 0; iter < maxiter; ++iter ) {
02952       done = 1;
02953       for ( i = 0; i < icnt; ++i ) {
02954         int a = ial[i];  int b = ibl[i];
02955         Vector vab = vel[a] - vel[b];
02956         Vector &rab = refab[i];
02957         BigReal rabsqi = dsqi[i];
02958         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02959         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02960           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02961           wc += outer(dp,rab);
02962           vel[a] += rmass[a] * dp;
02963           vel[b] -= rmass[b] * dp;
02964           done = 0;
02965         }
02966       }
02967       if ( done ) break;
02968       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02969     }
02970     if ( ! done ) {
02971       if ( dieOnError ) {
02972         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02973       } else {
02974         iout << iWARN <<
02975           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02976       }
02977     }
02978     // store data back to patch
02979     for ( i = 0; i < hgs; ++i ) {
02980       ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
02981     }
02982   }
02983   //  CkPrintf("Leaving rattle2!\n");
02984   // check that there isn't a constant needed here!
02985   *virial += wc / ( 0.5 * dt );
02986 
02987 }

void HomePatch::mollyAverage (  ) 

Definition at line 3206 of file HomePatch.C.

References average(), endi(), iout, iWARN(), Node::molecule, NAMD_die(), Node::Object(), Patch::p_avg, ResizeArray< Elem >::resize(), Node::simParameters, and simParams.

Referenced by positionsReady().

03207 {
03208   Molecule *mol = Node::Object()->molecule;
03209   SimParameters *simParams = Node::Object()->simParameters;
03210   BigReal tol = simParams->mollyTol;
03211   int maxiter = simParams->mollyIter;
03212   int i, iter;
03213   HGArrayBigReal dsq;
03214   BigReal tmp;
03215   HGArrayInt ial, ibl;
03216   HGArrayVector ref;  // reference position
03217   HGArrayVector refab;  // reference vector
03218   HGArrayBigReal rmass;  // 1 / mass
03219   BigReal *lambda;  // Lagrange multipliers
03220   CompAtom *avg;  // averaged position
03221   int numLambdas = 0;
03222   HGArrayInt fixed;  // is atom fixed?
03223 
03224   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
03225   p_avg.resize(numAtoms);
03226   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
03227 
03228   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03229     int hgs = atom[ig].hydrogenGroupSize;
03230     if ( hgs == 1 ) continue;  // only one atom in group
03231         for ( i = 0; i < hgs; ++i ) {
03232           ref[i] = atom[ig+i].position;
03233           rmass[i] = 1. / atom[ig+i].mass;
03234           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03235           if ( fixed[i] ) rmass[i] = 0.;
03236         }
03237         avg = &(p_avg[ig]);
03238         int icnt = 0;
03239 
03240   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
03241           if ( hgs != 3 ) {
03242             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
03243           }
03244           if ( !(fixed[1] && fixed[2]) ) {
03245             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03246           }
03247         }
03248         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03249     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
03250             if ( !(fixed[0] && fixed[i]) ) {
03251               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03252             }
03253           }
03254         }
03255         if ( icnt == 0 ) continue;  // no constraints
03256         numLambdas += icnt;
03257         molly_lambda.resize(numLambdas);
03258         lambda = &(molly_lambda[numLambdas - icnt]);
03259         for ( i = 0; i < icnt; ++i ) {
03260           refab[i] = ref[ial[i]] - ref[ibl[i]];
03261         }
03262         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
03263         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
03264         if ( iter == maxiter ) {
03265           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
03266         }
03267   }
03268 
03269   // for ( i=0; i<numAtoms; ++i ) {
03270   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
03271   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
03272   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
03273   //    }
03274   // }
03275 
03276 }

void HomePatch::mollyMollify ( Tensor virial  ) 

Definition at line 3280 of file HomePatch.C.

References Patch::f, Node::molecule, mollify(), NAMD_die(), Node::Object(), outer(), Patch::p_avg, ResizeArray< Elem >::resize(), Node::simParameters, simParams, and Results::slow.

03281 {
03282   Molecule *mol = Node::Object()->molecule;
03283   SimParameters *simParams = Node::Object()->simParameters;
03284   Tensor wc;  // constraint virial
03285   int i;
03286   HGArrayInt ial, ibl;
03287   HGArrayVector ref;  // reference position
03288   CompAtom *avg;  // averaged position
03289   HGArrayVector refab;  // reference vector
03290   HGArrayVector force;  // new force
03291   HGArrayBigReal rmass;  // 1 / mass
03292   BigReal *lambda;  // Lagrange multipliers
03293   int numLambdas = 0;
03294   HGArrayInt fixed;  // is atom fixed?
03295 
03296   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03297     int hgs = atom[ig].hydrogenGroupSize;
03298     if (hgs == 1 ) continue;  // only one atom in group
03299         for ( i = 0; i < hgs; ++i ) {
03300           ref[i] = atom[ig+i].position;
03301           force[i] = f[Results::slow][ig+i];
03302           rmass[i] = 1. / atom[ig+i].mass;
03303           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03304           if ( fixed[i] ) rmass[i] = 0.;
03305         }
03306         int icnt = 0;
03307         // c-ji I'm only going to mollify water for now
03308   if ( atom[ig].rigidBondLength > 0 ) {  // for water
03309           if ( hgs != 3 ) {
03310             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
03311           }
03312           if ( !(fixed[1] && fixed[2]) ) {
03313             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03314           }
03315         }
03316         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03317     if ( atom[ig+i].rigidBondLength > 0 ) {
03318             if ( !(fixed[0] && fixed[i]) ) {
03319               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03320             }
03321           }
03322         }
03323 
03324         if ( icnt == 0 ) continue;  // no constraints
03325         lambda = &(molly_lambda[numLambdas]);
03326         numLambdas += icnt;
03327         for ( i = 0; i < icnt; ++i ) {
03328           refab[i] = ref[ial[i]] - ref[ibl[i]];
03329         }
03330         avg = &(p_avg[ig]);
03331         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
03332         // store data back to patch
03333         for ( i = 0; i < hgs; ++i ) {
03334           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
03335           f[Results::slow][ig+i] = force[i];
03336         }
03337   }
03338   // check that there isn't a constant needed here!
03339   *virial += wc;
03340   p_avg.resize(0);
03341 }

void HomePatch::positionsReady ( int  doMigration = 0  ) 

Reimplemented from Patch.

Definition at line 1015 of file HomePatch.C.

References ProxyDataMsg::avgPlLen, ProxyDataMsg::avgPositionList, ResizeArray< Elem >::begin(), CompAtom::charge, COULOMB, Patch::cudaAtomPtr, ComputeNonbondedUtil::dielectric_1, doAtomMigration(), Flags::doGBIS, doGroupSizeCheck(), Flags::doLCPO, Flags::doLoweAndersen, doMarginCheck(), Flags::doMolly, doPairlistCheck(), ResizeArray< Elem >::end(), ProxyDataMsg::flags, Patch::flags, Patch::intRad, ProxyDataMsg::intRadList, j, Patch::lattice, Patch::lcpoType, ProxyDataMsg::lcpoTypeList, loweAndersenVelocities(), mollyAverage(), PatchMap::numNodesWithPatches(), ProxyMgr::Object(), PatchMap::Object(), Node::Object(), Patch::p, Patch::p_avg, ProxyDataMsg::patch, PATCH_PRIORITY, Patch::patchID, Patch::pExt, ProxyDataMsg::plLen, CompAtom::position, ProxyDataMsg::positionList, PRIORITY_SIZE, PROXY_DATA_PRIORITY, proxySendSpanning, proxySpanDim, CudaAtom::q, qmSwapAtoms(), rattleListValid, ResizeArray< Elem >::resize(), ComputeNonbondedUtil::scaling, ProxyMgr::sendProxyAll(), ProxyMgr::sendProxyData(), Flags::sequence, SET_PRIORITY, setGBISIntrinsicRadii(), setLcpoType(), Node::simParameters, simParams, ResizeArray< Elem >::size(), sortAtomsForCUDA(), CompAtomExt::sortOrder, Lattice::unscale(), Patch::v, ProxyDataMsg::velocityList, ProxyDataMsg::vlLen, CudaAtom::x, Vector::x, CudaAtom::y, Vector::y, CudaAtom::z, and Vector::z.

01016 {    
01017   SimParameters *simParams = Node::Object()->simParameters;
01018 
01019   flags.sequence++;
01020 
01021   if (!patchMapRead) { readPatchMap(); }
01022       
01023   if (numNeighbors && ! simParams->staticAtomAssignment) {
01024     if (doMigration) {
01025       rattleListValid = false;
01026       doAtomMigration();
01027     } else {
01028       doMarginCheck();
01029     }
01030   }
01031   
01032   if (doMigration && simParams->qmLSSOn)
01033       qmSwapAtoms();
01034 
01035 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
01036   if ( doMigration ) {
01037     int n = numAtoms;
01038     FullAtom *a_i = atom.begin();
01039     #if defined(NAMD_CUDA) || (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
01040     int *ao = new int[n];
01041     int nfree;
01042     if ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces ) {
01043       int k = 0;
01044       int k2 = n;
01045       for ( int j=0; j<n; ++j ) {
01046         // put fixed atoms at end
01047         if ( a_i[j].atomFixed ) ao[--k2] = j;
01048         else ao[k++] = j;
01049       }
01050       nfree = k;
01051     } else {
01052       nfree = n;
01053       for ( int j=0; j<n; ++j ) {
01054         ao[j] = j;
01055       }
01056     }
01057 
01058     sortAtomsForCUDA(ao,a_i,nfree,n);
01059   
01060     for ( int i=0; i<n; ++i ) { 
01061       a_i[i].sortOrder = ao[i];
01062     }
01063     delete [] ao;
01064     #else
01065       for (int i = 0; i < n; ++i) {
01066         a_i[i].sortOrder = i;
01067       }
01068     #endif
01069   }
01070 
01071   { 
01072     const double charge_scaling = sqrt(COULOMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1);
01073     const Vector ucenter = lattice.unscale(center);
01074     const BigReal ucenter_x = ucenter.x;
01075     const BigReal ucenter_y = ucenter.y;
01076     const BigReal ucenter_z = ucenter.z;
01077     const int n = numAtoms;
01078     #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
01079       int n_16 = n;
01080       n_16 = (n + 15) & (~15);
01081       cudaAtomList.resize(n_16);
01082       CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
01083       mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
01084       mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
01085       mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
01086       mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
01087     #else
01088     cudaAtomList.resize(n);
01089     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
01090     #endif
01091     const FullAtom *a = atom.begin();
01092     for ( int k=0; k<n; ++k ) {
01093       #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
01094         int j = a[k].sortOrder;
01095         atom_x[k] = a[j].position.x - ucenter_x;
01096         atom_y[k] = a[j].position.y - ucenter_y;
01097         atom_z[k] = a[j].position.z - ucenter_z;
01098         atom_q[k] = charge_scaling * a[j].charge;
01099       #else
01100       int j = a[k].sortOrder;
01101       ac[k].x = a[j].position.x - ucenter_x;
01102       ac[k].y = a[j].position.y - ucenter_y;
01103       ac[k].z = a[j].position.z - ucenter_z;
01104       ac[k].q = charge_scaling * a[j].charge;
01105       #endif
01106     }
01107   }
01108 #else
01109   doMigration = doMigration && numNeighbors;
01110 #endif
01111   doMigration = doMigration || ! patchMapRead;
01112 
01113   doMigration = doMigration || doAtomUpdate;
01114   doAtomUpdate = false;
01115 
01116   // Workaround for oversize groups
01117   doGroupSizeCheck();
01118 
01119   // Copy information needed by computes and proxys to Patch::p.
01120   p.resize(numAtoms);
01121   CompAtom *p_i = p.begin();
01122   pExt.resize(numAtoms);
01123   CompAtomExt *pExt_i = pExt.begin();
01124   FullAtom *a_i = atom.begin();
01125   int i; int n = numAtoms;
01126   for ( i=0; i<n; ++i ) { 
01127     p_i[i] = a_i[i]; 
01128     pExt_i[i] = a_i[i];
01129   }
01130 
01131   // Measure atom movement to test pairlist validity
01132   doPairlistCheck();
01133 
01134   if (flags.doMolly) mollyAverage();
01135   // BEGIN LA
01136   if (flags.doLoweAndersen) loweAndersenVelocities();
01137   // END LA
01138 
01139     if (flags.doGBIS) {
01140       //reset for next time step
01141       numGBISP1Arrived = 0;
01142       phase1BoxClosedCalled = false;
01143       numGBISP2Arrived = 0;
01144       phase2BoxClosedCalled = false;
01145       numGBISP3Arrived = 0;
01146       phase3BoxClosedCalled = false;
01147       if (doMigration || isNewProxyAdded)
01148         setGBISIntrinsicRadii();
01149     }
01150 
01151   if (flags.doLCPO) {
01152     if (doMigration || isNewProxyAdded) {
01153       setLcpoType();
01154     }
01155   }
01156 
01157   // Must Add Proxy Changes when migration completed!
01158   NodeIDList::iterator pli;
01159   int *pids = NULL;
01160   int pidsPreAllocated = 1;
01161   int npid;
01162   if (proxySendSpanning == 0) {
01163     npid = proxy.size();
01164     pids = new int[npid];
01165     pidsPreAllocated = 0;
01166     int *pidi = pids;
01167     int *pide = pids + proxy.size();
01168     int patchNodesLast =
01169       ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
01170     for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
01171     {
01172       if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
01173         *(--pide) = *pli;
01174       } else {
01175         *(pidi++) = *pli;
01176       }
01177     }
01178   }
01179   else {
01180 #ifdef NODEAWARE_PROXY_SPANNINGTREE
01181     #ifdef USE_NODEPATCHMGR
01182     npid = numNodeChild;
01183     pids = nodeChildren;
01184     #else
01185     npid = nChild;
01186     pids = child;
01187     #endif
01188 #else
01189     npid = nChild;
01190     pidsPreAllocated = 0;
01191     pids = new int[proxySpanDim];
01192     for (int i=0; i<nChild; i++) pids[i] = child[i];
01193 #endif
01194   }
01195   if (npid) { //have proxies
01196     int seq = flags.sequence;
01197     int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
01198     //begin to prepare proxy msg and send it
01199     int pdMsgPLLen = p.size();
01200     int pdMsgAvgPLLen = 0;
01201     if(flags.doMolly) {
01202         pdMsgAvgPLLen = p_avg.size();
01203     }
01204     // BEGIN LA
01205     int pdMsgVLLen = 0;
01206     if (flags.doLoweAndersen) {
01207         pdMsgVLLen = v.size();
01208     }
01209     // END LA
01210 
01211     int intRadLen = 0;
01212     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
01213             intRadLen = numAtoms * 2;
01214     }
01215 
01216     //LCPO
01217     int lcpoTypeLen = 0;
01218     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
01219             lcpoTypeLen = numAtoms;
01220     }
01221 
01222     int pdMsgPLExtLen = 0;
01223     if(doMigration || isNewProxyAdded) {
01224         pdMsgPLExtLen = pExt.size();
01225     }
01226 
01227     int cudaAtomLen = 0;
01228 #ifdef NAMD_CUDA
01229     cudaAtomLen = numAtoms;
01230 #endif
01231 
01232     #ifdef NAMD_MIC
01233       #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
01234         cudaAtomLen = numAtoms;
01235       #else
01236         cudaAtomLen = (numAtoms + 15) & (~15);
01237       #endif
01238     #endif
01239 
01240     ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
01241       lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
01242 
01243     SET_PRIORITY(nmsg,seq,priority);
01244     nmsg->patch = patchID;
01245     nmsg->flags = flags;
01246     nmsg->plLen = pdMsgPLLen;                
01247     //copying data to the newly created msg
01248     memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
01249     nmsg->avgPlLen = pdMsgAvgPLLen;        
01250     if(flags.doMolly) {
01251         memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
01252     }
01253     // BEGIN LA
01254     nmsg->vlLen = pdMsgVLLen;
01255     if (flags.doLoweAndersen) {
01256         memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
01257     }
01258     // END LA
01259 
01260     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
01261       for (int i = 0; i < numAtoms * 2; i++) {
01262         nmsg->intRadList[i] = intRad[i];
01263       }
01264     }
01265 
01266     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
01267       for (int i = 0; i < numAtoms; i++) {
01268         nmsg->lcpoTypeList[i] = lcpoType[i];
01269       }
01270     }
01271 
01272     nmsg->plExtLen = pdMsgPLExtLen;
01273     if(doMigration || isNewProxyAdded){     
01274         memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
01275     }
01276 
01277 // DMK
01278 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
01279     memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
01280 #endif
01281     
01282 #if NAMD_SeparateWaters != 0
01283     //DMK - Atom Separation (water vs. non-water)
01284     nmsg->numWaterAtoms = numWaterAtoms;
01285 #endif
01286 
01287 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
01288     nmsg->isFromImmMsgCall = 0;
01289 #endif
01290     
01291     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
01292     DebugFileTrace *dft = DebugFileTrace::Object();
01293     dft->openTrace();
01294     dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
01295     for(int i=0; i<npid; i++) {
01296         dft->writeTrace("%d ", pids[i]);
01297     }
01298     dft->writeTrace("\n");
01299     dft->closeTrace();
01300     #endif
01301 
01302 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
01303     if (isProxyChanged || localphs == NULL)
01304     {
01305 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
01306      //CmiAssert(isProxyChanged);
01307      if (nphs) {
01308        for (int i=0; i<nphs; i++) {
01309          CmiDestoryPersistent(localphs[i]);
01310        }
01311        delete [] localphs;
01312      }
01313      localphs = new PersistentHandle[npid];
01314      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;
01315      for (int i=0; i<npid; i++) {
01316 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
01317        if (proxySendSpanning)
01318            localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
01319        else
01320 #endif
01321        localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
01322      }
01323      nphs = npid;
01324     }
01325     CmiAssert(nphs == npid && localphs != NULL);
01326     CmiUsePersistentHandle(localphs, nphs);
01327 #endif
01328     if(doMigration || isNewProxyAdded) {
01329         ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
01330     }else{
01331         ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
01332     }
01333 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
01334     CmiUsePersistentHandle(NULL, 0);
01335 #endif
01336     isNewProxyAdded = 0;
01337   }
01338   isProxyChanged = 0;
01339   if(!pidsPreAllocated) delete [] pids;
01340   DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
01341 
01342 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
01343   positionPtrBegin = p.begin();
01344   positionPtrEnd = p.end();
01345 #endif
01346 
01347   if(flags.doMolly) {
01348       avgPositionPtrBegin = p_avg.begin();
01349       avgPositionPtrEnd = p_avg.end();
01350   }
01351   // BEGIN LA
01352   if (flags.doLoweAndersen) {
01353       velocityPtrBegin = v.begin();
01354       velocityPtrEnd = v.end();
01355   }
01356   // END LA
01357 
01358   Patch::positionsReady(doMigration);
01359 
01360   patchMapRead = 1;
01361 
01362   // gzheng
01363   Sync::Object()->PatchReady();
01364 }

void HomePatch::qmSwapAtoms (  ) 

Definition at line 956 of file HomePatch.C.

References ResizeArray< Elem >::begin(), CompAtom::charge, Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), if(), lssSubs(), Node::molecule, Node::Object(), Node::simParameters, simParams, and CompAtom::vdwType.

Referenced by positionsReady().

00957 {
00958     // This is used for LSS in QM/MM simulations.
00959     // Changes atom labels so that we effectively exchange solvent
00960     // molecules between classical and quantum modes.
00961     
00962     SimParameters *simParams = Node::Object()->simParameters;
00963     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
00964     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
00965     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
00966     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
00967     
00968     ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
00969             CkpvAccess(BOCclass_group).computeQMMgr) ;
00970     
00971     FullAtom *a_i = atom.begin();
00972     
00973     for (int i=0; i<numAtoms; ++i ) { 
00974         
00975         LSSSubsDat *subP = lssSubs(mgrP).find( LSSSubsDat(a_i[i].id) ) ;
00976         
00977         if ( subP != NULL ) {
00978             a_i[i].id = subP->newID ;
00979             a_i[i].vdwType = subP->newVdWType ;
00980             
00981             // If we are swappign a classical atom with a QM one, the charge
00982             // will need extra handling.
00983             if (qmAtomGroup[subP->newID] > 0 && simParams->PMEOn) {
00984                 // We make sure that the last atom charge calculated for the 
00985                 // QM atom being transfered here is available for PME
00986                 // in the next step.
00987                 
00988                 // Loops over all QM atoms (in all QM groups) comparing their 
00989                 // global indices (the sequential atom ID from NAMD).
00990                 for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
00991                     
00992                     if (qmAtmIndx[qmIter] == subP->newID) {
00993                         qmAtmChrg[qmIter] = subP->newCharge;
00994                         break;
00995                     }
00996                     
00997                 }
00998                 
00999                 // For QM atoms, the charge in the full atom structure is zero.
01000                 // Electrostatic interactions between QM atoms and their 
01001                 // environment is handled in ComputeQM.
01002                 a_i[i].charge = 0;
01003             }
01004             else {
01005                 // If we are swappign a QM atom with a Classical one, only the charge
01006                 // in the full atomstructure needs updating, since it used to be zero.
01007                 a_i[i].charge = subP->newCharge ;
01008             }
01009         }
01010     }
01011     
01012     return;
01013 }

int HomePatch::rattle1 ( const   BigReal,
Tensor virial,
SubmitReduction  
)

Definition at line 2273 of file HomePatch.C.

References addRattleForce(), buildRattleList(), endi(), iERROR(), iout, iWARN(), j, noconstList, Node::Object(), posNew, rattle1old(), rattleList, rattleListValid, rattleN(), rattlePair< 1 >(), rattleParam, settle1_SIMD< 1 >(), settle1_SIMD< 2 >(), settleList, Node::simParameters, simParams, TIMEFACTOR, velNew, WAT_TIP3, x, Vector::x, Vector::y, and Vector::z.

Referenced by Sequencer::minimize(), Sequencer::minimizeMoveDownhill(), and Sequencer::rattle1().

02274                                 {
02275 
02276   SimParameters *simParams = Node::Object()->simParameters;
02277   if (simParams->watmodel != WAT_TIP3 || ppreduction) {
02278     // Call old rattle1 -method instead
02279     return rattle1old(timestep, virial, ppreduction);
02280   }
02281 
02282   if (!rattleListValid) {
02283     buildRattleList();
02284     rattleListValid = true;
02285   }
02286 
02287   const int fixedAtomsOn = simParams->fixedAtomsOn;
02288   const int useSettle = simParams->useSettle;
02289   const BigReal dt = timestep / TIMEFACTOR;
02290   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02291   const BigReal tol2 = 2.0 * simParams->rigidTol;
02292   int maxiter = simParams->rigidIter;
02293   int dieOnError = simParams->rigidDie;
02294 
02295   Vector ref[10];  // reference position
02296   Vector pos[10];  // new position
02297   Vector vel[10];  // new velocity
02298 
02299   // Manual un-roll
02300   int n = (settleList.size()/2)*2;
02301   for (int j=0;j < n;j+=2) {
02302     int ig;
02303     ig = settleList[j];
02304     for (int i = 0; i < 3; ++i ) {
02305       ref[i] = atom[ig+i].position;
02306       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02307     }
02308     ig = settleList[j+1];
02309     for (int i = 0; i < 3; ++i ) {
02310       ref[i+3] = atom[ig+i].position;
02311       pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
02312     }
02313     settle1_SIMD<2>(ref, pos,
02314       settle_mOrmT, settle_mHrmT, settle_ra,
02315       settle_rb, settle_rc, settle_rra);
02316 
02317     ig = settleList[j];
02318     for (int i = 0; i < 3; ++i ) {
02319       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02320       posNew[ig+i] = pos[i];
02321     }
02322     ig = settleList[j+1];
02323     for (int i = 0; i < 3; ++i ) {
02324       velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
02325       posNew[ig+i] = pos[i+3];
02326     }
02327 
02328   }
02329 
02330   if (settleList.size() % 2) {
02331     int ig = settleList[settleList.size()-1];
02332     for (int i = 0; i < 3; ++i ) {
02333       ref[i] = atom[ig+i].position;
02334       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02335     }
02336     settle1_SIMD<1>(ref, pos,
02337             settle_mOrmT, settle_mHrmT, settle_ra,
02338             settle_rb, settle_rc, settle_rra);        
02339     for (int i = 0; i < 3; ++i ) {
02340       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02341       posNew[ig+i] = pos[i];
02342     }
02343   }
02344 
02345   int posParam = 0;
02346   for (int j=0;j < rattleList.size();++j) {
02347 
02348     BigReal refx[10];
02349     BigReal refy[10];
02350     BigReal refz[10];
02351 
02352     BigReal posx[10];
02353     BigReal posy[10];
02354     BigReal posz[10];
02355 
02356     int ig = rattleList[j].ig;
02357     int icnt = rattleList[j].icnt;
02358     int hgs = atom[ig].hydrogenGroupSize;
02359     for (int i = 0; i < hgs; ++i ) {
02360       ref[i] = atom[ig+i].position;
02361       pos[i] = atom[ig+i].position;
02362       if (!(fixedAtomsOn && atom[ig+i].atomFixed))
02363         pos[i] += atom[ig+i].velocity * dt;
02364       refx[i] = ref[i].x;
02365       refy[i] = ref[i].y;
02366       refz[i] = ref[i].z;
02367       posx[i] = pos[i].x;
02368       posy[i] = pos[i].y;
02369       posz[i] = pos[i].z;
02370     }
02371 
02372 
02373     bool done;
02374     bool consFailure;
02375     if (icnt == 1) {
02376       rattlePair<1>(&rattleParam[posParam],
02377         refx, refy, refz,
02378         posx, posy, posz);
02379       done = true;
02380       consFailure = false;
02381     } else {
02382       rattleN(icnt, &rattleParam[posParam],
02383         refx, refy, refz,
02384         posx, posy, posz,
02385         tol2, maxiter,
02386         done, consFailure);
02387     }
02388 
02389     // Advance position in rattleParam
02390     posParam += icnt;
02391 
02392     for (int i = 0; i < hgs; ++i ) {
02393       pos[i].x = posx[i];
02394       pos[i].y = posy[i];
02395       pos[i].z = posz[i];
02396     }
02397 
02398     for (int i = 0; i < hgs; ++i ) {
02399       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02400       posNew[ig+i] = pos[i];
02401     }
02402 
02403     if ( consFailure ) {
02404       if ( dieOnError ) {
02405         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02406         << (atom[ig].id + 1) << "!\n" << endi;
02407         return -1;  // triggers early exit
02408       } else {
02409         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02410         << (atom[ig].id + 1) << "!\n" << endi;
02411       }
02412     } else if ( ! done ) {
02413       if ( dieOnError ) {
02414         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02415         << (atom[ig].id + 1) << "!\n" << endi;
02416         return -1;  // triggers early exit
02417       } else {
02418         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02419         << (atom[ig].id + 1) << "!\n" << endi;
02420       }
02421     }
02422 
02423   }
02424 
02425   // Finally, we have to go through atoms that are not involved in rattle just so that we have
02426   // their positions and velocities up-to-date in posNew and velNew
02427   for (int j=0;j < noconstList.size();++j) {
02428     int ig = noconstList[j];
02429     int hgs = atom[ig].hydrogenGroupSize;
02430     for (int i = 0; i < hgs; ++i ) {
02431       velNew[ig+i] = atom[ig+i].velocity;
02432       posNew[ig+i] = atom[ig+i].position;
02433     }
02434   }
02435 
02436   if ( invdt == 0 ) {
02437     for (int ig = 0; ig < numAtoms; ++ig )
02438       atom[ig].position = posNew[ig];
02439   } else if ( virial == 0 ) {
02440     for (int ig = 0; ig < numAtoms; ++ig )
02441       atom[ig].velocity = velNew[ig];
02442   } else {
02443     Tensor wc;  // constraint virial
02444     addRattleForce(invdt, wc);
02445     *virial += wc;
02446   }
02447 
02448   return 0;
02449 }

int HomePatch::rattle1old ( const   BigReal,
Tensor virial,
SubmitReduction  
)

Definition at line 2452 of file HomePatch.C.

References Lattice::c(), endi(), Patch::f, iERROR(), iout, SubmitReduction::item(), iWARN(), Patch::lattice, Node::molecule, NAMD_die(), Results::normal, Node::Object(), Lattice::origin(), outer(), partition(), settle1(), settle1init(), Node::simParameters, simParams, TIMEFACTOR, WAT_SWM4, WAT_TIP4, Vector::x, Tensor::xx, Vector::y, Tensor::yy, z, Vector::z, and Tensor::zz.

Referenced by rattle1().

02454 {
02455   Molecule *mol = Node::Object()->molecule;
02456   SimParameters *simParams = Node::Object()->simParameters;
02457   const int fixedAtomsOn = simParams->fixedAtomsOn;
02458   const int useSettle = simParams->useSettle;
02459   const BigReal dt = timestep / TIMEFACTOR;
02460   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02461   BigReal tol2 = 2.0 * simParams->rigidTol;
02462   int maxiter = simParams->rigidIter;
02463   int dieOnError = simParams->rigidDie;
02464   int i, iter;
02465   BigReal dsq[10], tmp;
02466   int ial[10], ibl[10];
02467   Vector ref[10];  // reference position
02468   Vector refab[10];  // reference vector
02469   Vector pos[10];  // new position
02470   Vector vel[10];  // new velocity
02471   Vector netdp[10];  // total momentum change from constraint
02472   BigReal rmass[10];  // 1 / mass
02473   int fixed[10];  // is atom fixed?
02474   Tensor wc;  // constraint virial
02475   BigReal idz, zmin;
02476   int nslabs;
02477 
02478   // Initialize the settle algorithm with water parameters
02479   // settle1() assumes all waters are identical,
02480   // and will generate bad results if they are not.
02481   // XXX this will move to Molecule::build_atom_status when that 
02482   // version is debugged
02483   if ( ! settle_initialized ) {
02484     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02485       // find a water
02486       if (atom[ig].rigidBondLength > 0) {
02487         int oatm;
02488         if (simParams->watmodel == WAT_SWM4) {
02489           oatm = ig+3;  // skip over Drude and Lonepair
02490           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02491           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02492         }
02493         else {
02494           oatm = ig+1;
02495           // Avoid using the Om site to set this by mistake
02496           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02497             oatm += 1;
02498           }
02499         }
02500 
02501         // initialize settle water parameters
02502         settle1init(atom[ig].mass, atom[oatm].mass, 
02503                     atom[ig].rigidBondLength,
02504                     atom[oatm].rigidBondLength,
02505                     settle_mOrmT, settle_mHrmT, settle_ra,
02506                     settle_rb, settle_rc, settle_rra);
02507         settle_initialized = 1;
02508         break; // done with init
02509       }
02510     }
02511   }
02512 
02513   if (ppreduction) {
02514     nslabs = simParams->pressureProfileSlabs;
02515     idz = nslabs/lattice.c().z;
02516     zmin = lattice.origin().z - 0.5*lattice.c().z;
02517   }
02518 
02519   // Size of a hydrogen group for water
02520   int wathgsize = 3;
02521   int watmodel = simParams->watmodel;
02522   if (watmodel == WAT_TIP4) wathgsize = 4;
02523   else if (watmodel == WAT_SWM4) wathgsize = 5;
02524   
02525   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02526     int hgs = atom[ig].hydrogenGroupSize;
02527     if ( hgs == 1 ) continue;  // only one atom in group
02528     // cache data in local arrays and integrate positions normally
02529     int anyfixed = 0;
02530     for ( i = 0; i < hgs; ++i ) {
02531       ref[i] = atom[ig+i].position;
02532       pos[i] = atom[ig+i].position;
02533       vel[i] = atom[ig+i].velocity;
02534       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02535       //printf("rmass of %i is %f\n", ig+i, rmass[i]);
02536       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02537       //printf("fixed status of %i is %i\n", i, fixed[i]);
02538       // undo addVelocityToPosition to get proper reference coordinates
02539       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
02540     }
02541     int icnt = 0;
02542     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02543       if (hgs != wathgsize) {
02544         char errmsg[256];
02545         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02546                          "but the specified water model requires %d atoms.\n",
02547                          atom[ig].id+1, hgs, wathgsize);
02548         NAMD_die(errmsg);
02549       }
02550       // Use SETTLE for water unless some of the water atoms are fixed,
02551       if (useSettle && !anyfixed) {
02552         if (simParams->watmodel == WAT_SWM4) {
02553           // SWM4 ordering:  O D LP H1 H2
02554           // do swap(O,LP) and call settle with subarray O H1 H2
02555           // swap back after we return
02556           Vector lp_ref = ref[2];
02557           Vector lp_pos = pos[2];
02558           Vector lp_vel = vel[2];
02559           ref[2] = ref[0];
02560           pos[2] = pos[0];
02561           vel[2] = vel[0];
02562           settle1(ref+2, pos+2, vel+2, invdt,
02563                   settle_mOrmT, settle_mHrmT, settle_ra,
02564                   settle_rb, settle_rc, settle_rra);
02565           ref[0] = ref[2];
02566           pos[0] = pos[2];
02567           vel[0] = vel[2];
02568           ref[2] = lp_ref;
02569           pos[2] = lp_pos;
02570           vel[2] = lp_vel;
02571           // determine for LP updated pos and vel
02572           swm4_omrepos(ref, pos, vel, invdt);
02573         }
02574         else {
02575             settle1(ref, pos, vel, invdt,
02576                     settle_mOrmT, settle_mHrmT, settle_ra,
02577                     settle_rb, settle_rc, settle_rra);            
02578           if (simParams->watmodel == WAT_TIP4) {
02579             tip4_omrepos(ref, pos, vel, invdt);
02580           }
02581         }
02582 
02583         // which slab the hydrogen group will belong to
02584         // for pprofile calculations.
02585         int ppoffset, partition;
02586         if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02587           atom[ig+i].position = pos[i];
02588         } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02589           atom[ig+i].velocity = vel[i];
02590         } else for ( i = 0; i < wathgsize; ++i ) {
02591           Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
02592           Tensor vir = outer(df, ref[i]);
02593           wc += vir;
02594           f[Results::normal][ig+i] += df;
02595           atom[ig+i].velocity = vel[i];
02596           if (ppreduction) {
02597             // put all the atoms from a water in the same slab.  Atom 0
02598             // should be the parent atom.
02599             if (!i) {
02600               BigReal z = pos[i].z;
02601               partition = atom[ig].partition;
02602               int slab = (int)floor((z-zmin)*idz);
02603               if (slab < 0) slab += nslabs;
02604               else if (slab >= nslabs) slab -= nslabs;
02605               ppoffset = 3*(slab + nslabs*partition);
02606             }
02607             ppreduction->item(ppoffset  ) += vir.xx;
02608             ppreduction->item(ppoffset+1) += vir.yy;
02609             ppreduction->item(ppoffset+2) += vir.zz;
02610           }
02611         }
02612         continue;
02613       }
02614       if ( !(fixed[1] && fixed[2]) ) {
02615         dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02616       }
02617     }
02618     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02619       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02620         if ( !(fixed[0] && fixed[i]) ) {
02621           dsq[icnt] = tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
02622         }
02623       }
02624     }
02625     if ( icnt == 0 ) continue;  // no constraints
02626     for ( i = 0; i < icnt; ++i ) {
02627       refab[i] = ref[ial[i]] - ref[ibl[i]];
02628     }
02629     for ( i = 0; i < hgs; ++i ) {
02630       netdp[i] = 0.;
02631     }
02632     int done;
02633     int consFailure;
02634     for ( iter = 0; iter < maxiter; ++iter ) {
02635 //if (iter > 0) CkPrintf("iteration %d\n", iter);
02636       done = 1;
02637       consFailure = 0;
02638       for ( i = 0; i < icnt; ++i ) {
02639         int a = ial[i];  int b = ibl[i];
02640         Vector pab = pos[a] - pos[b];
02641               BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
02642         BigReal rabsq = dsq[i];
02643         BigReal diffsq = rabsq - pabsq;
02644         if ( fabs(diffsq) > (rabsq * tol2) ) {
02645           Vector &rab = refab[i];
02646           BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
02647           if ( rpab < ( rabsq * 1.0e-6 ) ) {
02648             done = 0;
02649             consFailure = 1;
02650             continue;
02651           }
02652           BigReal rma = rmass[a];
02653           BigReal rmb = rmass[b];
02654           BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
02655           Vector dp = rab * gab;
02656           pos[a] += rma * dp;
02657           pos[b] -= rmb * dp;
02658           if ( invdt != 0. ) {
02659             dp *= invdt;
02660             netdp[a] += dp;
02661             netdp[b] -= dp;
02662           }
02663           done = 0;
02664         }
02665       }
02666       if ( done ) break;
02667     }
02668 
02669     if ( consFailure ) {
02670       if ( dieOnError ) {
02671         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02672         << (atom[ig].id + 1) << "!\n" << endi;
02673         return -1;  // triggers early exit
02674       } else {
02675         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02676         << (atom[ig].id + 1) << "!\n" << endi;
02677       }
02678     } else if ( ! done ) {
02679       if ( dieOnError ) {
02680         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02681         << (atom[ig].id + 1) << "!\n" << endi;
02682         return -1;  // triggers early exit
02683       } else {
02684         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02685         << (atom[ig].id + 1) << "!\n" << endi;
02686       }
02687     }
02688 
02689     // store data back to patch
02690     int ppoffset, partition;
02691     if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
02692       atom[ig+i].position = pos[i];
02693     } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
02694       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02695     } else for ( i = 0; i < hgs; ++i ) {
02696       Force df = netdp[i] * invdt;
02697       Tensor vir = outer(df, ref[i]);
02698       wc += vir;
02699       f[Results::normal][ig+i] += df;
02700       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02701       if (ppreduction) {
02702         if (!i) {
02703           BigReal z = pos[i].z;
02704           int partition = atom[ig].partition;
02705           int slab = (int)floor((z-zmin)*idz);
02706           if (slab < 0) slab += nslabs;
02707           else if (slab >= nslabs) slab -= nslabs;
02708           ppoffset = 3*(slab + nslabs*partition);
02709         }
02710         ppreduction->item(ppoffset  ) += vir.xx;
02711         ppreduction->item(ppoffset+1) += vir.yy;
02712         ppreduction->item(ppoffset+2) += vir.zz;
02713       }
02714     }
02715   }
02716   if ( dt && virial ) *virial += wc;
02717 
02718   return 0;
02719 }

void HomePatch::rattle2 ( const   BigReal,
Tensor virial 
)

Definition at line 2722 of file HomePatch.C.

References endi(), iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Node::Object(), outer(), settle2(), Node::simParameters, simParams, TIMEFACTOR, WAT_SWM4, WAT_TIP4, Vector::x, Vector::y, and Vector::z.

02723 {
02724   Molecule *mol = Node::Object()->molecule;
02725   SimParameters *simParams = Node::Object()->simParameters;
02726   const int fixedAtomsOn = simParams->fixedAtomsOn;
02727   const int useSettle = simParams->useSettle;
02728   const BigReal dt = timestep / TIMEFACTOR;
02729   Tensor wc;  // constraint virial
02730   BigReal tol = simParams->rigidTol;
02731   int maxiter = simParams->rigidIter;
02732   int dieOnError = simParams->rigidDie;
02733   int i, iter;
02734   BigReal dsqi[10], tmp;
02735   int ial[10], ibl[10];
02736   Vector ref[10];  // reference position
02737   Vector refab[10];  // reference vector
02738   Vector vel[10];  // new velocity
02739   BigReal rmass[10];  // 1 / mass
02740   BigReal redmass[10];  // reduced mass
02741   int fixed[10];  // is atom fixed?
02742   
02743   // Size of a hydrogen group for water
02744   int wathgsize = 3;
02745   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02746   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02747 
02748   //  CkPrintf("In rattle2!\n");
02749   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02750     //    CkPrintf("ig=%d\n",ig);
02751     int hgs = atom[ig].hydrogenGroupSize;
02752     if ( hgs == 1 ) continue;  // only one atom in group
02753     // cache data in local arrays and integrate positions normally
02754     int anyfixed = 0;
02755     for ( i = 0; i < hgs; ++i ) {
02756       ref[i] = atom[ig+i].position;
02757       vel[i] = atom[ig+i].velocity;
02758       rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
02759       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02760       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02761     }
02762     int icnt = 0;
02763     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02764       if (hgs != wathgsize) {
02765         NAMD_bug("Hydrogen group error caught in rattle2().");
02766       }
02767       // Use SETTLE for water unless some of the water atoms are fixed,
02768       if (useSettle && !anyfixed) {
02769         if (simParams->watmodel == WAT_SWM4) {
02770           // SWM4 ordering:  O D LP H1 H2
02771           // do swap(O,LP) and call settle with subarray O H1 H2
02772           // swap back after we return
02773           Vector lp_ref = ref[2];
02774           // Vector lp_vel = vel[2];
02775           ref[2] = ref[0];
02776           vel[2] = vel[0];
02777           settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
02778           ref[0] = ref[2];
02779           vel[0] = vel[2];
02780           ref[2] = lp_ref;
02781           // vel[2] = vel[0];  // set LP vel to O vel
02782         }
02783         else {
02784           settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
02785           if (simParams->watmodel == WAT_TIP4) {
02786             vel[3] = vel[0];
02787           }
02788         }
02789         for (i=0; i<hgs; i++) {
02790           atom[ig+i].velocity = vel[i];
02791         }
02792         continue;
02793       }
02794       if ( !(fixed[1] && fixed[2]) ) {
02795         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02796         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02797       }
02798     }
02799     //    CkPrintf("Loop 2\n");
02800     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02801       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02802         if ( !(fixed[0] && fixed[i]) ) {
02803           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02804           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02805           ibl[icnt] = i;  ++icnt;
02806         }
02807       }
02808     }
02809     if ( icnt == 0 ) continue;  // no constraints
02810     //    CkPrintf("Loop 3\n");
02811     for ( i = 0; i < icnt; ++i ) {
02812       refab[i] = ref[ial[i]] - ref[ibl[i]];
02813     }
02814     //    CkPrintf("Loop 4\n");
02815     int done;
02816     for ( iter = 0; iter < maxiter; ++iter ) {
02817       done = 1;
02818       for ( i = 0; i < icnt; ++i ) {
02819         int a = ial[i];  int b = ibl[i];
02820         Vector vab = vel[a] - vel[b];
02821         Vector &rab = refab[i];
02822         BigReal rabsqi = dsqi[i];
02823         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02824         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02825           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02826           wc += outer(dp,rab);
02827           vel[a] += rmass[a] * dp;
02828           vel[b] -= rmass[b] * dp;
02829           done = 0;
02830         }
02831       }
02832       if ( done ) break;
02833       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02834     }
02835     if ( ! done ) {
02836       if ( dieOnError ) {
02837         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02838       } else {
02839         iout << iWARN <<
02840           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02841       }
02842     }
02843     // store data back to patch
02844     for ( i = 0; i < hgs; ++i ) {
02845       atom[ig+i].velocity = vel[i];
02846     }
02847   }
02848   //  CkPrintf("Leaving rattle2!\n");
02849   // check that there isn't a constant needed here!
02850   *virial += wc / ( 0.5 * dt );
02851 
02852 }

void HomePatch::receiveResult ( ProxyGBISP2ResultMsg msg  ) 

Definition at line 3180 of file HomePatch.C.

References Sequencer::awaken(), ProxyGBISP2ResultMsg::dEdaSum, ProxyGBISP2ResultMsg::dEdaSumLen, Patch::dHdrPrefix, Flags::doNonbonded, Patch::flags, if(), and ResizeArray< Elem >::size().

03180                                                        {
03181   ++numGBISP2Arrived;
03182   //accumulate dEda
03183   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
03184     dHdrPrefix[i] += msg->dEdaSum[i];
03185   }
03186   delete msg;
03187 
03188   if (flags.doNonbonded) {
03189     //awaken if phase 2 done
03190     if (phase2BoxClosedCalled == true &&
03191         numGBISP2Arrived==proxy.size() ) {
03192       sequencer->awaken();
03193     }
03194   } else {
03195     //awaken if all phases done on noWork step
03196     if (boxesOpen == 0 &&
03197         numGBISP1Arrived == proxy.size() &&
03198         numGBISP2Arrived == proxy.size() &&
03199         numGBISP3Arrived == proxy.size()) {
03200       sequencer->awaken();
03201     }
03202   }
03203 }

void HomePatch::receiveResult ( ProxyGBISP1ResultMsg msg  ) 

Definition at line 3155 of file HomePatch.C.

References Sequencer::awaken(), Flags::doNonbonded, Patch::flags, if(), Patch::psiFin, ProxyGBISP1ResultMsg::psiSum, ProxyGBISP1ResultMsg::psiSumLen, and ResizeArray< Elem >::size().

Referenced by ProxyMgr::recvResult().

03155                                                        {
03156   ++numGBISP1Arrived;
03157     for ( int i = 0; i < msg->psiSumLen; ++i ) {
03158       psiFin[i] += msg->psiSum[i];
03159     }
03160   delete msg;
03161 
03162   if (flags.doNonbonded) {
03163     //awaken if phase 1 done
03164     if (phase1BoxClosedCalled == true &&
03165         numGBISP1Arrived==proxy.size() ) {
03166         sequencer->awaken();
03167     }
03168   } else {
03169     //awaken if all phases done on noWork step
03170     if (boxesOpen == 0 &&
03171         numGBISP1Arrived == proxy.size() &&
03172         numGBISP2Arrived == proxy.size() &&
03173         numGBISP3Arrived == proxy.size()) {
03174       sequencer->awaken();
03175     }
03176   }
03177 }

void HomePatch::receiveResults ( ProxyCombinedResultRawMsg msg  ) 

Definition at line 927 of file HomePatch.C.

References OwnerBox< Owner, Data >::clientClose(), OwnerBox< Owner, Data >::clientOpen(), DebugM, Results::f, Patch::f, ProxyCombinedResultRawMsg::flLen, ProxyCombinedResultRawMsg::forceArr, Patch::forceBox, ProxyCombinedResultRawMsg::isForceNonZero, Results::maxNumForces, ProxyCombinedResultRawMsg::nodeSize, Patch::patchID, and x.

00928 {
00929     numGBISP3Arrived++;
00930   DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodeSize<<")\n");
00931     Results *r = forceBox.clientOpen(msg->nodeSize);
00932       register char* isNonZero = msg->isForceNonZero;
00933       register Force* f_i = msg->forceArr;
00934       for ( int k = 0; k < Results::maxNumForces; ++k )
00935       {
00936         Force *f = r->f[k];
00937                 int nf = msg->flLen[k];
00938 #ifdef ARCH_POWERPC
00939 #pragma disjoint (*f_i, *f)
00940 #endif
00941         for (int count = 0; count < nf; count++) {
00942           if(*isNonZero){
00943                         f[count].x += f_i->x;
00944                         f[count].y += f_i->y;
00945                         f[count].z += f_i->z;
00946                         f_i++;
00947           }
00948           isNonZero++;
00949         }
00950       }
00951     forceBox.clientClose(msg->nodeSize);
00952 
00953     delete msg;
00954 }

void HomePatch::receiveResults ( ProxyResultMsg msg  ) 

Definition at line 910 of file HomePatch.C.

References ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::clientClose(), OwnerBox< Owner, Data >::clientOpen(), DebugM, ResizeArray< Elem >::end(), Results::f, Patch::f, Patch::forceBox, ProxyResultMsg::forceList, Results::maxNumForces, ProxyResultMsg::node, and Patch::patchID.

00910                                                   {
00911   numGBISP3Arrived++;
00912   DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00913   int n = msg->node;
00914   Results *r = forceBox.clientOpen();
00915   for ( int k = 0; k < Results::maxNumForces; ++k )
00916   {
00917     Force *f = r->f[k];
00918     register ForceList::iterator f_i, f_e;
00919     f_i = msg->forceList[k]->begin();
00920     f_e = msg->forceList[k]->end();
00921     for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
00922   }
00923   forceBox.clientClose();
00924   delete msg;
00925 }

void HomePatch::receiveResults ( ProxyResultVarsizeMsg msg  ) 

Definition at line 886 of file HomePatch.C.

References OwnerBox< Owner, Data >::clientOpen(), DebugM, Results::f, ProxyResultVarsizeMsg::forceArr, Patch::forceBox, ProxyResultVarsizeMsg::isZero, Results::maxNumForces, ProxyResultVarsizeMsg::node, and Patch::patchID.

Referenced by ProxyMgr::recvResults().

00886                                                          {
00887 
00888     numGBISP3Arrived++;
00889     DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00890     int n = msg->node;
00891     Results *r = forceBox.clientOpen();
00892 
00893     char *iszeroPtr = msg->isZero;
00894     Force *msgFPtr = msg->forceArr;
00895 
00896     for ( int k = 0; k < Results::maxNumForces; ++k )
00897     {
00898       Force *rfPtr = r->f[k];
00899       for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
00900           if((*iszeroPtr)!=1) {
00901               *rfPtr += *msgFPtr;
00902               msgFPtr++;
00903           }
00904       }      
00905     }
00906     forceBox.clientClose();
00907     delete msg;
00908 }

void HomePatch::recvCheckpointAck (  ) 

Definition at line 3462 of file HomePatch.C.

Referenced by PatchMgr::recvCheckpointAck(), and recvCheckpointLoad().

03462                                   {  // initiating replica
03463   CkpvAccess(_qd)->process();
03464 }

void HomePatch::recvCheckpointLoad ( CheckpointAtomsMsg msg  ) 

Definition at line 3409 of file HomePatch.C.

References Patch::atomMapper, CheckpointAtomsMsg::atoms, ResizeArray< Elem >::begin(), Sequencer::berendsenPressure_count, CheckpointAtomsMsg::berendsenPressure_count, checkpoint_task, ResizeArray< Elem >::end(), CheckpointAtomsMsg::key, Patch::lattice, CheckpointAtomsMsg::lattice, NAMD_bug(), CheckpointAtomsMsg::numAtoms, PatchMgr::Object(), Node::Object(), Patch::patchID, CheckpointAtomsMsg::pe, CheckpointAtomsMsg::pid, rattleListValid, recvCheckpointAck(), AtomMapper::registerIDsFullAtom(), CheckpointAtomsMsg::replica, ResizeArray< Elem >::resize(), SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_STORE, SCRIPT_CHECKPOINT_SWAP, PatchMgr::sendCheckpointStore(), Node::simParameters, simParams, and AtomMapper::unregisterIDsFullAtom().

Referenced by PatchMgr::recvCheckpointLoad().

03409                                                           {  // initiating replica
03410   SimParameters *simParams = Node::Object()->simParameters;
03411   const int remote = simParams->scriptIntArg1;
03412   const char *key = simParams->scriptStringArg1;
03413   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
03414     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
03415   }
03416   if ( msg->replica != remote ) {
03417     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
03418   }
03419   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03420     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
03421     strcpy(newmsg->key,key);
03422     newmsg->lattice = lattice;
03423     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
03424     newmsg->pid = patchID;
03425     newmsg->pe = CkMyPe();
03426     newmsg->replica = CmiMyPartition();
03427     newmsg->numAtoms = numAtoms;
03428     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03429     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
03430   }
03431   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03432     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03433     lattice = msg->lattice;
03434     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
03435     numAtoms = msg->numAtoms;
03436     atom.resize(numAtoms);
03437     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03438     doAtomUpdate = true;
03439     rattleListValid = false;
03440     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03441   }
03442   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
03443     recvCheckpointAck();
03444   }
03445   delete msg;
03446 }

void HomePatch::recvCheckpointReq ( int  task,
const char *  key,
int  replica,
int  pe 
)

Definition at line 3379 of file HomePatch.C.

References HomePatch::checkpoint_t::atoms, CheckpointAtomsMsg::atoms, ResizeArray< Elem >::begin(), HomePatch::checkpoint_t::berendsenPressure_count, CheckpointAtomsMsg::berendsenPressure_count, checkpoints, HomePatch::checkpoint_t::lattice, CheckpointAtomsMsg::lattice, NAMD_die(), CheckpointAtomsMsg::numAtoms, HomePatch::checkpoint_t::numAtoms, PatchMgr::Object(), Patch::patchID, CheckpointAtomsMsg::pe, CheckpointAtomsMsg::pid, CheckpointAtomsMsg::replica, SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_SWAP, PatchMgr::sendCheckpointAck(), and PatchMgr::sendCheckpointLoad().

Referenced by PatchMgr::recvCheckpointReq().

03379                                                                                 {  // responding replica
03380   if ( task == SCRIPT_CHECKPOINT_FREE ) {
03381     if ( ! checkpoints.count(key) ) {
03382       NAMD_die("Unable to free checkpoint, requested key was never stored.");
03383     }
03384     delete checkpoints[key];
03385     checkpoints.erase(key);
03386     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
03387     return;
03388   }
03389   CheckpointAtomsMsg *msg;
03390   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
03391     if ( ! checkpoints.count(key) ) {
03392       NAMD_die("Unable to load checkpoint, requested key was never stored.");
03393     }
03394     checkpoint_t &cp = *checkpoints[key];
03395     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
03396     msg->lattice = cp.lattice;
03397     msg->berendsenPressure_count = cp.berendsenPressure_count;
03398     msg->numAtoms = cp.numAtoms;
03399     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
03400   } else {
03401     msg = new (0,1,0) CheckpointAtomsMsg;
03402   }
03403   msg->pid = patchID;
03404   msg->replica = CmiMyPartition();
03405   msg->pe = CkMyPe();
03406   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
03407 }

void HomePatch::recvCheckpointStore ( CheckpointAtomsMsg msg  ) 

Definition at line 3448 of file HomePatch.C.

References CheckpointAtomsMsg::atoms, HomePatch::checkpoint_t::atoms, ResizeArray< Elem >::begin(), CheckpointAtomsMsg::berendsenPressure_count, HomePatch::checkpoint_t::berendsenPressure_count, checkpoints, CheckpointAtomsMsg::key, CheckpointAtomsMsg::lattice, HomePatch::checkpoint_t::lattice, CheckpointAtomsMsg::numAtoms, HomePatch::checkpoint_t::numAtoms, PatchMgr::Object(), Patch::patchID, CheckpointAtomsMsg::pe, CheckpointAtomsMsg::replica, ResizeArray< Elem >::resize(), and PatchMgr::sendCheckpointAck().

Referenced by PatchMgr::recvCheckpointStore().

03448                                                            {  // responding replica
03449   if ( ! checkpoints.count(msg->key) ) {
03450     checkpoints[msg->key] = new checkpoint_t;
03451   }
03452   checkpoint_t &cp = *checkpoints[msg->key];
03453   cp.lattice = msg->lattice;
03454   cp.berendsenPressure_count = msg->berendsenPressure_count;
03455   cp.numAtoms = msg->numAtoms;
03456   cp.atoms.resize(cp.numAtoms);
03457   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
03458   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
03459   delete msg;
03460 }

void HomePatch::recvExchangeMsg ( ExchangeAtomsMsg msg  ) 

Definition at line 3499 of file HomePatch.C.

References Patch::atomMapper, ExchangeAtomsMsg::atoms, ResizeArray< Elem >::begin(), ResizeArray< Elem >::end(), ExchangeAtomsMsg::lattice, Patch::lattice, ExchangeAtomsMsg::numAtoms, rattleListValid, AtomMapper::registerIDsFullAtom(), ResizeArray< Elem >::resize(), and AtomMapper::unregisterIDsFullAtom().

Referenced by PatchMgr::recvExchangeMsg().

03499                                                      {
03500   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
03501   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03502   lattice = msg->lattice;
03503   numAtoms = msg->numAtoms;
03504   atom.resize(numAtoms);
03505   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03506   delete msg;
03507   CkpvAccess(_qd)->process();
03508   doAtomUpdate = true;
03509   rattleListValid = false;
03510   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03511 }

void HomePatch::recvExchangeReq ( int  req  ) 

Definition at line 3488 of file HomePatch.C.

References exchange_dst, exchange_msg, exchange_req, PatchMgr::Object(), and PatchMgr::sendExchangeMsg().

Referenced by exchangeAtoms(), and PatchMgr::recvExchangeReq().

03488                                        {
03489   exchange_req = req;
03490   if ( exchange_msg ) {
03491     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
03492     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
03493     exchange_msg = 0;
03494     exchange_req = -1;
03495     CkpvAccess(_qd)->process();
03496   }
03497 }

void HomePatch::recvNodeAwareSpanningTree ( ProxyNodeAwareSpanningTreeMsg msg  ) 

Definition at line 720 of file HomePatch.C.

Referenced by ProxyMgr::recvNodeAwareSpanningTreeOnHomePatch().

00720 {}

void HomePatch::recvSpanningTree ( int *  t,
int  n 
)

Definition at line 726 of file HomePatch.C.

References proxySpanDim, ResizeArray< Elem >::resize(), sendSpanningTree(), and ResizeArray< Elem >::size().

Referenced by ProxyMgr::recvSpanningTreeOnHomePatch().

00727 {
00728   int i;
00729   nChild=0;
00730   tree.resize(n);
00731   for (i=0; i<n; i++) {
00732     tree[i] = t[i];
00733   }
00734 
00735   for (i=1; i<=proxySpanDim; i++) {
00736     if (tree.size() <= i) break;
00737     child[i-1] = tree[i];
00738     nChild++;
00739   }
00740 
00741 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00742   isProxyChanged = 1;
00743 #endif
00744 
00745   // send down to children
00746   sendSpanningTree();
00747 }

void HomePatch::registerProxy ( RegisterProxyMsg  ) 

Definition at line 492 of file HomePatch.C.

References ResizeArray< Elem >::add(), ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::clientAdd(), DebugM, Patch::forceBox, RegisterProxyMsg::node, Patch::patchID, and ResizeArray< Elem >::size().

Referenced by ProxyMgr::recvRegisterProxy().

00492                                                    {
00493   DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
00494   proxy.add(msg->node);
00495   forceBox.clientAdd();
00496 
00497   isNewProxyAdded = 1;
00498 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00499   isProxyChanged = 1;
00500 #endif
00501 
00502   Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
00503   delete msg;
00504 }

void HomePatch::replaceForces ( ExtForce f  ) 

Definition at line 1366 of file HomePatch.C.

References Patch::f.

01367 {
01368   replacementForces = f;
01369 }

void HomePatch::revert ( void   ) 

Definition at line 3353 of file HomePatch.C.

References Patch::atomMapper, ResizeArray< Elem >::begin(), ResizeArray< Elem >::copy(), ResizeArray< Elem >::end(), Patch::lattice, rattleListValid, AtomMapper::registerIDsFullAtom(), ResizeArray< Elem >::size(), and AtomMapper::unregisterIDsFullAtom().

Referenced by Sequencer::algorithm().

03353                            {
03354   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03355 
03356   atom.copy(checkpoint_atom);
03357   numAtoms = atom.size();
03358   lattice = checkpoint_lattice;
03359 
03360   doAtomUpdate = true;
03361   rattleListValid = false;
03362 
03363   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03364 
03365   // DMK - Atom Separation (water vs. non-water)
03366   #if NAMD_SeparateWaters != 0
03367     numWaterAtoms = checkpoint_numWaterAtoms;
03368   #endif
03369 }

void HomePatch::runSequencer ( void   ) 

Definition at line 359 of file HomePatch.C.

References Sequencer::run().

00360 { sequencer->run(); }

void HomePatch::saveForce ( const int  ftag = Results::normal  ) 

Definition at line 1371 of file HomePatch.C.

References Patch::f, and ResizeArray< Elem >::resize().

Referenced by Sequencer::saveForce().

01372 {
01373   f_saved[ftag].resize(numAtoms);
01374   for ( int i = 0; i < numAtoms; ++i )
01375   {
01376     f_saved[ftag][i] = f[ftag][i];
01377   }
01378 }

void HomePatch::sendNodeAwareSpanningTree (  ) 

Definition at line 721 of file HomePatch.C.

00721 {}

void HomePatch::sendProxies (  ) 

Definition at line 558 of file HomePatch.C.

References ResizeArray< Elem >::begin(), ProxyMgr::Object(), Patch::patchID, ProxyMgr::sendProxies(), NodeProxyMgr::sendProxyList(), and ResizeArray< Elem >::size().

00559 {
00560 #if USE_NODEPATCHMGR
00561         CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
00562         NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
00563         npm->sendProxyList(patchID, proxy.begin(), proxy.size());
00564 #else
00565         ProxyMgr::Object()->sendProxies(patchID, proxy.begin(), proxy.size());
00566 #endif
00567 }

void HomePatch::sendSpanningTree (  ) 

Definition at line 749 of file HomePatch.C.

References ResizeArray< Elem >::copy(), ProxySpanningTreeMsg::node, ProxyMgr::Object(), ProxySpanningTreeMsg::patch, Patch::patchID, ProxyMgr::sendSpanningTree(), ResizeArray< Elem >::size(), and ProxySpanningTreeMsg::tree.

Referenced by buildSpanningTree(), and recvSpanningTree().

00750 {
00751   if (tree.size() == 0) return;
00752   ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
00753   msg->patch = patchID;
00754   msg->node = CkMyPe();
00755   msg->tree.copy(tree);  // copy data for thread safety
00756   ProxyMgr::Object()->sendSpanningTree(msg);  
00757 }

void HomePatch::setGBISIntrinsicRadii (  ) 

Definition at line 3025 of file HomePatch.C.

References Patch::intRad, MassToRadius(), MassToScreen(), Node::molecule, Node::Object(), ResizeArray< Elem >::resize(), ResizeArray< Elem >::setall(), Node::simParameters, and simParams.

Referenced by positionsReady().

03025                                       {
03026   intRad.resize(numAtoms*2);
03027   intRad.setall(0);
03028   Molecule *mol = Node::Object()->molecule;
03029   SimParameters *simParams = Node::Object()->simParameters;
03030   Real offset = simParams->coulomb_radius_offset;
03031   for (int i = 0; i < numAtoms; i++) {
03032     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
03033     Real screen = MassToScreen(atom[i].mass);//same
03034     intRad[2*i+0] = rad - offset;//r0
03035     intRad[2*i+1] = screen*(rad - offset);//s0
03036   }
03037 }

void HomePatch::setLcpoType (  ) 

Definition at line 3014 of file HomePatch.C.

References Molecule::getLcpoParamType(), Patch::lcpoType, Node::molecule, Node::Object(), Patch::pExt, and ResizeArray< Elem >::resize().

Referenced by positionsReady().

03014                             {
03015   Molecule *mol = Node::Object()->molecule;
03016   const int *lcpoParamType = mol->getLcpoParamType();
03017 
03018   lcpoType.resize(numAtoms);
03019   for (int i = 0; i < numAtoms; i++) {
03020     lcpoType[i] = lcpoParamType[pExt[i].id];
03021   }
03022 }

void HomePatch::submitLoadStats ( int  timestep  ) 

Definition at line 3513 of file HomePatch.C.

References LdbCoordinator::Object(), Patch::patchID, and LdbCoordinator::patchLoad().

Referenced by Sequencer::rebalanceLoad().

03514 {
03515   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
03516 }

void HomePatch::unregisterProxy ( UnregisterProxyMsg  ) 

Definition at line 506 of file HomePatch.C.

References ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::clientRemove(), ResizeArray< Elem >::del(), Patch::forceBox, and UnregisterProxyMsg::node.

Referenced by ProxyMgr::recvUnregisterProxy().

00506                                                        {
00507 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00508   isProxyChanged = 1;
00509 #endif
00510   int n = msg->node;
00511   NodeID *pe = proxy.begin();
00512   for ( ; *pe != n; ++pe );
00513   forceBox.clientRemove();
00514   proxy.del(pe - proxy.begin());
00515   delete msg;
00516 }

void HomePatch::useSequencer ( Sequencer sequencerPtr  ) 

Definition at line 355 of file HomePatch.C.

00356 { sequencer=sequencerPtr; }


Friends And Related Function Documentation

friend class ComputeGlobal [friend]

Definition at line 49 of file HomePatch.h.

friend class PatchMgr [friend]

Definition at line 47 of file HomePatch.h.

friend class Sequencer [friend]

Definition at line 48 of file HomePatch.h.


Member Data Documentation

int HomePatch::checkpoint_task

Definition at line 178 of file HomePatch.h.

Referenced by exchangeCheckpoint(), and recvCheckpointLoad().

std::map<std::string,checkpoint_t*> HomePatch::checkpoints

Definition at line 185 of file HomePatch.h.

Referenced by recvCheckpointReq(), and recvCheckpointStore().

int HomePatch::exchange_dst

Definition at line 191 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

ExchangeAtomsMsg* HomePatch::exchange_msg

Definition at line 194 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

int HomePatch::exchange_req

Definition at line 193 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

int HomePatch::exchange_src

Definition at line 192 of file HomePatch.h.

Referenced by exchangeAtoms().

int HomePatch::inMigration [protected]

Definition at line 238 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

LDObjHandle HomePatch::ldObjHandle

Definition at line 229 of file HomePatch.h.

Referenced by LdbCoordinator::initialize(), Sequencer::Sequencer(), Sequencer::suspend(), and Sequencer::terminate().

int HomePatch::marginViolations

Definition at line 110 of file HomePatch.h.

Referenced by doAtomMigration(), doMarginCheck(), Sequencer::multigratorPressure(), and Sequencer::submitReductions().

MigrateAtomsMsg* HomePatch::msgbuf[PatchMap::MaxOneAway] [protected]

Definition at line 240 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<int> HomePatch::noconstList

Definition at line 133 of file HomePatch.h.

Referenced by rattle1().

int HomePatch::numMlBuf [protected]

Definition at line 239 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<Vector> HomePatch::posNew

Definition at line 139 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

std::vector<RattleList> HomePatch::rattleList

Definition at line 131 of file HomePatch.h.

Referenced by rattle1().

bool HomePatch::rattleListValid

Definition at line 135 of file HomePatch.h.

Referenced by positionsReady(), rattle1(), recvCheckpointLoad(), recvExchangeMsg(), and revert().

std::vector<RattleParam> HomePatch::rattleParam

Definition at line 132 of file HomePatch.h.

Referenced by rattle1().

std::vector<int> HomePatch::settleList

Definition at line 130 of file HomePatch.h.

Referenced by rattle1().

std::vector<Vector> HomePatch::velNew

Definition at line 138 of file HomePatch.h.

Referenced by addRattleForce(), buildRattleList(), and rattle1().


The documentation for this class was generated from the following files:
Generated on Wed Sep 20 01:17:18 2017 for NAMD by  doxygen 1.4.7