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

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

Referenced by Sequencer::addForceToMomentum().

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 dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01947         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01948         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01949         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01950       }
01951     }
01952   } else {
01953     for ( int i = 0; i < num_atoms; ++i ) {
01954       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01955       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01956       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01957       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01958     }
01959   }
01960 }

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

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

Referenced by Sequencer::addForceToMomentum3().

01971       {
01972   SimParameters *simParams = Node::Object()->simParameters;
01973 
01974   if ( simParams->fixedAtomsOn ) {
01975     for ( int i = 0; i < num_atoms; ++i ) {
01976       if ( atom_arr[i].atomFixed ) {
01977         atom_arr[i].velocity = 0;
01978       } else {
01979         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01980         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01981             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01982         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01983             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01984         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01985             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01986       }
01987     }
01988   } else {
01989     for ( int i = 0; i < num_atoms; ++i ) {
01990       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01991       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01992           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01993       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01994           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01995       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01996           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01997     }
01998   }
01999 }

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

Definition at line 2372 of file HomePatch.C.

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

Referenced by rattle1().

02372                                                               {
02373   for (int ig = 0; ig < numAtoms; ++ig ) {
02374     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
02375     Tensor vir = outer(df, atom[ig].position);
02376     wc += vir;
02377     f[Results::normal][ig] += df;
02378     atom[ig].velocity = velNew[ig];
02379   }
02380 }

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

Definition at line 2001 of file HomePatch.C.

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

Referenced by Sequencer::addVelocityToPosition().

02005       {
02006   SimParameters *simParams = Node::Object()->simParameters;
02007   if ( simParams->fixedAtomsOn ) {
02008     for ( int i = 0; i < num_atoms; ++i ) {
02009       if ( ! atom_arr[i].atomFixed ) {
02010         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
02011         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
02012         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
02013       }
02014     }
02015   } else {
02016     for ( int i = 0; i < num_atoms; ++i ) {
02017       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
02018       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
02019       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
02020     }
02021   }
02022 }

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

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

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

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

Referenced by Sequencer::algorithm().

03450                                {
03451   checkpoint_atom.copy(atom);
03452   checkpoint_lattice = lattice;
03453 
03454   // DMK - Atom Separation (water vs. non-water)
03455   #if NAMD_SeparateWaters != 0
03456     checkpoint_numWaterAtoms = numWaterAtoms;
03457   #endif
03458 }

void HomePatch::depositMigration ( MigrateAtomsMsg msg  ) 

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

04000 {
04001 
04002   if (!inMigration) { // We have to buffer changes due to migration
04003                       // until our patch is in migration mode
04004     msgbuf[numMlBuf++] = msg;
04005     return;
04006   } 
04007 
04008 
04009   // DMK - Atom Separation (water vs. non-water)
04010   #if NAMD_SeparateWaters != 0
04011 
04012 
04013     // Merge the incoming list of atoms with the current list of
04014     //   atoms.  Note that mergeSeparatedAtomList() will apply any
04015     //   required transformations to the incoming atoms as it is
04016     //   separating them.
04017     mergeAtomList(msg->migrationList);
04018 
04019 
04020   #else
04021 
04022     // Merge the incoming list of atoms with the current list of
04023     // atoms.  Apply transformations to the incoming atoms as they are
04024     // added to this patch's list.
04025     {
04026       MigrationList &migrationList = msg->migrationList;
04027       MigrationList::iterator mi;
04028       Transform mother_transform;
04029       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
04030         DebugM(1,"Migrating atom " << mi->id << " to patch "
04031                   << patchID << " with position " << mi->position << "\n"); 
04032         if ( mi->migrationGroupSize ) {
04033           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
04034             Position pos = mi->position;
04035             int mgs = mi->migrationGroupSize;
04036             int c = 1;
04037             for ( int j=mi->hydrogenGroupSize; j<mgs;
04038                                 j+=(mi+j)->hydrogenGroupSize ) {
04039               pos += (mi+j)->position;
04040               ++c;
04041             }
04042             pos *= 1./c;
04043             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
04044             mother_transform = mi->transform;
04045             pos = lattice.nearest(pos,center,&mother_transform);
04046             mi->position = lattice.reverse_transform(mi->position,mi->transform);
04047             mi->position = lattice.apply_transform(mi->position,mother_transform);
04048             mi->transform = mother_transform;
04049           } else {
04050             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
04051             mother_transform = mi->transform;
04052           }
04053         } else {
04054           mi->position = lattice.reverse_transform(mi->position,mi->transform);
04055           mi->position = lattice.apply_transform(mi->position,mother_transform);
04056           mi->transform = mother_transform;
04057         }
04058         atom.add(*mi);
04059       }
04060     }
04061 
04062 
04063   #endif // if (NAMD_SeparateWaters != 0)
04064 
04065 
04066   numAtoms = atom.size();
04067   delete msg;
04068 
04069   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
04070   if (!--patchMigrationCounter) {
04071     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
04072     allMigrationIn = true;
04073     patchMigrationCounter = numNeighbors;
04074     if (migrationSuspended) {
04075       DebugM(3,"patch "<<patchID<<" is being awakened\n");
04076       migrationSuspended = false;
04077       sequencer->awaken();
04078       return;
04079     }
04080     else {
04081        DebugM(3,"patch "<<patchID<<" was not suspended\n");
04082     }
04083   }
04084 }

void HomePatch::doAtomMigration (  )  [protected]

Definition at line 3839 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_START, 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().

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

void HomePatch::doGroupSizeCheck (  )  [protected]

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

03719 {
03720   if ( ! flags.doNonbonded ) return;
03721 
03722   SimParameters *simParams = Node::Object()->simParameters;
03723   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
03724   BigReal maxrad2 = 0.;
03725 
03726   FullAtomList::iterator p_i = atom.begin();
03727   FullAtomList::iterator p_e = atom.end();
03728 
03729   while ( p_i != p_e ) {
03730     const int hgs = p_i->hydrogenGroupSize;
03731     if ( ! hgs ) break;  // avoid infinite loop on bug
03732     int ngs = hgs;
03733     if ( ngs > 5 ) ngs = 5;  // XXX why? limit to at most 5 atoms per group
03734     BigReal x = p_i->position.x;
03735     BigReal y = p_i->position.y;
03736     BigReal z = p_i->position.z;
03737     int i;
03738     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
03739       p_i[i].nonbondedGroupSize = 0;
03740       BigReal dx = p_i[i].position.x - x;
03741       BigReal dy = p_i[i].position.y - y;
03742       BigReal dz = p_i[i].position.z - z;
03743       BigReal r2 = dx * dx + dy * dy + dz * dz;
03744       if ( r2 > hgcut ) break;
03745       else if ( r2 > maxrad2 ) maxrad2 = r2;
03746     }
03747     p_i->nonbondedGroupSize = i;
03748     for ( ; i < hgs; ++i ) {
03749       p_i[i].nonbondedGroupSize = 1;
03750     }
03751     p_i += hgs;
03752   }
03753 
03754   if ( p_i != p_e ) {
03755     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
03756   }
03757 
03758   flags.maxGroupRadius = sqrt(maxrad2);
03759 
03760 }

void HomePatch::doMarginCheck (  )  [protected]

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

03763 {
03764   SimParameters *simParams = Node::Object()->simParameters;
03765 
03766   BigReal sysdima = lattice.a_r().unit() * lattice.a();
03767   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
03768   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
03769 
03770   BigReal minSize = simParams->patchDimension - simParams->margin;
03771 
03772   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
03773        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
03774        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
03775 
03776     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03777       "Possible solutions are to restart from a recent checkpoint,\n"
03778       "increase margin, or disable useFlexibleCell for liquid simulation.");
03779   }
03780 
03781   BigReal cutoff = simParams->cutoff;
03782 
03783   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
03784   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
03785   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
03786 
03787   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
03788     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03789       "There are probably many margin violations already reported.\n"
03790       "Possible solutions are to restart from a recent checkpoint,\n"
03791       "increase margin, or disable useFlexibleCell for liquid simulation.");
03792   }
03793 
03794   BigReal minx = min.x - margina;
03795   BigReal miny = min.y - marginb;
03796   BigReal minz = min.z - marginc;
03797   BigReal maxx = max.x + margina;
03798   BigReal maxy = max.y + marginb;
03799   BigReal maxz = max.z + marginc;
03800 
03801   int xdev, ydev, zdev;
03802   int problemCount = 0;
03803 
03804   FullAtomList::iterator p_i = atom.begin();
03805   FullAtomList::iterator p_e = atom.end();
03806   for ( ; p_i != p_e; ++p_i ) {
03807 
03808     ScaledPosition s = lattice.scale(p_i->position);
03809 
03810     // check if atom is within bounds
03811     if (s.x < minx) xdev = 0;
03812     else if (maxx <= s.x) xdev = 2; 
03813     else xdev = 1;
03814 
03815     if (s.y < miny) ydev = 0;
03816     else if (maxy <= s.y) ydev = 2; 
03817     else ydev = 1;
03818 
03819     if (s.z < minz) zdev = 0;
03820     else if (maxz <= s.z) zdev = 2; 
03821     else zdev = 1;
03822 
03823     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
03824         ++problemCount;
03825     }
03826 
03827   }
03828 
03829   marginViolations = problemCount;
03830   // if ( problemCount ) {
03831   //     iout << iERROR <<
03832   //       "Found " << problemCount << " margin violations!\n" << endi;
03833   // } 
03834 
03835 }

void HomePatch::doPairlistCheck (  )  [protected]

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

03630 {
03631 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
03632   char dpcbuf[32];
03633   sprintf(dpcbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::DO_PAIRLIST_CHECK], this->getPatchID());
03634   NAMD_EVENT_START_EX(1, NamdProfileEvent::DO_PAIRLIST_CHECK, dpcbuf);
03635 #endif
03636 
03637   SimParameters *simParams = Node::Object()->simParameters;
03638 
03639   if ( numAtoms == 0 || ! flags.usePairlists ) {
03640     flags.pairlistTolerance = 0.;
03641     flags.maxAtomMovement = 99999.;
03642 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
03643     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03644 #endif
03645     return;
03646   }
03647 
03648   int i; int n = numAtoms;
03649   CompAtom *p_i = p.begin();
03650 
03651   if ( flags.savePairlists ) {
03652     flags.pairlistTolerance = doPairlistCheck_newTolerance;
03653     flags.maxAtomMovement = 0.;
03654     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
03655     doPairlistCheck_lattice = lattice;
03656     doPairlistCheck_positions.resize(numAtoms);
03657     CompAtom *psave_i = doPairlistCheck_positions.begin();
03658     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
03659 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
03660     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03661 #endif
03662     return;
03663   }
03664 
03665   Lattice &lattice_old = doPairlistCheck_lattice;
03666   Position center_cur = lattice.unscale(center);
03667   Position center_old = lattice_old.unscale(center);
03668   Vector center_delta = center_cur - center_old;
03669   
03670   // find max deviation to corner (any neighbor shares a corner)
03671   BigReal max_cd = 0.;
03672   for ( i=0; i<2; ++i ) {
03673     for ( int j=0; j<2; ++j ) {
03674       for ( int k=0; k<2; ++k ) {
03675         ScaledPosition corner(  i ? min.x : max.x ,
03676                                 j ? min.y : max.y ,
03677                                 k ? min.z : max.z );
03678         Vector corner_delta =
03679                 lattice.unscale(corner) - lattice_old.unscale(corner);
03680         corner_delta -= center_delta;
03681         BigReal cd = corner_delta.length2();
03682         if ( cd > max_cd ) max_cd = cd;
03683       }
03684     }
03685   }
03686   max_cd = sqrt(max_cd);
03687 
03688   // find max deviation of atoms relative to center
03689   BigReal max_pd = 0.;
03690   CompAtom *p_old_i = doPairlistCheck_positions.begin();
03691   for ( i=0; i<n; ++i ) {
03692     Vector p_delta = p_i[i].position - p_old_i[i].position;
03693     p_delta -= center_delta;
03694     BigReal pd = p_delta.length2();
03695     if ( pd > max_pd ) max_pd = pd;
03696   }
03697   max_pd = sqrt(max_pd);
03698 
03699   BigReal max_tol = max_pd + max_cd;
03700 
03701   flags.maxAtomMovement = max_tol;
03702 
03703   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
03704 
03705   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
03706                                 doPairlistCheck_newTolerance ) ) {
03707     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
03708   }
03709 
03710   if ( max_tol > doPairlistCheck_newTolerance ) {
03711     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
03712   }
03713 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
03714     NAMD_EVENT_STOP(1, NamdProfileEvent::DO_PAIRLIST_CHECK);
03715 #endif
03716 }

void HomePatch::exchangeAtoms ( int  scriptTask  ) 

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

03574                                             {
03575   SimParameters *simParams = Node::Object()->simParameters;
03576   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
03577   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
03578     exchange_dst = (int) simParams->scriptArg1;
03579     // create and save outgoing message
03580     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
03581     exchange_msg->lattice = lattice;
03582     exchange_msg->pid = patchID;
03583     exchange_msg->numAtoms = numAtoms;
03584     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03585     if ( exchange_req >= 0 ) {
03586       recvExchangeReq(exchange_req);
03587     }
03588   }
03589   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
03590     exchange_src = (int) simParams->scriptArg2;
03591     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
03592   }
03593 }

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

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

03478                                                            {  // initiating replica
03479   SimParameters *simParams = Node::Object()->simParameters;
03480   checkpoint_task = scriptTask;
03481   const int remote = simParams->scriptIntArg1;
03482   const char *key = simParams->scriptStringArg1;
03483   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
03484 }

void HomePatch::gbisComputeAfterP1 (  ) 

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

03147                                    {
03148 
03149   SimParameters *simParams = Node::Object()->simParameters;
03150   BigReal alphaMax = simParams->alpha_max;
03151   BigReal delta = simParams->gbis_delta;
03152   BigReal beta = simParams->gbis_beta;
03153   BigReal gamma = simParams->gbis_gamma;
03154   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03155 
03156   BigReal rhoi;
03157   BigReal rhoi0;
03158   //calculate bornRad from psiSum
03159   for (int i = 0; i < numAtoms; i++) {
03160     rhoi0 = intRad[2*i];
03161     rhoi = rhoi0+coulomb_radius_offset;
03162     psiFin[i] += psiSum[i];
03163     psiFin[i] *= rhoi0;
03164     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
03165     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
03166 #ifdef PRINT_COMP
03167     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
03168 #endif
03169   }
03170 
03171   gbisP2Ready();
03172 }

void HomePatch::gbisComputeAfterP2 (  ) 

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

03175                                    {
03176 
03177   SimParameters *simParams = Node::Object()->simParameters;
03178   BigReal delta = simParams->gbis_delta;
03179   BigReal beta = simParams->gbis_beta;
03180   BigReal gamma = simParams->gbis_gamma;
03181   BigReal epsilon_s = simParams->solvent_dielectric;
03182   BigReal epsilon_p = simParams->dielectric;
03183   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
03184   BigReal epsilon_p_i = 1/simParams->dielectric;
03185   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
03186   BigReal kappa = simParams->kappa;
03187   BigReal fij, expkappa, Dij, dEdai, dedasum;
03188   BigReal rhoi, rhoi0, psii, nbetapsi;
03189   BigReal gammapsi2, tanhi, daidr;
03190   for (int i = 0; i < numAtoms; i++) {
03191     //add diagonal dEda term
03192     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
03193     fij = bornRad[i];//inf
03194     expkappa = exp(-kappa*fij);//0
03195     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
03196     //calculate dHij prefix
03197     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
03198                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
03199     dHdrPrefix[i] += dEdai;
03200     dedasum = dHdrPrefix[i];
03201 
03202     rhoi0 = intRad[2*i];
03203     rhoi = rhoi0+coulomb_radius_offset;
03204     psii = psiFin[i];
03205     nbetapsi = -beta*psii;
03206     gammapsi2 = gamma*psii*psii;
03207     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
03208     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
03209            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
03210     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
03211 #ifdef PRINT_COMP
03212     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
03213 #endif
03214   }
03215   gbisP3Ready();
03216 }

void HomePatch::gbisP2Ready (  ) 

Reimplemented from Patch.

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

03219                             {
03220   if (proxy.size() > 0) {
03221     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03222     for (int i = 0; i < proxy.size(); i++) {
03223       int node = proxy[i];
03224       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
03225       msg->patch = patchID;
03226       msg->origPe = CkMyPe();
03227       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
03228       msg->destPe = node;
03229       int seq = flags.sequence;
03230       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03231       SET_PRIORITY(msg,seq,priority);
03232       cp[node].recvData(msg);
03233     }
03234   }
03235   Patch::gbisP2Ready();
03236 }

void HomePatch::gbisP3Ready (  ) 

Reimplemented from Patch.

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

03239                             {
03240   if (proxy.size() > 0) {
03241     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03242     //only nonzero message should be sent for doFullElec
03243     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
03244     for (int i = 0; i < proxy.size(); i++) {
03245       int node = proxy[i];
03246       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
03247       msg->patch = patchID;
03248       msg->dHdrPrefixLen = msgAtoms;
03249       msg->origPe = CkMyPe();
03250       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
03251       msg->destPe = node;
03252       int seq = flags.sequence;
03253       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03254       SET_PRIORITY(msg,seq,priority);
03255       cp[node].recvData(msg);
03256     }
03257   }
03258   Patch::gbisP3Ready();
03259 }

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

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

void HomePatch::loweAndersenFinish (  ) 

Definition at line 3113 of file HomePatch.C.

References DebugM.

03114 {
03115     DebugM(2, "loweAndersenFinish\n");
03116     v.resize(0);
03117 }

void HomePatch::loweAndersenVelocities (  ) 

Definition at line 3098 of file HomePatch.C.

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

Referenced by positionsReady().

03099 {
03100     DebugM(2, "loweAndersenVelocities\n");
03101     Molecule *mol = Node::Object()->molecule;
03102     SimParameters *simParams = Node::Object()->simParameters;
03103     v.resize(numAtoms);
03104     for (int i = 0; i < numAtoms; ++i) {
03105         //v[i] = p[i];
03106         // co-opt CompAtom structure to pass velocity and mass information
03107         v[i].position = atom[i].velocity;
03108         v[i].charge = atom[i].mass;
03109     }
03110     DebugM(2, "loweAndersenVelocities\n");
03111 }

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

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

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

void HomePatch::mollyAverage (  ) 

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

03314 {
03315   Molecule *mol = Node::Object()->molecule;
03316   SimParameters *simParams = Node::Object()->simParameters;
03317   BigReal tol = simParams->mollyTol;
03318   int maxiter = simParams->mollyIter;
03319   int i, iter;
03320   HGArrayBigReal dsq;
03321   BigReal tmp;
03322   HGArrayInt ial, ibl;
03323   HGArrayVector ref;  // reference position
03324   HGArrayVector refab;  // reference vector
03325   HGArrayBigReal rmass;  // 1 / mass
03326   BigReal *lambda;  // Lagrange multipliers
03327   CompAtom *avg;  // averaged position
03328   int numLambdas = 0;
03329   HGArrayInt fixed;  // is atom fixed?
03330 
03331   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
03332   p_avg.resize(numAtoms);
03333   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
03334 
03335   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03336     int hgs = atom[ig].hydrogenGroupSize;
03337     if ( hgs == 1 ) continue;  // only one atom in group
03338         for ( i = 0; i < hgs; ++i ) {
03339           ref[i] = atom[ig+i].position;
03340           rmass[i] = 1. / atom[ig+i].mass;
03341           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03342           if ( fixed[i] ) rmass[i] = 0.;
03343         }
03344         avg = &(p_avg[ig]);
03345         int icnt = 0;
03346 
03347   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
03348           if ( hgs != 3 ) {
03349             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
03350           }
03351           if ( !(fixed[1] && fixed[2]) ) {
03352             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03353           }
03354         }
03355         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03356     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
03357             if ( !(fixed[0] && fixed[i]) ) {
03358               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03359             }
03360           }
03361         }
03362         if ( icnt == 0 ) continue;  // no constraints
03363         numLambdas += icnt;
03364         molly_lambda.resize(numLambdas);
03365         lambda = &(molly_lambda[numLambdas - icnt]);
03366         for ( i = 0; i < icnt; ++i ) {
03367           refab[i] = ref[ial[i]] - ref[ibl[i]];
03368         }
03369         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
03370         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
03371         if ( iter == maxiter ) {
03372           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
03373         }
03374   }
03375 
03376   // for ( i=0; i<numAtoms; ++i ) {
03377   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
03378   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
03379   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
03380   //    }
03381   // }
03382 
03383 }

void HomePatch::mollyMollify ( Tensor virial  ) 

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

03388 {
03389   Molecule *mol = Node::Object()->molecule;
03390   SimParameters *simParams = Node::Object()->simParameters;
03391   Tensor wc;  // constraint virial
03392   int i;
03393   HGArrayInt ial, ibl;
03394   HGArrayVector ref;  // reference position
03395   CompAtom *avg;  // averaged position
03396   HGArrayVector refab;  // reference vector
03397   HGArrayVector force;  // new force
03398   HGArrayBigReal rmass;  // 1 / mass
03399   BigReal *lambda;  // Lagrange multipliers
03400   int numLambdas = 0;
03401   HGArrayInt fixed;  // is atom fixed?
03402 
03403   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03404     int hgs = atom[ig].hydrogenGroupSize;
03405     if (hgs == 1 ) continue;  // only one atom in group
03406         for ( i = 0; i < hgs; ++i ) {
03407           ref[i] = atom[ig+i].position;
03408           force[i] = f[Results::slow][ig+i];
03409           rmass[i] = 1. / atom[ig+i].mass;
03410           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03411           if ( fixed[i] ) rmass[i] = 0.;
03412         }
03413         int icnt = 0;
03414         // c-ji I'm only going to mollify water for now
03415   if ( atom[ig].rigidBondLength > 0 ) {  // for water
03416           if ( hgs != 3 ) {
03417             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
03418           }
03419           if ( !(fixed[1] && fixed[2]) ) {
03420             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03421           }
03422         }
03423         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03424     if ( atom[ig+i].rigidBondLength > 0 ) {
03425             if ( !(fixed[0] && fixed[i]) ) {
03426               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03427             }
03428           }
03429         }
03430 
03431         if ( icnt == 0 ) continue;  // no constraints
03432         lambda = &(molly_lambda[numLambdas]);
03433         numLambdas += icnt;
03434         for ( i = 0; i < icnt; ++i ) {
03435           refab[i] = ref[ial[i]] - ref[ibl[i]];
03436         }
03437         avg = &(p_avg[ig]);
03438         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
03439         // store data back to patch
03440         for ( i = 0; i < hgs; ++i ) {
03441           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
03442           f[Results::slow][ig+i] = force[i];
03443         }
03444   }
03445   // check that there isn't a constant needed here!
03446   *virial += wc;
03447   p_avg.resize(0);
03448 }

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

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

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

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

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

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

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

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

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

void HomePatch::receiveResult ( ProxyGBISP2ResultMsg msg  ) 

Definition at line 3287 of file HomePatch.C.

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

03287                                                        {
03288   ++numGBISP2Arrived;
03289   //accumulate dEda
03290   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
03291     dHdrPrefix[i] += msg->dEdaSum[i];
03292   }
03293   delete msg;
03294 
03295   if (flags.doNonbonded) {
03296     //awaken if phase 2 done
03297     if (phase2BoxClosedCalled == true &&
03298         numGBISP2Arrived==proxy.size() ) {
03299       sequencer->awaken();
03300     }
03301   } else {
03302     //awaken if all phases done on noWork step
03303     if (boxesOpen == 0 &&
03304         numGBISP1Arrived == proxy.size() &&
03305         numGBISP2Arrived == proxy.size() &&
03306         numGBISP3Arrived == proxy.size()) {
03307       sequencer->awaken();
03308     }
03309   }
03310 }

void HomePatch::receiveResult ( ProxyGBISP1ResultMsg msg  ) 

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

03262                                                        {
03263   ++numGBISP1Arrived;
03264     for ( int i = 0; i < msg->psiSumLen; ++i ) {
03265       psiFin[i] += msg->psiSum[i];
03266     }
03267   delete msg;
03268 
03269   if (flags.doNonbonded) {
03270     //awaken if phase 1 done
03271     if (phase1BoxClosedCalled == true &&
03272         numGBISP1Arrived==proxy.size() ) {
03273         sequencer->awaken();
03274     }
03275   } else {
03276     //awaken if all phases done on noWork step
03277     if (boxesOpen == 0 &&
03278         numGBISP1Arrived == proxy.size() &&
03279         numGBISP2Arrived == proxy.size() &&
03280         numGBISP3Arrived == proxy.size()) {
03281       sequencer->awaken();
03282     }
03283   }
03284 }

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

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

03569                                   {  // initiating replica
03570   CkpvAccess(_qd)->process();
03571 }

void HomePatch::recvCheckpointLoad ( CheckpointAtomsMsg msg  ) 

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

03516                                                           {  // initiating replica
03517   SimParameters *simParams = Node::Object()->simParameters;
03518   const int remote = simParams->scriptIntArg1;
03519   const char *key = simParams->scriptStringArg1;
03520   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
03521     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
03522   }
03523   if ( msg->replica != remote ) {
03524     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
03525   }
03526   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03527     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
03528     strcpy(newmsg->key,key);
03529     newmsg->lattice = lattice;
03530     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
03531     newmsg->pid = patchID;
03532     newmsg->pe = CkMyPe();
03533     newmsg->replica = CmiMyPartition();
03534     newmsg->numAtoms = numAtoms;
03535     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03536     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
03537   }
03538   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03539     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03540     lattice = msg->lattice;
03541     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
03542     numAtoms = msg->numAtoms;
03543     atom.resize(numAtoms);
03544     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03545     doAtomUpdate = true;
03546     rattleListValid = false;
03547     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03548   }
03549   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
03550     recvCheckpointAck();
03551   }
03552   delete msg;
03553 }

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

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

03486                                                                                 {  // responding replica
03487   if ( task == SCRIPT_CHECKPOINT_FREE ) {
03488     if ( ! checkpoints.count(key) ) {
03489       NAMD_die("Unable to free checkpoint, requested key was never stored.");
03490     }
03491     delete checkpoints[key];
03492     checkpoints.erase(key);
03493     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
03494     return;
03495   }
03496   CheckpointAtomsMsg *msg;
03497   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
03498     if ( ! checkpoints.count(key) ) {
03499       NAMD_die("Unable to load checkpoint, requested key was never stored.");
03500     }
03501     checkpoint_t &cp = *checkpoints[key];
03502     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
03503     msg->lattice = cp.lattice;
03504     msg->berendsenPressure_count = cp.berendsenPressure_count;
03505     msg->numAtoms = cp.numAtoms;
03506     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
03507   } else {
03508     msg = new (0,1,0) CheckpointAtomsMsg;
03509   }
03510   msg->pid = patchID;
03511   msg->replica = CmiMyPartition();
03512   msg->pe = CkMyPe();
03513   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
03514 }

void HomePatch::recvCheckpointStore ( CheckpointAtomsMsg msg  ) 

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

03555                                                            {  // responding replica
03556   if ( ! checkpoints.count(msg->key) ) {
03557     checkpoints[msg->key] = new checkpoint_t;
03558   }
03559   checkpoint_t &cp = *checkpoints[msg->key];
03560   cp.lattice = msg->lattice;
03561   cp.berendsenPressure_count = msg->berendsenPressure_count;
03562   cp.numAtoms = msg->numAtoms;
03563   cp.atoms.resize(cp.numAtoms);
03564   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
03565   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
03566   delete msg;
03567 }

void HomePatch::recvExchangeMsg ( ExchangeAtomsMsg msg  ) 

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

03606                                                      {
03607   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
03608   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03609   lattice = msg->lattice;
03610   numAtoms = msg->numAtoms;
03611   atom.resize(numAtoms);
03612   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03613   delete msg;
03614   CkpvAccess(_qd)->process();
03615   doAtomUpdate = true;
03616   rattleListValid = false;
03617   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03618 }

void HomePatch::recvExchangeReq ( int  req  ) 

Definition at line 3595 of file HomePatch.C.

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

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

03595                                        {
03596   exchange_req = req;
03597   if ( exchange_msg ) {
03598     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
03599     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
03600     exchange_msg = 0;
03601     exchange_req = -1;
03602     CkpvAccess(_qd)->process();
03603   }
03604 }

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

01332 {
01333   replacementForces = f;
01334 }

void HomePatch::revert ( void   ) 

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

03460                            {
03461   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03462 
03463   atom.copy(checkpoint_atom);
03464   numAtoms = atom.size();
03465   lattice = checkpoint_lattice;
03466 
03467   doAtomUpdate = true;
03468   rattleListValid = false;
03469 
03470   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03471 
03472   // DMK - Atom Separation (water vs. non-water)
03473   #if NAMD_SeparateWaters != 0
03474     numWaterAtoms = checkpoint_numWaterAtoms;
03475   #endif
03476 }

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

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

Referenced by Sequencer::saveForce().

01337 {
01338   f_saved[ftag].resize(numAtoms);
01339   for ( int i = 0; i < numAtoms; ++i )
01340   {
01341     f_saved[ftag][i] = f[ftag][i];
01342   }
01343 }

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

03132                                       {
03133   intRad.resize(numAtoms*2);
03134   intRad.setall(0);
03135   Molecule *mol = Node::Object()->molecule;
03136   SimParameters *simParams = Node::Object()->simParameters;
03137   Real offset = simParams->coulomb_radius_offset;
03138   for (int i = 0; i < numAtoms; i++) {
03139     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
03140     Real screen = MassToScreen(atom[i].mass);//same
03141     intRad[2*i+0] = rad - offset;//r0
03142     intRad[2*i+1] = screen*(rad - offset);//s0
03143   }
03144 }

void HomePatch::setLcpoType (  ) 

Definition at line 3121 of file HomePatch.C.

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

Referenced by positionsReady().

03121                             {
03122   Molecule *mol = Node::Object()->molecule;
03123   const int *lcpoParamType = mol->getLcpoParamType();
03124 
03125   lcpoType.resize(numAtoms);
03126   for (int i = 0; i < numAtoms; i++) {
03127     lcpoType[i] = lcpoParamType[pExt[i].id];
03128   }
03129 }

void HomePatch::submitLoadStats ( int  timestep  ) 

Definition at line 3620 of file HomePatch.C.

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

Referenced by Sequencer::rebalanceLoad().

03621 {
03622   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
03623 }

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 16 Jun 2022 for NAMD by  doxygen 1.6.1