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 1900 of file HomePatch.C.

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

Referenced by Sequencer::addForceToMomentum().

01905       {
01906   SimParameters *simParams = Node::Object()->simParameters;
01907 
01908   if ( simParams->fixedAtomsOn ) {
01909     for ( int i = 0; i < num_atoms; ++i ) {
01910       if ( atom_arr[i].atomFixed ) {
01911         atom_arr[i].velocity = 0;
01912       } else {
01913         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01914         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01915         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01916         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01917       }
01918     }
01919   } else {
01920     for ( int i = 0; i < num_atoms; ++i ) {
01921       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01922       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01923       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01924       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01925     }
01926   }
01927 }

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 1929 of file HomePatch.C.

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

Referenced by Sequencer::addForceToMomentum3().

01938       {
01939   SimParameters *simParams = Node::Object()->simParameters;
01940 
01941   if ( simParams->fixedAtomsOn ) {
01942     for ( int i = 0; i < num_atoms; ++i ) {
01943       if ( atom_arr[i].atomFixed ) {
01944         atom_arr[i].velocity = 0;
01945       } else {
01946         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01947         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01948             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01949         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01950             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01951         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01952             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01953       }
01954     }
01955   } else {
01956     for ( int i = 0; i < num_atoms; ++i ) {
01957       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01958       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01959           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01960       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01961           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01962       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01963           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01964     }
01965   }
01966 }

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

Definition at line 2339 of file HomePatch.C.

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

Referenced by rattle1().

02339                                                               {
02340   for (int ig = 0; ig < numAtoms; ++ig ) {
02341     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
02342     Tensor vir = outer(df, atom[ig].position);
02343     wc += vir;
02344     f[Results::normal][ig] += df;
02345     atom[ig].velocity = velNew[ig];
02346   }
02347 }

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

Definition at line 1968 of file HomePatch.C.

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

Referenced by Sequencer::addVelocityToPosition().

01972       {
01973   SimParameters *simParams = Node::Object()->simParameters;
01974   if ( simParams->fixedAtomsOn ) {
01975     for ( int i = 0; i < num_atoms; ++i ) {
01976       if ( ! atom_arr[i].atomFixed ) {
01977         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01978         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01979         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01980       }
01981     }
01982   } else {
01983     for ( int i = 0; i < num_atoms; ++i ) {
01984       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01985       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01986       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01987     }
01988   }
01989 }

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 2197 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().

02197                                 {
02198 
02199   SimParameters *simParams = Node::Object()->simParameters;
02200   const int fixedAtomsOn = simParams->fixedAtomsOn;
02201   const int useSettle = simParams->useSettle;
02202 
02203   // Re-size to containt numAtoms elements
02204   velNew.resize(numAtoms);
02205   posNew.resize(numAtoms);
02206 
02207   // Size of a hydrogen group for water
02208   int wathgsize = 3;
02209   int watmodel = simParams->watmodel;
02210   if (watmodel == WAT_TIP4) wathgsize = 4;
02211   else if (watmodel == WAT_SWM4) wathgsize = 5;
02212 
02213   // Initialize the settle algorithm with water parameters
02214   // settle1() assumes all waters are identical,
02215   // and will generate bad results if they are not.
02216   // XXX this will move to Molecule::build_atom_status when that 
02217   // version is debugged
02218   if ( ! settle_initialized ) {
02219     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02220       // find a water
02221       if (atom[ig].rigidBondLength > 0) {
02222         int oatm;
02223         if (simParams->watmodel == WAT_SWM4) {
02224           oatm = ig+3;  // skip over Drude and Lonepair
02225           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02226           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02227         }
02228         else {
02229           oatm = ig+1;
02230           // Avoid using the Om site to set this by mistake
02231           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02232             oatm += 1;
02233           }
02234         }
02235 
02236         // initialize settle water parameters
02237         settle1init(atom[ig].mass, atom[oatm].mass, 
02238                     atom[ig].rigidBondLength,
02239                     atom[oatm].rigidBondLength,
02240                     settle_mOrmT, settle_mHrmT, settle_ra,
02241                     settle_rb, settle_rc, settle_rra);
02242         settle_initialized = 1;
02243         break; // done with init
02244       }
02245     }
02246   }
02247   
02248   Vector ref[10];
02249   BigReal rmass[10];
02250   BigReal dsq[10];
02251   int fixed[10];
02252   int ial[10];
02253   int ibl[10];
02254 
02255   int numSettle = 0;
02256   int numRattle = 0;
02257   int posRattleParam = 0;
02258 
02259   settleList.clear();
02260   rattleList.clear();
02261   noconstList.clear();
02262   rattleParam.clear();
02263 
02264   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02265     int hgs = atom[ig].hydrogenGroupSize;
02266     if ( hgs == 1 ) {
02267       // only one atom in group
02268       noconstList.push_back(ig);
02269       continue;
02270     }
02271     int anyfixed = 0;
02272     for (int i = 0; i < hgs; ++i ) {
02273       ref[i] = atom[ig+i].position;
02274       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02275       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02276       if ( fixed[i] ) {
02277         anyfixed = 1;
02278         rmass[i] = 0.;
02279       }
02280     }
02281     int icnt = 0;
02282     BigReal tmp = atom[ig].rigidBondLength;
02283     if (tmp > 0.0) {  // for water
02284       if (hgs != wathgsize) {
02285         char errmsg[256];
02286         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02287                          "but the specified water model requires %d atoms.\n",
02288                          atom[ig].id+1, hgs, wathgsize);
02289         NAMD_die(errmsg);
02290       }
02291       // Use SETTLE for water unless some of the water atoms are fixed,
02292       if (useSettle && !anyfixed) {
02293         // Store to Settle -list
02294         settleList.push_back(ig);
02295         continue;
02296       }
02297       if ( !(fixed[1] && fixed[2]) ) {
02298         dsq[icnt] = tmp * tmp;
02299         ial[icnt] = 1;
02300         ibl[icnt] = 2;
02301         ++icnt;
02302       }
02303     } // if (tmp > 0.0)
02304     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02305       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02306         if ( !(fixed[0] && fixed[i]) ) {
02307           dsq[icnt] = tmp * tmp;
02308           ial[icnt] = 0;
02309           ibl[icnt] = i;
02310           ++icnt;
02311         }
02312       }
02313     }
02314     if ( icnt == 0 ) {
02315       // no constraints
02316       noconstList.push_back(ig);
02317       continue;  
02318     }
02319     // Store to Rattle -list
02320     RattleList rattleListElem;
02321     rattleListElem.ig  = ig;
02322     rattleListElem.icnt = icnt;
02323     rattleList.push_back(rattleListElem);
02324     for (int i = 0; i < icnt; ++i ) {
02325       int a = ial[i];
02326       int b = ibl[i];
02327       RattleParam rattleParamElem;
02328       rattleParamElem.ia = a;
02329       rattleParamElem.ib = b;
02330       rattleParamElem.dsq = dsq[i];
02331       rattleParamElem.rma = rmass[a];
02332       rattleParamElem.rmb = rmass[b];
02333       rattleParam.push_back(rattleParamElem);
02334     }
02335   }
02336 
02337 }

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 3417 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

03417                                {
03418   checkpoint_atom.copy(atom);
03419   checkpoint_lattice = lattice;
03420 
03421   // DMK - Atom Separation (water vs. non-water)
03422   #if NAMD_SeparateWaters != 0
03423     checkpoint_numWaterAtoms = numWaterAtoms;
03424   #endif
03425 }

void HomePatch::depositMigration ( MigrateAtomsMsg msg  ) 

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

03965 {
03966 
03967   if (!inMigration) { // We have to buffer changes due to migration
03968                       // until our patch is in migration mode
03969     msgbuf[numMlBuf++] = msg;
03970     return;
03971   } 
03972 
03973 
03974   // DMK - Atom Separation (water vs. non-water)
03975   #if NAMD_SeparateWaters != 0
03976 
03977 
03978     // Merge the incoming list of atoms with the current list of
03979     //   atoms.  Note that mergeSeparatedAtomList() will apply any
03980     //   required transformations to the incoming atoms as it is
03981     //   separating them.
03982     mergeAtomList(msg->migrationList);
03983 
03984 
03985   #else
03986 
03987     // Merge the incoming list of atoms with the current list of
03988     // atoms.  Apply transformations to the incoming atoms as they are
03989     // added to this patch's list.
03990     {
03991       MigrationList &migrationList = msg->migrationList;
03992       MigrationList::iterator mi;
03993       Transform mother_transform;
03994       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
03995         DebugM(1,"Migrating atom " << mi->id << " to patch "
03996                   << patchID << " with position " << mi->position << "\n"); 
03997         if ( mi->migrationGroupSize ) {
03998           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
03999             Position pos = mi->position;
04000             int mgs = mi->migrationGroupSize;
04001             int c = 1;
04002             for ( int j=mi->hydrogenGroupSize; j<mgs;
04003                                 j+=(mi+j)->hydrogenGroupSize ) {
04004               pos += (mi+j)->position;
04005               ++c;
04006             }
04007             pos *= 1./c;
04008             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
04009             mother_transform = mi->transform;
04010             pos = lattice.nearest(pos,center,&mother_transform);
04011             mi->position = lattice.reverse_transform(mi->position,mi->transform);
04012             mi->position = lattice.apply_transform(mi->position,mother_transform);
04013             mi->transform = mother_transform;
04014           } else {
04015             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
04016             mother_transform = mi->transform;
04017           }
04018         } else {
04019           mi->position = lattice.reverse_transform(mi->position,mi->transform);
04020           mi->position = lattice.apply_transform(mi->position,mother_transform);
04021           mi->transform = mother_transform;
04022         }
04023         atom.add(*mi);
04024       }
04025     }
04026 
04027 
04028   #endif // if (NAMD_SeparateWaters != 0)
04029 
04030 
04031   numAtoms = atom.size();
04032   delete msg;
04033 
04034   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
04035   if (!--patchMigrationCounter) {
04036     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
04037     allMigrationIn = true;
04038     patchMigrationCounter = numNeighbors;
04039     if (migrationSuspended) {
04040       DebugM(3,"patch "<<patchID<<" is being awakened\n");
04041       migrationSuspended = false;
04042       sequencer->awaken();
04043       return;
04044     }
04045     else {
04046        DebugM(3,"patch "<<patchID<<" was not suspended\n");
04047     }
04048   }
04049 }

void HomePatch::doAtomMigration (  )  [protected]

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

03809 {
03810   int i;
03811 
03812   for (i=0; i<numNeighbors; i++) {
03813     realInfo[i].mList.resize(0);
03814   }
03815 
03816   // Purge the AtomMap
03817   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03818 
03819   // Determine atoms that need to migrate
03820 
03821   BigReal minx = min.x;
03822   BigReal miny = min.y;
03823   BigReal minz = min.z;
03824   BigReal maxx = max.x;
03825   BigReal maxy = max.y;
03826   BigReal maxz = max.z;
03827 
03828   int xdev, ydev, zdev;
03829   int delnum = 0;
03830 
03831   FullAtomList::iterator atom_i = atom.begin();
03832   FullAtomList::iterator atom_e = atom.end();
03833 
03834   // DMK - Atom Separation (water vs. non-water)
03835   #if NAMD_SeparateWaters != 0
03836     FullAtomList::iterator atom_first = atom_i;
03837     int numLostWaterAtoms = 0;
03838   #endif
03839 
03840   while ( atom_i != atom_e ) {
03841     // Even though this code iterates through all atoms successively
03842     // it moves entire hydrogen/migration groups as follows:
03843     // Only the parent atom of the hydrogen/migration group has
03844     // nonzero migrationGroupSize.  Values determined for xdev,ydev,zdev
03845     // will persist through the remaining group members so that each
03846     // following atom will again be added to the same mList.
03847     if ( atom_i->migrationGroupSize ) {
03848       Position pos = atom_i->position;
03849       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
03850         // If there are multiple hydrogen groups in a migration group
03851         // (e.g. for supporting lone pairs)
03852         // the following code takes the average position (midpoint)
03853         // of their parents.
03854         int mgs = atom_i->migrationGroupSize;
03855         int c = 1;
03856         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
03857                                 j+=(atom_i+j)->hydrogenGroupSize ) {
03858           pos += (atom_i+j)->position;
03859           ++c;
03860         }
03861         pos *= 1./c;
03862         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
03863       }
03864 
03865       // Scaling the position below transforms space within patch from
03866       // what could have been a rotated parallelepiped into
03867       // orthogonal coordinates, where we can use minmax comparison
03868       // to detect which of our nearest neighbors this
03869       // parent atom might have entered.
03870       ScaledPosition s = lattice.scale(pos);
03871 
03872       // check if atom is within bounds
03873       if (s.x < minx) xdev = 0;
03874       else if (maxx <= s.x) xdev = 2;
03875       else xdev = 1;
03876 
03877       if (s.y < miny) ydev = 0;
03878       else if (maxy <= s.y) ydev = 2;
03879       else ydev = 1;
03880 
03881       if (s.z < minz) zdev = 0;
03882       else if (maxz <= s.z) zdev = 2;
03883       else zdev = 1;
03884 
03885     }
03886 
03887     if (mInfo[xdev][ydev][zdev]) { // process atom for migration
03888                                     // Don't migrate if destination is myself
03889 
03890       // See if we have a migration list already
03891       MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
03892       DebugM(3,"Migrating atom " << atom_i->id << " from patch "
03893                 << patchID << " with position " << atom_i->position << "\n");
03894       mCur.add(*atom_i);
03895 
03896       ++delnum;
03897 
03898 
03899       // DMK - Atom Separation (water vs. non-water)
03900       #if NAMD_SeparateWaters != 0
03901         // Check to see if this atom is part of a water molecule.  If
03902         //   so, numWaterAtoms needs to be adjusted to refect the lost of
03903         //   this atom.
03904         // NOTE: The atom separation code assumes that if the oxygen
03905         //   atom of the hydrogen group making up the water molecule is
03906         //   migrated to another HomePatch, the hydrogens will also
03907         //   move!!!
03908         int atomIndex = atom_i - atom_first;
03909         if (atomIndex < numWaterAtoms)
03910           numLostWaterAtoms++;
03911       #endif
03912 
03913 
03914     } else {
03915       // By keeping track of delnum total being deleted from FullAtomList
03916       // the else clause allows us to fill holes as we visit each atom.
03917 
03918       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
03919 
03920     }
03921 
03922     ++atom_i;
03923   }
03924 
03925   // DMK - Atom Separation (water vs. non-water)
03926   #if NAMD_SeparateWaters != 0
03927     numWaterAtoms -= numLostWaterAtoms;
03928   #endif
03929 
03930 
03931   int delpos = numAtoms - delnum;
03932   DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
03933   atom.del(delpos,delnum);
03934 
03935   numAtoms = atom.size();
03936 
03937   PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
03938 
03939   // signal depositMigration() that we are inMigration mode
03940   inMigration = true;
03941 
03942   // Drain the migration message buffer
03943   for (i=0; i<numMlBuf; i++) {
03944      DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
03945      depositMigration(msgbuf[i]);
03946   }
03947   numMlBuf = 0;
03948 
03949   NAMD_EVENT_STOP(1, NamdProfileEvent::ATOM_MIGRATIONS);
03950 
03951   if (!allMigrationIn) {
03952     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
03953     migrationSuspended = true;
03954     sequencer->suspend();
03955     migrationSuspended = false;
03956   }
03957   allMigrationIn = false;
03958 
03959   inMigration = false;
03960   marginViolations = 0;
03961 }

void HomePatch::doGroupSizeCheck (  )  [protected]

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

03688 {
03689   if ( ! flags.doNonbonded ) return;
03690 
03691   SimParameters *simParams = Node::Object()->simParameters;
03692   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
03693   BigReal maxrad2 = 0.;
03694 
03695   FullAtomList::iterator p_i = atom.begin();
03696   FullAtomList::iterator p_e = atom.end();
03697 
03698   while ( p_i != p_e ) {
03699     const int hgs = p_i->hydrogenGroupSize;
03700     if ( ! hgs ) break;  // avoid infinite loop on bug
03701     int ngs = hgs;
03702     if ( ngs > 5 ) ngs = 5;  // XXX why? limit to at most 5 atoms per group
03703     BigReal x = p_i->position.x;
03704     BigReal y = p_i->position.y;
03705     BigReal z = p_i->position.z;
03706     int i;
03707     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
03708       p_i[i].nonbondedGroupSize = 0;
03709       BigReal dx = p_i[i].position.x - x;
03710       BigReal dy = p_i[i].position.y - y;
03711       BigReal dz = p_i[i].position.z - z;
03712       BigReal r2 = dx * dx + dy * dy + dz * dz;
03713       if ( r2 > hgcut ) break;
03714       else if ( r2 > maxrad2 ) maxrad2 = r2;
03715     }
03716     p_i->nonbondedGroupSize = i;
03717     for ( ; i < hgs; ++i ) {
03718       p_i[i].nonbondedGroupSize = 1;
03719     }
03720     p_i += hgs;
03721   }
03722 
03723   if ( p_i != p_e ) {
03724     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
03725   }
03726 
03727   flags.maxGroupRadius = sqrt(maxrad2);
03728 
03729 }

void HomePatch::doMarginCheck (  )  [protected]

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

03732 {
03733   SimParameters *simParams = Node::Object()->simParameters;
03734 
03735   BigReal sysdima = lattice.a_r().unit() * lattice.a();
03736   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
03737   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
03738 
03739   BigReal minSize = simParams->patchDimension - simParams->margin;
03740 
03741   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
03742        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
03743        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
03744 
03745     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03746       "Possible solutions are to restart from a recent checkpoint,\n"
03747       "increase margin, or disable useFlexibleCell for liquid simulation.");
03748   }
03749 
03750   BigReal cutoff = simParams->cutoff;
03751 
03752   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
03753   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
03754   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
03755 
03756   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
03757     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03758       "There are probably many margin violations already reported.\n"
03759       "Possible solutions are to restart from a recent checkpoint,\n"
03760       "increase margin, or disable useFlexibleCell for liquid simulation.");
03761   }
03762 
03763   BigReal minx = min.x - margina;
03764   BigReal miny = min.y - marginb;
03765   BigReal minz = min.z - marginc;
03766   BigReal maxx = max.x + margina;
03767   BigReal maxy = max.y + marginb;
03768   BigReal maxz = max.z + marginc;
03769 
03770   int xdev, ydev, zdev;
03771   int problemCount = 0;
03772 
03773   FullAtomList::iterator p_i = atom.begin();
03774   FullAtomList::iterator p_e = atom.end();
03775   for ( ; p_i != p_e; ++p_i ) {
03776 
03777     ScaledPosition s = lattice.scale(p_i->position);
03778 
03779     // check if atom is within bounds
03780     if (s.x < minx) xdev = 0;
03781     else if (maxx <= s.x) xdev = 2; 
03782     else xdev = 1;
03783 
03784     if (s.y < miny) ydev = 0;
03785     else if (maxy <= s.y) ydev = 2; 
03786     else ydev = 1;
03787 
03788     if (s.z < minz) zdev = 0;
03789     else if (maxz <= s.z) zdev = 2; 
03790     else zdev = 1;
03791 
03792     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
03793         ++problemCount;
03794     }
03795 
03796   }
03797 
03798   marginViolations = problemCount;
03799   // if ( problemCount ) {
03800   //     iout << iERROR <<
03801   //       "Found " << problemCount << " margin violations!\n" << endi;
03802   // } 
03803 
03804 }

void HomePatch::doPairlistCheck (  )  [protected]

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

03597 {
03598 #if 0
03599 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
03600   char dpcbuf[32];
03601   sprintf(dpcbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::DO_PAIRLIST_CHECK], this->getPatchID());
03602   NAMD_EVENT_START_EX(1, NamdProfileEvent::DO_PAIRLIST_CHECK, dpcbuf);
03603 #endif
03604 #endif
03605 
03606   SimParameters *simParams = Node::Object()->simParameters;
03607 
03608   if ( numAtoms == 0 || ! flags.usePairlists ) {
03609     flags.pairlistTolerance = 0.;
03610     flags.maxAtomMovement = 99999.;
03611 #if 0
03612     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03613 #endif
03614     return;
03615   }
03616 
03617   int i; int n = numAtoms;
03618   CompAtom *p_i = p.begin();
03619 
03620   if ( flags.savePairlists ) {
03621     flags.pairlistTolerance = doPairlistCheck_newTolerance;
03622     flags.maxAtomMovement = 0.;
03623     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
03624     doPairlistCheck_lattice = lattice;
03625     doPairlistCheck_positions.resize(numAtoms);
03626     CompAtom *psave_i = doPairlistCheck_positions.begin();
03627     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
03628 #if 0
03629     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03630 #endif
03631     return;
03632   }
03633 
03634   Lattice &lattice_old = doPairlistCheck_lattice;
03635   Position center_cur = lattice.unscale(center);
03636   Position center_old = lattice_old.unscale(center);
03637   Vector center_delta = center_cur - center_old;
03638   
03639   // find max deviation to corner (any neighbor shares a corner)
03640   BigReal max_cd = 0.;
03641   for ( i=0; i<2; ++i ) {
03642     for ( int j=0; j<2; ++j ) {
03643       for ( int k=0; k<2; ++k ) {
03644         ScaledPosition corner(  i ? min.x : max.x ,
03645                                 j ? min.y : max.y ,
03646                                 k ? min.z : max.z );
03647         Vector corner_delta =
03648                 lattice.unscale(corner) - lattice_old.unscale(corner);
03649         corner_delta -= center_delta;
03650         BigReal cd = corner_delta.length2();
03651         if ( cd > max_cd ) max_cd = cd;
03652       }
03653     }
03654   }
03655   max_cd = sqrt(max_cd);
03656 
03657   // find max deviation of atoms relative to center
03658   BigReal max_pd = 0.;
03659   CompAtom *p_old_i = doPairlistCheck_positions.begin();
03660   for ( i=0; i<n; ++i ) {
03661     Vector p_delta = p_i[i].position - p_old_i[i].position;
03662     p_delta -= center_delta;
03663     BigReal pd = p_delta.length2();
03664     if ( pd > max_pd ) max_pd = pd;
03665   }
03666   max_pd = sqrt(max_pd);
03667 
03668   BigReal max_tol = max_pd + max_cd;
03669 
03670   flags.maxAtomMovement = max_tol;
03671 
03672   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
03673 
03674   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
03675                                 doPairlistCheck_newTolerance ) ) {
03676     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
03677   }
03678 
03679   if ( max_tol > doPairlistCheck_newTolerance ) {
03680     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
03681   }
03682 #if 0
03683   NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03684 #endif
03685 }

void HomePatch::exchangeAtoms ( int  scriptTask  ) 

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

03541                                             {
03542   SimParameters *simParams = Node::Object()->simParameters;
03543   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
03544   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
03545     exchange_dst = (int) simParams->scriptArg1;
03546     // create and save outgoing message
03547     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
03548     exchange_msg->lattice = lattice;
03549     exchange_msg->pid = patchID;
03550     exchange_msg->numAtoms = numAtoms;
03551     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03552     if ( exchange_req >= 0 ) {
03553       recvExchangeReq(exchange_req);
03554     }
03555   }
03556   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
03557     exchange_src = (int) simParams->scriptArg2;
03558     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
03559   }
03560 }

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

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

03445                                                            {  // initiating replica
03446   SimParameters *simParams = Node::Object()->simParameters;
03447   checkpoint_task = scriptTask;
03448   const int remote = simParams->scriptIntArg1;
03449   const char *key = simParams->scriptStringArg1;
03450   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
03451 }

void HomePatch::gbisComputeAfterP1 (  ) 

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

03114                                    {
03115 
03116   SimParameters *simParams = Node::Object()->simParameters;
03117   BigReal alphaMax = simParams->alpha_max;
03118   BigReal delta = simParams->gbis_delta;
03119   BigReal beta = simParams->gbis_beta;
03120   BigReal gamma = simParams->gbis_gamma;
03121   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03122 
03123   BigReal rhoi;
03124   BigReal rhoi0;
03125   //calculate bornRad from psiSum
03126   for (int i = 0; i < numAtoms; i++) {
03127     rhoi0 = intRad[2*i];
03128     rhoi = rhoi0+coulomb_radius_offset;
03129     psiFin[i] += psiSum[i];
03130     psiFin[i] *= rhoi0;
03131     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
03132     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
03133 #ifdef PRINT_COMP
03134     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
03135 #endif
03136   }
03137 
03138   gbisP2Ready();
03139 }

void HomePatch::gbisComputeAfterP2 (  ) 

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

03142                                    {
03143 
03144   SimParameters *simParams = Node::Object()->simParameters;
03145   BigReal delta = simParams->gbis_delta;
03146   BigReal beta = simParams->gbis_beta;
03147   BigReal gamma = simParams->gbis_gamma;
03148   BigReal epsilon_s = simParams->solvent_dielectric;
03149   BigReal epsilon_p = simParams->dielectric;
03150   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
03151   BigReal epsilon_p_i = 1/simParams->dielectric;
03152   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03153   BigReal kappa = simParams->kappa;
03154   BigReal fij, expkappa, Dij, dEdai, dedasum;
03155   BigReal rhoi, rhoi0, psii, nbetapsi;
03156   BigReal gammapsi2, tanhi, daidr;
03157   for (int i = 0; i < numAtoms; i++) {
03158     //add diagonal dEda term
03159     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
03160     fij = bornRad[i];//inf
03161     expkappa = exp(-kappa*fij);//0
03162     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
03163     //calculate dHij prefix
03164     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
03165                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
03166     dHdrPrefix[i] += dEdai;
03167     dedasum = dHdrPrefix[i];
03168 
03169     rhoi0 = intRad[2*i];
03170     rhoi = rhoi0+coulomb_radius_offset;
03171     psii = psiFin[i];
03172     nbetapsi = -beta*psii;
03173     gammapsi2 = gamma*psii*psii;
03174     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
03175     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
03176            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
03177     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
03178 #ifdef PRINT_COMP
03179     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
03180 #endif
03181   }
03182   gbisP3Ready();
03183 }

void HomePatch::gbisP2Ready (  ) 

Reimplemented from Patch.

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

03186                             {
03187   if (proxy.size() > 0) {
03188     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03189     for (int i = 0; i < proxy.size(); i++) {
03190       int node = proxy[i];
03191       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
03192       msg->patch = patchID;
03193       msg->origPe = CkMyPe();
03194       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
03195       msg->destPe = node;
03196       int seq = flags.sequence;
03197       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03198       SET_PRIORITY(msg,seq,priority);
03199       cp[node].recvData(msg);
03200     }
03201   }
03202   Patch::gbisP2Ready();
03203 }

void HomePatch::gbisP3Ready (  ) 

Reimplemented from Patch.

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

03206                             {
03207   if (proxy.size() > 0) {
03208     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03209     //only nonzero message should be sent for doFullElec
03210     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
03211     for (int i = 0; i < proxy.size(); i++) {
03212       int node = proxy[i];
03213       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
03214       msg->patch = patchID;
03215       msg->dHdrPrefixLen = msgAtoms;
03216       msg->origPe = CkMyPe();
03217       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
03218       msg->destPe = node;
03219       int seq = flags.sequence;
03220       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03221       SET_PRIORITY(msg,seq,priority);
03222       cp[node].recvData(msg);
03223     }
03224   }
03225   Patch::gbisP3Ready();
03226 }

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 1991 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().

01993 {
01994   Molecule *mol = Node::Object()->molecule;
01995   SimParameters *simParams = Node::Object()->simParameters;
01996   const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
01997   const int fixedAtomsOn = simParams->fixedAtomsOn;
01998   const BigReal dt = timestep / TIMEFACTOR;
01999   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02000   int i, ia, ib, j;
02001   int dieOnError = simParams->rigidDie;
02002   Tensor wc;  // constraint virial
02003   BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
02004   int nslabs;
02005 
02006   // start data for hard wall boundary between drude and its host atom
02007   // static int Count=0;
02008   int Idx;
02009   double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
02010   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;
02011   double dot_v_r_1, dot_v_r_2;
02012   double vb_cm, dr_a, dr_b;
02013   // end data for hard wall boundary between drude and its host atom
02014 
02015   // start calculation of hard wall boundary between drude and its host atom
02016   if (simParams->drudeHardWallOn) {
02017     if (ppreduction) {
02018       nslabs = simParams->pressureProfileSlabs;
02019       idz = nslabs/lattice.c().z;
02020       zmin = lattice.origin().z - 0.5*lattice.c().z;
02021     }
02022 
02023     r_wall = simParams->drudeBondLen;
02024     r_wall_SQ = r_wall*r_wall;
02025     // Count++;
02026     for (i=1; i<numAtoms; i++)  {
02027       if ( (0.05 < atom[i].mass) && ((atom[i].mass < 1.0)) ) { // drude particle
02028         ia = i-1;
02029         ib = i;
02030 
02031         v_ab = atom[ib].position - atom[ia].position;
02032         rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
02033 
02034         if (rab_SQ > r_wall_SQ) {  // to impose the hard wall constraint
02035           rab = sqrt(rab_SQ);
02036           if ( (rab > (2.0*r_wall)) && dieOnError ) {  // unexpected situation
02037             iout << iERROR << "HardWallDrude> "
02038               << "The drude is too far away from atom "
02039               << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
02040             return -1;  // triggers early exit
02041           }
02042 
02043           v_ab.x /= rab;
02044           v_ab.y /= rab;
02045           v_ab.z /= rab;
02046 
02047           if ( fixedAtomsOn && atom[ia].atomFixed ) {  // the heavy atom is fixed
02048             if (atom[ib].atomFixed) {  // the drude is fixed too
02049               continue;
02050             }
02051             else {  // only the heavy atom is fixed
02052               dot_v_r_2 = atom[ib].velocity.x*v_ab.x
02053                 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
02054               vb_2 = dot_v_r_2 * v_ab;
02055               vp_2 = atom[ib].velocity - vb_2;
02056 
02057               dr = rab - r_wall;
02058               if(dot_v_r_2 == 0.0) {
02059                 delta_T = maxtime;
02060               }
02061               else {
02062                 delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
02063                 if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
02064               }
02065 
02066               dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
02067 
02068               vb_2 = dot_v_r_2 * v_ab;
02069 
02070               new_vel_a = atom[ia].velocity;
02071               new_vel_b = vp_2 + vb_2;
02072 
02073               dr_b = -dr + delta_T*dot_v_r_2;  // L = L_0 + dT *v_new, v was flipped
02074 
02075               new_pos_a = atom[ia].position;
02076               new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
02077             }
02078           }
02079           else {
02080             mass_a = atom[ia].mass;
02081             mass_b = atom[ib].mass;
02082             mass_sum = mass_a+mass_b;
02083 
02084             dot_v_r_1 = atom[ia].velocity.x*v_ab.x
02085               + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
02086             vb_1 = dot_v_r_1 * v_ab;
02087             vp_1 = atom[ia].velocity - vb_1;
02088 
02089             dot_v_r_2 = atom[ib].velocity.x*v_ab.x
02090               + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
02091             vb_2 = dot_v_r_2 * v_ab;
02092             vp_2 = atom[ib].velocity - vb_2;
02093 
02094             vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
02095 
02096             dot_v_r_1 -= vb_cm;
02097             dot_v_r_2 -= vb_cm;
02098 
02099             dr = rab - r_wall;
02100 
02101             if(dot_v_r_2 == dot_v_r_1) {
02102                 delta_T = maxtime;
02103             }
02104             else {
02105               delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);  // the time since the collision occurs
02106               if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
02107             }
02108             
02109             // the relative velocity between ia and ib. Drawn according to T_Drude
02110             v_Bond = sqrt(kbt/mass_b);
02111 
02112             // reflect the velocity along bond vector and scale down
02113             dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
02114             dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
02115 
02116             dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
02117             dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
02118 
02119             new_pos_a = atom[ia].position + dr_a*v_ab;  // correct the position
02120             new_pos_b = atom[ib].position + dr_b*v_ab;
02121             // atom[ia].position += (dr_a*v_ab);  // correct the position
02122             // atom[ib].position += (dr_b*v_ab);
02123             
02124             dot_v_r_1 += vb_cm;
02125             dot_v_r_2 += vb_cm;
02126 
02127             vb_1 = dot_v_r_1 * v_ab;
02128             vb_2 = dot_v_r_2 * v_ab;
02129 
02130             new_vel_a = vp_1 + vb_1;
02131             new_vel_b = vp_2 + vb_2;
02132           }
02133 
02134           int ppoffset, partition;
02135           if ( invdt == 0 ) {
02136             atom[ia].position = new_pos_a;
02137             atom[ib].position = new_pos_b;
02138           }
02139           else if ( virial == 0 ) {
02140             atom[ia].velocity = new_vel_a;
02141             atom[ib].velocity = new_vel_b;
02142           }
02143           else {
02144             for ( j = 0; j < 2; j++ ) {
02145               if (j==0) {  // atom ia, heavy atom
02146                 Idx = ia;
02147                 new_pos = &new_pos_a;
02148                 new_vel = &new_vel_a;
02149               }
02150               else if (j==1) {  // atom ib, drude
02151                 Idx = ib;
02152                 new_pos = &new_pos_b;
02153                 new_vel = &new_vel_b;
02154               }
02155               Force df = (*new_vel - atom[Idx].velocity) *
02156                 ( atom[Idx].mass * invdt );
02157               Tensor vir = outer(df, atom[Idx].position);
02158               wc += vir;
02159               atom[Idx].velocity = *new_vel;
02160               atom[Idx].position = *new_pos;
02161 
02162               if (ppreduction) {
02163                 if (!j) {
02164                   BigReal z = new_pos->z;
02165                   int partition = atom[Idx].partition;
02166                   int slab = (int)floor((z-zmin)*idz);
02167                   if (slab < 0) slab += nslabs;
02168                   else if (slab >= nslabs) slab -= nslabs;
02169                   ppoffset = 3*(slab + nslabs*partition);
02170                 }
02171                 ppreduction->item(ppoffset  ) += vir.xx;
02172                 ppreduction->item(ppoffset+1) += vir.yy;
02173                 ppreduction->item(ppoffset+2) += vir.zz;
02174               }
02175 
02176             }
02177           }
02178         }                               
02179       }
02180     }
02181 
02182     // if ( (Count>10000) && (Count%10==0) ) {
02183     //   v_ab = atom[1].position - atom[0].position;
02184     //   rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
02185     //   iout << "DBG_R: " << Count << "  " << sqrt(rab_SQ) << "\n" << endi;
02186     // }
02187 
02188   }
02189 
02190   // end calculation of hard wall boundary between drude and its host atom
02191 
02192   if ( dt && virial ) *virial += wc;
02193 
02194   return 0;
02195 }

void HomePatch::loweAndersenFinish (  ) 

Definition at line 3080 of file HomePatch.C.

References DebugM.

03081 {
03082     DebugM(2, "loweAndersenFinish\n");
03083     v.resize(0);
03084 }

void HomePatch::loweAndersenVelocities (  ) 

Definition at line 3065 of file HomePatch.C.

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

Referenced by positionsReady().

03066 {
03067     DebugM(2, "loweAndersenVelocities\n");
03068     Molecule *mol = Node::Object()->molecule;
03069     SimParameters *simParams = Node::Object()->simParameters;
03070     v.resize(numAtoms);
03071     for (int i = 0; i < numAtoms; ++i) {
03072         //v[i] = p[i];
03073         // co-opt CompAtom structure to pass velocity and mass information
03074         v[i].position = atom[i].velocity;
03075         v[i].charge = atom[i].mass;
03076     }
03077     DebugM(2, "loweAndersenVelocities\n");
03078 }

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

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

02931 {
02932   Molecule *mol = Node::Object()->molecule;
02933   SimParameters *simParams = Node::Object()->simParameters;
02934   Force *f1 = f[Results::normal].begin();
02935   const int fixedAtomsOn = simParams->fixedAtomsOn;
02936   const int useSettle = simParams->useSettle;
02937   const BigReal dt = timestep / TIMEFACTOR;
02938   Tensor wc;  // constraint virial
02939   BigReal tol = simParams->rigidTol;
02940   int maxiter = simParams->rigidIter;
02941   int dieOnError = simParams->rigidDie;
02942   int i, iter;
02943   BigReal dsqi[10], tmp;
02944   int ial[10], ibl[10];
02945   Vector ref[10];  // reference position
02946   Vector refab[10];  // reference vector
02947   Vector vel[10];  // new velocity
02948   BigReal rmass[10];  // 1 / mass
02949   BigReal redmass[10];  // reduced mass
02950   int fixed[10];  // is atom fixed?
02951   
02952   // Size of a hydrogen group for water
02953   int wathgsize = 3;
02954   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02955   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02956 
02957   //  CkPrintf("In rattle2!\n");
02958   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02959     //    CkPrintf("ig=%d\n",ig);
02960     int hgs = atom[ig].hydrogenGroupSize;
02961     if ( hgs == 1 ) continue;  // only one atom in group
02962     // cache data in local arrays and integrate positions normally
02963     int anyfixed = 0;
02964     for ( i = 0; i < hgs; ++i ) {
02965       ref[i] = atom[ig+i].position;
02966       vel[i] = ( forces ? f1[ig+i] : atom[ig+i].velocity );
02967       rmass[i] = 1.0;
02968       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02969       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02970     }
02971     int icnt = 0;
02972     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02973       if (hgs != wathgsize) {
02974         NAMD_bug("Hydrogen group error caught in rattle2().");
02975       }
02976       // Use SETTLE for water unless some of the water atoms are fixed,
02977       if (useSettle && !anyfixed) {
02978         if (simParams->watmodel == WAT_SWM4) {
02979           // SWM4 ordering:  O D LP H1 H2
02980           // do swap(O,LP) and call settle with subarray O H1 H2
02981           // swap back after we return
02982           Vector lp_ref = ref[2];
02983           // Vector lp_vel = vel[2];
02984           ref[2] = ref[0];
02985           vel[2] = vel[0];
02986           settle2(1.0, 1.0, ref+2, vel+2, dt, virial);
02987           ref[0] = ref[2];
02988           vel[0] = vel[2];
02989           ref[2] = lp_ref;
02990           // vel[2] = vel[0];  // set LP vel to O vel
02991         }
02992         else {
02993           settle2(1.0, 1.0, ref, vel, dt, virial);
02994           if (simParams->watmodel == WAT_TIP4) {
02995             vel[3] = vel[0];
02996           }
02997         }
02998         for (i=0; i<hgs; i++) {
02999           ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
03000         }
03001         continue;
03002       }
03003       if ( !(fixed[1] && fixed[2]) ) {
03004         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
03005         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03006       }
03007     }
03008     //    CkPrintf("Loop 2\n");
03009     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03010       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
03011         if ( !(fixed[0] && fixed[i]) ) {
03012           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
03013           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
03014           ibl[icnt] = i;  ++icnt;
03015         }
03016       }
03017     }
03018     if ( icnt == 0 ) continue;  // no constraints
03019     //    CkPrintf("Loop 3\n");
03020     for ( i = 0; i < icnt; ++i ) {
03021       refab[i] = ref[ial[i]] - ref[ibl[i]];
03022     }
03023     //    CkPrintf("Loop 4\n");
03024     int done;
03025     for ( iter = 0; iter < maxiter; ++iter ) {
03026       done = 1;
03027       for ( i = 0; i < icnt; ++i ) {
03028         int a = ial[i];  int b = ibl[i];
03029         Vector vab = vel[a] - vel[b];
03030         Vector &rab = refab[i];
03031         BigReal rabsqi = dsqi[i];
03032         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
03033         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
03034           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
03035           wc += outer(dp,rab);
03036           vel[a] += rmass[a] * dp;
03037           vel[b] -= rmass[b] * dp;
03038           done = 0;
03039         }
03040       }
03041       if ( done ) break;
03042       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
03043     }
03044     if ( ! done ) {
03045       if ( dieOnError ) {
03046         NAMD_die("Exceeded maximum number of iterations in rattle2().");
03047       } else {
03048         iout << iWARN <<
03049           "Exceeded maximum number of iterations in rattle2().\n" << endi;
03050       }
03051     }
03052     // store data back to patch
03053     for ( i = 0; i < hgs; ++i ) {
03054       ( forces ? f1[ig+i] : atom[ig+i].velocity ) = vel[i];
03055     }
03056   }
03057   //  CkPrintf("Leaving rattle2!\n");
03058   // check that there isn't a constant needed here!
03059   *virial += wc / ( 0.5 * dt );
03060 
03061 }

void HomePatch::mollyAverage (  ) 

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

03281 {
03282   Molecule *mol = Node::Object()->molecule;
03283   SimParameters *simParams = Node::Object()->simParameters;
03284   BigReal tol = simParams->mollyTol;
03285   int maxiter = simParams->mollyIter;
03286   int i, iter;
03287   HGArrayBigReal dsq;
03288   BigReal tmp;
03289   HGArrayInt ial, ibl;
03290   HGArrayVector ref;  // reference position
03291   HGArrayVector refab;  // reference vector
03292   HGArrayBigReal rmass;  // 1 / mass
03293   BigReal *lambda;  // Lagrange multipliers
03294   CompAtom *avg;  // averaged position
03295   int numLambdas = 0;
03296   HGArrayInt fixed;  // is atom fixed?
03297 
03298   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
03299   p_avg.resize(numAtoms);
03300   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
03301 
03302   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03303     int hgs = atom[ig].hydrogenGroupSize;
03304     if ( hgs == 1 ) continue;  // only one atom in group
03305         for ( i = 0; i < hgs; ++i ) {
03306           ref[i] = atom[ig+i].position;
03307           rmass[i] = 1. / atom[ig+i].mass;
03308           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03309           if ( fixed[i] ) rmass[i] = 0.;
03310         }
03311         avg = &(p_avg[ig]);
03312         int icnt = 0;
03313 
03314   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
03315           if ( hgs != 3 ) {
03316             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
03317           }
03318           if ( !(fixed[1] && fixed[2]) ) {
03319             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03320           }
03321         }
03322         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03323     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
03324             if ( !(fixed[0] && fixed[i]) ) {
03325               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03326             }
03327           }
03328         }
03329         if ( icnt == 0 ) continue;  // no constraints
03330         numLambdas += icnt;
03331         molly_lambda.resize(numLambdas);
03332         lambda = &(molly_lambda[numLambdas - icnt]);
03333         for ( i = 0; i < icnt; ++i ) {
03334           refab[i] = ref[ial[i]] - ref[ibl[i]];
03335         }
03336         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
03337         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
03338         if ( iter == maxiter ) {
03339           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
03340         }
03341   }
03342 
03343   // for ( i=0; i<numAtoms; ++i ) {
03344   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
03345   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
03346   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
03347   //    }
03348   // }
03349 
03350 }

void HomePatch::mollyMollify ( Tensor virial  ) 

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

03355 {
03356   Molecule *mol = Node::Object()->molecule;
03357   SimParameters *simParams = Node::Object()->simParameters;
03358   Tensor wc;  // constraint virial
03359   int i;
03360   HGArrayInt ial, ibl;
03361   HGArrayVector ref;  // reference position
03362   CompAtom *avg;  // averaged position
03363   HGArrayVector refab;  // reference vector
03364   HGArrayVector force;  // new force
03365   HGArrayBigReal rmass;  // 1 / mass
03366   BigReal *lambda;  // Lagrange multipliers
03367   int numLambdas = 0;
03368   HGArrayInt fixed;  // is atom fixed?
03369 
03370   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03371     int hgs = atom[ig].hydrogenGroupSize;
03372     if (hgs == 1 ) continue;  // only one atom in group
03373         for ( i = 0; i < hgs; ++i ) {
03374           ref[i] = atom[ig+i].position;
03375           force[i] = f[Results::slow][ig+i];
03376           rmass[i] = 1. / atom[ig+i].mass;
03377           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03378           if ( fixed[i] ) rmass[i] = 0.;
03379         }
03380         int icnt = 0;
03381         // c-ji I'm only going to mollify water for now
03382   if ( atom[ig].rigidBondLength > 0 ) {  // for water
03383           if ( hgs != 3 ) {
03384             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
03385           }
03386           if ( !(fixed[1] && fixed[2]) ) {
03387             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03388           }
03389         }
03390         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03391     if ( atom[ig+i].rigidBondLength > 0 ) {
03392             if ( !(fixed[0] && fixed[i]) ) {
03393               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03394             }
03395           }
03396         }
03397 
03398         if ( icnt == 0 ) continue;  // no constraints
03399         lambda = &(molly_lambda[numLambdas]);
03400         numLambdas += icnt;
03401         for ( i = 0; i < icnt; ++i ) {
03402           refab[i] = ref[ial[i]] - ref[ibl[i]];
03403         }
03404         avg = &(p_avg[ig]);
03405         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
03406         // store data back to patch
03407         for ( i = 0; i < hgs; ++i ) {
03408           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
03409           f[Results::slow][ig+i] = force[i];
03410         }
03411   }
03412   // check that there isn't a constant needed here!
03413   *virial += wc;
03414   p_avg.resize(0);
03415 }

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 2349 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().

02350                                 {
02351 
02352   SimParameters *simParams = Node::Object()->simParameters;
02353   if (simParams->watmodel != WAT_TIP3 || ppreduction) {
02354     // Call old rattle1 -method instead
02355     return rattle1old(timestep, virial, ppreduction);
02356   }
02357 
02358   if (!rattleListValid) {
02359     buildRattleList();
02360     rattleListValid = true;
02361   }
02362 
02363   const int fixedAtomsOn = simParams->fixedAtomsOn;
02364   const int useSettle = simParams->useSettle;
02365   const BigReal dt = timestep / TIMEFACTOR;
02366   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02367   const BigReal tol2 = 2.0 * simParams->rigidTol;
02368   int maxiter = simParams->rigidIter;
02369   int dieOnError = simParams->rigidDie;
02370 
02371   Vector ref[10];  // reference position
02372   Vector pos[10];  // new position
02373   Vector vel[10];  // new velocity
02374 
02375   // Manual un-roll
02376   int n = (settleList.size()/2)*2;
02377   for (int j=0;j < n;j+=2) {
02378     int ig;
02379     ig = settleList[j];
02380     for (int i = 0; i < 3; ++i ) {
02381       ref[i] = atom[ig+i].position;
02382       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02383     }
02384     ig = settleList[j+1];
02385     for (int i = 0; i < 3; ++i ) {
02386       ref[i+3] = atom[ig+i].position;
02387       pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
02388     }
02389     settle1_SIMD<2>(ref, pos,
02390       settle_mOrmT, settle_mHrmT, settle_ra,
02391       settle_rb, settle_rc, settle_rra);
02392 
02393     ig = settleList[j];
02394     for (int i = 0; i < 3; ++i ) {
02395       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02396       posNew[ig+i] = pos[i];
02397     }
02398     ig = settleList[j+1];
02399     for (int i = 0; i < 3; ++i ) {
02400       velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
02401       posNew[ig+i] = pos[i+3];
02402     }
02403 
02404   }
02405 
02406   if (settleList.size() % 2) {
02407     int ig = settleList[settleList.size()-1];
02408     for (int i = 0; i < 3; ++i ) {
02409       ref[i] = atom[ig+i].position;
02410       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02411     }
02412     settle1_SIMD<1>(ref, pos,
02413             settle_mOrmT, settle_mHrmT, settle_ra,
02414             settle_rb, settle_rc, settle_rra);        
02415     for (int i = 0; i < 3; ++i ) {
02416       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02417       posNew[ig+i] = pos[i];
02418     }
02419   }
02420 
02421   int posParam = 0;
02422   for (int j=0;j < rattleList.size();++j) {
02423 
02424     BigReal refx[10];
02425     BigReal refy[10];
02426     BigReal refz[10];
02427 
02428     BigReal posx[10];
02429     BigReal posy[10];
02430     BigReal posz[10];
02431 
02432     int ig = rattleList[j].ig;
02433     int icnt = rattleList[j].icnt;
02434     int hgs = atom[ig].hydrogenGroupSize;
02435     for (int i = 0; i < hgs; ++i ) {
02436       ref[i] = atom[ig+i].position;
02437       pos[i] = atom[ig+i].position;
02438       if (!(fixedAtomsOn && atom[ig+i].atomFixed)) {
02439         pos[i] += atom[ig+i].velocity * dt;
02440       }
02441       refx[i] = ref[i].x;
02442       refy[i] = ref[i].y;
02443       refz[i] = ref[i].z;
02444       posx[i] = pos[i].x;
02445       posy[i] = pos[i].y;
02446       posz[i] = pos[i].z;
02447     }
02448 
02449     bool done;
02450     bool consFailure;
02451     if (icnt == 1) {
02452       rattlePair<1>(&rattleParam[posParam],
02453         refx, refy, refz,
02454         posx, posy, posz,
02455         consFailure);
02456       done = true;
02457     } else {
02458       rattleN(icnt, &rattleParam[posParam],
02459         refx, refy, refz,
02460         posx, posy, posz,
02461         tol2, maxiter,
02462         done, consFailure);
02463     }
02464 
02465     // Advance position in rattleParam
02466     posParam += icnt;
02467 
02468     for (int i = 0; i < hgs; ++i ) {
02469       pos[i].x = posx[i];
02470       pos[i].y = posy[i];
02471       pos[i].z = posz[i];
02472     }
02473 
02474     for (int i = 0; i < hgs; ++i ) {
02475       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02476       posNew[ig+i] = pos[i];
02477     }
02478 
02479     if ( consFailure ) {
02480       if ( dieOnError ) {
02481         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02482         << (atom[ig].id + 1) << "!\n" << endi;
02483         return -1;  // triggers early exit
02484       } else {
02485         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02486         << (atom[ig].id + 1) << "!\n" << endi;
02487       }
02488     } else if ( ! done ) {
02489       if ( dieOnError ) {
02490         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02491         << (atom[ig].id + 1) << "!\n" << endi;
02492         return -1;  // triggers early exit
02493       } else {
02494         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02495         << (atom[ig].id + 1) << "!\n" << endi;
02496       }
02497     }
02498   }
02499   // Finally, we have to go through atoms that are not involved in rattle just so that we have
02500   // their positions and velocities up-to-date in posNew and velNew
02501   for (int j=0;j < noconstList.size();++j) {
02502     int ig = noconstList[j];
02503     int hgs = atom[ig].hydrogenGroupSize;
02504     for (int i = 0; i < hgs; ++i ) {
02505       velNew[ig+i] = atom[ig+i].velocity;
02506       posNew[ig+i] = atom[ig+i].position;
02507     }
02508   }
02509 
02510   if ( invdt == 0 ) {
02511     for (int ig = 0; ig < numAtoms; ++ig )
02512       atom[ig].position = posNew[ig];
02513   } else if ( virial == 0 ) {
02514     for (int ig = 0; ig < numAtoms; ++ig )
02515       atom[ig].velocity = velNew[ig];
02516   } else {
02517     Tensor wc;  // constraint virial
02518     addRattleForce(invdt, wc);
02519     *virial += wc;
02520   }
02521 
02522   return 0;
02523 }

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

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

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

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

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

02797 {
02798   Molecule *mol = Node::Object()->molecule;
02799   SimParameters *simParams = Node::Object()->simParameters;
02800   const int fixedAtomsOn = simParams->fixedAtomsOn;
02801   const int useSettle = simParams->useSettle;
02802   const BigReal dt = timestep / TIMEFACTOR;
02803   Tensor wc;  // constraint virial
02804   BigReal tol = simParams->rigidTol;
02805   int maxiter = simParams->rigidIter;
02806   int dieOnError = simParams->rigidDie;
02807   int i, iter;
02808   BigReal dsqi[10], tmp;
02809   int ial[10], ibl[10];
02810   Vector ref[10];  // reference position
02811   Vector refab[10];  // reference vector
02812   Vector vel[10];  // new velocity
02813   BigReal rmass[10];  // 1 / mass
02814   BigReal redmass[10];  // reduced mass
02815   int fixed[10];  // is atom fixed?
02816   
02817   // Size of a hydrogen group for water
02818   int wathgsize = 3;
02819   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02820   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02821 
02822   //  CkPrintf("In rattle2!\n");
02823   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02824     //    CkPrintf("ig=%d\n",ig);
02825     int hgs = atom[ig].hydrogenGroupSize;
02826     if ( hgs == 1 ) continue;  // only one atom in group
02827     // cache data in local arrays and integrate positions normally
02828     int anyfixed = 0;
02829     for ( i = 0; i < hgs; ++i ) {
02830       ref[i] = atom[ig+i].position;
02831       vel[i] = atom[ig+i].velocity;
02832       rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
02833       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02834       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02835     }
02836     int icnt = 0;
02837     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02838       if (hgs != wathgsize) {
02839         NAMD_bug("Hydrogen group error caught in rattle2().");
02840       }
02841       // Use SETTLE for water unless some of the water atoms are fixed,
02842       if (useSettle && !anyfixed) {
02843         if (simParams->watmodel == WAT_SWM4) {
02844           // SWM4 ordering:  O D LP H1 H2
02845           // do swap(O,LP) and call settle with subarray O H1 H2
02846           // swap back after we return
02847           Vector lp_ref = ref[2];
02848           // Vector lp_vel = vel[2];
02849           ref[2] = ref[0];
02850           vel[2] = vel[0];
02851           settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
02852           ref[0] = ref[2];
02853           vel[0] = vel[2];
02854           ref[2] = lp_ref;
02855           // vel[2] = vel[0];  // set LP vel to O vel
02856         }
02857         else {
02858           settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
02859           if (simParams->watmodel == WAT_TIP4) {
02860             vel[3] = vel[0];
02861           }
02862         }
02863         for (i=0; i<hgs; i++) {
02864           atom[ig+i].velocity = vel[i];
02865         }
02866         continue;
02867       }
02868       if ( !(fixed[1] && fixed[2]) ) {
02869         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02870         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02871       }
02872     }
02873     //    CkPrintf("Loop 2\n");
02874     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02875       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02876         if ( !(fixed[0] && fixed[i]) ) {
02877           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02878           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02879           ibl[icnt] = i;  ++icnt;
02880         }
02881       }
02882     }
02883     if ( icnt == 0 ) continue;  // no constraints
02884     //    CkPrintf("Loop 3\n");
02885     for ( i = 0; i < icnt; ++i ) {
02886       refab[i] = ref[ial[i]] - ref[ibl[i]];
02887     }
02888     //    CkPrintf("Loop 4\n");
02889     int done;
02890     for ( iter = 0; iter < maxiter; ++iter ) {
02891       done = 1;
02892       for ( i = 0; i < icnt; ++i ) {
02893         int a = ial[i];  int b = ibl[i];
02894         Vector vab = vel[a] - vel[b];
02895         Vector &rab = refab[i];
02896         BigReal rabsqi = dsqi[i];
02897         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02898         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02899           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02900           wc += outer(dp,rab);
02901           vel[a] += rmass[a] * dp;
02902           vel[b] -= rmass[b] * dp;
02903           done = 0;
02904         }
02905       }
02906       if ( done ) break;
02907       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02908     }
02909     if ( ! done ) {
02910       if ( dieOnError ) {
02911         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02912       } else {
02913         iout << iWARN <<
02914           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02915       }
02916     }
02917     // store data back to patch
02918     for ( i = 0; i < hgs; ++i ) {
02919       atom[ig+i].velocity = vel[i];
02920     }
02921   }
02922   //  CkPrintf("Leaving rattle2!\n");
02923   // check that there isn't a constant needed here!
02924   *virial += wc / ( 0.5 * dt );
02925 
02926 }

void HomePatch::receiveResult ( ProxyGBISP2ResultMsg msg  ) 

Definition at line 3254 of file HomePatch.C.

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

03254                                                        {
03255   ++numGBISP2Arrived;
03256   //accumulate dEda
03257   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
03258     dHdrPrefix[i] += msg->dEdaSum[i];
03259   }
03260   delete msg;
03261 
03262   if (flags.doNonbonded) {
03263     //awaken if phase 2 done
03264     if (phase2BoxClosedCalled == true &&
03265         numGBISP2Arrived==proxy.size() ) {
03266       sequencer->awaken();
03267     }
03268   } else {
03269     //awaken if all phases done on noWork step
03270     if (boxesOpen == 0 &&
03271         numGBISP1Arrived == proxy.size() &&
03272         numGBISP2Arrived == proxy.size() &&
03273         numGBISP3Arrived == proxy.size()) {
03274       sequencer->awaken();
03275     }
03276   }
03277 }

void HomePatch::receiveResult ( ProxyGBISP1ResultMsg msg  ) 

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

03229                                                        {
03230   ++numGBISP1Arrived;
03231     for ( int i = 0; i < msg->psiSumLen; ++i ) {
03232       psiFin[i] += msg->psiSum[i];
03233     }
03234   delete msg;
03235 
03236   if (flags.doNonbonded) {
03237     //awaken if phase 1 done
03238     if (phase1BoxClosedCalled == true &&
03239         numGBISP1Arrived==proxy.size() ) {
03240         sequencer->awaken();
03241     }
03242   } else {
03243     //awaken if all phases done on noWork step
03244     if (boxesOpen == 0 &&
03245         numGBISP1Arrived == proxy.size() &&
03246         numGBISP2Arrived == proxy.size() &&
03247         numGBISP3Arrived == proxy.size()) {
03248       sequencer->awaken();
03249     }
03250   }
03251 }

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 3536 of file HomePatch.C.

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

03536                                   {  // initiating replica
03537   CkpvAccess(_qd)->process();
03538 }

void HomePatch::recvCheckpointLoad ( CheckpointAtomsMsg msg  ) 

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

03483                                                           {  // initiating replica
03484   SimParameters *simParams = Node::Object()->simParameters;
03485   const int remote = simParams->scriptIntArg1;
03486   const char *key = simParams->scriptStringArg1;
03487   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
03488     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
03489   }
03490   if ( msg->replica != remote ) {
03491     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
03492   }
03493   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03494     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
03495     strcpy(newmsg->key,key);
03496     newmsg->lattice = lattice;
03497     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
03498     newmsg->pid = patchID;
03499     newmsg->pe = CkMyPe();
03500     newmsg->replica = CmiMyPartition();
03501     newmsg->numAtoms = numAtoms;
03502     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03503     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
03504   }
03505   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03506     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03507     lattice = msg->lattice;
03508     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
03509     numAtoms = msg->numAtoms;
03510     atom.resize(numAtoms);
03511     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03512     doAtomUpdate = true;
03513     rattleListValid = false;
03514     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03515   }
03516   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
03517     recvCheckpointAck();
03518   }
03519   delete msg;
03520 }

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

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

03453                                                                                 {  // responding replica
03454   if ( task == SCRIPT_CHECKPOINT_FREE ) {
03455     if ( ! checkpoints.count(key) ) {
03456       NAMD_die("Unable to free checkpoint, requested key was never stored.");
03457     }
03458     delete checkpoints[key];
03459     checkpoints.erase(key);
03460     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
03461     return;
03462   }
03463   CheckpointAtomsMsg *msg;
03464   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
03465     if ( ! checkpoints.count(key) ) {
03466       NAMD_die("Unable to load checkpoint, requested key was never stored.");
03467     }
03468     checkpoint_t &cp = *checkpoints[key];
03469     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
03470     msg->lattice = cp.lattice;
03471     msg->berendsenPressure_count = cp.berendsenPressure_count;
03472     msg->numAtoms = cp.numAtoms;
03473     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
03474   } else {
03475     msg = new (0,1,0) CheckpointAtomsMsg;
03476   }
03477   msg->pid = patchID;
03478   msg->replica = CmiMyPartition();
03479   msg->pe = CkMyPe();
03480   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
03481 }

void HomePatch::recvCheckpointStore ( CheckpointAtomsMsg msg  ) 

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

03522                                                            {  // responding replica
03523   if ( ! checkpoints.count(msg->key) ) {
03524     checkpoints[msg->key] = new checkpoint_t;
03525   }
03526   checkpoint_t &cp = *checkpoints[msg->key];
03527   cp.lattice = msg->lattice;
03528   cp.berendsenPressure_count = msg->berendsenPressure_count;
03529   cp.numAtoms = msg->numAtoms;
03530   cp.atoms.resize(cp.numAtoms);
03531   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
03532   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
03533   delete msg;
03534 }

void HomePatch::recvExchangeMsg ( ExchangeAtomsMsg msg  ) 

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

03573                                                      {
03574   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
03575   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03576   lattice = msg->lattice;
03577   numAtoms = msg->numAtoms;
03578   atom.resize(numAtoms);
03579   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03580   delete msg;
03581   CkpvAccess(_qd)->process();
03582   doAtomUpdate = true;
03583   rattleListValid = false;
03584   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03585 }

void HomePatch::recvExchangeReq ( int  req  ) 

Definition at line 3562 of file HomePatch.C.

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

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

03562                                        {
03563   exchange_req = req;
03564   if ( exchange_msg ) {
03565     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
03566     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
03567     exchange_msg = 0;
03568     exchange_req = -1;
03569     CkpvAccess(_qd)->process();
03570   }
03571 }

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 3427 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().

03427                            {
03428   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03429 
03430   atom.copy(checkpoint_atom);
03431   numAtoms = atom.size();
03432   lattice = checkpoint_lattice;
03433 
03434   doAtomUpdate = true;
03435   rattleListValid = false;
03436 
03437   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03438 
03439   // DMK - Atom Separation (water vs. non-water)
03440   #if NAMD_SeparateWaters != 0
03441     numWaterAtoms = checkpoint_numWaterAtoms;
03442   #endif
03443 }

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 3099 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().

03099                                       {
03100   intRad.resize(numAtoms*2);
03101   intRad.setall(0);
03102   Molecule *mol = Node::Object()->molecule;
03103   SimParameters *simParams = Node::Object()->simParameters;
03104   Real offset = simParams->coulomb_radius_offset;
03105   for (int i = 0; i < numAtoms; i++) {
03106     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
03107     Real screen = MassToScreen(atom[i].mass);//same
03108     intRad[2*i+0] = rad - offset;//r0
03109     intRad[2*i+1] = screen*(rad - offset);//s0
03110   }
03111 }

void HomePatch::setLcpoType (  ) 

Definition at line 3088 of file HomePatch.C.

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

Referenced by positionsReady().

03088                             {
03089   Molecule *mol = Node::Object()->molecule;
03090   const int *lcpoParamType = mol->getLcpoParamType();
03091 
03092   lcpoType.resize(numAtoms);
03093   for (int i = 0; i < numAtoms; i++) {
03094     lcpoType[i] = lcpoParamType[pExt[i].id];
03095   }
03096 }

void HomePatch::submitLoadStats ( int  timestep  ) 

Definition at line 3587 of file HomePatch.C.

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

Referenced by Sequencer::rebalanceLoad().

03588 {
03589   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
03590 }

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 5 Aug 2020 for NAMD by  doxygen 1.6.1