HomePatch Class Reference

#include <HomePatch.h>

Inheritance diagram for HomePatch:
Patch

List of all members.

Classes

struct  checkpoint_t
struct  RattleList

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 (FullAtom *__restrict atom_arr, const Force *__restrict force_arr, const BigReal dt, int num_atoms) __attribute__((__noinline__))
void addForceToMomentum3 (FullAtom *__restrict atom_arr, const Force *__restrict force_arr1, const Force *__restrict force_arr2, const Force *__restrict force_arr3, const BigReal dt1, const BigReal dt2, const BigReal dt3, int num_atoms) __attribute__((__noinline__))
void addVelocityToPosition (FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms) __attribute__((__noinline__))
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

Detailed Description

Definition at line 273 of file HomePatch.h.


Constructor & Destructor Documentation

HomePatch::~HomePatch (  ) 

Definition at line 321 of file HomePatch.C.

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

00322 {
00323     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
00324 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00325     ptnTree.resize(0);
00326     #ifdef USE_NODEPATCHMGR
00327     delete [] nodeChildren;    
00328     #endif
00329 #endif
00330   delete [] child;
00331 }


Member Function Documentation

void HomePatch::addForceToMomentum ( FullAtom *__restrict  atom_arr,
const Force *__restrict  force_arr,
const BigReal  dt,
int  num_atoms 
)

Definition at line 1844 of file HomePatch.C.

References SimParameters::fixedAtomsOn, Node::Object(), and Node::simParameters.

Referenced by Sequencer::addForceToMomentum().

01849       {
01850   SimParameters *simParams = Node::Object()->simParameters;
01851 
01852   if ( simParams->fixedAtomsOn ) {
01853     for ( int i = 0; i < num_atoms; ++i ) {
01854       if ( atom_arr[i].atomFixed ) {
01855         atom_arr[i].velocity = 0;
01856       } else {
01857         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01858         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01859         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01860         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01861       }
01862     }
01863   } else {
01864     for ( int i = 0; i < num_atoms; ++i ) {
01865       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01866       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01867       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01868       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01869     }
01870   }
01871 }

void HomePatch::addForceToMomentum3 ( FullAtom *__restrict  atom_arr,
const Force *__restrict  force_arr1,
const Force *__restrict  force_arr2,
const Force *__restrict  force_arr3,
const BigReal  dt1,
const BigReal  dt2,
const BigReal  dt3,
int  num_atoms 
)

Definition at line 1873 of file HomePatch.C.

References SimParameters::fixedAtomsOn, Node::Object(), and Node::simParameters.

Referenced by Sequencer::addForceToMomentum3().

01882       {
01883   SimParameters *simParams = Node::Object()->simParameters;
01884 
01885   if ( simParams->fixedAtomsOn ) {
01886     for ( int i = 0; i < num_atoms; ++i ) {
01887       if ( atom_arr[i].atomFixed ) {
01888         atom_arr[i].velocity = 0;
01889       } else {
01890         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01891         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01892             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01893         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01894             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01895         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01896             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01897       }
01898     }
01899   } else {
01900     for ( int i = 0; i < num_atoms; ++i ) {
01901       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01902       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01903           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01904       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01905           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01906       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01907           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01908     }
01909   }
01910 }

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

Definition at line 2283 of file HomePatch.C.

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

Referenced by rattle1().

02283                                                               {
02284   for (int ig = 0; ig < numAtoms; ++ig ) {
02285     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
02286     Tensor vir = outer(df, atom[ig].position);
02287     wc += vir;
02288     f[Results::normal][ig] += df;
02289     atom[ig].velocity = velNew[ig];
02290   }
02291 }

void HomePatch::addVelocityToPosition ( FullAtom *__restrict  atom_arr,
const BigReal  dt,
int  num_atoms 
)

Definition at line 1912 of file HomePatch.C.

References SimParameters::fixedAtomsOn, Node::Object(), and Node::simParameters.

Referenced by Sequencer::addVelocityToPosition().

01916       {
01917   SimParameters *simParams = Node::Object()->simParameters;
01918   if ( simParams->fixedAtomsOn ) {
01919     for ( int i = 0; i < num_atoms; ++i ) {
01920       if ( ! atom_arr[i].atomFixed ) {
01921         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01922         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01923         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01924       }
01925     }
01926   } else {
01927     for ( int i = 0; i < num_atoms; ++i ) {
01928       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01929       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01930       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01931     }
01932   }
01933 }

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

Implements Patch.

Definition at line 334 of file HomePatch.C.

References Sequencer::awaken(), Patch::boxesOpen, DebugM, Patch::dEdaSumBox, Flags::doGBIS, Flags::doNonbonded, Patch::f, Patch::flags, ExtForce::force, OwnerBox< Owner, Data >::isOpen(), j, Results::maxNumForces, Results::normal, Patch::numAtoms, Patch::patchID, Patch::psiSumBox, and ResizeArray< Elem >::size().

00334                                  {
00335   // begin gbis
00336   if (box == 5) {// end of phase 1
00337     phase1BoxClosedCalled = true;
00338     if (!psiSumBox.isOpen() && numGBISP1Arrived == proxy.size()) {
00339       if (flags.doGBIS && flags.doNonbonded) {
00340         sequencer->awaken();
00341       }
00342     } else {
00343       //need to wait until proxies arrive before awakening
00344     }
00345   } else if (box == 6) {// intRad
00346     //do nothing
00347   } else if (box == 7) {// bornRad
00348     //do nothing
00349   } else if (box == 8) {// end of phase 2
00350     phase2BoxClosedCalled = true;
00351     //if no proxies, AfterP1 can't be called from receive
00352     //so it will be called from here
00353     if (!dEdaSumBox.isOpen() && numGBISP2Arrived == proxy.size()) {
00354       if (flags.doGBIS && flags.doNonbonded) {
00355         sequencer->awaken();
00356       }
00357     } else {
00358       //need to wait until proxies arrive before awakening
00359     }
00360   } else if (box == 9) {
00361     //do nothing
00362   } else if (box == 10) {
00363     //lcpoType Box closed: do nothing
00364   } else {
00365     //do nothing
00366   }
00367   // end gbis
00368 
00369   if ( ! --boxesOpen ) {
00370     if ( replacementForces ) {
00371       for ( int i = 0; i < numAtoms; ++i ) {
00372         if ( replacementForces[i].replace ) {
00373           for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
00374           f[Results::normal][i] = replacementForces[i].force;
00375         }
00376       }
00377       replacementForces = 0;
00378     }
00379     DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
00380         << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
00381     // only awaken suspended threads.  Then say it is suspended.
00382 
00383     phase3BoxClosedCalled = true;
00384     if (flags.doGBIS) {
00385       if (flags.doNonbonded) {
00386         sequencer->awaken();
00387       } else {
00388         if (numGBISP1Arrived == proxy.size() &&
00389           numGBISP2Arrived == proxy.size() &&
00390           numGBISP3Arrived == proxy.size()) {
00391           sequencer->awaken();//all boxes closed and all proxies arrived
00392         }
00393       }
00394     } else {//non-gbis awaken
00395       sequencer->awaken();
00396     }
00397   } else {
00398     DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
00399   }
00400 }

void HomePatch::buildRattleList (  ) 

Definition at line 2141 of file HomePatch.C.

References RattleParam::dsq, SimParameters::fixedAtomsOn, RattleParam::ia, RattleParam::ib, HomePatch::RattleList::icnt, HomePatch::RattleList::ig, NAMD_die(), noconstList, Node::Object(), posNew, rattleList, rattleParam, RattleParam::rma, RattleParam::rmb, settle1init(), settleList, Node::simParameters, SimParameters::useSettle, velNew, WAT_SWM4, WAT_TIP4, and SimParameters::watmodel.

Referenced by rattle1().

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

void HomePatch::buildSpanningTree ( void   ) 

Definition at line 674 of file HomePatch.C.

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

Referenced by ProxyMgr::buildProxySpanningTree().

00675 {
00676   nChild = 0;
00677   int psize = proxy.size();
00678   if (psize == 0) return;
00679   NodeIDList oldtree;  oldtree.swap(tree);
00680   int oldsize = oldtree.size();
00681   tree.resize(psize + 1);
00682   tree.setall(-1);
00683   tree[0] = CkMyPe();
00684   int s=1, e=psize+1;
00685   NodeIDList::iterator pli;
00686   int patchNodesLast =
00687     ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00688   int nNonPatch = 0;
00689 #if 1
00690     // try to put it to the same old tree
00691   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
00692   {
00693     int oldindex = oldtree.find(*pli);
00694     if (oldindex != -1 && oldindex < psize) {
00695       tree[oldindex] = *pli;
00696     }
00697   }
00698   s=1; e=psize;
00699   for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
00700   {
00701     if (tree.find(*pli) != -1) continue;    // already assigned
00702     if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
00703       while (tree[e] != -1) { e--; if (e==-1) e = psize; }
00704       tree[e] = *pli;
00705     } else {
00706       while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
00707       tree[s] = *pli;
00708       nNonPatch++;
00709     }
00710   }
00711 #if 1
00712   if (oldsize==0 && nNonPatch) {
00713     // first time, sort by distance
00714     qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
00715   }
00716 #endif
00717 
00718   //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
00719 
00720 #if USE_TOPOMAP && USE_SPANNING_TREE
00721   
00722   //Right now only works for spanning trees with two levels
00723   int *treecopy = new int [psize];
00724   int subroots[proxySpanDim];
00725   int subsizes[proxySpanDim];
00726   int subtrees[proxySpanDim][proxySpanDim];
00727   int idxes[proxySpanDim];
00728   int i = 0;
00729 
00730   for(i=0;i<proxySpanDim;i++){
00731     subsizes[i] = 0;
00732     idxes[i] = i;
00733   }
00734   
00735   for(i=0;i<psize;i++){
00736     treecopy[i] = tree[i+1];
00737   }
00738   
00739   TopoManager tmgr;
00740   tmgr.sortRanksByHops(treecopy,nNonPatch);
00741   tmgr.sortRanksByHops(treecopy+nNonPatch,
00742                                                 psize-nNonPatch);  
00743   
00744   /* build tree and subtrees */
00745   nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
00746   delete [] treecopy;
00747   
00748   for(int i=1;i<psize+1;i++){
00749     int isSubroot=0;
00750     for(int j=0;j<nChild;j++)
00751       if(tree[i]==subroots[j]){
00752         isSubroot=1;
00753         break;
00754       }
00755     if(isSubroot) continue;
00756     
00757     int bAdded = 0;
00758     tmgr.sortIndexByHops(tree[i], subroots,
00759                                                   idxes, proxySpanDim);
00760     for(int j=0;j<proxySpanDim;j++){
00761       if(subsizes[idxes[j]]<proxySpanDim){
00762         subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
00763         bAdded = 1; 
00764         break;
00765       }
00766     }
00767     if( psize > proxySpanDim && ! bAdded ) {
00768       NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
00769     }
00770   }
00771 
00772 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
00773   
00774   for (int i=1; i<=proxySpanDim; i++) {
00775     if (tree.size() <= i) break;
00776     child[i-1] = tree[i];
00777     nChild++;
00778   }
00779 #endif
00780 #endif
00781   
00782 #if 0
00783   // for debugging
00784   CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
00785   for (int i=0; i<psize+1; i++) {
00786     CkPrintf("%d ", tree[i]);
00787   }
00788   CkPrintf("\n");
00789 #endif
00790   // send to children nodes
00791   sendSpanningTree();
00792 }

void HomePatch::checkpoint ( void   ) 

Definition at line 3361 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

03361                                {
03362   checkpoint_atom.copy(atom);
03363   checkpoint_lattice = lattice;
03364 
03365   // DMK - Atom Separation (water vs. non-water)
03366   #if NAMD_SeparateWaters != 0
03367     checkpoint_numWaterAtoms = numWaterAtoms;
03368   #endif
03369 }

void HomePatch::depositMigration ( MigrateAtomsMsg msg  ) 

Definition at line 3908 of file HomePatch.C.

References ResizeArray< Elem >::add(), Lattice::apply_transform(), Sequencer::awaken(), ResizeArray< Elem >::begin(), DebugM, ResizeArray< Elem >::end(), CompAtom::hydrogenGroupSize, CompAtomExt::id, inMigration, Patch::lattice, FullAtom::migrationGroupSize, MigrateAtomsMsg::migrationList, msgbuf, Lattice::nearest(), numMlBuf, Patch::patchID, CompAtom::position, Lattice::reverse_transform(), ResizeArray< Elem >::size(), and FullAtom::transform.

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

03909 {
03910 
03911   if (!inMigration) { // We have to buffer changes due to migration
03912                       // until our patch is in migration mode
03913     msgbuf[numMlBuf++] = msg;
03914     return;
03915   } 
03916 
03917 
03918   // DMK - Atom Separation (water vs. non-water)
03919   #if NAMD_SeparateWaters != 0
03920 
03921 
03922     // Merge the incoming list of atoms with the current list of
03923     //   atoms.  Note that mergeSeparatedAtomList() will apply any
03924     //   required transformations to the incoming atoms as it is
03925     //   separating them.
03926     mergeAtomList(msg->migrationList);
03927 
03928 
03929   #else
03930 
03931     // Merge the incoming list of atoms with the current list of
03932     // atoms.  Apply transformations to the incoming atoms as they are
03933     // added to this patch's list.
03934     {
03935       MigrationList &migrationList = msg->migrationList;
03936       MigrationList::iterator mi;
03937       Transform mother_transform;
03938       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
03939         DebugM(1,"Migrating atom " << mi->id << " to patch "
03940                   << patchID << " with position " << mi->position << "\n"); 
03941         if ( mi->migrationGroupSize ) {
03942           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
03943             Position pos = mi->position;
03944             int mgs = mi->migrationGroupSize;
03945             int c = 1;
03946             for ( int j=mi->hydrogenGroupSize; j<mgs;
03947                                 j+=(mi+j)->hydrogenGroupSize ) {
03948               pos += (mi+j)->position;
03949               ++c;
03950             }
03951             pos *= 1./c;
03952             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
03953             mother_transform = mi->transform;
03954             pos = lattice.nearest(pos,center,&mother_transform);
03955             mi->position = lattice.reverse_transform(mi->position,mi->transform);
03956             mi->position = lattice.apply_transform(mi->position,mother_transform);
03957             mi->transform = mother_transform;
03958           } else {
03959             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
03960             mother_transform = mi->transform;
03961           }
03962         } else {
03963           mi->position = lattice.reverse_transform(mi->position,mi->transform);
03964           mi->position = lattice.apply_transform(mi->position,mother_transform);
03965           mi->transform = mother_transform;
03966         }
03967         atom.add(*mi);
03968       }
03969     }
03970 
03971 
03972   #endif // if (NAMD_SeparateWaters != 0)
03973 
03974 
03975   numAtoms = atom.size();
03976   delete msg;
03977 
03978   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
03979   if (!--patchMigrationCounter) {
03980     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
03981     allMigrationIn = true;
03982     patchMigrationCounter = numNeighbors;
03983     if (migrationSuspended) {
03984       DebugM(3,"patch "<<patchID<<" is being awakened\n");
03985       migrationSuspended = false;
03986       sequencer->awaken();
03987       return;
03988     }
03989     else {
03990        DebugM(3,"patch "<<patchID<<" was not suspended\n");
03991     }
03992   }
03993 }

void HomePatch::doAtomMigration (  )  [protected]

Definition at line 3752 of file HomePatch.C.

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

Referenced by positionsReady().

03753 {
03754   int i;
03755 
03756   for (i=0; i<numNeighbors; i++) {
03757     realInfo[i].mList.resize(0);
03758   }
03759 
03760   // Purge the AtomMap
03761   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03762 
03763   // Determine atoms that need to migrate
03764 
03765   BigReal minx = min.x;
03766   BigReal miny = min.y;
03767   BigReal minz = min.z;
03768   BigReal maxx = max.x;
03769   BigReal maxy = max.y;
03770   BigReal maxz = max.z;
03771 
03772   int xdev, ydev, zdev;
03773   int delnum = 0;
03774 
03775   FullAtomList::iterator atom_i = atom.begin();
03776   FullAtomList::iterator atom_e = atom.end();
03777 
03778   // DMK - Atom Separation (water vs. non-water)
03779   #if NAMD_SeparateWaters != 0
03780     FullAtomList::iterator atom_first = atom_i;
03781     int numLostWaterAtoms = 0;
03782   #endif
03783 
03784   while ( atom_i != atom_e ) {
03785     // Even though this code iterates through all atoms successively
03786     // it moves entire hydrogen/migration groups as follows:
03787     // Only the parent atom of the hydrogen/migration group has
03788     // nonzero migrationGroupSize.  Values determined for xdev,ydev,zdev
03789     // will persist through the remaining group members so that each
03790     // following atom will again be added to the same mList.
03791     if ( atom_i->migrationGroupSize ) {
03792       Position pos = atom_i->position;
03793       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
03794         // If there are multiple hydrogen groups in a migration group
03795         // (e.g. for supporting lone pairs)
03796         // the following code takes the average position (midpoint)
03797         // of their parents.
03798         int mgs = atom_i->migrationGroupSize;
03799         int c = 1;
03800         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
03801                                 j+=(atom_i+j)->hydrogenGroupSize ) {
03802           pos += (atom_i+j)->position;
03803           ++c;
03804         }
03805         pos *= 1./c;
03806         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
03807       }
03808 
03809       // Scaling the position below transforms space within patch from
03810       // what could have been a rotated parallelepiped into
03811       // orthogonal coordinates, where we can use minmax comparison
03812       // to detect which of our nearest neighbors this
03813       // parent atom might have entered.
03814       ScaledPosition s = lattice.scale(pos);
03815 
03816       // check if atom is within bounds
03817       if (s.x < minx) xdev = 0;
03818       else if (maxx <= s.x) xdev = 2;
03819       else xdev = 1;
03820 
03821       if (s.y < miny) ydev = 0;
03822       else if (maxy <= s.y) ydev = 2;
03823       else ydev = 1;
03824 
03825       if (s.z < minz) zdev = 0;
03826       else if (maxz <= s.z) zdev = 2;
03827       else zdev = 1;
03828 
03829     }
03830 
03831     if (mInfo[xdev][ydev][zdev]) { // process atom for migration
03832                                     // Don't migrate if destination is myself
03833 
03834       // See if we have a migration list already
03835       MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
03836       DebugM(3,"Migrating atom " << atom_i->id << " from patch "
03837                 << patchID << " with position " << atom_i->position << "\n");
03838       mCur.add(*atom_i);
03839 
03840       ++delnum;
03841 
03842 
03843       // DMK - Atom Separation (water vs. non-water)
03844       #if NAMD_SeparateWaters != 0
03845         // Check to see if this atom is part of a water molecule.  If
03846         //   so, numWaterAtoms needs to be adjusted to refect the lost of
03847         //   this atom.
03848         // NOTE: The atom separation code assumes that if the oxygen
03849         //   atom of the hydrogen group making up the water molecule is
03850         //   migrated to another HomePatch, the hydrogens will also
03851         //   move!!!
03852         int atomIndex = atom_i - atom_first;
03853         if (atomIndex < numWaterAtoms)
03854           numLostWaterAtoms++;
03855       #endif
03856 
03857 
03858     } else {
03859       // By keeping track of delnum total being deleted from FullAtomList
03860       // the else clause allows us to fill holes as we visit each atom.
03861 
03862       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
03863 
03864     }
03865 
03866     ++atom_i;
03867   }
03868 
03869   // DMK - Atom Separation (water vs. non-water)
03870   #if NAMD_SeparateWaters != 0
03871     numWaterAtoms -= numLostWaterAtoms;
03872   #endif
03873 
03874 
03875   int delpos = numAtoms - delnum;
03876   DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
03877   atom.del(delpos,delnum);
03878 
03879   numAtoms = atom.size();
03880 
03881   PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
03882 
03883   // signal depositMigration() that we are inMigration mode
03884   inMigration = true;
03885 
03886   // Drain the migration message buffer
03887   for (i=0; i<numMlBuf; i++) {
03888      DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
03889      depositMigration(msgbuf[i]);
03890   }
03891   numMlBuf = 0;
03892 
03893   NAMD_EVENT_STOP(1, NamdProfileEvent::ATOM_MIGRATIONS);
03894 
03895   if (!allMigrationIn) {
03896     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
03897     migrationSuspended = true;
03898     sequencer->suspend();
03899     migrationSuspended = false;
03900   }
03901   allMigrationIn = false;
03902 
03903   inMigration = false;
03904   marginViolations = 0;
03905 }

void HomePatch::doGroupSizeCheck (  )  [protected]

Definition at line 3631 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Flags::doNonbonded, ResizeArray< Elem >::end(), Patch::flags, SimParameters::hgroupCutoff, CompAtom::hydrogenGroupSize, Flags::maxGroupRadius, NAMD_bug(), CompAtom::nonbondedGroupSize, Node::Object(), CompAtom::position, Node::simParameters, Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

03632 {
03633   if ( ! flags.doNonbonded ) return;
03634 
03635   SimParameters *simParams = Node::Object()->simParameters;
03636   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
03637   BigReal maxrad2 = 0.;
03638 
03639   FullAtomList::iterator p_i = atom.begin();
03640   FullAtomList::iterator p_e = atom.end();
03641 
03642   while ( p_i != p_e ) {
03643     const int hgs = p_i->hydrogenGroupSize;
03644     if ( ! hgs ) break;  // avoid infinite loop on bug
03645     int ngs = hgs;
03646     if ( ngs > 5 ) ngs = 5;  // XXX why? limit to at most 5 atoms per group
03647     BigReal x = p_i->position.x;
03648     BigReal y = p_i->position.y;
03649     BigReal z = p_i->position.z;
03650     int i;
03651     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
03652       p_i[i].nonbondedGroupSize = 0;
03653       BigReal dx = p_i[i].position.x - x;
03654       BigReal dy = p_i[i].position.y - y;
03655       BigReal dz = p_i[i].position.z - z;
03656       BigReal r2 = dx * dx + dy * dy + dz * dz;
03657       if ( r2 > hgcut ) break;
03658       else if ( r2 > maxrad2 ) maxrad2 = r2;
03659     }
03660     p_i->nonbondedGroupSize = i;
03661     for ( ; i < hgs; ++i ) {
03662       p_i[i].nonbondedGroupSize = 1;
03663     }
03664     p_i += hgs;
03665   }
03666 
03667   if ( p_i != p_e ) {
03668     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
03669   }
03670 
03671   flags.maxGroupRadius = sqrt(maxrad2);
03672 
03673 }

void HomePatch::doMarginCheck (  )  [protected]

Definition at line 3675 of file HomePatch.C.

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

Referenced by positionsReady().

03676 {
03677   SimParameters *simParams = Node::Object()->simParameters;
03678 
03679   BigReal sysdima = lattice.a_r().unit() * lattice.a();
03680   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
03681   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
03682 
03683   BigReal minSize = simParams->patchDimension - simParams->margin;
03684 
03685   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
03686        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
03687        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
03688 
03689     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03690       "Possible solutions are to restart from a recent checkpoint,\n"
03691       "increase margin, or disable useFlexibleCell for liquid simulation.");
03692   }
03693 
03694   BigReal cutoff = simParams->cutoff;
03695 
03696   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
03697   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
03698   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
03699 
03700   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
03701     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03702       "There are probably many margin violations already reported.\n"
03703       "Possible solutions are to restart from a recent checkpoint,\n"
03704       "increase margin, or disable useFlexibleCell for liquid simulation.");
03705   }
03706 
03707   BigReal minx = min.x - margina;
03708   BigReal miny = min.y - marginb;
03709   BigReal minz = min.z - marginc;
03710   BigReal maxx = max.x + margina;
03711   BigReal maxy = max.y + marginb;
03712   BigReal maxz = max.z + marginc;
03713 
03714   int xdev, ydev, zdev;
03715   int problemCount = 0;
03716 
03717   FullAtomList::iterator p_i = atom.begin();
03718   FullAtomList::iterator p_e = atom.end();
03719   for ( ; p_i != p_e; ++p_i ) {
03720 
03721     ScaledPosition s = lattice.scale(p_i->position);
03722 
03723     // check if atom is within bounds
03724     if (s.x < minx) xdev = 0;
03725     else if (maxx <= s.x) xdev = 2; 
03726     else xdev = 1;
03727 
03728     if (s.y < miny) ydev = 0;
03729     else if (maxy <= s.y) ydev = 2; 
03730     else ydev = 1;
03731 
03732     if (s.z < minz) zdev = 0;
03733     else if (maxz <= s.z) zdev = 2; 
03734     else zdev = 1;
03735 
03736     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
03737         ++problemCount;
03738     }
03739 
03740   }
03741 
03742   marginViolations = problemCount;
03743   // if ( problemCount ) {
03744   //     iout << iERROR <<
03745   //       "Found " << problemCount << " margin violations!\n" << endi;
03746   // } 
03747 
03748 }

void HomePatch::doPairlistCheck (  )  [protected]

Definition at line 3540 of file HomePatch.C.

References ResizeArray< Elem >::begin(), Patch::flags, Patch::getPatchID(), Patch::lattice, Vector::length2(), Flags::maxAtomMovement, NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NamdProfileEventStr, Node::Object(), SimParameters::pairlistGrow, SimParameters::pairlistShrink, Flags::pairlistTolerance, SimParameters::pairlistTrigger, CompAtom::position, ResizeArray< Elem >::resize(), Flags::savePairlists, Node::simParameters, Lattice::unscale(), Flags::usePairlists, Vector::x, Vector::y, and Vector::z.

Referenced by positionsReady().

03541 {
03542 #if 0
03543 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
03544   char dpcbuf[32];
03545   sprintf(dpcbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::DO_PAIRLIST_CHECK], this->getPatchID());
03546   NAMD_EVENT_START_EX(1, NamdProfileEvent::DO_PAIRLIST_CHECK, dpcbuf);
03547 #endif
03548 #endif
03549 
03550   SimParameters *simParams = Node::Object()->simParameters;
03551 
03552   if ( numAtoms == 0 || ! flags.usePairlists ) {
03553     flags.pairlistTolerance = 0.;
03554     flags.maxAtomMovement = 99999.;
03555 #if 0
03556     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03557 #endif
03558     return;
03559   }
03560 
03561   int i; int n = numAtoms;
03562   CompAtom *p_i = p.begin();
03563 
03564   if ( flags.savePairlists ) {
03565     flags.pairlistTolerance = doPairlistCheck_newTolerance;
03566     flags.maxAtomMovement = 0.;
03567     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
03568     doPairlistCheck_lattice = lattice;
03569     doPairlistCheck_positions.resize(numAtoms);
03570     CompAtom *psave_i = doPairlistCheck_positions.begin();
03571     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
03572 #if 0
03573     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03574 #endif
03575     return;
03576   }
03577 
03578   Lattice &lattice_old = doPairlistCheck_lattice;
03579   Position center_cur = lattice.unscale(center);
03580   Position center_old = lattice_old.unscale(center);
03581   Vector center_delta = center_cur - center_old;
03582   
03583   // find max deviation to corner (any neighbor shares a corner)
03584   BigReal max_cd = 0.;
03585   for ( i=0; i<2; ++i ) {
03586     for ( int j=0; j<2; ++j ) {
03587       for ( int k=0; k<2; ++k ) {
03588         ScaledPosition corner(  i ? min.x : max.x ,
03589                                 j ? min.y : max.y ,
03590                                 k ? min.z : max.z );
03591         Vector corner_delta =
03592                 lattice.unscale(corner) - lattice_old.unscale(corner);
03593         corner_delta -= center_delta;
03594         BigReal cd = corner_delta.length2();
03595         if ( cd > max_cd ) max_cd = cd;
03596       }
03597     }
03598   }
03599   max_cd = sqrt(max_cd);
03600 
03601   // find max deviation of atoms relative to center
03602   BigReal max_pd = 0.;
03603   CompAtom *p_old_i = doPairlistCheck_positions.begin();
03604   for ( i=0; i<n; ++i ) {
03605     Vector p_delta = p_i[i].position - p_old_i[i].position;
03606     p_delta -= center_delta;
03607     BigReal pd = p_delta.length2();
03608     if ( pd > max_pd ) max_pd = pd;
03609   }
03610   max_pd = sqrt(max_pd);
03611 
03612   BigReal max_tol = max_pd + max_cd;
03613 
03614   flags.maxAtomMovement = max_tol;
03615 
03616   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
03617 
03618   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
03619                                 doPairlistCheck_newTolerance ) ) {
03620     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
03621   }
03622 
03623   if ( max_tol > doPairlistCheck_newTolerance ) {
03624     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
03625   }
03626 #if 0
03627   NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03628 #endif
03629 }

void HomePatch::exchangeAtoms ( int  scriptTask  ) 

Definition at line 3485 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, SimParameters::scriptArg1, SimParameters::scriptArg2, PatchMgr::sendExchangeReq(), and Node::simParameters.

Referenced by Sequencer::algorithm().

03485                                             {
03486   SimParameters *simParams = Node::Object()->simParameters;
03487   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
03488   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
03489     exchange_dst = (int) simParams->scriptArg1;
03490     // create and save outgoing message
03491     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
03492     exchange_msg->lattice = lattice;
03493     exchange_msg->pid = patchID;
03494     exchange_msg->numAtoms = numAtoms;
03495     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03496     if ( exchange_req >= 0 ) {
03497       recvExchangeReq(exchange_req);
03498     }
03499   }
03500   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
03501     exchange_src = (int) simParams->scriptArg2;
03502     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
03503   }
03504 }

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

Definition at line 3389 of file HomePatch.C.

References checkpoint_task, PatchMgr::Object(), Node::Object(), Patch::patchID, SimParameters::scriptIntArg1, SimParameters::scriptStringArg1, PatchMgr::sendCheckpointReq(), and Node::simParameters.

Referenced by Sequencer::algorithm().

03389                                                            {  // initiating replica
03390   SimParameters *simParams = Node::Object()->simParameters;
03391   checkpoint_task = scriptTask;
03392   const int remote = simParams->scriptIntArg1;
03393   const char *key = simParams->scriptStringArg1;
03394   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
03395 }

void HomePatch::gbisComputeAfterP1 (  ) 

Definition at line 3058 of file HomePatch.C.

References SimParameters::alpha_max, Patch::bornRad, SimParameters::coulomb_radius_offset, Patch::flags, SimParameters::gbis_beta, SimParameters::gbis_delta, SimParameters::gbis_gamma, gbisP2Ready(), Patch::intRad, Node::Object(), Patch::pExt, Patch::psiFin, Patch::psiSum, Flags::sequence, and Node::simParameters.

03058                                    {
03059 
03060   SimParameters *simParams = Node::Object()->simParameters;
03061   BigReal alphaMax = simParams->alpha_max;
03062   BigReal delta = simParams->gbis_delta;
03063   BigReal beta = simParams->gbis_beta;
03064   BigReal gamma = simParams->gbis_gamma;
03065   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03066 
03067   BigReal rhoi;
03068   BigReal rhoi0;
03069   //calculate bornRad from psiSum
03070   for (int i = 0; i < numAtoms; i++) {
03071     rhoi0 = intRad[2*i];
03072     rhoi = rhoi0+coulomb_radius_offset;
03073     psiFin[i] += psiSum[i];
03074     psiFin[i] *= rhoi0;
03075     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
03076     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
03077 #ifdef PRINT_COMP
03078     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
03079 #endif
03080   }
03081 
03082   gbisP2Ready();
03083 }

void HomePatch::gbisComputeAfterP2 (  ) 

Definition at line 3086 of file HomePatch.C.

References Patch::bornRad, COULOMB, SimParameters::coulomb_radius_offset, Patch::dEdaSum, Patch::dHdrPrefix, SimParameters::dielectric, Patch::flags, SimParameters::gbis_beta, SimParameters::gbis_delta, SimParameters::gbis_gamma, gbisP3Ready(), Patch::intRad, SimParameters::kappa, Node::Object(), Patch::pExt, Patch::psiFin, Flags::sequence, Node::simParameters, and SimParameters::solvent_dielectric.

03086                                    {
03087 
03088   SimParameters *simParams = Node::Object()->simParameters;
03089   BigReal delta = simParams->gbis_delta;
03090   BigReal beta = simParams->gbis_beta;
03091   BigReal gamma = simParams->gbis_gamma;
03092   BigReal epsilon_s = simParams->solvent_dielectric;
03093   BigReal epsilon_p = simParams->dielectric;
03094   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
03095   BigReal epsilon_p_i = 1/simParams->dielectric;
03096   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03097   BigReal kappa = simParams->kappa;
03098   BigReal fij, expkappa, Dij, dEdai, dedasum;
03099   BigReal rhoi, rhoi0, psii, nbetapsi;
03100   BigReal gammapsi2, tanhi, daidr;
03101   for (int i = 0; i < numAtoms; i++) {
03102     //add diagonal dEda term
03103     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
03104     fij = bornRad[i];//inf
03105     expkappa = exp(-kappa*fij);//0
03106     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
03107     //calculate dHij prefix
03108     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
03109                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
03110     dHdrPrefix[i] += dEdai;
03111     dedasum = dHdrPrefix[i];
03112 
03113     rhoi0 = intRad[2*i];
03114     rhoi = rhoi0+coulomb_radius_offset;
03115     psii = psiFin[i];
03116     nbetapsi = -beta*psii;
03117     gammapsi2 = gamma*psii*psii;
03118     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
03119     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
03120            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
03121     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
03122 #ifdef PRINT_COMP
03123     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
03124 #endif
03125   }
03126   gbisP3Ready();
03127 }

void HomePatch::gbisP2Ready (  ) 

Reimplemented from Patch.

Definition at line 3130 of file HomePatch.C.

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

Referenced by gbisComputeAfterP1().

03130                             {
03131   if (proxy.size() > 0) {
03132     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03133     for (int i = 0; i < proxy.size(); i++) {
03134       int node = proxy[i];
03135       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
03136       msg->patch = patchID;
03137       msg->origPe = CkMyPe();
03138       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
03139       msg->destPe = node;
03140       int seq = flags.sequence;
03141       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03142       SET_PRIORITY(msg,seq,priority);
03143       cp[node].recvData(msg);
03144     }
03145   }
03146   Patch::gbisP2Ready();
03147 }

void HomePatch::gbisP3Ready (  ) 

Reimplemented from Patch.

Definition at line 3150 of file HomePatch.C.

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

Referenced by gbisComputeAfterP2().

03150                             {
03151   if (proxy.size() > 0) {
03152     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03153     //only nonzero message should be sent for doFullElec
03154     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
03155     for (int i = 0; i < proxy.size(); i++) {
03156       int node = proxy[i];
03157       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
03158       msg->patch = patchID;
03159       msg->dHdrPrefixLen = msgAtoms;
03160       msg->origPe = CkMyPe();
03161       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
03162       msg->destPe = node;
03163       int seq = flags.sequence;
03164       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03165       SET_PRIORITY(msg,seq,priority);
03166       cp[node].recvData(msg);
03167     }
03168   }
03169   Patch::gbisP3Ready();
03170 }

FullAtomList& HomePatch::getAtomList (  )  [inline]

Definition at line 455 of file HomePatch.h.

Referenced by ComputeHomePatch::doWork(), and dumpbench().

00455 { return (atom); }

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

Definition at line 1935 of file HomePatch.C.

References BOLTZMANN, Lattice::c(), SimParameters::drudeBondLen, SimParameters::drudeHardWallOn, SimParameters::drudeTemp, endi(), SimParameters::fixedAtomsOn, iERROR(), iout, SubmitReduction::item(), Patch::lattice, Node::molecule, Node::Object(), Lattice::origin(), outer(), partition(), SimParameters::pressureProfileSlabs, SimParameters::rigidDie, Node::simParameters, TIMEFACTOR, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz.

Referenced by Sequencer::hardWallDrude().

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

void HomePatch::loweAndersenFinish (  ) 

Definition at line 3024 of file HomePatch.C.

References DebugM.

03025 {
03026     DebugM(2, "loweAndersenFinish\n");
03027     v.resize(0);
03028 }

void HomePatch::loweAndersenVelocities (  ) 

Definition at line 3009 of file HomePatch.C.

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

Referenced by positionsReady().

03010 {
03011     DebugM(2, "loweAndersenVelocities\n");
03012     Molecule *mol = Node::Object()->molecule;
03013     SimParameters *simParams = Node::Object()->simParameters;
03014     v.resize(numAtoms);
03015     for (int i = 0; i < numAtoms; ++i) {
03016         //v[i] = p[i];
03017         // co-opt CompAtom structure to pass velocity and mass information
03018         v[i].position = atom[i].velocity;
03019         v[i].charge = atom[i].mass;
03020     }
03021     DebugM(2, "loweAndersenVelocities\n");
03022 }

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

Definition at line 2874 of file HomePatch.C.

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

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

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

void HomePatch::mollyAverage (  ) 

Definition at line 3224 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::mollyMollify ( Tensor virial  ) 

Definition at line 3298 of file HomePatch.C.

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

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

void HomePatch::positionsReady ( int  doMigration = 0  ) 

Reimplemented from Patch.

Definition at line 925 of file HomePatch.C.

References ProxyDataMsg::avgPlLen, ProxyDataMsg::avgPositionList, Patch::avgPositionPtrBegin, Patch::avgPositionPtrEnd, ResizeArray< Elem >::begin(), CompAtom::charge, COULOMB, ProxyDataMsg::cudaAtomList, Patch::cudaAtomPtr, DebugM, ComputeNonbondedUtil::dielectric_1, doAtomMigration(), Flags::doGBIS, doGroupSizeCheck(), Flags::doLCPO, Flags::doLoweAndersen, doMarginCheck(), Flags::doMolly, doPairlistCheck(), ResizeArray< Elem >::end(), SimParameters::fixedAtomsForces, SimParameters::fixedAtomsOn, ProxyDataMsg::flags, Patch::flags, Patch::getPatchID(), Patch::intRad, ProxyDataMsg::intRadList, j, Patch::lattice, Patch::lcpoType, ProxyDataMsg::lcpoTypeList, loweAndersenVelocities(), mollyAverage(), NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NamdProfileEventStr, Patch::numAtoms, PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), Sync::Object(), ProxyMgr::Object(), PatchMap::Object(), Node::Object(), Patch::p_avg, ProxyDataMsg::patch, PATCH_PRIORITY, Patch::patchID, Sync::PatchReady(), Patch::pExt, ProxyDataMsg::plExtLen, ProxyDataMsg::plLen, CompAtom::position, ProxyDataMsg::positionExtList, ProxyDataMsg::positionList, PRIORITY_SIZE, PROXY_DATA_PRIORITY, proxySendSpanning, proxySpanDim, CudaAtom::q, SimParameters::qmLSSOn, qmSwapAtoms(), rattleListValid, ResizeArray< Elem >::resize(), ComputeNonbondedUtil::scaling, ProxyMgr::sendProxyAll(), ProxyMgr::sendProxyData(), Flags::sequence, SET_PRIORITY, setGBISIntrinsicRadii(), setLcpoType(), Node::simParameters, ResizeArray< Elem >::size(), sortAtomsForCUDA(), CompAtomExt::sortOrder, SimParameters::staticAtomAssignment, Lattice::unscale(), Patch::v, ProxyDataMsg::velocityList, Patch::velocityPtrBegin, Patch::velocityPtrEnd, ProxyDataMsg::vlLen, CudaAtom::x, Vector::x, CudaAtom::y, Vector::y, CudaAtom::z, and Vector::z.

00926 {    
00927   SimParameters *simParams = Node::Object()->simParameters;
00928 
00929   flags.sequence++;
00930 
00931   if (!patchMapRead) { readPatchMap(); }
00932       
00933   if (numNeighbors && ! simParams->staticAtomAssignment) {
00934     if (doMigration) {
00935       rattleListValid = false;
00936       doAtomMigration();
00937     } else {
00938       doMarginCheck();
00939     }
00940   }
00941 
00942 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
00943   char prbuf[32];
00944   sprintf(prbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::POSITIONS_READY], this->getPatchID());
00945   NAMD_EVENT_START_EX(1, NamdProfileEvent::POSITIONS_READY, prbuf);
00946 #endif
00947 
00948   if (doMigration && simParams->qmLSSOn)
00949       qmSwapAtoms();
00950 
00951 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
00952   if ( doMigration ) {
00953     int n = numAtoms;
00954     FullAtom *a_i = atom.begin();
00955     #if defined(NAMD_CUDA) || (defined(NAMD_MIC) && MIC_SORT_ATOMS != 0)
00956     int *ao = new int[n];
00957     int nfree;
00958     if ( simParams->fixedAtomsOn && ! simParams->fixedAtomsForces ) {
00959       int k = 0;
00960       int k2 = n;
00961       for ( int j=0; j<n; ++j ) {
00962         // put fixed atoms at end
00963         if ( a_i[j].atomFixed ) ao[--k2] = j;
00964         else ao[k++] = j;
00965       }
00966       nfree = k;
00967     } else {
00968       nfree = n;
00969       for ( int j=0; j<n; ++j ) {
00970         ao[j] = j;
00971       }
00972     }
00973 
00974     sortAtomsForCUDA(ao,a_i,nfree,n);
00975   
00976     for ( int i=0; i<n; ++i ) { 
00977       a_i[i].sortOrder = ao[i];
00978     }
00979     delete [] ao;
00980     #else
00981       for (int i = 0; i < n; ++i) {
00982         a_i[i].sortOrder = i;
00983       }
00984     #endif
00985   }
00986 
00987   { 
00988     const double charge_scaling = sqrt(COULOMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1);
00989     const Vector ucenter = lattice.unscale(center);
00990     const BigReal ucenter_x = ucenter.x;
00991     const BigReal ucenter_y = ucenter.y;
00992     const BigReal ucenter_z = ucenter.z;
00993     const int n = numAtoms;
00994     #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
00995       int n_16 = n;
00996       n_16 = (n + 15) & (~15);
00997       cudaAtomList.resize(n_16);
00998       CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
00999       mic_position_t *atom_x = ((mic_position_t*)ac) + (0 * n_16);
01000       mic_position_t *atom_y = ((mic_position_t*)ac) + (1 * n_16);
01001       mic_position_t *atom_z = ((mic_position_t*)ac) + (2 * n_16);
01002       mic_position_t *atom_q = ((mic_position_t*)ac) + (3 * n_16);
01003     #else
01004     cudaAtomList.resize(n);
01005     CudaAtom *ac = cudaAtomPtr = cudaAtomList.begin();
01006     #endif
01007     const FullAtom *a = atom.begin();
01008     for ( int k=0; k<n; ++k ) {
01009       #if defined(NAMD_MIC) && (MIC_HANDCODE_FORCE_SOA_VS_AOS == 0)
01010         int j = a[k].sortOrder;
01011         atom_x[k] = a[j].position.x - ucenter_x;
01012         atom_y[k] = a[j].position.y - ucenter_y;
01013         atom_z[k] = a[j].position.z - ucenter_z;
01014         atom_q[k] = charge_scaling * a[j].charge;
01015       #else
01016       int j = a[k].sortOrder;
01017       ac[k].x = a[j].position.x - ucenter_x;
01018       ac[k].y = a[j].position.y - ucenter_y;
01019       ac[k].z = a[j].position.z - ucenter_z;
01020       ac[k].q = charge_scaling * a[j].charge;
01021       #endif
01022     }
01023   }
01024 #else
01025   doMigration = doMigration && numNeighbors;
01026 #endif
01027   doMigration = doMigration || ! patchMapRead;
01028 
01029   doMigration = doMigration || doAtomUpdate;
01030   doAtomUpdate = false;
01031 
01032   // Workaround for oversize groups:
01033   // reset nonbondedGroupSize (ngs) before force calculation,
01034   // making sure that subset of hydrogen group starting with 
01035   // parent atom are all within 0.5 * hgroupCutoff.
01036   // XXX hydrogentGroupSize remains constant but is checked for nonzero
01037   // XXX should be skipped for CUDA, ngs not used by CUDA kernels
01038   // XXX should this also be skipped for KNL kernels?
01039   // ngs used by ComputeNonbondedBase.h - CPU nonbonded kernels
01040   // ngs used by ComputeGBIS.C - CPU GB nonbonded kernels
01041 #if ! defined(NAMD_CUDA)
01042   doGroupSizeCheck();
01043 #endif
01044 
01045   // Copy information needed by computes and proxys to Patch::p.
01046   // Resize only if atoms were migrated
01047   if (doMigration) {
01048     p.resize(numAtoms);
01049     pExt.resize(numAtoms);
01050   }
01051   CompAtom *p_i = p.begin();
01052   CompAtomExt *pExt_i = pExt.begin();
01053   FullAtom *a_i = atom.begin();
01054   int i; int n = numAtoms;
01055   for ( i=0; i<n; ++i ) { 
01056     p_i[i] = a_i[i]; 
01057     pExt_i[i] = a_i[i];
01058   }
01059 
01060   // Measure atom movement to test pairlist validity
01061   doPairlistCheck();
01062 
01063   if (flags.doMolly) mollyAverage();
01064   // BEGIN LA
01065   if (flags.doLoweAndersen) loweAndersenVelocities();
01066   // END LA
01067 
01068     if (flags.doGBIS) {
01069       //reset for next time step
01070       numGBISP1Arrived = 0;
01071       phase1BoxClosedCalled = false;
01072       numGBISP2Arrived = 0;
01073       phase2BoxClosedCalled = false;
01074       numGBISP3Arrived = 0;
01075       phase3BoxClosedCalled = false;
01076       if (doMigration || isNewProxyAdded)
01077         setGBISIntrinsicRadii();
01078     }
01079 
01080   if (flags.doLCPO) {
01081     if (doMigration || isNewProxyAdded) {
01082       setLcpoType();
01083     }
01084   }
01085 
01086   // Must Add Proxy Changes when migration completed!
01087   NodeIDList::iterator pli;
01088   int *pids = NULL;
01089   int pidsPreAllocated = 1;
01090   int npid;
01091   if (proxySendSpanning == 0) {
01092     npid = proxy.size();
01093     pids = new int[npid];
01094     pidsPreAllocated = 0;
01095     int *pidi = pids;
01096     int *pide = pids + proxy.size();
01097     int patchNodesLast =
01098       ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
01099     for ( pli = proxy.begin(); pli != proxy.end(); ++pli )
01100     {
01101       if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(*pli) ) {
01102         *(--pide) = *pli;
01103       } else {
01104         *(pidi++) = *pli;
01105       }
01106     }
01107   }
01108   else {
01109 #ifdef NODEAWARE_PROXY_SPANNINGTREE
01110     #ifdef USE_NODEPATCHMGR
01111     npid = numNodeChild;
01112     pids = nodeChildren;
01113     #else
01114     npid = nChild;
01115     pids = child;
01116     #endif
01117 #else
01118     npid = nChild;
01119     pidsPreAllocated = 0;
01120     pids = new int[proxySpanDim];
01121     for (int i=0; i<nChild; i++) pids[i] = child[i];
01122 #endif
01123   }
01124   if (npid) { //have proxies
01125     int seq = flags.sequence;
01126     int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
01127     //begin to prepare proxy msg and send it
01128     int pdMsgPLLen = p.size();
01129     int pdMsgAvgPLLen = 0;
01130     if(flags.doMolly) {
01131         pdMsgAvgPLLen = p_avg.size();
01132     }
01133     // BEGIN LA
01134     int pdMsgVLLen = 0;
01135     if (flags.doLoweAndersen) {
01136         pdMsgVLLen = v.size();
01137     }
01138     // END LA
01139 
01140     int intRadLen = 0;
01141     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
01142             intRadLen = numAtoms * 2;
01143     }
01144 
01145     //LCPO
01146     int lcpoTypeLen = 0;
01147     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
01148             lcpoTypeLen = numAtoms;
01149     }
01150 
01151     int pdMsgPLExtLen = 0;
01152     if(doMigration || isNewProxyAdded) {
01153         pdMsgPLExtLen = pExt.size();
01154     }
01155 
01156     int cudaAtomLen = 0;
01157 #ifdef NAMD_CUDA
01158     cudaAtomLen = numAtoms;
01159 #endif
01160 
01161     #ifdef NAMD_MIC
01162       #if MIC_HANDCODE_FORCE_SOA_VS_AOS != 0
01163         cudaAtomLen = numAtoms;
01164       #else
01165         cudaAtomLen = (numAtoms + 15) & (~15);
01166       #endif
01167     #endif
01168 
01169     ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgVLLen, intRadLen,
01170       lcpoTypeLen, pdMsgPLExtLen, cudaAtomLen, PRIORITY_SIZE) ProxyDataMsg; // BEGIN LA, END LA
01171 
01172     SET_PRIORITY(nmsg,seq,priority);
01173     nmsg->patch = patchID;
01174     nmsg->flags = flags;
01175     nmsg->plLen = pdMsgPLLen;                
01176     //copying data to the newly created msg
01177     memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
01178     nmsg->avgPlLen = pdMsgAvgPLLen;        
01179     if(flags.doMolly) {
01180         memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
01181     }
01182     // BEGIN LA
01183     nmsg->vlLen = pdMsgVLLen;
01184     if (flags.doLoweAndersen) {
01185         memcpy(nmsg->velocityList, v.begin(), sizeof(CompAtom)*pdMsgVLLen);
01186     }
01187     // END LA
01188 
01189     if (flags.doGBIS && (doMigration || isNewProxyAdded)) {
01190       for (int i = 0; i < numAtoms * 2; i++) {
01191         nmsg->intRadList[i] = intRad[i];
01192       }
01193     }
01194 
01195     if (flags.doLCPO && (doMigration || isNewProxyAdded)) {
01196       for (int i = 0; i < numAtoms; i++) {
01197         nmsg->lcpoTypeList[i] = lcpoType[i];
01198       }
01199     }
01200 
01201     nmsg->plExtLen = pdMsgPLExtLen;
01202     if(doMigration || isNewProxyAdded){     
01203         memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
01204     }
01205 
01206 // DMK
01207 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
01208     memcpy(nmsg->cudaAtomList, cudaAtomPtr, sizeof(CudaAtom)*cudaAtomLen);
01209 #endif
01210     
01211 #if NAMD_SeparateWaters != 0
01212     //DMK - Atom Separation (water vs. non-water)
01213     nmsg->numWaterAtoms = numWaterAtoms;
01214 #endif
01215 
01216 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
01217     nmsg->isFromImmMsgCall = 0;
01218 #endif
01219     
01220     #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
01221     DebugFileTrace *dft = DebugFileTrace::Object();
01222     dft->openTrace();
01223     dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
01224     for(int i=0; i<npid; i++) {
01225         dft->writeTrace("%d ", pids[i]);
01226     }
01227     dft->writeTrace("\n");
01228     dft->closeTrace();
01229     #endif
01230 
01231 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
01232     if (isProxyChanged || localphs == NULL)
01233     {
01234 //CmiPrintf("[%d] Build persistent: isProxyChanged: %d %p\n", CkMyPe(), isProxyChanged, localphs);
01235      //CmiAssert(isProxyChanged);
01236      if (nphs) {
01237        for (int i=0; i<nphs; i++) {
01238          CmiDestoryPersistent(localphs[i]);
01239        }
01240        delete [] localphs;
01241      }
01242      localphs = new PersistentHandle[npid];
01243      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;
01244      for (int i=0; i<npid; i++) {
01245 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR)
01246        if (proxySendSpanning)
01247            localphs[i] = CmiCreateNodePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
01248        else
01249 #endif
01250        localphs[i] = CmiCreatePersistent(pids[i], persist_size, sizeof(envelope)+sizeof(ProxyDataMsg));
01251      }
01252      nphs = npid;
01253     }
01254     CmiAssert(nphs == npid && localphs != NULL);
01255     CmiUsePersistentHandle(localphs, nphs);
01256 #endif
01257     if(doMigration || isNewProxyAdded) {
01258         ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
01259     }else{
01260         ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
01261     }
01262 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
01263     CmiUsePersistentHandle(NULL, 0);
01264 #endif
01265     isNewProxyAdded = 0;
01266   }
01267   isProxyChanged = 0;
01268   if(!pidsPreAllocated) delete [] pids;
01269   DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
01270 
01271 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
01272   positionPtrBegin = p.begin();
01273   positionPtrEnd = p.end();
01274 #endif
01275 
01276   if(flags.doMolly) {
01277       avgPositionPtrBegin = p_avg.begin();
01278       avgPositionPtrEnd = p_avg.end();
01279   }
01280   // BEGIN LA
01281   if (flags.doLoweAndersen) {
01282       velocityPtrBegin = v.begin();
01283       velocityPtrEnd = v.end();
01284   }
01285   // END LA
01286 
01287   Patch::positionsReady(doMigration);
01288 
01289   patchMapRead = 1;
01290 
01291   // gzheng
01292   Sync::Object()->PatchReady();
01293 
01294   NAMD_EVENT_STOP(1, NamdProfileEvent::POSITIONS_READY);
01295 
01296 }

void HomePatch::qmSwapAtoms (  ) 

Definition at line 866 of file HomePatch.C.

References ResizeArray< Elem >::begin(), CompAtom::charge, SortedArray< Elem >::find(), Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), CompAtomExt::id, lssSubs(), Node::molecule, LSSSubsDat::newCharge, LSSSubsDat::newID, LSSSubsDat::newVdWType, Patch::numAtoms, Node::Object(), SimParameters::PMEOn, Node::simParameters, and CompAtom::vdwType.

Referenced by positionsReady().

00867 {
00868     // This is used for LSS in QM/MM simulations.
00869     // Changes atom labels so that we effectively exchange solvent
00870     // molecules between classical and quantum modes.
00871     
00872     SimParameters *simParams = Node::Object()->simParameters;
00873     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
00874     const Real * const qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
00875     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
00876     Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
00877     
00878     ComputeQMMgr *mgrP = CProxy_ComputeQMMgr::ckLocalBranch(
00879             CkpvAccess(BOCclass_group).computeQMMgr) ;
00880     
00881     FullAtom *a_i = atom.begin();
00882     
00883     for (int i=0; i<numAtoms; ++i ) { 
00884         
00885         LSSSubsDat *subP = lssSubs(mgrP).find( LSSSubsDat(a_i[i].id) ) ;
00886         
00887         if ( subP != NULL ) {
00888             a_i[i].id = subP->newID ;
00889             a_i[i].vdwType = subP->newVdWType ;
00890             
00891             // If we are swappign a classical atom with a QM one, the charge
00892             // will need extra handling.
00893             if (qmAtomGroup[subP->newID] > 0 && simParams->PMEOn) {
00894                 // We make sure that the last atom charge calculated for the 
00895                 // QM atom being transfered here is available for PME
00896                 // in the next step.
00897                 
00898                 // Loops over all QM atoms (in all QM groups) comparing their 
00899                 // global indices (the sequential atom ID from NAMD).
00900                 for (int qmIter=0; qmIter<numQMAtms; qmIter++) {
00901                     
00902                     if (qmAtmIndx[qmIter] == subP->newID) {
00903                         qmAtmChrg[qmIter] = subP->newCharge;
00904                         break;
00905                     }
00906                     
00907                 }
00908                 
00909                 // For QM atoms, the charge in the full atom structure is zero.
00910                 // Electrostatic interactions between QM atoms and their 
00911                 // environment is handled in ComputeQM.
00912                 a_i[i].charge = 0;
00913             }
00914             else {
00915                 // If we are swappign a QM atom with a Classical one, only the charge
00916                 // in the full atomstructure needs updating, since it used to be zero.
00917                 a_i[i].charge = subP->newCharge ;
00918             }
00919         }
00920     }
00921     
00922     return;
00923 }

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

Definition at line 2293 of file HomePatch.C.

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

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

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

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

Definition at line 2470 of file HomePatch.C.

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

Referenced by rattle1().

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

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

Definition at line 2740 of file HomePatch.C.

References endi(), SimParameters::fixedAtomsOn, iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Node::Object(), outer(), SimParameters::rigidDie, SimParameters::rigidIter, SimParameters::rigidTol, settle2(), Node::simParameters, TIMEFACTOR, SimParameters::useSettle, WAT_SWM4, WAT_TIP4, SimParameters::watmodel, Vector::x, Vector::y, and Vector::z.

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

void HomePatch::receiveResult ( ProxyGBISP2ResultMsg msg  ) 

Definition at line 3198 of file HomePatch.C.

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

03198                                                        {
03199   ++numGBISP2Arrived;
03200   //accumulate dEda
03201   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
03202     dHdrPrefix[i] += msg->dEdaSum[i];
03203   }
03204   delete msg;
03205 
03206   if (flags.doNonbonded) {
03207     //awaken if phase 2 done
03208     if (phase2BoxClosedCalled == true &&
03209         numGBISP2Arrived==proxy.size() ) {
03210       sequencer->awaken();
03211     }
03212   } else {
03213     //awaken if all phases done on noWork step
03214     if (boxesOpen == 0 &&
03215         numGBISP1Arrived == proxy.size() &&
03216         numGBISP2Arrived == proxy.size() &&
03217         numGBISP3Arrived == proxy.size()) {
03218       sequencer->awaken();
03219     }
03220   }
03221 }

void HomePatch::receiveResult ( ProxyGBISP1ResultMsg msg  ) 

Definition at line 3173 of file HomePatch.C.

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

Referenced by ProxyMgr::recvResult().

03173                                                        {
03174   ++numGBISP1Arrived;
03175     for ( int i = 0; i < msg->psiSumLen; ++i ) {
03176       psiFin[i] += msg->psiSum[i];
03177     }
03178   delete msg;
03179 
03180   if (flags.doNonbonded) {
03181     //awaken if phase 1 done
03182     if (phase1BoxClosedCalled == true &&
03183         numGBISP1Arrived==proxy.size() ) {
03184         sequencer->awaken();
03185     }
03186   } else {
03187     //awaken if all phases done on noWork step
03188     if (boxesOpen == 0 &&
03189         numGBISP1Arrived == proxy.size() &&
03190         numGBISP2Arrived == proxy.size() &&
03191         numGBISP3Arrived == proxy.size()) {
03192       sequencer->awaken();
03193     }
03194   }
03195 }

void HomePatch::receiveResults ( ProxyCombinedResultRawMsg msg  ) 

Definition at line 837 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, Vector::x, Vector::y, and Vector::z.

00838 {
00839     numGBISP3Arrived++;
00840   DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodeSize<<")\n");
00841     Results *r = forceBox.clientOpen(msg->nodeSize);
00842       register char* isNonZero = msg->isForceNonZero;
00843       register Force* f_i = msg->forceArr;
00844       for ( int k = 0; k < Results::maxNumForces; ++k )
00845       {
00846         Force *f = r->f[k];
00847                 int nf = msg->flLen[k];
00848 #ifdef ARCH_POWERPC
00849 #pragma disjoint (*f_i, *f)
00850 #endif
00851         for (int count = 0; count < nf; count++) {
00852           if(*isNonZero){
00853                         f[count].x += f_i->x;
00854                         f[count].y += f_i->y;
00855                         f[count].z += f_i->z;
00856                         f_i++;
00857           }
00858           isNonZero++;
00859         }
00860       }
00861     forceBox.clientClose(msg->nodeSize);
00862 
00863     delete msg;
00864 }

void HomePatch::receiveResults ( ProxyResultMsg msg  ) 

Definition at line 820 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.

00820                                                   {
00821   numGBISP3Arrived++;
00822   DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00823   int n = msg->node;
00824   Results *r = forceBox.clientOpen();
00825   for ( int k = 0; k < Results::maxNumForces; ++k )
00826   {
00827     Force *f = r->f[k];
00828     register ForceList::iterator f_i, f_e;
00829     f_i = msg->forceList[k]->begin();
00830     f_e = msg->forceList[k]->end();
00831     for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
00832   }
00833   forceBox.clientClose();
00834   delete msg;
00835 }

void HomePatch::receiveResults ( ProxyResultVarsizeMsg msg  ) 

Definition at line 796 of file HomePatch.C.

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

Referenced by ProxyMgr::recvResults().

00796                                                          {
00797 
00798     numGBISP3Arrived++;
00799     DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00800     int n = msg->node;
00801     Results *r = forceBox.clientOpen();
00802 
00803     char *iszeroPtr = msg->isZero;
00804     Force *msgFPtr = msg->forceArr;
00805 
00806     for ( int k = 0; k < Results::maxNumForces; ++k )
00807     {
00808       Force *rfPtr = r->f[k];
00809       for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
00810           if((*iszeroPtr)!=1) {
00811               *rfPtr += *msgFPtr;
00812               msgFPtr++;
00813           }
00814       }      
00815     }
00816     forceBox.clientClose();
00817     delete msg;
00818 }

void HomePatch::recvCheckpointAck (  ) 

Definition at line 3480 of file HomePatch.C.

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

03480                                   {  // initiating replica
03481   CkpvAccess(_qd)->process();
03482 }

void HomePatch::recvCheckpointLoad ( CheckpointAtomsMsg msg  ) 

Definition at line 3427 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, SimParameters::scriptIntArg1, SimParameters::scriptStringArg1, PatchMgr::sendCheckpointStore(), Node::simParameters, and AtomMapper::unregisterIDsFullAtom().

Referenced by PatchMgr::recvCheckpointLoad().

03427                                                           {  // initiating replica
03428   SimParameters *simParams = Node::Object()->simParameters;
03429   const int remote = simParams->scriptIntArg1;
03430   const char *key = simParams->scriptStringArg1;
03431   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
03432     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
03433   }
03434   if ( msg->replica != remote ) {
03435     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
03436   }
03437   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03438     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
03439     strcpy(newmsg->key,key);
03440     newmsg->lattice = lattice;
03441     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
03442     newmsg->pid = patchID;
03443     newmsg->pe = CkMyPe();
03444     newmsg->replica = CmiMyPartition();
03445     newmsg->numAtoms = numAtoms;
03446     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03447     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
03448   }
03449   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03450     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03451     lattice = msg->lattice;
03452     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
03453     numAtoms = msg->numAtoms;
03454     atom.resize(numAtoms);
03455     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03456     doAtomUpdate = true;
03457     rattleListValid = false;
03458     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03459   }
03460   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
03461     recvCheckpointAck();
03462   }
03463   delete msg;
03464 }

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

Definition at line 3397 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().

03397                                                                                 {  // responding replica
03398   if ( task == SCRIPT_CHECKPOINT_FREE ) {
03399     if ( ! checkpoints.count(key) ) {
03400       NAMD_die("Unable to free checkpoint, requested key was never stored.");
03401     }
03402     delete checkpoints[key];
03403     checkpoints.erase(key);
03404     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
03405     return;
03406   }
03407   CheckpointAtomsMsg *msg;
03408   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
03409     if ( ! checkpoints.count(key) ) {
03410       NAMD_die("Unable to load checkpoint, requested key was never stored.");
03411     }
03412     checkpoint_t &cp = *checkpoints[key];
03413     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
03414     msg->lattice = cp.lattice;
03415     msg->berendsenPressure_count = cp.berendsenPressure_count;
03416     msg->numAtoms = cp.numAtoms;
03417     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
03418   } else {
03419     msg = new (0,1,0) CheckpointAtomsMsg;
03420   }
03421   msg->pid = patchID;
03422   msg->replica = CmiMyPartition();
03423   msg->pe = CkMyPe();
03424   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
03425 }

void HomePatch::recvCheckpointStore ( CheckpointAtomsMsg msg  ) 

Definition at line 3466 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().

03466                                                            {  // responding replica
03467   if ( ! checkpoints.count(msg->key) ) {
03468     checkpoints[msg->key] = new checkpoint_t;
03469   }
03470   checkpoint_t &cp = *checkpoints[msg->key];
03471   cp.lattice = msg->lattice;
03472   cp.berendsenPressure_count = msg->berendsenPressure_count;
03473   cp.numAtoms = msg->numAtoms;
03474   cp.atoms.resize(cp.numAtoms);
03475   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
03476   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
03477   delete msg;
03478 }

void HomePatch::recvExchangeMsg ( ExchangeAtomsMsg msg  ) 

Definition at line 3517 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().

03517                                                      {
03518   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
03519   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03520   lattice = msg->lattice;
03521   numAtoms = msg->numAtoms;
03522   atom.resize(numAtoms);
03523   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03524   delete msg;
03525   CkpvAccess(_qd)->process();
03526   doAtomUpdate = true;
03527   rattleListValid = false;
03528   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03529 }

void HomePatch::recvExchangeReq ( int  req  ) 

Definition at line 3506 of file HomePatch.C.

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

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

03506                                        {
03507   exchange_req = req;
03508   if ( exchange_msg ) {
03509     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
03510     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
03511     exchange_msg = 0;
03512     exchange_req = -1;
03513     CkpvAccess(_qd)->process();
03514   }
03515 }

void HomePatch::recvNodeAwareSpanningTree ( ProxyNodeAwareSpanningTreeMsg msg  ) 

Definition at line 630 of file HomePatch.C.

Referenced by ProxyMgr::recvNodeAwareSpanningTreeOnHomePatch().

00630 {}

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

Definition at line 636 of file HomePatch.C.

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

Referenced by ProxyMgr::recvSpanningTreeOnHomePatch().

00637 {
00638   int i;
00639   nChild=0;
00640   tree.resize(n);
00641   for (i=0; i<n; i++) {
00642     tree[i] = t[i];
00643   }
00644 
00645   for (i=1; i<=proxySpanDim; i++) {
00646     if (tree.size() <= i) break;
00647     child[i-1] = tree[i];
00648     nChild++;
00649   }
00650 
00651 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00652   isProxyChanged = 1;
00653 #endif
00654 
00655   // send down to children
00656   sendSpanningTree();
00657 }

void HomePatch::registerProxy ( RegisterProxyMsg msg  ) 

Definition at line 402 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().

00402                                                    {
00403   DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
00404   proxy.add(msg->node);
00405   forceBox.clientAdd();
00406 
00407   isNewProxyAdded = 1;
00408 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00409   isProxyChanged = 1;
00410 #endif
00411 
00412   Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
00413   delete msg;
00414 }

void HomePatch::replaceForces ( ExtForce f  ) 

Definition at line 1298 of file HomePatch.C.

01299 {
01300   replacementForces = f;
01301 }

void HomePatch::revert ( void   ) 

Definition at line 3371 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().

03371                            {
03372   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03373 
03374   atom.copy(checkpoint_atom);
03375   numAtoms = atom.size();
03376   lattice = checkpoint_lattice;
03377 
03378   doAtomUpdate = true;
03379   rattleListValid = false;
03380 
03381   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03382 
03383   // DMK - Atom Separation (water vs. non-water)
03384   #if NAMD_SeparateWaters != 0
03385     numWaterAtoms = checkpoint_numWaterAtoms;
03386   #endif
03387 }

void HomePatch::runSequencer ( void   ) 

Definition at line 269 of file HomePatch.C.

References Sequencer::run().

Referenced by Node::run().

00270 { sequencer->run(); }

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

Definition at line 1303 of file HomePatch.C.

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

Referenced by Sequencer::saveForce().

01304 {
01305   f_saved[ftag].resize(numAtoms);
01306   for ( int i = 0; i < numAtoms; ++i )
01307   {
01308     f_saved[ftag][i] = f[ftag][i];
01309   }
01310 }

void HomePatch::sendNodeAwareSpanningTree (  ) 

Definition at line 631 of file HomePatch.C.

00631 {}

void HomePatch::sendProxies (  ) 

Definition at line 468 of file HomePatch.C.

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

Referenced by ProxyMgr::buildProxySpanningTree2().

00469 {
00470 #if USE_NODEPATCHMGR
00471         CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
00472         NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
00473         npm->sendProxyList(patchID, proxy.begin(), proxy.size());
00474 #else
00475         ProxyMgr::Object()->sendProxies(patchID, proxy.begin(), proxy.size());
00476 #endif
00477 }

void HomePatch::sendSpanningTree (  ) 

Definition at line 659 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().

00660 {
00661   if (tree.size() == 0) return;
00662   ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
00663   msg->patch = patchID;
00664   msg->node = CkMyPe();
00665   msg->tree.copy(tree);  // copy data for thread safety
00666   ProxyMgr::Object()->sendSpanningTree(msg);  
00667 }

void HomePatch::setGBISIntrinsicRadii (  ) 

Definition at line 3043 of file HomePatch.C.

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

Referenced by positionsReady().

03043                                       {
03044   intRad.resize(numAtoms*2);
03045   intRad.setall(0);
03046   Molecule *mol = Node::Object()->molecule;
03047   SimParameters *simParams = Node::Object()->simParameters;
03048   Real offset = simParams->coulomb_radius_offset;
03049   for (int i = 0; i < numAtoms; i++) {
03050     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
03051     Real screen = MassToScreen(atom[i].mass);//same
03052     intRad[2*i+0] = rad - offset;//r0
03053     intRad[2*i+1] = screen*(rad - offset);//s0
03054   }
03055 }

void HomePatch::setLcpoType (  ) 

Definition at line 3032 of file HomePatch.C.

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

Referenced by positionsReady().

03032                             {
03033   Molecule *mol = Node::Object()->molecule;
03034   const int *lcpoParamType = mol->getLcpoParamType();
03035 
03036   lcpoType.resize(numAtoms);
03037   for (int i = 0; i < numAtoms; i++) {
03038     lcpoType[i] = lcpoParamType[pExt[i].id];
03039   }
03040 }

void HomePatch::submitLoadStats ( int  timestep  ) 

Definition at line 3531 of file HomePatch.C.

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

Referenced by Sequencer::rebalanceLoad().

03532 {
03533   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
03534 }

void HomePatch::unregisterProxy ( UnregisterProxyMsg msg  ) 

Definition at line 416 of file HomePatch.C.

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

Referenced by ProxyMgr::recvUnregisterProxy().

00416                                                        {
00417 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00418   isProxyChanged = 1;
00419 #endif
00420   int n = msg->node;
00421   NodeID *pe = proxy.begin();
00422   for ( ; *pe != n; ++pe );
00423   forceBox.clientRemove();
00424   proxy.del(pe - proxy.begin());
00425   delete msg;
00426 }

void HomePatch::useSequencer ( Sequencer sequencerPtr  ) 

Definition at line 265 of file HomePatch.C.

00266 { sequencer=sequencerPtr; }


Friends And Related Function Documentation

friend class ComputeGlobal [friend]

Definition at line 276 of file HomePatch.h.

friend class PatchMgr [friend]

Definition at line 274 of file HomePatch.h.

friend class Sequencer [friend]

Definition at line 275 of file HomePatch.h.


Member Data Documentation

Definition at line 428 of file HomePatch.h.

Referenced by exchangeCheckpoint(), and recvCheckpointLoad().

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

Definition at line 435 of file HomePatch.h.

Referenced by recvCheckpointReq(), and recvCheckpointStore().

Definition at line 441 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

Definition at line 444 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

Definition at line 443 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

Definition at line 442 of file HomePatch.h.

Referenced by exchangeAtoms().

int HomePatch::inMigration [protected]

Definition at line 492 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

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

Definition at line 494 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<int> HomePatch::noconstList

Definition at line 383 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

int HomePatch::numMlBuf [protected]

Definition at line 493 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<Vector> HomePatch::posNew

Definition at line 389 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

Definition at line 381 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

Definition at line 385 of file HomePatch.h.

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

Definition at line 382 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

std::vector<int> HomePatch::settleList

Definition at line 380 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

std::vector<Vector> HomePatch::velNew

Definition at line 388 of file HomePatch.h.

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


The documentation for this class was generated from the following files:

Generated on 19 Oct 2019 for NAMD by  doxygen 1.6.1