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

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

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


Member Function Documentation

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

Definition at line 1824 of file HomePatch.C.

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

Referenced by Sequencer::addForceToMomentum().

01829       {
01830   SimParameters *simParams = Node::Object()->simParameters;
01831 
01832   if ( simParams->fixedAtomsOn ) {
01833     for ( int i = 0; i < num_atoms; ++i ) {
01834       if ( atom_arr[i].atomFixed ) {
01835         atom_arr[i].velocity = 0;
01836       } else {
01837         BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01838         atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01839         atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01840         atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01841       }
01842     }
01843   } else {
01844     for ( int i = 0; i < num_atoms; ++i ) {
01845       BigReal dt_mass = dt * atom_arr[i].recipMass;  // dt/mass
01846       atom_arr[i].velocity.x += force_arr[i].x * dt_mass;
01847       atom_arr[i].velocity.y += force_arr[i].y * dt_mass;
01848       atom_arr[i].velocity.z += force_arr[i].z * dt_mass;
01849     }
01850   }
01851 }

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

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

Referenced by Sequencer::addForceToMomentum3().

01862       {
01863   SimParameters *simParams = Node::Object()->simParameters;
01864 
01865   if ( simParams->fixedAtomsOn ) {
01866     for ( int i = 0; i < num_atoms; ++i ) {
01867       if ( atom_arr[i].atomFixed ) {
01868         atom_arr[i].velocity = 0;
01869       } else {
01870         BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01871         atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01872             + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01873         atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01874             + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01875         atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01876             + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01877       }
01878     }
01879   } else {
01880     for ( int i = 0; i < num_atoms; ++i ) {
01881       BigReal rmass = atom_arr[i].recipMass;  // 1/mass
01882       atom_arr[i].velocity.x += (force_arr1[i].x*dt1 
01883           + force_arr2[i].x*dt2 + force_arr3[i].x*dt3) * rmass;
01884       atom_arr[i].velocity.y += (force_arr1[i].y*dt1 
01885           + force_arr2[i].y*dt2 + force_arr3[i].y*dt3) * rmass;
01886       atom_arr[i].velocity.z += (force_arr1[i].z*dt1 
01887           + force_arr2[i].z*dt2 + force_arr3[i].z*dt3) * rmass;
01888     }
01889   }
01890 }

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

Definition at line 2263 of file HomePatch.C.

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

Referenced by rattle1().

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

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

Definition at line 1892 of file HomePatch.C.

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

Referenced by Sequencer::addVelocityToPosition().

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

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

Implements Patch.

Definition at line 426 of file HomePatch.C.

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

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

void HomePatch::buildRattleList (  ) 

Definition at line 2121 of file HomePatch.C.

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

Referenced by rattle1().

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

void HomePatch::buildSpanningTree ( void   ) 

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

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

void HomePatch::checkpoint ( void   ) 

Definition at line 3343 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

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

void HomePatch::depositMigration ( MigrateAtomsMsg  ) 

Definition at line 3852 of file HomePatch.C.

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

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

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

void HomePatch::doAtomMigration (  )  [protected]

Definition at line 3715 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::doGroupSizeCheck (  )  [protected]

Definition at line 3594 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::doMarginCheck (  )  [protected]

Definition at line 3638 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::doPairlistCheck (  )  [protected]

Definition at line 3519 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::exchangeAtoms ( int  scriptTask  ) 

Definition at line 3467 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

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

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

Definition at line 3371 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

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

void HomePatch::gbisComputeAfterP1 (  ) 

Definition at line 3040 of file HomePatch.C.

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

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

void HomePatch::gbisComputeAfterP2 (  ) 

Definition at line 3068 of file HomePatch.C.

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

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

void HomePatch::gbisP2Ready (  ) 

Reimplemented from Patch.

Definition at line 3112 of file HomePatch.C.

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

Referenced by gbisComputeAfterP1().

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

void HomePatch::gbisP3Ready (  ) 

Reimplemented from Patch.

Definition at line 3132 of file HomePatch.C.

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

Referenced by gbisComputeAfterP2().

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

FullAtomList& HomePatch::getAtomList (  )  [inline]

Definition at line 274 of file HomePatch.h.

Referenced by ComputeHomePatch::doWork().

00274 { return (atom); }

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

Definition at line 1915 of file HomePatch.C.

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

Referenced by Sequencer::hardWallDrude().

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

void HomePatch::loweAndersenFinish (  ) 

Definition at line 3006 of file HomePatch.C.

References DebugM.

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

void HomePatch::loweAndersenVelocities (  ) 

Definition at line 2991 of file HomePatch.C.

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

Referenced by positionsReady().

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

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

Definition at line 2856 of file HomePatch.C.

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

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

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

void HomePatch::mollyAverage (  ) 

Definition at line 3206 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::mollyMollify ( Tensor virial  ) 

Definition at line 3280 of file HomePatch.C.

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

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

void HomePatch::positionsReady ( int  doMigration = 0  ) 

Reimplemented from Patch.

Definition at line 1017 of file HomePatch.C.

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

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

void HomePatch::qmSwapAtoms (  ) 

Definition at line 958 of file HomePatch.C.

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

Referenced by positionsReady().

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

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

Definition at line 2273 of file HomePatch.C.

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

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

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

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

Definition at line 2452 of file HomePatch.C.

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

Referenced by rattle1().

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

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

Definition at line 2722 of file HomePatch.C.

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

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

void HomePatch::receiveResult ( ProxyGBISP2ResultMsg msg  ) 

Definition at line 3180 of file HomePatch.C.

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

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

void HomePatch::receiveResult ( ProxyGBISP1ResultMsg msg  ) 

Definition at line 3155 of file HomePatch.C.

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

Referenced by ProxyMgr::recvResult().

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

void HomePatch::receiveResults ( ProxyCombinedResultRawMsg msg  ) 

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

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

void HomePatch::receiveResults ( ProxyResultMsg msg  ) 

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

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

void HomePatch::receiveResults ( ProxyResultVarsizeMsg msg  ) 

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

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

void HomePatch::recvCheckpointAck (  ) 

Definition at line 3462 of file HomePatch.C.

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

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

void HomePatch::recvCheckpointLoad ( CheckpointAtomsMsg msg  ) 

Definition at line 3409 of file HomePatch.C.

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

Referenced by PatchMgr::recvCheckpointLoad().

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

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

Definition at line 3379 of file HomePatch.C.

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

Referenced by PatchMgr::recvCheckpointReq().

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

void HomePatch::recvCheckpointStore ( CheckpointAtomsMsg msg  ) 

Definition at line 3448 of file HomePatch.C.

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

Referenced by PatchMgr::recvCheckpointStore().

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

void HomePatch::recvExchangeMsg ( ExchangeAtomsMsg msg  ) 

Definition at line 3499 of file HomePatch.C.

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

Referenced by PatchMgr::recvExchangeMsg().

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

void HomePatch::recvExchangeReq ( int  req  ) 

Definition at line 3488 of file HomePatch.C.

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

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

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

void HomePatch::recvNodeAwareSpanningTree ( ProxyNodeAwareSpanningTreeMsg msg  ) 

Definition at line 722 of file HomePatch.C.

Referenced by ProxyMgr::recvNodeAwareSpanningTreeOnHomePatch().

00722 {}

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

Definition at line 728 of file HomePatch.C.

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

Referenced by ProxyMgr::recvSpanningTreeOnHomePatch().

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

void HomePatch::registerProxy ( RegisterProxyMsg  ) 

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

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

void HomePatch::replaceForces ( ExtForce f  ) 

Definition at line 1368 of file HomePatch.C.

References Patch::f.

01369 {
01370   replacementForces = f;
01371 }

void HomePatch::revert ( void   ) 

Definition at line 3353 of file HomePatch.C.

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

Referenced by Sequencer::algorithm().

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

void HomePatch::runSequencer ( void   ) 

Definition at line 361 of file HomePatch.C.

References Sequencer::run().

00362 { sequencer->run(); }

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

Definition at line 1373 of file HomePatch.C.

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

Referenced by Sequencer::saveForce().

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

void HomePatch::sendNodeAwareSpanningTree (  ) 

Definition at line 723 of file HomePatch.C.

00723 {}

void HomePatch::sendProxies (  ) 

Definition at line 560 of file HomePatch.C.

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

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

void HomePatch::sendSpanningTree (  ) 

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

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

void HomePatch::setGBISIntrinsicRadii (  ) 

Definition at line 3025 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::setLcpoType (  ) 

Definition at line 3014 of file HomePatch.C.

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

Referenced by positionsReady().

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

void HomePatch::submitLoadStats ( int  timestep  ) 

Definition at line 3513 of file HomePatch.C.

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

Referenced by Sequencer::rebalanceLoad().

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

void HomePatch::unregisterProxy ( UnregisterProxyMsg  ) 

Definition at line 508 of file HomePatch.C.

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

Referenced by ProxyMgr::recvUnregisterProxy().

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

void HomePatch::useSequencer ( Sequencer sequencerPtr  ) 

Definition at line 357 of file HomePatch.C.

00358 { 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 247 of file HomePatch.h.

Referenced by exchangeCheckpoint(), and recvCheckpointLoad().

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

Definition at line 254 of file HomePatch.h.

Referenced by recvCheckpointReq(), and recvCheckpointStore().

int HomePatch::exchange_dst

Definition at line 260 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

ExchangeAtomsMsg* HomePatch::exchange_msg

Definition at line 263 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

int HomePatch::exchange_req

Definition at line 262 of file HomePatch.h.

Referenced by exchangeAtoms(), and recvExchangeReq().

int HomePatch::exchange_src

Definition at line 261 of file HomePatch.h.

Referenced by exchangeAtoms().

int HomePatch::inMigration [protected]

Definition at line 307 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

LDObjHandle HomePatch::ldObjHandle

Definition at line 298 of file HomePatch.h.

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

int HomePatch::marginViolations

Definition at line 152 of file HomePatch.h.

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

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

Definition at line 309 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<int> HomePatch::noconstList

Definition at line 202 of file HomePatch.h.

Referenced by rattle1().

int HomePatch::numMlBuf [protected]

Definition at line 308 of file HomePatch.h.

Referenced by depositMigration(), and doAtomMigration().

std::vector<Vector> HomePatch::posNew

Definition at line 208 of file HomePatch.h.

Referenced by buildRattleList(), and rattle1().

std::vector<RattleList> HomePatch::rattleList

Definition at line 200 of file HomePatch.h.

Referenced by rattle1().

bool HomePatch::rattleListValid

Definition at line 204 of file HomePatch.h.

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

std::vector<RattleParam> HomePatch::rattleParam

Definition at line 201 of file HomePatch.h.

Referenced by rattle1().

std::vector<int> HomePatch::settleList

Definition at line 199 of file HomePatch.h.

Referenced by rattle1().

std::vector<Vector> HomePatch::velNew

Definition at line 207 of file HomePatch.h.

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


The documentation for this class was generated from the following files:
Generated on Wed Nov 22 01:17:21 2017 for NAMD by  doxygen 1.4.7