HomePatch Class Reference

#include <HomePatch.h>

Inheritance diagram for HomePatch:

Patch List of all members.

Public Member Functions

 ~HomePatch ()
void registerProxy (RegisterProxyMsg *)
void unregisterProxy (UnregisterProxyMsg *)
void receiveResults (ProxyResultVarsizeMsg *msg)
void receiveResults (ProxyResultMsg *msg)
void receiveResult (ProxyGBISP1ResultMsg *msg)
void receiveResult (ProxyGBISP2ResultMsg *msg)
void receiveResults (ProxyCombinedResultRawMsg *msg)
void depositMigration (MigrateAtomsMsg *)
void useSequencer (Sequencer *sequencerPtr)
void runSequencer (void)
void positionsReady (int doMigration=0)
void saveForce (const int ftag=Results::normal)
void addForceToMomentum (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

Classes

struct  checkpoint_t
struct  RattleList

Detailed Description

Definition at line 88 of file HomePatch.h.


Constructor & Destructor Documentation

HomePatch::~HomePatch (  ) 

Definition at line 304 of file HomePatch.C.

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

00305 {
00306     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
00307 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00308     ptnTree.resize(0);
00309     #ifdef USE_NODEPATCHMGR
00310     delete [] nodeChildren;    
00311     #endif
00312 #endif
00313   delete [] child;
00314 }


Member Function Documentation

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

Definition at line 1715 of file HomePatch.C.

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

Referenced by Sequencer::addForceToMomentum().

01720       {
01721   SimParameters *simParams = Node::Object()->simParameters;
01722 
01723   if ( simParams->fixedAtomsOn ) {
01724     for ( int i = 0; i < num_atoms; ++i ) {
01725       if ( atom_arr[i].atomFixed ) {
01726         atom_arr[i].velocity = 0;
01727       } else {
01728         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01729         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01730         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01731         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01732       }
01733     }
01734   } else {
01735     for ( int i = 0; i < num_atoms; ++i ) {
01736       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01737       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01738       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01739       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01740     }
01741   }
01742 }

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

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

Referenced by Sequencer::addForceToMomentum3().

01753       {
01754   SimParameters *simParams = Node::Object()->simParameters;
01755 
01756   if ( simParams->fixedAtomsOn ) {
01757     for ( int i = 0; i < num_atoms; ++i ) {
01758       if ( atom_arr[i].atomFixed ) {
01759         atom_arr[i].velocity = 0;
01760       } else {
01761         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01762         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01763             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01764         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01765             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01766         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01767             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01768       }
01769     }
01770   } else {
01771     for ( int i = 0; i < num_atoms; ++i ) {
01772       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01773       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01774           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01775       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01776           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01777       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01778           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01779     }
01780   }
01781 }

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

Definition at line 2154 of file HomePatch.C.

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

Referenced by rattle1().

02154                                                               {
02155   for (int ig = 0; ig < numAtoms; ++ig ) {
02156     Force df = (velNew[ig] - atom[ig].velocity) * ( atom[ig].mass * invdt );
02157     Tensor vir = outer(df, atom[ig].position);
02158     wc += vir;
02159     f[Results::normal][ig] += df;
02160     atom[ig].velocity = velNew[ig];
02161   }
02162 }

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

Definition at line 1783 of file HomePatch.C.

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

Referenced by Sequencer::addVelocityToPosition().

01787       {
01788   SimParameters *simParams = Node::Object()->simParameters;
01789   if ( simParams->fixedAtomsOn ) {
01790     for ( int i = 0; i < num_atoms; ++i ) {
01791       if ( ! atom_arr[i].atomFixed ) {
01792         atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01793         atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01794         atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01795       }
01796     }
01797   } else {
01798     for ( int i = 0; i < num_atoms; ++i ) {
01799       atom_arr[i].position.x  +=  atom_arr[i].velocity.x * dt;
01800       atom_arr[i].position.y  +=  atom_arr[i].velocity.y * dt;
01801       atom_arr[i].position.z  +=  atom_arr[i].velocity.z * dt;
01802     }
01803   }
01804 }

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

Implements Patch.

Definition at line 317 of file HomePatch.C.

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

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

void HomePatch::buildRattleList (  ) 

Definition at line 2012 of file HomePatch.C.

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

Referenced by rattle1().

02012                                 {
02013 
02014   SimParameters *simParams = Node::Object()->simParameters;
02015   const int fixedAtomsOn = simParams->fixedAtomsOn;
02016   const int useSettle = simParams->useSettle;
02017 
02018   // Re-size to containt numAtoms elements
02019   velNew.resize(numAtoms);
02020   posNew.resize(numAtoms);
02021 
02022   // Size of a hydrogen group for water
02023   int wathgsize = 3;
02024   int watmodel = simParams->watmodel;
02025   if (watmodel == WAT_TIP4) wathgsize = 4;
02026   else if (watmodel == WAT_SWM4) wathgsize = 5;
02027 
02028   // Initialize the settle algorithm with water parameters
02029   // settle1() assumes all waters are identical,
02030   // and will generate bad results if they are not.
02031   // XXX this will move to Molecule::build_atom_status when that 
02032   // version is debugged
02033   if ( ! settle_initialized ) {
02034     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02035       // find a water
02036       if (atom[ig].rigidBondLength > 0) {
02037         int oatm;
02038         if (simParams->watmodel == WAT_SWM4) {
02039           oatm = ig+3;  // skip over Drude and Lonepair
02040           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02041           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02042         }
02043         else {
02044           oatm = ig+1;
02045           // Avoid using the Om site to set this by mistake
02046           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02047             oatm += 1;
02048           }
02049         }
02050 
02051         // initialize settle water parameters
02052         settle1init(atom[ig].mass, atom[oatm].mass, 
02053                     atom[ig].rigidBondLength,
02054                     atom[oatm].rigidBondLength,
02055                     settle_mOrmT, settle_mHrmT, settle_ra,
02056                     settle_rb, settle_rc, settle_rra);
02057         settle_initialized = 1;
02058         break; // done with init
02059       }
02060     }
02061   }
02062   
02063   Vector ref[10];
02064   BigReal rmass[10];
02065   BigReal dsq[10];
02066   int fixed[10];
02067   int ial[10];
02068   int ibl[10];
02069 
02070   int numSettle = 0;
02071   int numRattle = 0;
02072   int posRattleParam = 0;
02073 
02074   settleList.clear();
02075   rattleList.clear();
02076   noconstList.clear();
02077   rattleParam.clear();
02078 
02079   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02080     int hgs = atom[ig].hydrogenGroupSize;
02081     if ( hgs == 1 ) {
02082       // only one atom in group
02083       noconstList.push_back(ig);
02084       continue;
02085     }
02086     int anyfixed = 0;
02087     for (int i = 0; i < hgs; ++i ) {
02088       ref[i] = atom[ig+i].position;
02089       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02090       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02091       if ( fixed[i] ) {
02092         anyfixed = 1;
02093         rmass[i] = 0.;
02094       }
02095     }
02096     int icnt = 0;
02097     BigReal tmp = atom[ig].rigidBondLength;
02098     if (tmp > 0.0) {  // for water
02099       if (hgs != wathgsize) {
02100         char errmsg[256];
02101         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02102                          "but the specified water model requires %d atoms.\n",
02103                          atom[ig].id+1, hgs, wathgsize);
02104         NAMD_die(errmsg);
02105       }
02106       // Use SETTLE for water unless some of the water atoms are fixed,
02107       if (useSettle && !anyfixed) {
02108         // Store to Settle -list
02109         settleList.push_back(ig);
02110         continue;
02111       }
02112       if ( !(fixed[1] && fixed[2]) ) {
02113         dsq[icnt] = tmp * tmp;
02114         ial[icnt] = 1;
02115         ibl[icnt] = 2;
02116         ++icnt;
02117       }
02118     } // if (tmp > 0.0)
02119     for (int i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02120       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02121         if ( !(fixed[0] && fixed[i]) ) {
02122           dsq[icnt] = tmp * tmp;
02123           ial[icnt] = 0;
02124           ibl[icnt] = i;
02125           ++icnt;
02126         }
02127       }
02128     }
02129     if ( icnt == 0 ) {
02130       // no constraints
02131       noconstList.push_back(ig);
02132       continue;  
02133     }
02134     // Store to Rattle -list
02135     RattleList rattleListElem;
02136     rattleListElem.ig  = ig;
02137     rattleListElem.icnt = icnt;
02138     rattleList.push_back(rattleListElem);
02139     for (int i = 0; i < icnt; ++i ) {
02140       int a = ial[i];
02141       int b = ibl[i];
02142       RattleParam rattleParamElem;
02143       rattleParamElem.ia = a;
02144       rattleParamElem.ib = b;
02145       rattleParamElem.dsq = dsq[i];
02146       rattleParamElem.rma = rmass[a];
02147       rattleParamElem.rmb = rmass[b];
02148       rattleParam.push_back(rattleParamElem);
02149     }
02150   }
02151 
02152 }

void HomePatch::buildSpanningTree ( void   ) 

Definition at line 657 of file HomePatch.C.

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

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

void HomePatch::checkpoint ( void   ) 

Definition at line 3234 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

03234                                {
03235   checkpoint_atom.copy(atom);
03236   checkpoint_lattice = lattice;
03237 
03238   // DMK - Atom Separation (water vs. non-water)
03239   #if NAMD_SeparateWaters != 0
03240     checkpoint_numWaterAtoms = numWaterAtoms;
03241   #endif
03242 }

void HomePatch::depositMigration ( MigrateAtomsMsg  ) 

Definition at line 3743 of file HomePatch.C.

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

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

03744 {
03745 
03746   if (!inMigration) { // We have to buffer changes due to migration
03747                       // until our patch is in migration mode
03748     msgbuf[numMlBuf++] = msg;
03749     return;
03750   } 
03751 
03752 
03753   // DMK - Atom Separation (water vs. non-water)
03754   #if NAMD_SeparateWaters != 0
03755 
03756 
03757     // Merge the incoming list of atoms with the current list of
03758     //   atoms.  Note that mergeSeparatedAtomList() will apply any
03759     //   required transformations to the incoming atoms as it is
03760     //   separating them.
03761     mergeAtomList(msg->migrationList);
03762 
03763 
03764   #else
03765 
03766     // Merge the incoming list of atoms with the current list of
03767     // atoms.  Apply transformations to the incoming atoms as they are
03768     // added to this patch's list.
03769     {
03770       MigrationList &migrationList = msg->migrationList;
03771       MigrationList::iterator mi;
03772       Transform mother_transform;
03773       for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
03774         DebugM(1,"Migrating atom " << mi->id << " to patch "
03775                   << patchID << " with position " << mi->position << "\n"); 
03776         if ( mi->migrationGroupSize ) {
03777           if ( mi->migrationGroupSize != mi->hydrogenGroupSize ) {
03778             Position pos = mi->position;
03779             int mgs = mi->migrationGroupSize;
03780             int c = 1;
03781             for ( int j=mi->hydrogenGroupSize; j<mgs;
03782                                 j+=(mi+j)->hydrogenGroupSize ) {
03783               pos += (mi+j)->position;
03784               ++c;
03785             }
03786             pos *= 1./c;
03787             // iout << "mgroup " << mi->id << " at " << pos << "\n" << endi;
03788             mother_transform = mi->transform;
03789             pos = lattice.nearest(pos,center,&mother_transform);
03790             mi->position = lattice.reverse_transform(mi->position,mi->transform);
03791             mi->position = lattice.apply_transform(mi->position,mother_transform);
03792             mi->transform = mother_transform;
03793           } else {
03794             mi->position = lattice.nearest(mi->position,center,&(mi->transform));
03795             mother_transform = mi->transform;
03796           }
03797         } else {
03798           mi->position = lattice.reverse_transform(mi->position,mi->transform);
03799           mi->position = lattice.apply_transform(mi->position,mother_transform);
03800           mi->transform = mother_transform;
03801         }
03802         atom.add(*mi);
03803       }
03804     }
03805 
03806 
03807   #endif // if (NAMD_SeparateWaters != 0)
03808 
03809 
03810   numAtoms = atom.size();
03811   delete msg;
03812 
03813   DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
03814   if (!--patchMigrationCounter) {
03815     DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
03816     allMigrationIn = true;
03817     patchMigrationCounter = numNeighbors;
03818     if (migrationSuspended) {
03819       DebugM(3,"patch "<<patchID<<" is being awakened\n");
03820       migrationSuspended = false;
03821       sequencer->awaken();
03822       return;
03823     }
03824     else {
03825        DebugM(3,"patch "<<patchID<<" was not suspended\n");
03826     }
03827   }
03828 }

void HomePatch::doAtomMigration (  )  [protected]

Definition at line 3606 of file HomePatch.C.

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

Referenced by positionsReady().

03607 {
03608   int i;
03609 
03610   for (i=0; i<numNeighbors; i++) {
03611     realInfo[i].mList.resize(0);
03612   }
03613 
03614   // Purge the AtomMap
03615   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03616 
03617   // Determine atoms that need to migrate
03618 
03619   BigReal minx = min.x;
03620   BigReal miny = min.y;
03621   BigReal minz = min.z;
03622   BigReal maxx = max.x;
03623   BigReal maxy = max.y;
03624   BigReal maxz = max.z;
03625 
03626   int xdev, ydev, zdev;
03627   int delnum = 0;
03628 
03629   FullAtomList::iterator atom_i = atom.begin();
03630   FullAtomList::iterator atom_e = atom.end();
03631 
03632   // DMK - Atom Separation (water vs. non-water)
03633   #if NAMD_SeparateWaters != 0
03634     FullAtomList::iterator atom_first = atom_i;
03635     int numLostWaterAtoms = 0;
03636   #endif
03637 
03638   while ( atom_i != atom_e ) {
03639     if ( atom_i->migrationGroupSize ) {
03640       Position pos = atom_i->position;
03641       if ( atom_i->migrationGroupSize != atom_i->hydrogenGroupSize ) {
03642         int mgs = atom_i->migrationGroupSize;
03643         int c = 1;
03644         for ( int j=atom_i->hydrogenGroupSize; j<mgs;
03645                                 j+=(atom_i+j)->hydrogenGroupSize ) {
03646           pos += (atom_i+j)->position;
03647           ++c;
03648         }
03649         pos *= 1./c;
03650         // iout << "mgroup " << atom_i->id << " at " << pos << "\n" << endi;
03651       }
03652 
03653       ScaledPosition s = lattice.scale(pos);
03654 
03655       // check if atom is within bounds
03656       if (s.x < minx) xdev = 0;
03657       else if (maxx <= s.x) xdev = 2;
03658       else xdev = 1;
03659 
03660       if (s.y < miny) ydev = 0;
03661       else if (maxy <= s.y) ydev = 2;
03662       else ydev = 1;
03663 
03664       if (s.z < minz) zdev = 0;
03665       else if (maxz <= s.z) zdev = 2;
03666       else zdev = 1;
03667 
03668     }
03669 
03670     if (mInfo[xdev][ydev][zdev]) { // process atom for migration
03671                                     // Don't migrate if destination is myself
03672 
03673       // See if we have a migration list already
03674       MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
03675       DebugM(3,"Migrating atom " << atom_i->id << " from patch "
03676                 << patchID << " with position " << atom_i->position << "\n");
03677       mCur.add(*atom_i);
03678 
03679       ++delnum;
03680 
03681 
03682       // DMK - Atom Separation (water vs. non-water)
03683       #if NAMD_SeparateWaters != 0
03684         // Check to see if this atom is part of a water molecule.  If
03685         //   so, numWaterAtoms needs to be adjusted to refect the lost of
03686         //   this atom.
03687         // NOTE: The atom separation code assumes that if the oxygen
03688         //   atom of the hydrogen group making up the water molecule is
03689         //   migrated to another HomePatch, the hydrogens will also
03690         //   move!!!
03691         int atomIndex = atom_i - atom_first;
03692         if (atomIndex < numWaterAtoms)
03693           numLostWaterAtoms++;
03694       #endif
03695 
03696 
03697     } else {
03698 
03699       if ( delnum ) { *(atom_i-delnum) = *atom_i; }
03700 
03701     }
03702 
03703     ++atom_i;
03704   }
03705 
03706   // DMK - Atom Separation (water vs. non-water)
03707   #if NAMD_SeparateWaters != 0
03708     numWaterAtoms -= numLostWaterAtoms;
03709   #endif
03710 
03711 
03712   int delpos = numAtoms - delnum;
03713   DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
03714   atom.del(delpos,delnum);
03715 
03716   numAtoms = atom.size();
03717 
03718   PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
03719 
03720   // signal depositMigration() that we are inMigration mode
03721   inMigration = true;
03722 
03723   // Drain the migration message buffer
03724   for (i=0; i<numMlBuf; i++) {
03725      DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
03726      depositMigration(msgbuf[i]);
03727   }
03728   numMlBuf = 0;
03729      
03730   if (!allMigrationIn) {
03731     DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
03732     migrationSuspended = true;
03733     sequencer->suspend();
03734     migrationSuspended = false;
03735   }
03736   allMigrationIn = false;
03737 
03738   inMigration = false;
03739   marginViolations = 0;
03740 }

void HomePatch::doGroupSizeCheck (  )  [protected]

Definition at line 3485 of file HomePatch.C.

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

Referenced by positionsReady().

03486 {
03487   if ( ! flags.doNonbonded ) return;
03488 
03489   SimParameters *simParams = Node::Object()->simParameters;
03490   BigReal hgcut = 0.5 * simParams->hgroupCutoff;  hgcut *= hgcut;
03491   BigReal maxrad2 = 0.;
03492 
03493   FullAtomList::iterator p_i = atom.begin();
03494   FullAtomList::iterator p_e = atom.end();
03495 
03496   while ( p_i != p_e ) {
03497     const int hgs = p_i->hydrogenGroupSize;
03498     if ( ! hgs ) break;  // avoid infinite loop on bug
03499     int ngs = hgs;
03500     if ( ngs > 5 ) ngs = 5;  // limit to at most 5 atoms per group
03501     BigReal x = p_i->position.x;
03502     BigReal y = p_i->position.y;
03503     BigReal z = p_i->position.z;
03504     int i;
03505     for ( i = 1; i < ngs; ++i ) {  // limit spatial extent
03506       p_i[i].nonbondedGroupSize = 0;
03507       BigReal dx = p_i[i].position.x - x;
03508       BigReal dy = p_i[i].position.y - y;
03509       BigReal dz = p_i[i].position.z - z;
03510       BigReal r2 = dx * dx + dy * dy + dz * dz;
03511       if ( r2 > hgcut ) break;
03512       else if ( r2 > maxrad2 ) maxrad2 = r2;
03513     }
03514     p_i->nonbondedGroupSize = i;
03515     for ( ; i < hgs; ++i ) {
03516       p_i[i].nonbondedGroupSize = 1;
03517     }
03518     p_i += hgs;
03519   }
03520 
03521   if ( p_i != p_e ) {
03522     NAMD_bug("hydrogenGroupSize is zero in HomePatch::doGroupSizeCheck");
03523   }
03524 
03525   flags.maxGroupRadius = sqrt(maxrad2);
03526 
03527 }

void HomePatch::doMarginCheck (  )  [protected]

Definition at line 3529 of file HomePatch.C.

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

Referenced by positionsReady().

03530 {
03531   SimParameters *simParams = Node::Object()->simParameters;
03532 
03533   BigReal sysdima = lattice.a_r().unit() * lattice.a();
03534   BigReal sysdimb = lattice.b_r().unit() * lattice.b();
03535   BigReal sysdimc = lattice.c_r().unit() * lattice.c();
03536 
03537   BigReal minSize = simParams->patchDimension - simParams->margin;
03538 
03539   if ( ( aAwayDist*sysdima < minSize*0.9999 ) ||
03540        ( bAwayDist*sysdimb < minSize*0.9999 ) ||
03541        ( cAwayDist*sysdimc < minSize*0.9999 ) ) {
03542 
03543     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03544       "Possible solutions are to restart from a recent checkpoint,\n"
03545       "increase margin, or disable useFlexibleCell for liquid simulation.");
03546   }
03547 
03548   BigReal cutoff = simParams->cutoff;
03549 
03550   BigReal margina = 0.5 * ( aAwayDist - cutoff / sysdima );
03551   BigReal marginb = 0.5 * ( bAwayDist - cutoff / sysdimb );
03552   BigReal marginc = 0.5 * ( cAwayDist - cutoff / sysdimc );
03553 
03554   if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
03555     NAMD_die("Periodic cell has become too small for original patch grid!\n"
03556       "There are probably many margin violations already reported.\n"
03557       "Possible solutions are to restart from a recent checkpoint,\n"
03558       "increase margin, or disable useFlexibleCell for liquid simulation.");
03559   }
03560 
03561   BigReal minx = min.x - margina;
03562   BigReal miny = min.y - marginb;
03563   BigReal minz = min.z - marginc;
03564   BigReal maxx = max.x + margina;
03565   BigReal maxy = max.y + marginb;
03566   BigReal maxz = max.z + marginc;
03567 
03568   int xdev, ydev, zdev;
03569   int problemCount = 0;
03570 
03571   FullAtomList::iterator p_i = atom.begin();
03572   FullAtomList::iterator p_e = atom.end();
03573   for ( ; p_i != p_e; ++p_i ) {
03574 
03575     ScaledPosition s = lattice.scale(p_i->position);
03576 
03577     // check if atom is within bounds
03578     if (s.x < minx) xdev = 0;
03579     else if (maxx <= s.x) xdev = 2; 
03580     else xdev = 1;
03581 
03582     if (s.y < miny) ydev = 0;
03583     else if (maxy <= s.y) ydev = 2; 
03584     else ydev = 1;
03585 
03586     if (s.z < minz) zdev = 0;
03587     else if (maxz <= s.z) zdev = 2; 
03588     else zdev = 1;
03589 
03590     if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
03591         ++problemCount;
03592     }
03593 
03594   }
03595 
03596   marginViolations = problemCount;
03597   // if ( problemCount ) {
03598   //     iout << iERROR <<
03599   //       "Found " << problemCount << " margin violations!\n" << endi;
03600   // } 
03601 
03602 }

void HomePatch::doPairlistCheck (  )  [protected]

Definition at line 3410 of file HomePatch.C.

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

Referenced by positionsReady().

03411 {
03412   SimParameters *simParams = Node::Object()->simParameters;
03413 
03414   if ( numAtoms == 0 || ! flags.usePairlists ) {
03415     flags.pairlistTolerance = 0.;
03416     flags.maxAtomMovement = 99999.;
03417     return;
03418   }
03419 
03420   int i; int n = numAtoms;
03421   CompAtom *p_i = p.begin();
03422 
03423   if ( flags.savePairlists ) {
03424     flags.pairlistTolerance = doPairlistCheck_newTolerance;
03425     flags.maxAtomMovement = 0.;
03426     doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
03427     doPairlistCheck_lattice = lattice;
03428     doPairlistCheck_positions.resize(numAtoms);
03429     CompAtom *psave_i = doPairlistCheck_positions.begin();
03430     for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
03431     return;
03432   }
03433 
03434   Lattice &lattice_old = doPairlistCheck_lattice;
03435   Position center_cur = lattice.unscale(center);
03436   Position center_old = lattice_old.unscale(center);
03437   Vector center_delta = center_cur - center_old;
03438   
03439   // find max deviation to corner (any neighbor shares a corner)
03440   BigReal max_cd = 0.;
03441   for ( i=0; i<2; ++i ) {
03442     for ( int j=0; j<2; ++j ) {
03443       for ( int k=0; k<2; ++k ) {
03444         ScaledPosition corner(  i ? min.x : max.x ,
03445                                 j ? min.y : max.y ,
03446                                 k ? min.z : max.z );
03447         Vector corner_delta =
03448                 lattice.unscale(corner) - lattice_old.unscale(corner);
03449         corner_delta -= center_delta;
03450         BigReal cd = corner_delta.length2();
03451         if ( cd > max_cd ) max_cd = cd;
03452       }
03453     }
03454   }
03455   max_cd = sqrt(max_cd);
03456 
03457   // find max deviation of atoms relative to center
03458   BigReal max_pd = 0.;
03459   CompAtom *p_old_i = doPairlistCheck_positions.begin();
03460   for ( i=0; i<n; ++i ) {
03461     Vector p_delta = p_i[i].position - p_old_i[i].position;
03462     p_delta -= center_delta;
03463     BigReal pd = p_delta.length2();
03464     if ( pd > max_pd ) max_pd = pd;
03465   }
03466   max_pd = sqrt(max_pd);
03467 
03468   BigReal max_tol = max_pd + max_cd;
03469 
03470   flags.maxAtomMovement = max_tol;
03471 
03472   // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
03473 
03474   if ( max_tol > ( (1. - simParams->pairlistTrigger) *
03475                                 doPairlistCheck_newTolerance ) ) {
03476     doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
03477   }
03478 
03479   if ( max_tol > doPairlistCheck_newTolerance ) {
03480     doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
03481   }
03482 
03483 }

void HomePatch::exchangeAtoms ( int  scriptTask  ) 

Definition at line 3358 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

03358                                             {
03359   SimParameters *simParams = Node::Object()->simParameters;
03360   // CkPrintf("exchangeAtoms %d %d %d %d\n", CmiMyPartition(), scriptTask, (int)(simParams->scriptArg1), (int)(simParams->scriptArg2));
03361   if ( scriptTask == SCRIPT_ATOMSEND || scriptTask == SCRIPT_ATOMSENDRECV ) {
03362     exchange_dst = (int) simParams->scriptArg1;
03363     // create and save outgoing message
03364     exchange_msg = new (numAtoms,0) ExchangeAtomsMsg;
03365     exchange_msg->lattice = lattice;
03366     exchange_msg->pid = patchID;
03367     exchange_msg->numAtoms = numAtoms;
03368     memcpy(exchange_msg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03369     if ( exchange_req >= 0 ) {
03370       recvExchangeReq(exchange_req);
03371     }
03372   }
03373   if ( scriptTask == SCRIPT_ATOMRECV || scriptTask == SCRIPT_ATOMSENDRECV ) {
03374     exchange_src = (int) simParams->scriptArg2;
03375     PatchMgr::Object()->sendExchangeReq(patchID, exchange_src);
03376   }
03377 }

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

Definition at line 3262 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

03262                                                            {  // initiating replica
03263   SimParameters *simParams = Node::Object()->simParameters;
03264   checkpoint_task = scriptTask;
03265   const int remote = simParams->scriptIntArg1;
03266   const char *key = simParams->scriptStringArg1;
03267   PatchMgr::Object()->sendCheckpointReq(patchID, remote, key, scriptTask);
03268 }

void HomePatch::gbisComputeAfterP1 (  ) 

Definition at line 2931 of file HomePatch.C.

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

02931                                    {
02932 
02933   SimParameters *simParams = Node::Object()->simParameters;
02934   BigReal alphaMax = simParams->alpha_max;
02935   BigReal delta = simParams->gbis_delta;
02936   BigReal beta = simParams->gbis_beta;
02937   BigReal gamma = simParams->gbis_gamma;
02938   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
02939 
02940   BigReal rhoi;
02941   BigReal rhoi0;
02942   //calculate bornRad from psiSum
02943   for (int i = 0; i < numAtoms; i++) {
02944     rhoi0 = intRad[2*i];
02945     rhoi = rhoi0+coulomb_radius_offset;
02946     psiFin[i] += psiSum[i];
02947     psiFin[i] *= rhoi0;
02948     bornRad[i]=1/(1/rhoi0-1/rhoi*tanh(psiFin[i]*(delta+psiFin[i]*(-beta+gamma*psiFin[i]))));
02949     bornRad[i] = (bornRad[i] > alphaMax) ? alphaMax : bornRad[i];
02950 #ifdef PRINT_COMP
02951     CkPrintf("BORNRAD(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,bornRad[i]);
02952 #endif
02953   }
02954 
02955   gbisP2Ready();
02956 }

void HomePatch::gbisComputeAfterP2 (  ) 

Definition at line 2959 of file HomePatch.C.

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

02959                                    {
02960 
02961   SimParameters *simParams = Node::Object()->simParameters;
02962   BigReal delta = simParams->gbis_delta;
02963   BigReal beta = simParams->gbis_beta;
02964   BigReal gamma = simParams->gbis_gamma;
02965   BigReal epsilon_s = simParams->solvent_dielectric;
02966   BigReal epsilon_p = simParams->dielectric;
02967   BigReal epsilon_s_i = 1/simParams->solvent_dielectric;
02968   BigReal epsilon_p_i = 1/simParams->dielectric;
02969   BigReal coulomb_radius_offset = simParams->coulomb_radius_offset;
02970   BigReal kappa = simParams->kappa;
02971   BigReal fij, expkappa, Dij, dEdai, dedasum;
02972   BigReal rhoi, rhoi0, psii, nbetapsi;
02973   BigReal gammapsi2, tanhi, daidr;
02974   for (int i = 0; i < numAtoms; i++) {
02975     //add diagonal dEda term
02976     dHdrPrefix[i] += dEdaSum[i];//accumulated from proxies
02977     fij = bornRad[i];//inf
02978     expkappa = exp(-kappa*fij);//0
02979     Dij = epsilon_p_i - expkappa*epsilon_s_i;//dielectric term
02980     //calculate dHij prefix
02981     dEdai = -0.5*COULOMB*atom[i].charge*atom[i].charge
02982                   *(kappa*epsilon_s_i*expkappa-Dij/fij)/bornRad[i];
02983     dHdrPrefix[i] += dEdai;
02984     dedasum = dHdrPrefix[i];
02985 
02986     rhoi0 = intRad[2*i];
02987     rhoi = rhoi0+coulomb_radius_offset;
02988     psii = psiFin[i];
02989     nbetapsi = -beta*psii;
02990     gammapsi2 = gamma*psii*psii;
02991     tanhi = tanh(psii*(delta+nbetapsi+gammapsi2));
02992     daidr = bornRad[i]*bornRad[i]*rhoi0/rhoi*(1-tanhi*tanhi)
02993            * (delta+nbetapsi+nbetapsi+gammapsi2+gammapsi2+gammapsi2);
02994     dHdrPrefix[i] *= daidr;//dHdrPrefix previously equaled dEda
02995 #ifdef PRINT_COMP
02996     CkPrintf("DHDR(%04i)[%04i] = % .4e\n",flags.sequence,pExt[i].id,dHdrPrefix[i]);
02997 #endif
02998   }
02999   gbisP3Ready();
03000 }

void HomePatch::gbisP2Ready (  ) 

Reimplemented from Patch.

Definition at line 3003 of file HomePatch.C.

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

Referenced by gbisComputeAfterP1().

03003                             {
03004   if (proxy.size() > 0) {
03005     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03006     for (int i = 0; i < proxy.size(); i++) {
03007       int node = proxy[i];
03008       ProxyGBISP2DataMsg *msg=new(numAtoms,PRIORITY_SIZE) ProxyGBISP2DataMsg;
03009       msg->patch = patchID;
03010       msg->origPe = CkMyPe();
03011       memcpy(msg->bornRad,bornRad.begin(),numAtoms*sizeof(Real));
03012       msg->destPe = node;
03013       int seq = flags.sequence;
03014       int priority = GB2_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03015       SET_PRIORITY(msg,seq,priority);
03016       cp[node].recvData(msg);
03017     }
03018   }
03019   Patch::gbisP2Ready();
03020 }

void HomePatch::gbisP3Ready (  ) 

Reimplemented from Patch.

Definition at line 3023 of file HomePatch.C.

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

Referenced by gbisComputeAfterP2().

03023                             {
03024   if (proxy.size() > 0) {
03025     CProxy_ProxyMgr cp(CkpvAccess(BOCclass_group).proxyMgr);
03026     //only nonzero message should be sent for doFullElec
03027     int msgAtoms = (flags.doFullElectrostatics) ? numAtoms : 0;
03028     for (int i = 0; i < proxy.size(); i++) {
03029       int node = proxy[i];
03030       ProxyGBISP3DataMsg *msg = new(msgAtoms,PRIORITY_SIZE) ProxyGBISP3DataMsg;
03031       msg->patch = patchID;
03032       msg->dHdrPrefixLen = msgAtoms;
03033       msg->origPe = CkMyPe();
03034       memcpy(msg->dHdrPrefix, dHdrPrefix.begin(), msgAtoms*sizeof(Real));
03035       msg->destPe = node;
03036       int seq = flags.sequence;
03037       int priority = GB3_PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
03038       SET_PRIORITY(msg,seq,priority);
03039       cp[node].recvData(msg);
03040     }
03041   }
03042   Patch::gbisP3Ready();
03043 }

FullAtomList& HomePatch::getAtomList (  )  [inline]

Definition at line 270 of file HomePatch.h.

Referenced by ComputeHomePatch::doWork().

00270 { return (atom); }

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

Definition at line 1806 of file HomePatch.C.

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

Referenced by Sequencer::hardWallDrude().

01808 {
01809   Molecule *mol = Node::Object()->molecule;
01810   SimParameters *simParams = Node::Object()->simParameters;
01811   const BigReal kbt=BOLTZMANN*simParams->drudeTemp;
01812   const int fixedAtomsOn = simParams->fixedAtomsOn;
01813   const BigReal dt = timestep / TIMEFACTOR;
01814   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
01815   int i, ia, ib, j;
01816   int dieOnError = simParams->rigidDie;
01817   Tensor wc;  // constraint virial
01818   BigReal idz, zmin, delta_T, maxtime=timestep,v_Bond;
01819   int nslabs;
01820 
01821   // start data for hard wall boundary between drude and its host atom
01822   // static int Count=0;
01823   int Idx;
01824   double r_wall, r_wall_SQ, rab, rab_SQ, dr, mass_a, mass_b, mass_sum;
01825   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;
01826   double dot_v_r_1, dot_v_r_2;
01827   double vb_cm, dr_a, dr_b;
01828   // end data for hard wall boundary between drude and its host atom
01829 
01830   // start calculation of hard wall boundary between drude and its host atom
01831   if (simParams->drudeHardWallOn) {
01832     if (ppreduction) {
01833       nslabs = simParams->pressureProfileSlabs;
01834       idz = nslabs/lattice.c().z;
01835       zmin = lattice.origin().z - 0.5*lattice.c().z;
01836     }
01837 
01838     r_wall = simParams->drudeBondLen;
01839     r_wall_SQ = r_wall*r_wall;
01840     // Count++;
01841     for (i=1; i<numAtoms; i++)  {
01842       if ( (atom[i].mass > 0.01) && ((atom[i].mass < 1.0)) ) { // drude particle
01843         ia = i-1;
01844         ib = i;
01845 
01846         v_ab = atom[ib].position - atom[ia].position;
01847         rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
01848 
01849         if (rab_SQ > r_wall_SQ) {  // to impose the hard wall constraint
01850           rab = sqrt(rab_SQ);
01851           if ( (rab > (2.0*r_wall)) && dieOnError ) {  // unexpected situation
01852             iout << iERROR << "HardWallDrude> "
01853               << "The drude is too far away from atom "
01854               << (atom[ia].id + 1) << " d = " << rab << "!\n" << endi;
01855             return -1;  // triggers early exit
01856           }
01857 
01858           v_ab.x /= rab;
01859           v_ab.y /= rab;
01860           v_ab.z /= rab;
01861 
01862           if ( fixedAtomsOn && atom[ia].atomFixed ) {  // the heavy atom is fixed
01863             if (atom[ib].atomFixed) {  // the drude is fixed too
01864               continue;
01865             }
01866             else {  // only the heavy atom is fixed
01867               dot_v_r_2 = atom[ib].velocity.x*v_ab.x
01868                 + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
01869               vb_2 = dot_v_r_2 * v_ab;
01870               vp_2 = atom[ib].velocity - vb_2;
01871 
01872               dr = rab - r_wall;
01873               if(dot_v_r_2 == 0.0) {
01874                 delta_T = maxtime;
01875               }
01876               else {
01877                 delta_T = dr/fabs(dot_v_r_2); // the time since the collision occurs
01878                 if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
01879               }
01880 
01881               dot_v_r_2 = -dot_v_r_2*sqrt(kbt/atom[ib].mass)/fabs(dot_v_r_2);
01882 
01883               vb_2 = dot_v_r_2 * v_ab;
01884 
01885               new_vel_a = atom[ia].velocity;
01886               new_vel_b = vp_2 + vb_2;
01887 
01888               dr_b = -dr + delta_T*dot_v_r_2;  // L = L_0 + dT *v_new, v was flipped
01889 
01890               new_pos_a = atom[ia].position;
01891               new_pos_b = atom[ib].position + dr_b*v_ab; // correct the position
01892             }
01893           }
01894           else {
01895             mass_a = atom[ia].mass;
01896             mass_b = atom[ib].mass;
01897             mass_sum = mass_a+mass_b;
01898 
01899             dot_v_r_1 = atom[ia].velocity.x*v_ab.x
01900               + atom[ia].velocity.y*v_ab.y + atom[ia].velocity.z*v_ab.z;
01901             vb_1 = dot_v_r_1 * v_ab;
01902             vp_1 = atom[ia].velocity - vb_1;
01903 
01904             dot_v_r_2 = atom[ib].velocity.x*v_ab.x
01905               + atom[ib].velocity.y*v_ab.y + atom[ib].velocity.z*v_ab.z;
01906             vb_2 = dot_v_r_2 * v_ab;
01907             vp_2 = atom[ib].velocity - vb_2;
01908 
01909             vb_cm = (mass_a*dot_v_r_1 + mass_b*dot_v_r_2)/mass_sum;
01910 
01911             dot_v_r_1 -= vb_cm;
01912             dot_v_r_2 -= vb_cm;
01913 
01914             dr = rab - r_wall;
01915 
01916             if(dot_v_r_2 == dot_v_r_1) {
01917                 delta_T = maxtime;
01918             }
01919             else {
01920               delta_T = dr/fabs(dot_v_r_2 - dot_v_r_1);  // the time since the collision occurs
01921               if(delta_T > maxtime ) delta_T = maxtime; // make sure it is not crazy
01922             }
01923             
01924             // the relative velocity between ia and ib. Drawn according to T_Drude
01925             v_Bond = sqrt(kbt/mass_b);
01926 
01927             // reflect the velocity along bond vector and scale down
01928             dot_v_r_1 = -dot_v_r_1*v_Bond*mass_b/(fabs(dot_v_r_1)*mass_sum);
01929             dot_v_r_2 = -dot_v_r_2*v_Bond*mass_a/(fabs(dot_v_r_2)*mass_sum);
01930 
01931             dr_a = dr*mass_b/mass_sum + delta_T*dot_v_r_1;
01932             dr_b = -dr*mass_a/mass_sum + delta_T*dot_v_r_2;
01933 
01934             new_pos_a = atom[ia].position + dr_a*v_ab;  // correct the position
01935             new_pos_b = atom[ib].position + dr_b*v_ab;
01936             // atom[ia].position += (dr_a*v_ab);  // correct the position
01937             // atom[ib].position += (dr_b*v_ab);
01938             
01939             dot_v_r_1 += vb_cm;
01940             dot_v_r_2 += vb_cm;
01941 
01942             vb_1 = dot_v_r_1 * v_ab;
01943             vb_2 = dot_v_r_2 * v_ab;
01944 
01945             new_vel_a = vp_1 + vb_1;
01946             new_vel_b = vp_2 + vb_2;
01947           }
01948 
01949           int ppoffset, partition;
01950           if ( invdt == 0 ) {
01951             atom[ia].position = new_pos_a;
01952             atom[ib].position = new_pos_b;
01953           }
01954           else if ( virial == 0 ) {
01955             atom[ia].velocity = new_vel_a;
01956             atom[ib].velocity = new_vel_b;
01957           }
01958           else {
01959             for ( j = 0; j < 2; j++ ) {
01960               if (j==0) {  // atom ia, heavy atom
01961                 Idx = ia;
01962                 new_pos = &new_pos_a;
01963                 new_vel = &new_vel_a;
01964               }
01965               else if (j==1) {  // atom ib, drude
01966                 Idx = ib;
01967                 new_pos = &new_pos_b;
01968                 new_vel = &new_vel_b;
01969               }
01970               Force df = (*new_vel - atom[Idx].velocity) *
01971                 ( atom[Idx].mass * invdt );
01972               Tensor vir = outer(df, atom[Idx].position);
01973               wc += vir;
01974               atom[Idx].velocity = *new_vel;
01975               atom[Idx].position = *new_pos;
01976 
01977               if (ppreduction) {
01978                 if (!j) {
01979                   BigReal z = new_pos->z;
01980                   int partition = atom[Idx].partition;
01981                   int slab = (int)floor((z-zmin)*idz);
01982                   if (slab < 0) slab += nslabs;
01983                   else if (slab >= nslabs) slab -= nslabs;
01984                   ppoffset = 3*(slab + nslabs*partition);
01985                 }
01986                 ppreduction->item(ppoffset  ) += vir.xx;
01987                 ppreduction->item(ppoffset+1) += vir.yy;
01988                 ppreduction->item(ppoffset+2) += vir.zz;
01989               }
01990 
01991             }
01992           }
01993         }                               
01994       }
01995     }
01996 
01997     // if ( (Count>10000) && (Count%10==0) ) {
01998     //   v_ab = atom[1].position - atom[0].position;
01999     //   rab_SQ = v_ab.x*v_ab.x + v_ab.y*v_ab.y + v_ab.z*v_ab.z;
02000     //   iout << "DBG_R: " << Count << "  " << sqrt(rab_SQ) << "\n" << endi;
02001     // }
02002 
02003   }
02004 
02005   // end calculation of hard wall boundary between drude and its host atom
02006 
02007   if ( dt && virial ) *virial += wc;
02008 
02009   return 0;
02010 }

void HomePatch::loweAndersenFinish (  ) 

Definition at line 2897 of file HomePatch.C.

References DebugM.

02898 {
02899     DebugM(2, "loweAndersenFinish\n");
02900     v.resize(0);
02901 }

void HomePatch::loweAndersenVelocities (  ) 

Definition at line 2882 of file HomePatch.C.

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

Referenced by positionsReady().

02883 {
02884     DebugM(2, "loweAndersenVelocities\n");
02885     Molecule *mol = Node::Object()->molecule;
02886     SimParameters *simParams = Node::Object()->simParameters;
02887     v.resize(numAtoms);
02888     for (int i = 0; i < numAtoms; ++i) {
02889         //v[i] = p[i];
02890         // co-opt CompAtom structure to pass velocity and mass information
02891         v[i].position = atom[i].velocity;
02892         v[i].charge = atom[i].mass;
02893     }
02894     DebugM(2, "loweAndersenVelocities\n");
02895 }

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

Definition at line 2747 of file HomePatch.C.

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

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

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

void HomePatch::mollyAverage (  ) 

Definition at line 3097 of file HomePatch.C.

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

Referenced by positionsReady().

03098 {
03099   Molecule *mol = Node::Object()->molecule;
03100   SimParameters *simParams = Node::Object()->simParameters;
03101   BigReal tol = simParams->mollyTol;
03102   int maxiter = simParams->mollyIter;
03103   int i, iter;
03104   HGArrayBigReal dsq;
03105   BigReal tmp;
03106   HGArrayInt ial, ibl;
03107   HGArrayVector ref;  // reference position
03108   HGArrayVector refab;  // reference vector
03109   HGArrayBigReal rmass;  // 1 / mass
03110   BigReal *lambda;  // Lagrange multipliers
03111   CompAtom *avg;  // averaged position
03112   int numLambdas = 0;
03113   HGArrayInt fixed;  // is atom fixed?
03114 
03115   //  iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
03116   p_avg.resize(numAtoms);
03117   for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
03118 
03119   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03120     int hgs = atom[ig].hydrogenGroupSize;
03121     if ( hgs == 1 ) continue;  // only one atom in group
03122         for ( i = 0; i < hgs; ++i ) {
03123           ref[i] = atom[ig+i].position;
03124           rmass[i] = 1. / atom[ig+i].mass;
03125           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03126           if ( fixed[i] ) rmass[i] = 0.;
03127         }
03128         avg = &(p_avg[ig]);
03129         int icnt = 0;
03130 
03131   if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
03132           if ( hgs != 3 ) {
03133             NAMD_die("Hydrogen group error caught in mollyAverage().  It's a bug!\n");
03134           }
03135           if ( !(fixed[1] && fixed[2]) ) {
03136             dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03137           }
03138         }
03139         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03140     if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
03141             if ( !(fixed[0] && fixed[i]) ) {
03142               dsq[icnt] =  tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03143             }
03144           }
03145         }
03146         if ( icnt == 0 ) continue;  // no constraints
03147         numLambdas += icnt;
03148         molly_lambda.resize(numLambdas);
03149         lambda = &(molly_lambda[numLambdas - icnt]);
03150         for ( i = 0; i < icnt; ++i ) {
03151           refab[i] = ref[ial[i]] - ref[ibl[i]];
03152         }
03153         //      iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
03154         iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
03155         if ( iter == maxiter ) {
03156           iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
03157         }
03158   }
03159 
03160   // for ( i=0; i<numAtoms; ++i ) {
03161   //    if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
03162   //      iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
03163   //        << p[i].position << " to " << p_avg[i].position << "\n" << endi;
03164   //    }
03165   // }
03166 
03167 }

void HomePatch::mollyMollify ( Tensor virial  ) 

Definition at line 3171 of file HomePatch.C.

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

03172 {
03173   Molecule *mol = Node::Object()->molecule;
03174   SimParameters *simParams = Node::Object()->simParameters;
03175   Tensor wc;  // constraint virial
03176   int i;
03177   HGArrayInt ial, ibl;
03178   HGArrayVector ref;  // reference position
03179   CompAtom *avg;  // averaged position
03180   HGArrayVector refab;  // reference vector
03181   HGArrayVector force;  // new force
03182   HGArrayBigReal rmass;  // 1 / mass
03183   BigReal *lambda;  // Lagrange multipliers
03184   int numLambdas = 0;
03185   HGArrayInt fixed;  // is atom fixed?
03186 
03187   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
03188     int hgs = atom[ig].hydrogenGroupSize;
03189     if (hgs == 1 ) continue;  // only one atom in group
03190         for ( i = 0; i < hgs; ++i ) {
03191           ref[i] = atom[ig+i].position;
03192           force[i] = f[Results::slow][ig+i];
03193           rmass[i] = 1. / atom[ig+i].mass;
03194           fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
03195           if ( fixed[i] ) rmass[i] = 0.;
03196         }
03197         int icnt = 0;
03198         // c-ji I'm only going to mollify water for now
03199   if ( atom[ig].rigidBondLength > 0 ) {  // for water
03200           if ( hgs != 3 ) {
03201             NAMD_die("Hydrogen group error caught in mollyMollify().  It's a bug!\n");
03202           }
03203           if ( !(fixed[1] && fixed[2]) ) {
03204             ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
03205           }
03206         }
03207         for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
03208     if ( atom[ig+i].rigidBondLength > 0 ) {
03209             if ( !(fixed[0] && fixed[i]) ) {
03210               ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
03211             }
03212           }
03213         }
03214 
03215         if ( icnt == 0 ) continue;  // no constraints
03216         lambda = &(molly_lambda[numLambdas]);
03217         numLambdas += icnt;
03218         for ( i = 0; i < icnt; ++i ) {
03219           refab[i] = ref[ial[i]] - ref[ibl[i]];
03220         }
03221         avg = &(p_avg[ig]);
03222         mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
03223         // store data back to patch
03224         for ( i = 0; i < hgs; ++i ) {
03225           wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
03226           f[Results::slow][ig+i] = force[i];
03227         }
03228   }
03229   // check that there isn't a constant needed here!
03230   *virial += wc;
03231   p_avg.resize(0);
03232 }

void HomePatch::positionsReady ( int  doMigration = 0  ) 

Reimplemented from Patch.

Definition at line 908 of file HomePatch.C.

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

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

void HomePatch::qmSwapAtoms (  ) 

Definition at line 849 of file HomePatch.C.

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

Referenced by positionsReady().

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

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

Definition at line 2164 of file HomePatch.C.

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

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

02165                                 {
02166 
02167   SimParameters *simParams = Node::Object()->simParameters;
02168   if (simParams->watmodel != WAT_TIP3 || ppreduction) {
02169     // Call old rattle1 -method instead
02170     return rattle1old(timestep, virial, ppreduction);
02171   }
02172 
02173   if (!rattleListValid) {
02174     buildRattleList();
02175     rattleListValid = true;
02176   }
02177 
02178   const int fixedAtomsOn = simParams->fixedAtomsOn;
02179   const int useSettle = simParams->useSettle;
02180   const BigReal dt = timestep / TIMEFACTOR;
02181   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02182   const BigReal tol2 = 2.0 * simParams->rigidTol;
02183   int maxiter = simParams->rigidIter;
02184   int dieOnError = simParams->rigidDie;
02185 
02186   Vector ref[10];  // reference position
02187   Vector pos[10];  // new position
02188   Vector vel[10];  // new velocity
02189 
02190   // Manual un-roll
02191   int n = (settleList.size()/2)*2;
02192   for (int j=0;j < n;j+=2) {
02193     int ig;
02194     ig = settleList[j];
02195     for (int i = 0; i < 3; ++i ) {
02196       ref[i] = atom[ig+i].position;
02197       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02198     }
02199     ig = settleList[j+1];
02200     for (int i = 0; i < 3; ++i ) {
02201       ref[i+3] = atom[ig+i].position;
02202       pos[i+3] = atom[ig+i].position + atom[ig+i].velocity * dt;
02203     }
02204     settle1_SIMD<2>(ref, pos,
02205       settle_mOrmT, settle_mHrmT, settle_ra,
02206       settle_rb, settle_rc, settle_rra);
02207 
02208     ig = settleList[j];
02209     for (int i = 0; i < 3; ++i ) {
02210       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02211       posNew[ig+i] = pos[i];
02212     }
02213     ig = settleList[j+1];
02214     for (int i = 0; i < 3; ++i ) {
02215       velNew[ig+i] = (pos[i+3] - ref[i+3])*invdt;
02216       posNew[ig+i] = pos[i+3];
02217     }
02218 
02219   }
02220 
02221   if (settleList.size() % 2) {
02222     int ig = settleList[settleList.size()-1];
02223     for (int i = 0; i < 3; ++i ) {
02224       ref[i] = atom[ig+i].position;
02225       pos[i] = atom[ig+i].position + atom[ig+i].velocity * dt;
02226     }
02227     settle1_SIMD<1>(ref, pos,
02228             settle_mOrmT, settle_mHrmT, settle_ra,
02229             settle_rb, settle_rc, settle_rra);        
02230     for (int i = 0; i < 3; ++i ) {
02231       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02232       posNew[ig+i] = pos[i];
02233     }
02234   }
02235 
02236   int posParam = 0;
02237   for (int j=0;j < rattleList.size();++j) {
02238 
02239     BigReal refx[10];
02240     BigReal refy[10];
02241     BigReal refz[10];
02242 
02243     BigReal posx[10];
02244     BigReal posy[10];
02245     BigReal posz[10];
02246 
02247     int ig = rattleList[j].ig;
02248     int icnt = rattleList[j].icnt;
02249     int hgs = atom[ig].hydrogenGroupSize;
02250     for (int i = 0; i < hgs; ++i ) {
02251       ref[i] = atom[ig+i].position;
02252       pos[i] = atom[ig+i].position;
02253       if (!(fixedAtomsOn && atom[ig+i].atomFixed))
02254         pos[i] += atom[ig+i].velocity * dt;
02255       refx[i] = ref[i].x;
02256       refy[i] = ref[i].y;
02257       refz[i] = ref[i].z;
02258       posx[i] = pos[i].x;
02259       posy[i] = pos[i].y;
02260       posz[i] = pos[i].z;
02261     }
02262 
02263 
02264     bool done;
02265     bool consFailure;
02266     if (icnt == 1) {
02267       rattlePair<1>(&rattleParam[posParam],
02268         refx, refy, refz,
02269         posx, posy, posz);
02270       done = true;
02271       consFailure = false;
02272     } else {
02273       rattleN(icnt, &rattleParam[posParam],
02274         refx, refy, refz,
02275         posx, posy, posz,
02276         tol2, maxiter,
02277         done, consFailure);
02278     }
02279 
02280     // Advance position in rattleParam
02281     posParam += icnt;
02282 
02283     for (int i = 0; i < hgs; ++i ) {
02284       pos[i].x = posx[i];
02285       pos[i].y = posy[i];
02286       pos[i].z = posz[i];
02287     }
02288 
02289     for (int i = 0; i < hgs; ++i ) {
02290       velNew[ig+i] = (pos[i] - ref[i])*invdt;
02291       posNew[ig+i] = pos[i];
02292     }
02293 
02294     if ( consFailure ) {
02295       if ( dieOnError ) {
02296         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02297         << (atom[ig].id + 1) << "!\n" << endi;
02298         return -1;  // triggers early exit
02299       } else {
02300         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02301         << (atom[ig].id + 1) << "!\n" << endi;
02302       }
02303     } else if ( ! done ) {
02304       if ( dieOnError ) {
02305         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02306         << (atom[ig].id + 1) << "!\n" << endi;
02307         return -1;  // triggers early exit
02308       } else {
02309         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02310         << (atom[ig].id + 1) << "!\n" << endi;
02311       }
02312     }
02313 
02314   }
02315 
02316   // Finally, we have to go through atoms that are not involved in rattle just so that we have
02317   // their positions and velocities up-to-date in posNew and velNew
02318   for (int j=0;j < noconstList.size();++j) {
02319     int ig = noconstList[j];
02320     int hgs = atom[ig].hydrogenGroupSize;
02321     for (int i = 0; i < hgs; ++i ) {
02322       velNew[ig+i] = atom[ig+i].velocity;
02323       posNew[ig+i] = atom[ig+i].position;
02324     }
02325   }
02326 
02327   if ( invdt == 0 ) {
02328     for (int ig = 0; ig < numAtoms; ++ig )
02329       atom[ig].position = posNew[ig];
02330   } else if ( virial == 0 ) {
02331     for (int ig = 0; ig < numAtoms; ++ig )
02332       atom[ig].velocity = velNew[ig];
02333   } else {
02334     Tensor wc;  // constraint virial
02335     addRattleForce(invdt, wc);
02336     *virial += wc;
02337   }
02338 
02339   return 0;
02340 }

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

Definition at line 2343 of file HomePatch.C.

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

Referenced by rattle1().

02345 {
02346   Molecule *mol = Node::Object()->molecule;
02347   SimParameters *simParams = Node::Object()->simParameters;
02348   const int fixedAtomsOn = simParams->fixedAtomsOn;
02349   const int useSettle = simParams->useSettle;
02350   const BigReal dt = timestep / TIMEFACTOR;
02351   const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
02352   BigReal tol2 = 2.0 * simParams->rigidTol;
02353   int maxiter = simParams->rigidIter;
02354   int dieOnError = simParams->rigidDie;
02355   int i, iter;
02356   BigReal dsq[10], tmp;
02357   int ial[10], ibl[10];
02358   Vector ref[10];  // reference position
02359   Vector refab[10];  // reference vector
02360   Vector pos[10];  // new position
02361   Vector vel[10];  // new velocity
02362   Vector netdp[10];  // total momentum change from constraint
02363   BigReal rmass[10];  // 1 / mass
02364   int fixed[10];  // is atom fixed?
02365   Tensor wc;  // constraint virial
02366   BigReal idz, zmin;
02367   int nslabs;
02368 
02369   // Initialize the settle algorithm with water parameters
02370   // settle1() assumes all waters are identical,
02371   // and will generate bad results if they are not.
02372   // XXX this will move to Molecule::build_atom_status when that 
02373   // version is debugged
02374   if ( ! settle_initialized ) {
02375     for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02376       // find a water
02377       if (atom[ig].rigidBondLength > 0) {
02378         int oatm;
02379         if (simParams->watmodel == WAT_SWM4) {
02380           oatm = ig+3;  // skip over Drude and Lonepair
02381           //printf("ig=%d  mass_ig=%g  oatm=%d  mass_oatm=%g\n",
02382           //    ig, atom[ig].mass, oatm, atom[oatm].mass);
02383         }
02384         else {
02385           oatm = ig+1;
02386           // Avoid using the Om site to set this by mistake
02387           if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
02388             oatm += 1;
02389           }
02390         }
02391 
02392         // initialize settle water parameters
02393         settle1init(atom[ig].mass, atom[oatm].mass, 
02394                     atom[ig].rigidBondLength,
02395                     atom[oatm].rigidBondLength,
02396                     settle_mOrmT, settle_mHrmT, settle_ra,
02397                     settle_rb, settle_rc, settle_rra);
02398         settle_initialized = 1;
02399         break; // done with init
02400       }
02401     }
02402   }
02403 
02404   if (ppreduction) {
02405     nslabs = simParams->pressureProfileSlabs;
02406     idz = nslabs/lattice.c().z;
02407     zmin = lattice.origin().z - 0.5*lattice.c().z;
02408   }
02409 
02410   // Size of a hydrogen group for water
02411   int wathgsize = 3;
02412   int watmodel = simParams->watmodel;
02413   if (watmodel == WAT_TIP4) wathgsize = 4;
02414   else if (watmodel == WAT_SWM4) wathgsize = 5;
02415   
02416   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02417     int hgs = atom[ig].hydrogenGroupSize;
02418     if ( hgs == 1 ) continue;  // only one atom in group
02419     // cache data in local arrays and integrate positions normally
02420     int anyfixed = 0;
02421     for ( i = 0; i < hgs; ++i ) {
02422       ref[i] = atom[ig+i].position;
02423       pos[i] = atom[ig+i].position;
02424       vel[i] = atom[ig+i].velocity;
02425       rmass[i] = (atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.);
02426       //printf("rmass of %i is %f\n", ig+i, rmass[i]);
02427       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02428       //printf("fixed status of %i is %i\n", i, fixed[i]);
02429       // undo addVelocityToPosition to get proper reference coordinates
02430       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; } else pos[i] += vel[i] * dt;
02431     }
02432     int icnt = 0;
02433     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02434       if (hgs != wathgsize) {
02435         char errmsg[256];
02436         sprintf(errmsg, "Water molecule starting with atom %d contains %d atoms "
02437                          "but the specified water model requires %d atoms.\n",
02438                          atom[ig].id+1, hgs, wathgsize);
02439         NAMD_die(errmsg);
02440       }
02441       // Use SETTLE for water unless some of the water atoms are fixed,
02442       if (useSettle && !anyfixed) {
02443         if (simParams->watmodel == WAT_SWM4) {
02444           // SWM4 ordering:  O D LP H1 H2
02445           // do swap(O,LP) and call settle with subarray O H1 H2
02446           // swap back after we return
02447           Vector lp_ref = ref[2];
02448           Vector lp_pos = pos[2];
02449           Vector lp_vel = vel[2];
02450           ref[2] = ref[0];
02451           pos[2] = pos[0];
02452           vel[2] = vel[0];
02453           settle1(ref+2, pos+2, vel+2, invdt,
02454                   settle_mOrmT, settle_mHrmT, settle_ra,
02455                   settle_rb, settle_rc, settle_rra);
02456           ref[0] = ref[2];
02457           pos[0] = pos[2];
02458           vel[0] = vel[2];
02459           ref[2] = lp_ref;
02460           pos[2] = lp_pos;
02461           vel[2] = lp_vel;
02462           // determine for LP updated pos and vel
02463           swm4_omrepos(ref, pos, vel, invdt);
02464         }
02465         else {
02466             settle1(ref, pos, vel, invdt,
02467                     settle_mOrmT, settle_mHrmT, settle_ra,
02468                     settle_rb, settle_rc, settle_rra);            
02469           if (simParams->watmodel == WAT_TIP4) {
02470             tip4_omrepos(ref, pos, vel, invdt);
02471           }
02472         }
02473 
02474         // which slab the hydrogen group will belong to
02475         // for pprofile calculations.
02476         int ppoffset, partition;
02477         if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02478           atom[ig+i].position = pos[i];
02479         } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
02480           atom[ig+i].velocity = vel[i];
02481         } else for ( i = 0; i < wathgsize; ++i ) {
02482           Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
02483           Tensor vir = outer(df, ref[i]);
02484           wc += vir;
02485           f[Results::normal][ig+i] += df;
02486           atom[ig+i].velocity = vel[i];
02487           if (ppreduction) {
02488             // put all the atoms from a water in the same slab.  Atom 0
02489             // should be the parent atom.
02490             if (!i) {
02491               BigReal z = pos[i].z;
02492               partition = atom[ig].partition;
02493               int slab = (int)floor((z-zmin)*idz);
02494               if (slab < 0) slab += nslabs;
02495               else if (slab >= nslabs) slab -= nslabs;
02496               ppoffset = 3*(slab + nslabs*partition);
02497             }
02498             ppreduction->item(ppoffset  ) += vir.xx;
02499             ppreduction->item(ppoffset+1) += vir.yy;
02500             ppreduction->item(ppoffset+2) += vir.zz;
02501           }
02502         }
02503         continue;
02504       }
02505       if ( !(fixed[1] && fixed[2]) ) {
02506         dsq[icnt] = tmp * tmp;  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02507       }
02508     }
02509     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02510       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02511         if ( !(fixed[0] && fixed[i]) ) {
02512           dsq[icnt] = tmp * tmp;  ial[icnt] = 0;  ibl[icnt] = i;  ++icnt;
02513         }
02514       }
02515     }
02516     if ( icnt == 0 ) continue;  // no constraints
02517     for ( i = 0; i < icnt; ++i ) {
02518       refab[i] = ref[ial[i]] - ref[ibl[i]];
02519     }
02520     for ( i = 0; i < hgs; ++i ) {
02521       netdp[i] = 0.;
02522     }
02523     int done;
02524     int consFailure;
02525     for ( iter = 0; iter < maxiter; ++iter ) {
02526 //if (iter > 0) CkPrintf("iteration %d\n", iter);
02527       done = 1;
02528       consFailure = 0;
02529       for ( i = 0; i < icnt; ++i ) {
02530         int a = ial[i];  int b = ibl[i];
02531         Vector pab = pos[a] - pos[b];
02532               BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
02533         BigReal rabsq = dsq[i];
02534         BigReal diffsq = rabsq - pabsq;
02535         if ( fabs(diffsq) > (rabsq * tol2) ) {
02536           Vector &rab = refab[i];
02537           BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
02538           if ( rpab < ( rabsq * 1.0e-6 ) ) {
02539             done = 0;
02540             consFailure = 1;
02541             continue;
02542           }
02543           BigReal rma = rmass[a];
02544           BigReal rmb = rmass[b];
02545           BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
02546           Vector dp = rab * gab;
02547           pos[a] += rma * dp;
02548           pos[b] -= rmb * dp;
02549           if ( invdt != 0. ) {
02550             dp *= invdt;
02551             netdp[a] += dp;
02552             netdp[b] -= dp;
02553           }
02554           done = 0;
02555         }
02556       }
02557       if ( done ) break;
02558     }
02559 
02560     if ( consFailure ) {
02561       if ( dieOnError ) {
02562         iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
02563         << (atom[ig].id + 1) << "!\n" << endi;
02564         return -1;  // triggers early exit
02565       } else {
02566         iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
02567         << (atom[ig].id + 1) << "!\n" << endi;
02568       }
02569     } else if ( ! done ) {
02570       if ( dieOnError ) {
02571         iout << iERROR << "Exceeded RATTLE iteration limit for atom "
02572         << (atom[ig].id + 1) << "!\n" << endi;
02573         return -1;  // triggers early exit
02574       } else {
02575         iout << iWARN << "Exceeded RATTLE iteration limit for atom "
02576         << (atom[ig].id + 1) << "!\n" << endi;
02577       }
02578     }
02579 
02580     // store data back to patch
02581     int ppoffset, partition;
02582     if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
02583       atom[ig+i].position = pos[i];
02584     } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
02585       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02586     } else for ( i = 0; i < hgs; ++i ) {
02587       Force df = netdp[i] * invdt;
02588       Tensor vir = outer(df, ref[i]);
02589       wc += vir;
02590       f[Results::normal][ig+i] += df;
02591       atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
02592       if (ppreduction) {
02593         if (!i) {
02594           BigReal z = pos[i].z;
02595           int partition = atom[ig].partition;
02596           int slab = (int)floor((z-zmin)*idz);
02597           if (slab < 0) slab += nslabs;
02598           else if (slab >= nslabs) slab -= nslabs;
02599           ppoffset = 3*(slab + nslabs*partition);
02600         }
02601         ppreduction->item(ppoffset  ) += vir.xx;
02602         ppreduction->item(ppoffset+1) += vir.yy;
02603         ppreduction->item(ppoffset+2) += vir.zz;
02604       }
02605     }
02606   }
02607   if ( dt && virial ) *virial += wc;
02608 
02609   return 0;
02610 }

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

Definition at line 2613 of file HomePatch.C.

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

02614 {
02615   Molecule *mol = Node::Object()->molecule;
02616   SimParameters *simParams = Node::Object()->simParameters;
02617   const int fixedAtomsOn = simParams->fixedAtomsOn;
02618   const int useSettle = simParams->useSettle;
02619   const BigReal dt = timestep / TIMEFACTOR;
02620   Tensor wc;  // constraint virial
02621   BigReal tol = simParams->rigidTol;
02622   int maxiter = simParams->rigidIter;
02623   int dieOnError = simParams->rigidDie;
02624   int i, iter;
02625   BigReal dsqi[10], tmp;
02626   int ial[10], ibl[10];
02627   Vector ref[10];  // reference position
02628   Vector refab[10];  // reference vector
02629   Vector vel[10];  // new velocity
02630   BigReal rmass[10];  // 1 / mass
02631   BigReal redmass[10];  // reduced mass
02632   int fixed[10];  // is atom fixed?
02633   
02634   // Size of a hydrogen group for water
02635   int wathgsize = 3;
02636   if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02637   else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02638 
02639   //  CkPrintf("In rattle2!\n");
02640   for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
02641     //    CkPrintf("ig=%d\n",ig);
02642     int hgs = atom[ig].hydrogenGroupSize;
02643     if ( hgs == 1 ) continue;  // only one atom in group
02644     // cache data in local arrays and integrate positions normally
02645     int anyfixed = 0;
02646     for ( i = 0; i < hgs; ++i ) {
02647       ref[i] = atom[ig+i].position;
02648       vel[i] = atom[ig+i].velocity;
02649       rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
02650       fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
02651       if ( fixed[i] ) { anyfixed = 1; rmass[i] = 0.; }
02652     }
02653     int icnt = 0;
02654     if ( ( tmp = atom[ig].rigidBondLength ) > 0 ) {  // for water
02655       if (hgs != wathgsize) {
02656         NAMD_bug("Hydrogen group error caught in rattle2().");
02657       }
02658       // Use SETTLE for water unless some of the water atoms are fixed,
02659       if (useSettle && !anyfixed) {
02660         if (simParams->watmodel == WAT_SWM4) {
02661           // SWM4 ordering:  O D LP H1 H2
02662           // do swap(O,LP) and call settle with subarray O H1 H2
02663           // swap back after we return
02664           Vector lp_ref = ref[2];
02665           // Vector lp_vel = vel[2];
02666           ref[2] = ref[0];
02667           vel[2] = vel[0];
02668           settle2(atom[ig].mass, atom[ig+3].mass, ref+2, vel+2, dt, virial);
02669           ref[0] = ref[2];
02670           vel[0] = vel[2];
02671           ref[2] = lp_ref;
02672           // vel[2] = vel[0];  // set LP vel to O vel
02673         }
02674         else {
02675           settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
02676           if (simParams->watmodel == WAT_TIP4) {
02677             vel[3] = vel[0];
02678           }
02679         }
02680         for (i=0; i<hgs; i++) {
02681           atom[ig+i].velocity = vel[i];
02682         }
02683         continue;
02684       }
02685       if ( !(fixed[1] && fixed[2]) ) {
02686         redmass[icnt] = 1. / (rmass[1] + rmass[2]);
02687         dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 1;  ibl[icnt] = 2;  ++icnt;
02688       }
02689     }
02690     //    CkPrintf("Loop 2\n");
02691     for ( i = 1; i < hgs; ++i ) {  // normal bonds to mother atom
02692       if ( ( tmp = atom[ig+i].rigidBondLength ) > 0 ) {
02693         if ( !(fixed[0] && fixed[i]) ) {
02694           redmass[icnt] = 1. / (rmass[0] + rmass[i]);
02695           dsqi[icnt] = 1. / (tmp * tmp);  ial[icnt] = 0;
02696           ibl[icnt] = i;  ++icnt;
02697         }
02698       }
02699     }
02700     if ( icnt == 0 ) continue;  // no constraints
02701     //    CkPrintf("Loop 3\n");
02702     for ( i = 0; i < icnt; ++i ) {
02703       refab[i] = ref[ial[i]] - ref[ibl[i]];
02704     }
02705     //    CkPrintf("Loop 4\n");
02706     int done;
02707     for ( iter = 0; iter < maxiter; ++iter ) {
02708       done = 1;
02709       for ( i = 0; i < icnt; ++i ) {
02710         int a = ial[i];  int b = ibl[i];
02711         Vector vab = vel[a] - vel[b];
02712         Vector &rab = refab[i];
02713         BigReal rabsqi = dsqi[i];
02714         BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
02715         if ( (fabs(rvab) * dt * rabsqi) > tol ) {
02716           Vector dp = rab * (-rvab * redmass[i] * rabsqi);
02717           wc += outer(dp,rab);
02718           vel[a] += rmass[a] * dp;
02719           vel[b] -= rmass[b] * dp;
02720           done = 0;
02721         }
02722       }
02723       if ( done ) break;
02724       //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
02725     }
02726     if ( ! done ) {
02727       if ( dieOnError ) {
02728         NAMD_die("Exceeded maximum number of iterations in rattle2().");
02729       } else {
02730         iout << iWARN <<
02731           "Exceeded maximum number of iterations in rattle2().\n" << endi;
02732       }
02733     }
02734     // store data back to patch
02735     for ( i = 0; i < hgs; ++i ) {
02736       atom[ig+i].velocity = vel[i];
02737     }
02738   }
02739   //  CkPrintf("Leaving rattle2!\n");
02740   // check that there isn't a constant needed here!
02741   *virial += wc / ( 0.5 * dt );
02742 
02743 }

void HomePatch::receiveResult ( ProxyGBISP2ResultMsg msg  ) 

Definition at line 3071 of file HomePatch.C.

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

03071                                                        {
03072   ++numGBISP2Arrived;
03073   //accumulate dEda
03074   for ( int i = 0; i < msg->dEdaSumLen; ++i ) {
03075     dHdrPrefix[i] += msg->dEdaSum[i];
03076   }
03077   delete msg;
03078 
03079   if (flags.doNonbonded) {
03080     //awaken if phase 2 done
03081     if (phase2BoxClosedCalled == true &&
03082         numGBISP2Arrived==proxy.size() ) {
03083       sequencer->awaken();
03084     }
03085   } else {
03086     //awaken if all phases done on noWork step
03087     if (boxesOpen == 0 &&
03088         numGBISP1Arrived == proxy.size() &&
03089         numGBISP2Arrived == proxy.size() &&
03090         numGBISP3Arrived == proxy.size()) {
03091       sequencer->awaken();
03092     }
03093   }
03094 }

void HomePatch::receiveResult ( ProxyGBISP1ResultMsg msg  ) 

Definition at line 3046 of file HomePatch.C.

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

Referenced by ProxyMgr::recvResult().

03046                                                        {
03047   ++numGBISP1Arrived;
03048     for ( int i = 0; i < msg->psiSumLen; ++i ) {
03049       psiFin[i] += msg->psiSum[i];
03050     }
03051   delete msg;
03052 
03053   if (flags.doNonbonded) {
03054     //awaken if phase 1 done
03055     if (phase1BoxClosedCalled == true &&
03056         numGBISP1Arrived==proxy.size() ) {
03057         sequencer->awaken();
03058     }
03059   } else {
03060     //awaken if all phases done on noWork step
03061     if (boxesOpen == 0 &&
03062         numGBISP1Arrived == proxy.size() &&
03063         numGBISP2Arrived == proxy.size() &&
03064         numGBISP3Arrived == proxy.size()) {
03065       sequencer->awaken();
03066     }
03067   }
03068 }

void HomePatch::receiveResults ( ProxyCombinedResultRawMsg msg  ) 

Definition at line 820 of file HomePatch.C.

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

00821 {
00822     numGBISP3Arrived++;
00823   DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodeSize<<")\n");
00824     Results *r = forceBox.clientOpen(msg->nodeSize);
00825       register char* isNonZero = msg->isForceNonZero;
00826       register Force* f_i = msg->forceArr;
00827       for ( int k = 0; k < Results::maxNumForces; ++k )
00828       {
00829         Force *f = r->f[k];
00830                 int nf = msg->flLen[k];
00831 #ifdef ARCH_POWERPC
00832 #pragma disjoint (*f_i, *f)
00833 #endif
00834         for (int count = 0; count < nf; count++) {
00835           if(*isNonZero){
00836                         f[count].x += f_i->x;
00837                         f[count].y += f_i->y;
00838                         f[count].z += f_i->z;
00839                         f_i++;
00840           }
00841           isNonZero++;
00842         }
00843       }
00844     forceBox.clientClose(msg->nodeSize);
00845 
00846     delete msg;
00847 }

void HomePatch::receiveResults ( ProxyResultMsg msg  ) 

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

00803                                                   {
00804   numGBISP3Arrived++;
00805   DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00806   int n = msg->node;
00807   Results *r = forceBox.clientOpen();
00808   for ( int k = 0; k < Results::maxNumForces; ++k )
00809   {
00810     Force *f = r->f[k];
00811     register ForceList::iterator f_i, f_e;
00812     f_i = msg->forceList[k]->begin();
00813     f_e = msg->forceList[k]->end();
00814     for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
00815   }
00816   forceBox.clientClose();
00817   delete msg;
00818 }

void HomePatch::receiveResults ( ProxyResultVarsizeMsg msg  ) 

Definition at line 779 of file HomePatch.C.

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

Referenced by ProxyMgr::recvResults().

00779                                                          {
00780 
00781     numGBISP3Arrived++;
00782     DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00783     int n = msg->node;
00784     Results *r = forceBox.clientOpen();
00785 
00786     char *iszeroPtr = msg->isZero;
00787     Force *msgFPtr = msg->forceArr;
00788 
00789     for ( int k = 0; k < Results::maxNumForces; ++k )
00790     {
00791       Force *rfPtr = r->f[k];
00792       for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
00793           if((*iszeroPtr)!=1) {
00794               *rfPtr += *msgFPtr;
00795               msgFPtr++;
00796           }
00797       }      
00798     }
00799     forceBox.clientClose();
00800     delete msg;
00801 }

void HomePatch::recvCheckpointAck (  ) 

Definition at line 3353 of file HomePatch.C.

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

03353                                   {  // initiating replica
03354   CkpvAccess(_qd)->process();
03355 }

void HomePatch::recvCheckpointLoad ( CheckpointAtomsMsg msg  ) 

Definition at line 3300 of file HomePatch.C.

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

Referenced by PatchMgr::recvCheckpointLoad().

03300                                                           {  // initiating replica
03301   SimParameters *simParams = Node::Object()->simParameters;
03302   const int remote = simParams->scriptIntArg1;
03303   const char *key = simParams->scriptStringArg1;
03304   if ( checkpoint_task == SCRIPT_CHECKPOINT_FREE ) {
03305     NAMD_bug("HomePatch::recvCheckpointLoad called during checkpointFree.");
03306   }
03307   if ( msg->replica != remote ) {
03308     NAMD_bug("HomePatch::recvCheckpointLoad message from wrong replica.");
03309   }
03310   if ( checkpoint_task == SCRIPT_CHECKPOINT_STORE || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03311     CheckpointAtomsMsg *newmsg = new (numAtoms,1+strlen(key),0) CheckpointAtomsMsg;
03312     strcpy(newmsg->key,key);
03313     newmsg->lattice = lattice;
03314     newmsg->berendsenPressure_count = sequencer->berendsenPressure_count;
03315     newmsg->pid = patchID;
03316     newmsg->pe = CkMyPe();
03317     newmsg->replica = CmiMyPartition();
03318     newmsg->numAtoms = numAtoms;
03319     memcpy(newmsg->atoms,atom.begin(),numAtoms*sizeof(FullAtom));
03320     PatchMgr::Object()->sendCheckpointStore(newmsg, remote, msg->pe);
03321   }
03322   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD || checkpoint_task == SCRIPT_CHECKPOINT_SWAP ) {
03323     atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03324     lattice = msg->lattice;
03325     sequencer->berendsenPressure_count = msg->berendsenPressure_count;
03326     numAtoms = msg->numAtoms;
03327     atom.resize(numAtoms);
03328     memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03329     doAtomUpdate = true;
03330     rattleListValid = false;
03331     if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03332   }
03333   if ( checkpoint_task == SCRIPT_CHECKPOINT_LOAD ) {
03334     recvCheckpointAck();
03335   }
03336   delete msg;
03337 }

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

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

03270                                                                                 {  // responding replica
03271   if ( task == SCRIPT_CHECKPOINT_FREE ) {
03272     if ( ! checkpoints.count(key) ) {
03273       NAMD_die("Unable to free checkpoint, requested key was never stored.");
03274     }
03275     delete checkpoints[key];
03276     checkpoints.erase(key);
03277     PatchMgr::Object()->sendCheckpointAck(patchID, replica, pe);
03278     return;
03279   }
03280   CheckpointAtomsMsg *msg;
03281   if ( task == SCRIPT_CHECKPOINT_LOAD || task == SCRIPT_CHECKPOINT_SWAP ) {
03282     if ( ! checkpoints.count(key) ) {
03283       NAMD_die("Unable to load checkpoint, requested key was never stored.");
03284     }
03285     checkpoint_t &cp = *checkpoints[key];
03286     msg = new (cp.numAtoms,1,0) CheckpointAtomsMsg;
03287     msg->lattice = cp.lattice;
03288     msg->berendsenPressure_count = cp.berendsenPressure_count;
03289     msg->numAtoms = cp.numAtoms;
03290     memcpy(msg->atoms,cp.atoms.begin(),cp.numAtoms*sizeof(FullAtom));
03291   } else {
03292     msg = new (0,1,0) CheckpointAtomsMsg;
03293   }
03294   msg->pid = patchID;
03295   msg->replica = CmiMyPartition();
03296   msg->pe = CkMyPe();
03297   PatchMgr::Object()->sendCheckpointLoad(msg, replica, pe);
03298 }

void HomePatch::recvCheckpointStore ( CheckpointAtomsMsg msg  ) 

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

03339                                                            {  // responding replica
03340   if ( ! checkpoints.count(msg->key) ) {
03341     checkpoints[msg->key] = new checkpoint_t;
03342   }
03343   checkpoint_t &cp = *checkpoints[msg->key];
03344   cp.lattice = msg->lattice;
03345   cp.berendsenPressure_count = msg->berendsenPressure_count;
03346   cp.numAtoms = msg->numAtoms;
03347   cp.atoms.resize(cp.numAtoms);
03348   memcpy(cp.atoms.begin(),msg->atoms,cp.numAtoms*sizeof(FullAtom));
03349   PatchMgr::Object()->sendCheckpointAck(patchID, msg->replica, msg->pe);
03350   delete msg;
03351 }

void HomePatch::recvExchangeMsg ( ExchangeAtomsMsg msg  ) 

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

03390                                                      {
03391   // CkPrintf("recvExchangeMsg %d %d\n", CmiMyPartition(), exchange_src);
03392   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03393   lattice = msg->lattice;
03394   numAtoms = msg->numAtoms;
03395   atom.resize(numAtoms);
03396   memcpy(atom.begin(),msg->atoms,numAtoms*sizeof(FullAtom));
03397   delete msg;
03398   CkpvAccess(_qd)->process();
03399   doAtomUpdate = true;
03400   rattleListValid = false;
03401   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03402 }

void HomePatch::recvExchangeReq ( int  req  ) 

Definition at line 3379 of file HomePatch.C.

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

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

03379                                        {
03380   exchange_req = req;
03381   if ( exchange_msg ) {
03382     // CkPrintf("recvExchangeReq %d %d\n", CmiMyPartition(), exchange_dst);
03383     PatchMgr::Object()->sendExchangeMsg(exchange_msg, exchange_dst, exchange_req);
03384     exchange_msg = 0;
03385     exchange_req = -1;
03386     CkpvAccess(_qd)->process();
03387   }
03388 }

void HomePatch::recvNodeAwareSpanningTree ( ProxyNodeAwareSpanningTreeMsg msg  ) 

Definition at line 613 of file HomePatch.C.

Referenced by ProxyMgr::recvNodeAwareSpanningTreeOnHomePatch().

00613 {}

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

Definition at line 619 of file HomePatch.C.

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

Referenced by ProxyMgr::recvSpanningTreeOnHomePatch().

00620 {
00621   int i;
00622   nChild=0;
00623   tree.resize(n);
00624   for (i=0; i<n; i++) {
00625     tree[i] = t[i];
00626   }
00627 
00628   for (i=1; i<=proxySpanDim; i++) {
00629     if (tree.size() <= i) break;
00630     child[i-1] = tree[i];
00631     nChild++;
00632   }
00633 
00634 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00635   isProxyChanged = 1;
00636 #endif
00637 
00638   // send down to children
00639   sendSpanningTree();
00640 }

void HomePatch::registerProxy ( RegisterProxyMsg  ) 

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

00385                                                    {
00386   DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
00387   proxy.add(msg->node);
00388   forceBox.clientAdd();
00389 
00390   isNewProxyAdded = 1;
00391 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00392   isProxyChanged = 1;
00393 #endif
00394 
00395   Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
00396   delete msg;
00397 }

void HomePatch::replaceForces ( ExtForce f  ) 

Definition at line 1259 of file HomePatch.C.

References Patch::f.

01260 {
01261   replacementForces = f;
01262 }

void HomePatch::revert ( void   ) 

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

03244                            {
03245   atomMapper->unregisterIDsFullAtom(atom.begin(),atom.end());
03246 
03247   atom.copy(checkpoint_atom);
03248   numAtoms = atom.size();
03249   lattice = checkpoint_lattice;
03250 
03251   doAtomUpdate = true;
03252   rattleListValid = false;
03253 
03254   if ( ! numNeighbors ) atomMapper->registerIDsFullAtom(atom.begin(),atom.end());
03255 
03256   // DMK - Atom Separation (water vs. non-water)
03257   #if NAMD_SeparateWaters != 0
03258     numWaterAtoms = checkpoint_numWaterAtoms;
03259   #endif
03260 }

void HomePatch::runSequencer ( void   ) 

Definition at line 252 of file HomePatch.C.

References Sequencer::run().

00253 { sequencer->run(); }

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

Definition at line 1264 of file HomePatch.C.

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

Referenced by Sequencer::saveForce().

01265 {
01266   f_saved[ftag].resize(numAtoms);
01267   for ( int i = 0; i < numAtoms; ++i )
01268   {
01269     f_saved[ftag][i] = f[ftag][i];
01270   }
01271 }

void HomePatch::sendNodeAwareSpanningTree (  ) 

Definition at line 614 of file HomePatch.C.

00614 {}

void HomePatch::sendProxies (  ) 

Definition at line 451 of file HomePatch.C.

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

00452 {
00453 #if USE_NODEPATCHMGR
00454         CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
00455         NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
00456         npm->sendProxyList(patchID, proxy.begin(), proxy.size());
00457 #else
00458         ProxyMgr::Object()->sendProxies(patchID, proxy.begin(), proxy.size());
00459 #endif
00460 }

void HomePatch::sendSpanningTree (  ) 

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

00643 {
00644   if (tree.size() == 0) return;
00645   ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
00646   msg->patch = patchID;
00647   msg->node = CkMyPe();
00648   msg->tree.copy(tree);  // copy data for thread safety
00649   ProxyMgr::Object()->sendSpanningTree(msg);  
00650 }

void HomePatch::setGBISIntrinsicRadii (  ) 

Definition at line 2916 of file HomePatch.C.

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

Referenced by positionsReady().

02916                                       {
02917   intRad.resize(numAtoms*2);
02918   intRad.setall(0);
02919   Molecule *mol = Node::Object()->molecule;
02920   SimParameters *simParams = Node::Object()->simParameters;
02921   Real offset = simParams->coulomb_radius_offset;
02922   for (int i = 0; i < numAtoms; i++) {
02923     Real rad = MassToRadius(atom[i].mass);//in ComputeGBIS.inl
02924     Real screen = MassToScreen(atom[i].mass);//same
02925     intRad[2*i+0] = rad - offset;//r0
02926     intRad[2*i+1] = screen*(rad - offset);//s0
02927   }
02928 }

void HomePatch::setLcpoType (  ) 

Definition at line 2905 of file HomePatch.C.

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

Referenced by positionsReady().

02905                             {
02906   Molecule *mol = Node::Object()->molecule;
02907   const int *lcpoParamType = mol->getLcpoParamType();
02908 
02909   lcpoType.resize(numAtoms);
02910   for (int i = 0; i < numAtoms; i++) {
02911     lcpoType[i] = lcpoParamType[pExt[i].id];
02912   }
02913 }

void HomePatch::submitLoadStats ( int  timestep  ) 

Definition at line 3404 of file HomePatch.C.

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

Referenced by Sequencer::rebalanceLoad().

03405 {
03406   LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
03407 }

void HomePatch::unregisterProxy ( UnregisterProxyMsg  ) 

Definition at line 399 of file HomePatch.C.

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

Referenced by ProxyMgr::recvUnregisterProxy().

00399                                                        {
00400 #if CMK_PERSISTENT_COMM && USE_PERSISTENT_TREE
00401   isProxyChanged = 1;
00402 #endif
00403   int n = msg->node;
00404   NodeID *pe = proxy.begin();
00405   for ( ; *pe != n; ++pe );
00406   forceBox.clientRemove();
00407   proxy.del(pe - proxy.begin());
00408   delete msg;
00409 }

void HomePatch::useSequencer ( Sequencer sequencerPtr  ) 

Definition at line 248 of file HomePatch.C.

00249 { sequencer=sequencerPtr; }


Friends And Related Function Documentation

friend class ComputeGlobal [friend]

Definition at line 91 of file HomePatch.h.

friend class PatchMgr [friend]

Definition at line 89 of file HomePatch.h.

friend class Sequencer [friend]

Definition at line 90 of file HomePatch.h.


Member Data Documentation

int HomePatch::checkpoint_task

Definition at line 243 of file HomePatch.h.

Referenced by exchangeCheckpoint(), and recvCheckpointLoad().

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

Definition at line 250 of file HomePatch.h.

Referenced by recvCheckpointReq(), and recvCheckpointStore().

int HomePatch::exchange_dst

Definition at line 256 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

ExchangeAtomsMsg* HomePatch::exchange_msg

Definition at line 259 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

int HomePatch::exchange_req

Definition at line 258 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

int HomePatch::exchange_src

Definition at line 257 of file HomePatch.h.

Referenced by exchangeAtoms().

int HomePatch::inMigration [protected]

Definition at line 303 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

LDObjHandle HomePatch::ldObjHandle

Definition at line 294 of file HomePatch.h.

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

int HomePatch::marginViolations

Definition at line 148 of file HomePatch.h.

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

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

Definition at line 305 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<int> HomePatch::noconstList

Definition at line 198 of file HomePatch.h.

Referenced by rattle1().

int HomePatch::numMlBuf [protected]

Definition at line 304 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<Vector> HomePatch::posNew

Definition at line 204 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

std::vector<RattleList> HomePatch::rattleList

Definition at line 196 of file HomePatch.h.

Referenced by rattle1().

bool HomePatch::rattleListValid

Definition at line 200 of file HomePatch.h.

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

std::vector<RattleParam> HomePatch::rattleParam

Definition at line 197 of file HomePatch.h.

Referenced by rattle1().

std::vector<int> HomePatch::settleList

Definition at line 195 of file HomePatch.h.

Referenced by rattle1().

std::vector<Vector> HomePatch::velNew

Definition at line 203 of file HomePatch.h.

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


The documentation for this class was generated from the following files:
Generated on Sat Jun 23 01:17:20 2018 for NAMD by  doxygen 1.4.7