#include <HomePatch.h>
Inheritance diagram for HomePatch:
|
|
Definition at line 334 of file HomePatch.C. 00335 {
00336 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00337 ptnTree.resize(0);
00338 delete [] children;
00339 #ifdef USE_NODEPATCHMGR
00340 delete [] nodeChildren;
00341 #endif
00342 #else
00343 delete [] child;
00344 #endif
00345 }
|
|
||||||||||||||||
|
Definition at line 1237 of file HomePatch.C. References ResizeArray< Elem >::begin(), BigReal, ResizeArray< Elem >::const_begin(), SimParameters::fixedAtomsOn, Force, ForceList, namd_reciprocal, Node::Object(), Node::simParameters, simParams, FullAtom::velocity, Vector::x, Vector::y, and Vector::z. Referenced by Sequencer::addForceToMomentum(). 01239 {
01240 SimParameters *simParams = Node::Object()->simParameters;
01241 const BigReal dt = timestep / TIMEFACTOR;
01242 ForceList *f_use = (useSaved ? f_saved : f);
01243
01244 if ( simParams->fixedAtomsOn ) {
01245 for ( int i = 0; i < numAtoms; ++i ) {
01246 if ( atom[i].atomFixed ) {
01247 atom[i].velocity = 0;
01248 } else {
01249 BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.);
01250 atom[i].velocity += f_use[ftag][i] * recip_val;
01251 }
01252 }
01253 } else {
01254 FullAtom *atom_arr = atom.begin();
01255 const Force *force_arr = f_use[ftag].const_begin();
01256 #ifdef ARCH_POWERPC
01257 #pragma disjoint (*force_arr, *atom_arr)
01258 #endif
01259 for ( int i = 0; i < numAtoms; ++i ) {
01260 if (atom[i].mass == 0.) continue;
01261 BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.);
01262 //printf("Taking reciprocal of mass %f\n", atom[i].mass);
01263 atom_arr[i].velocity.x += force_arr[i].x * recip_val;
01264 atom_arr[i].velocity.y += force_arr[i].y * recip_val;
01265 atom_arr[i].velocity.z += force_arr[i].z * recip_val;
01266 }
01267 }
01268 }
|
|
|
Definition at line 1270 of file HomePatch.C. References ResizeArray< Elem >::begin(), BigReal, SimParameters::fixedAtomsOn, Node::Object(), CompAtom::position, Node::simParameters, simParams, FullAtom::velocity, Vector::x, Vector::y, and Vector::z. Referenced by Sequencer::addVelocityToPosition(). 01271 {
01272 SimParameters *simParams = Node::Object()->simParameters;
01273 const BigReal dt = timestep / TIMEFACTOR;
01274 if ( simParams->fixedAtomsOn ) {
01275 for ( int i = 0; i < numAtoms; ++i ) {
01276 if ( ! atom[i].atomFixed ) atom[i].position += atom[i].velocity * dt;
01277 }
01278 } else {
01279 FullAtom *atom_arr = atom.begin();
01280 for ( int i = 0; i < numAtoms; ++i ) {
01281 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
01282 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
01283 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
01284 }
01285 }
01286 }
|
|
|
Implements Patch. Definition at line 348 of file HomePatch.C. References Sequencer::awaken(), DebugM, ExtForce::force, j, ExtForce::replace, and Sequencer::thread. 00349 {
00350 if ( ! --boxesOpen )
00351 {
00352 if ( replacementForces ) {
00353 for ( int i = 0; i < numAtoms; ++i ) {
00354 if ( replacementForces[i].replace ) {
00355 for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
00356 f[Results::normal][i] = replacementForces[i].force;
00357 }
00358 }
00359 replacementForces = 0;
00360 }
00361 DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
00362 << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
00363 // only awaken suspended threads. Then say it is suspended.
00364 sequencer->awaken();
00365 return;
00366 }
00367 else
00368 {
00369 DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
00370 }
00371 }
|
|
|
Definition at line 645 of file HomePatch.C. References ResizeArrayIter< T >::begin(), compDistance(), ResizeArrayIter< T >::end(), ResizeArray< Elem >::find(), j, NAMD_bug(), NodeIDList, PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), PatchMap::Object(), ProxyListIter, proxySpanDim, ResizeArray< Elem >::resize(), sendSpanningTree(), ResizeArray< Elem >::setall(), and ResizeArray< Elem >::size(). Referenced by ProxyMgr::buildProxySpanningTree(). 00646 {
00647 nChild = 0;
00648 int psize = proxy.size();
00649 if (psize == 0) return;
00650 NodeIDList oldtree = tree;
00651 int oldsize = oldtree.size();
00652 tree.resize(psize + 1);
00653 tree.setall(-1);
00654 tree[0] = CkMyPe();
00655 int s=1, e=psize+1;
00656 ProxyListIter pli(proxy);
00657 int patchNodesLast =
00658 ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00659 int nNonPatch = 0;
00660 #if 0
00661 int* pelists = new int[psize];
00662 for (int i=0; i<psize; i++) pelists[i] = -1;
00663 for ( pli = pli.begin(); pli != pli.end(); ++pli ) {
00664 int idx = rand()%psize;
00665 while (pelists[idx] != -1) { idx++; if (idx == psize) idx = 0; }
00666 pelists[idx] = pli->node;
00667 }
00668 for ( int i=0; i<psize; i++)
00669 {
00670 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pelists[i]) ) {
00671 tree[--e] = pelists[i];
00672 } else {
00673 tree[s++] = pelists[i];
00674 nNonPatch++;
00675 }
00676 }
00677 delete [] pelists;
00678 #else
00679 // try to put it to the same old tree
00680 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00681 {
00682 int oldindex = oldtree.find(pli->node);
00683 if (oldindex != -1 && oldindex < psize) {
00684 tree[oldindex] = pli->node;
00685 }
00686 }
00687 s=1; e=psize;
00688 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00689 {
00690 if (tree.find(pli->node) != -1) continue; // already assigned
00691 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
00692 while (tree[e] != -1) { e--; if (e==-1) e = psize; }
00693 tree[e] = pli->node;
00694 } else {
00695 while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
00696 tree[s] = pli->node;
00697 nNonPatch++;
00698 }
00699 }
00700 #if 1
00701 if (oldsize==0 && nNonPatch) {
00702 // first time, sort by distance
00703 qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
00704 }
00705 #endif
00706
00707 //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
00708
00709 #if USE_TOPOMAP && USE_SPANNING_TREE
00710
00711 //Right now only works for spanning trees with two levels
00712 int *treecopy = new int [psize];
00713 int subroots[proxySpanDim];
00714 int subsizes[proxySpanDim];
00715 int subtrees[proxySpanDim][proxySpanDim];
00716 int idxes[proxySpanDim];
00717 int i = 0;
00718
00719 for(i=0;i<proxySpanDim;i++){
00720 subsizes[i] = 0;
00721 idxes[i] = i;
00722 }
00723
00724 for(i=0;i<psize;i++){
00725 treecopy[i] = tree[i+1];
00726 }
00727
00728 TopoManager tmgr;
00729 tmgr.sortRanksByHops(treecopy,nNonPatch);
00730 tmgr.sortRanksByHops(treecopy+nNonPatch,
00731 psize-nNonPatch);
00732
00733 /* build tree and subtrees */
00734 nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
00735 delete [] treecopy;
00736
00737 for(int i=1;i<psize+1;i++){
00738 int isSubroot=0;
00739 for(int j=0;j<nChild;j++)
00740 if(tree[i]==subroots[j]){
00741 isSubroot=1;
00742 break;
00743 }
00744 if(isSubroot) continue;
00745
00746 int bAdded = 0;
00747 tmgr.sortIndexByHops(tree[i], subroots,
00748 idxes, proxySpanDim);
00749 for(int j=0;j<proxySpanDim;j++){
00750 if(subsizes[idxes[j]]<proxySpanDim){
00751 subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
00752 bAdded = 1;
00753 break;
00754 }
00755 }
00756 if( psize > proxySpanDim && ! bAdded ) {
00757 NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
00758 }
00759 }
00760
00761 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
00762
00763 for (int i=1; i<=proxySpanDim; i++) {
00764 if (tree.size() <= i) break;
00765 child[i-1] = tree[i];
00766 nChild++;
00767 }
00768 #endif
00769 #endif
00770
00771 #if 0
00772 // for debugging
00773 CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
00774 for (int i=0; i<psize+1; i++) {
00775 CkPrintf("%d ", tree[i]);
00776 }
00777 CkPrintf("\n");
00778 #endif
00779 // send to children nodes
00780 sendSpanningTree();
00781 }
|
|
|
Definition at line 1799 of file HomePatch.C. References FullAtomList. Referenced by Sequencer::algorithm(). 01799 {
01800 FullAtomList tmp_a(&atom); checkpoint_atom = tmp_a;
01801 checkpoint_lattice = lattice;
01802
01803 // DMK - Atom Separation (water vs. non-water)
01804 #if NAMD_SeparateWaters != 0
01805 checkpoint_numWaterAtoms = numWaterAtoms;
01806 #endif
01807 }
|
|
|
Definition at line 2147 of file HomePatch.C. References ResizeArray< Elem >::add(), Lattice::apply_transform(), Sequencer::awaken(), ResizeArray< Elem >::begin(), DebugM, ResizeArray< Elem >::end(), MigrationList, MigrateAtomsMsg::migrationList, msgbuf, Lattice::nearest(), numMlBuf, Lattice::reverse_transform(), and ResizeArray< Elem >::size(). Referenced by MigrateAtomsCombinedMsg::distribute(), doAtomMigration(), and PatchMgr::recvMigrateAtoms(). 02148 {
02149
02150 if (!inMigration) { // We have to buffer changes due to migration
02151 // until our patch is in migration mode
02152 msgbuf[numMlBuf++] = msg;
02153 return;
02154 }
02155
02156
02157 // DMK - Atom Separation (water vs. non-water)
02158 #if NAMD_SeparateWaters != 0
02159
02160
02161 // Merge the incoming list of atoms with the current list of
02162 // atoms. Note that mergeSeparatedAtomList() will apply any
02163 // required transformations to the incoming atoms as it is
02164 // separating them.
02165 mergeAtomList(msg->migrationList);
02166
02167
02168 #else
02169
02170 // Merge the incoming list of atoms with the current list of
02171 // atoms. Apply transformations to the incoming atoms as they are
02172 // added to this patch's list.
02173 {
02174 MigrationList &migrationList = msg->migrationList;
02175 MigrationList::iterator mi;
02176 Transform mother_transform;
02177 for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
02178 DebugM(1,"Migrating atom " << mi->id << " to patch "
02179 << patchID << " with position " << mi->position << "\n");
02180 if ( mi->hydrogenGroupSize ) {
02181 mi->position = lattice.nearest(mi->position,center,&(mi->transform));
02182 mother_transform = mi->transform;
02183 } else {
02184 mi->position = lattice.reverse_transform(mi->position,mi->transform);
02185 mi->position = lattice.apply_transform(mi->position,mother_transform);
02186 mi->transform = mother_transform;
02187 }
02188 atom.add(*mi);
02189 }
02190 }
02191
02192
02193 #endif // if (NAMD_SeparateWaters != 0)
02194
02195
02196 numAtoms = atom.size();
02197 delete msg;
02198
02199 DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
02200 if (!--patchMigrationCounter) {
02201 DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
02202 allMigrationIn = true;
02203 patchMigrationCounter = numNeighbors;
02204 if (migrationSuspended) {
02205 DebugM(3,"patch "<<patchID<<" is being awakened\n");
02206 migrationSuspended = false;
02207 sequencer->awaken();
02208 return;
02209 }
02210 else {
02211 DebugM(3,"patch "<<patchID<<" was not suspended\n");
02212 }
02213 }
02214 }
|
|
|
Definition at line 2022 of file HomePatch.C. References ResizeArray< Elem >::add(), ResizeArray< Elem >::begin(), BigReal, DebugM, ResizeArray< Elem >::del(), depositMigration(), ResizeArray< Elem >::end(), inMigration, marginViolations, MigrationList, MigrationInfo::mList, msgbuf, PatchMgr::Object(), AtomMap::Object(), ResizeArray< Elem >::resize(), Lattice::scale(), ScaledPosition, PatchMgr::sendMigrationMsgs(), ResizeArray< Elem >::size(), Sequencer::suspend(), AtomMap::unregisterIDs(), Vector::x, Vector::y, and Vector::z. Referenced by positionsReady(). 02023 {
02024 int i;
02025
02026 for (i=0; i<numNeighbors; i++) {
02027 realInfo[i].mList.resize(0);
02028 }
02029
02030 // Purge the AtomMap
02031 AtomMap::Object()->unregisterIDs(patchID,pExt.begin(),pExt.end());
02032
02033 // Determine atoms that need to migrate
02034
02035 BigReal minx = min.x;
02036 BigReal miny = min.y;
02037 BigReal minz = min.z;
02038 BigReal maxx = max.x;
02039 BigReal maxy = max.y;
02040 BigReal maxz = max.z;
02041
02042 int xdev, ydev, zdev;
02043 int delnum = 0;
02044
02045 FullAtomList::iterator atom_i = atom.begin();
02046 FullAtomList::iterator atom_e = atom.end();
02047
02048 // DMK - Atom Separation (water vs. non-water)
02049 #if NAMD_SeparateWaters != 0
02050 FullAtomList::iterator atom_first = atom_i;
02051 int numLostWaterAtoms = 0;
02052 #endif
02053
02054 while ( atom_i != atom_e ) {
02055 if ( atom_i->hydrogenGroupSize ) {
02056
02057 ScaledPosition s = lattice.scale(atom_i->position);
02058
02059 // check if atom is within bounds
02060 if (s.x < minx) xdev = 0;
02061 else if (maxx <= s.x) xdev = 2;
02062 else xdev = 1;
02063
02064 if (s.y < miny) ydev = 0;
02065 else if (maxy <= s.y) ydev = 2;
02066 else ydev = 1;
02067
02068 if (s.z < minz) zdev = 0;
02069 else if (maxz <= s.z) zdev = 2;
02070 else zdev = 1;
02071
02072 }
02073
02074 if (mInfo[xdev][ydev][zdev]) { // process atom for migration
02075 // Don't migrate if destination is myself
02076
02077 // See if we have a migration list already
02078 MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
02079 DebugM(3,"Migrating atom " << atom_i->id << " from patch "
02080 << patchID << " with position " << atom_i->position << "\n");
02081 mCur.add(*atom_i);
02082
02083 ++delnum;
02084
02085
02086 // DMK - Atom Separation (water vs. non-water)
02087 #if NAMD_SeparateWaters != 0
02088 // Check to see if this atom is part of a water molecule. If
02089 // so, numWaterAtoms needs to be adjusted to refect the lost of
02090 // this atom.
02091 // NOTE: The atom separation code assumes that if the oxygen
02092 // atom of the hydrogen group making up the water molecule is
02093 // migrated to another HomePatch, the hydrogens will also
02094 // move!!!
02095 int atomIndex = atom_i - atom_first;
02096 if (atomIndex < numWaterAtoms)
02097 numLostWaterAtoms++;
02098 #endif
02099
02100
02101 } else {
02102
02103 if ( delnum ) { *(atom_i-delnum) = *atom_i; }
02104
02105 }
02106
02107 ++atom_i;
02108 }
02109
02110 // DMK - Atom Separation (water vs. non-water)
02111 #if NAMD_SeparateWaters != 0
02112 numWaterAtoms -= numLostWaterAtoms;
02113 #endif
02114
02115
02116 int delpos = numAtoms - delnum;
02117 DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
02118 atom.del(delpos,delnum);
02119
02120 numAtoms = atom.size();
02121
02122 PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
02123
02124 // signal depositMigration() that we are inMigration mode
02125 inMigration = true;
02126
02127 // Drain the migration message buffer
02128 for (i=0; i<numMlBuf; i++) {
02129 DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
02130 depositMigration(msgbuf[i]);
02131 }
02132 numMlBuf = 0;
02133
02134 if (!allMigrationIn) {
02135 DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
02136 migrationSuspended = true;
02137 sequencer->suspend();
02138 migrationSuspended = false;
02139 }
02140 allMigrationIn = false;
02141
02142 inMigration = false;
02143 marginViolations = 0;
02144 }
|
|
|
Definition at line 1901 of file HomePatch.C. References ResizeArray< Elem >::begin(), BigReal, Flags::doNonbonded, ResizeArray< Elem >::end(), SimParameters::hgroupCutoff, Flags::maxGroupRadius, Node::Object(), Node::simParameters, and simParams. Referenced by positionsReady(). 01902 {
01903 if ( ! flags.doNonbonded ) return;
01904
01905 SimParameters *simParams = Node::Object()->simParameters;
01906 BigReal hgcut = 0.5 * simParams->hgroupCutoff; hgcut *= hgcut;
01907 BigReal maxrad2 = 0.;
01908
01909 FullAtomList::iterator p_i = atom.begin();
01910 FullAtomList::iterator p_e = atom.end();
01911
01912 while ( p_i != p_e ) {
01913 int hgs = p_i->hydrogenGroupSize;
01914 p_i->nonbondedGroupSize = hgs;
01915 BigReal x = p_i->position.x;
01916 BigReal y = p_i->position.y;
01917 BigReal z = p_i->position.z;
01918 ++p_i;
01919 int oversize = 0;
01920 // limit spatial extent
01921 for ( int i = 1; i < hgs; ++i ) {
01922 p_i->nonbondedGroupSize = 0;
01923 BigReal dx = p_i->position.x - x;
01924 BigReal dy = p_i->position.y - y;
01925 BigReal dz = p_i->position.z - z;
01926 BigReal r2 = dx * dx + dy * dy + dz * dz;
01927 ++p_i;
01928 if ( r2 > hgcut ) oversize = 1;
01929 else if ( r2 > maxrad2 ) maxrad2 = r2;
01930 }
01931 // also limit to at most 4 atoms per group
01932 if ( oversize || hgs > 4 ) {
01933 p_i -= hgs;
01934 for ( int i = 0; i < hgs; ++i ) {
01935 p_i->nonbondedGroupSize = 1;
01936 ++p_i;
01937 }
01938 }
01939 }
01940
01941 flags.maxGroupRadius = sqrt(maxrad2);
01942
01943 }
|
|
|
Definition at line 1945 of file HomePatch.C. References Lattice::a(), Lattice::a_r(), Lattice::b(), Lattice::b_r(), ResizeArray< Elem >::begin(), BigReal, Lattice::c(), Lattice::c_r(), SimParameters::cutoff, ResizeArray< Elem >::end(), SimParameters::margin, marginViolations, NAMD_die(), Node::Object(), SimParameters::patchDimension, Lattice::scale(), ScaledPosition, Node::simParameters, simParams, Vector::unit(), Vector::x, Vector::y, and Vector::z. Referenced by positionsReady(). 01946 {
01947 SimParameters *simParams = Node::Object()->simParameters;
01948
01949 BigReal sysdima = lattice.a_r().unit() * lattice.a();
01950 BigReal sysdimb = lattice.b_r().unit() * lattice.b();
01951 BigReal sysdimc = lattice.c_r().unit() * lattice.c();
01952
01953 BigReal minSize = simParams->patchDimension - simParams->margin;
01954
01955 if ( ( (max.x - min.x)*aAway*sysdima < minSize*0.9999 ) ||
01956 ( (max.y - min.y)*bAway*sysdimb < minSize*0.9999 ) ||
01957 ( (max.z - min.z)*cAway*sysdimc < minSize*0.9999 ) ) {
01958
01959 NAMD_die("Periodic cell has become too small for original patch grid!\n"
01960 "Possible solutions are to restart from a recent checkpoint,\n"
01961 "increase margin, or disable useFlexibleCell for liquid simulation.");
01962 }
01963
01964 BigReal cutoff = simParams->cutoff;
01965
01966 BigReal margina = 0.5 * ( (max.x - min.x) * aAway - cutoff / sysdima );
01967 BigReal marginb = 0.5 * ( (max.y - min.y) * bAway - cutoff / sysdimb );
01968 BigReal marginc = 0.5 * ( (max.z - min.z) * cAway - cutoff / sysdimc );
01969
01970 if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
01971 NAMD_die("Periodic cell has become too small for original patch grid!\n"
01972 "There are probably many margin violations already reported.\n"
01973 "Possible solutions are to restart from a recent checkpoint,\n"
01974 "increase margin, or disable useFlexibleCell for liquid simulation.");
01975 }
01976
01977 BigReal minx = min.x - margina;
01978 BigReal miny = min.y - marginb;
01979 BigReal minz = min.z - marginc;
01980 BigReal maxx = max.x + margina;
01981 BigReal maxy = max.y + marginb;
01982 BigReal maxz = max.z + marginc;
01983
01984 int xdev, ydev, zdev;
01985 int problemCount = 0;
01986
01987 FullAtomList::iterator p_i = atom.begin();
01988 FullAtomList::iterator p_e = atom.end();
01989 for ( ; p_i != p_e; ++p_i ) {
01990
01991 ScaledPosition s = lattice.scale(p_i->position);
01992
01993 // check if atom is within bounds
01994 if (s.x < minx) xdev = 0;
01995 else if (maxx <= s.x) xdev = 2;
01996 else xdev = 1;
01997
01998 if (s.y < miny) ydev = 0;
01999 else if (maxy <= s.y) ydev = 2;
02000 else ydev = 1;
02001
02002 if (s.z < minz) zdev = 0;
02003 else if (maxz <= s.z) zdev = 2;
02004 else zdev = 1;
02005
02006 if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
02007 ++problemCount;
02008 }
02009
02010 }
02011
02012 marginViolations = problemCount;
02013 // if ( problemCount ) {
02014 // iout << iERROR <<
02015 // "Found " << problemCount << " margin violations!\n" << endi;
02016 // }
02017
02018 }
|
|
|
Definition at line 1826 of file HomePatch.C. References ResizeArray< Elem >::begin(), BigReal, j, Vector::length2(), Flags::maxAtomMovement, Node::Object(), SimParameters::pairlistGrow, SimParameters::pairlistShrink, Flags::pairlistTolerance, SimParameters::pairlistTrigger, CompAtom::position, Position, ResizeArray< Elem >::resize(), Flags::savePairlists, ScaledPosition, Node::simParameters, simParams, Lattice::unscale(), Flags::usePairlists, Vector::x, Vector::y, and Vector::z. Referenced by positionsReady(). 01827 {
01828 SimParameters *simParams = Node::Object()->simParameters;
01829
01830 if ( ! flags.usePairlists ) {
01831 flags.pairlistTolerance = 0.;
01832 flags.maxAtomMovement = 99999.;
01833 return;
01834 }
01835
01836 int i; int n = numAtoms;
01837 CompAtom *p_i = p.begin();
01838
01839 if ( flags.savePairlists ) {
01840 flags.pairlistTolerance = doPairlistCheck_newTolerance;
01841 flags.maxAtomMovement = 0.;
01842 doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
01843 doPairlistCheck_lattice = lattice;
01844 doPairlistCheck_positions.resize(numAtoms);
01845 CompAtom *psave_i = doPairlistCheck_positions.begin();
01846 for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
01847 return;
01848 }
01849
01850 Lattice &lattice_old = doPairlistCheck_lattice;
01851 Position center_cur = lattice.unscale(center);
01852 Position center_old = lattice_old.unscale(center);
01853 Vector center_delta = center_cur - center_old;
01854
01855 // find max deviation to corner (any neighbor shares a corner)
01856 BigReal max_cd = 0.;
01857 for ( i=0; i<2; ++i ) {
01858 for ( int j=0; j<2; ++j ) {
01859 for ( int k=0; k<2; ++k ) {
01860 ScaledPosition corner( i ? min.x : max.x ,
01861 j ? min.y : max.y ,
01862 k ? min.z : max.z );
01863 Vector corner_delta =
01864 lattice.unscale(corner) - lattice_old.unscale(corner);
01865 corner_delta -= center_delta;
01866 BigReal cd = corner_delta.length2();
01867 if ( cd > max_cd ) max_cd = cd;
01868 }
01869 }
01870 }
01871 max_cd = sqrt(max_cd);
01872
01873 // find max deviation of atoms relative to center
01874 BigReal max_pd = 0.;
01875 CompAtom *p_old_i = doPairlistCheck_positions.begin();
01876 for ( i=0; i<n; ++i ) {
01877 Vector p_delta = p_i[i].position - p_old_i[i].position;
01878 p_delta -= center_delta;
01879 BigReal pd = p_delta.length2();
01880 if ( pd > max_pd ) max_pd = pd;
01881 }
01882 max_pd = sqrt(max_pd);
01883
01884 BigReal max_tol = max_pd + max_cd;
01885
01886 flags.maxAtomMovement = max_tol;
01887
01888 // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
01889
01890 if ( max_tol > ( (1. - simParams->pairlistTrigger) *
01891 doPairlistCheck_newTolerance ) ) {
01892 doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
01893 }
01894
01895 if ( max_tol > doPairlistCheck_newTolerance ) {
01896 doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
01897 }
01898
01899 }
|
|
|
Definition at line 118 of file HomePatch.h. References FullAtomList. Referenced by ComputeHomePatch::doWork(), ComputeGridForce::doWork(), dumpbench(), and ComputeGridForce::finishWork(). 00118 { return (atom); }
|
|
|
Definition at line 1662 of file HomePatch.C. References average(), BigReal, SimParameters::fixedAtomsOn, HGArrayBigReal, HGArrayInt, HGArrayVector, iout, iWARN(), Node::molecule, SimParameters::mollyIter, SimParameters::mollyTol, NAMD_die(), Node::Object(), ResizeArray< Elem >::resize(), Molecule::rigid_bond_length(), Node::simParameters, and simParams. Referenced by positionsReady(). 01663 {
01664 Molecule *mol = Node::Object()->molecule;
01665 SimParameters *simParams = Node::Object()->simParameters;
01666 BigReal tol = simParams->mollyTol;
01667 int maxiter = simParams->mollyIter;
01668 int i, iter;
01669 HGArrayBigReal dsq;
01670 BigReal tmp;
01671 HGArrayInt ial, ibl;
01672 HGArrayVector ref; // reference position
01673 HGArrayVector refab; // reference vector
01674 HGArrayBigReal rmass; // 1 / mass
01675 BigReal *lambda; // Lagrange multipliers
01676 CompAtom *avg; // averaged position
01677 int numLambdas = 0;
01678 HGArrayInt fixed; // is atom fixed?
01679
01680 // iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
01681 p_avg.resize(numAtoms);
01682 for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
01683
01684 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01685 int hgs = atom[ig].hydrogenGroupSize;
01686 if ( hgs == 1 ) continue; // only one atom in group
01687 for ( i = 0; i < hgs; ++i ) {
01688 ref[i] = atom[ig+i].position;
01689 rmass[i] = 1. / atom[ig+i].mass;
01690 fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
01691 if ( fixed[i] ) rmass[i] = 0.;
01692 }
01693 avg = &(p_avg[ig]);
01694 int icnt = 0;
01695
01696 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) ) { // for water
01697 if ( hgs != 3 ) {
01698 NAMD_die("Hydrogen group error caught in mollyAverage(). It's a bug!\n");
01699 }
01700 if ( !(fixed[1] && fixed[2]) ) {
01701 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01702 }
01703 }
01704 for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
01705 if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) ) {
01706 if ( !(fixed[0] && fixed[i]) ) {
01707 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01708 }
01709 }
01710 }
01711 if ( icnt == 0 ) continue; // no constraints
01712 numLambdas += icnt;
01713 molly_lambda.resize(numLambdas);
01714 lambda = &(molly_lambda[numLambdas - icnt]);
01715 for ( i = 0; i < icnt; ++i ) {
01716 refab[i] = ref[ial[i]] - ref[ibl[i]];
01717 }
01718 // iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
01719 iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
01720 if ( iter == maxiter ) {
01721 iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
01722 }
01723 }
01724
01725 // for ( i=0; i<numAtoms; ++i ) {
01726 // if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
01727 // iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
01728 // << p[i].position << " to " << p_avg[i].position << "\n" << endi;
01729 // }
01730 // }
01731
01732 }
|
|
|
Definition at line 1736 of file HomePatch.C. References BigReal, SimParameters::fixedAtomsOn, HGArrayBigReal, HGArrayInt, HGArrayVector, Node::molecule, mollify(), NAMD_die(), Node::Object(), outer(), ResizeArray< Elem >::resize(), Molecule::rigid_bond_length(), Node::simParameters, and simParams. Referenced by Sequencer::runComputeObjects(). 01737 {
01738 Molecule *mol = Node::Object()->molecule;
01739 SimParameters *simParams = Node::Object()->simParameters;
01740 Tensor wc; // constraint virial
01741 int i;
01742 HGArrayInt ial, ibl;
01743 HGArrayVector ref; // reference position
01744 CompAtom *avg; // averaged position
01745 HGArrayVector refab; // reference vector
01746 HGArrayVector force; // new force
01747 HGArrayBigReal rmass; // 1 / mass
01748 BigReal *lambda; // Lagrange multipliers
01749 int numLambdas = 0;
01750 HGArrayInt fixed; // is atom fixed?
01751
01752 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01753 int hgs = atom[ig].hydrogenGroupSize;
01754 if (hgs == 1 ) continue; // only one atom in group
01755 for ( i = 0; i < hgs; ++i ) {
01756 ref[i] = atom[ig+i].position;
01757 force[i] = f[Results::slow][ig+i];
01758 rmass[i] = 1. / atom[ig+i].mass;
01759 fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
01760 if ( fixed[i] ) rmass[i] = 0.;
01761 }
01762 int icnt = 0;
01763 // c-ji I'm only going to mollify water for now
01764 if ( ( mol->rigid_bond_length(atom[ig].id) ) ) { // for water
01765 if ( hgs != 3 ) {
01766 NAMD_die("Hydrogen group error caught in mollyMollify(). It's a bug!\n");
01767 }
01768 if ( !(fixed[1] && fixed[2]) ) {
01769 ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01770 }
01771 }
01772 for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
01773 if ( ( mol->rigid_bond_length(atom[ig+i].id) ) ) {
01774 if ( !(fixed[0] && fixed[i]) ) {
01775 ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01776 }
01777 }
01778 }
01779
01780 if ( icnt == 0 ) continue; // no constraints
01781 lambda = &(molly_lambda[numLambdas]);
01782 numLambdas += icnt;
01783 for ( i = 0; i < icnt; ++i ) {
01784 refab[i] = ref[ial[i]] - ref[ibl[i]];
01785 }
01786 avg = &(p_avg[ig]);
01787 mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
01788 // store data back to patch
01789 for ( i = 0; i < hgs; ++i ) {
01790 wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
01791 f[Results::slow][ig+i] = force[i];
01792 }
01793 }
01794 // check that there isn't a constant needed here!
01795 *virial += wc;
01796 p_avg.resize(0);
01797 }
|
|
|
Reimplemented from Patch. Definition at line 864 of file HomePatch.C. References ProxyDataMsg::avgPlLen, ProxyDataMsg::avgPositionList, ResizeArrayIter< T >::begin(), ResizeArray< Elem >::begin(), DebugM, doAtomMigration(), doGroupSizeCheck(), doMarginCheck(), Flags::doMolly, doPairlistCheck(), ResizeArray< Elem >::end(), ResizeArrayIter< T >::end(), ProxyDataMsg::flags, mollyAverage(), PatchMap::numNodesWithPatches(), PatchMap::numPatchesOnNode(), Sync::Object(), ProxyMgr::Object(), PatchMap::Object(), ProxyDataMsg::patch, PATCH_PRIORITY, Sync::PatchReady(), ProxyDataMsg::plExtLen, ProxyDataMsg::plLen, ProxyDataMsg::positionExtList, ProxyDataMsg::positionList, Patch::positionsReady(), PROXY_DATA_PRIORITY, ProxyListIter, proxySendSpanning, ResizeArray< Elem >::resize(), ProxyMgr::sendProxyAll(), ProxyMgr::sendProxyData(), Flags::sequence, SET_PRIORITY, and ResizeArray< Elem >::size(). Referenced by Sequencer::runComputeObjects(). 00865 {
00866 flags.sequence++;
00867
00868 if (!patchMapRead) { readPatchMap(); }
00869
00870 if (numNeighbors) {
00871 if (doMigration) {
00872 doAtomMigration();
00873 } else {
00874 doMarginCheck();
00875 }
00876 }
00877 doMigration = (doMigration && numNeighbors) || ! patchMapRead;
00878
00879 // Workaround for oversize groups
00880 doGroupSizeCheck();
00881
00882 // Copy information needed by computes and proxys to Patch::p.
00883 p.resize(numAtoms);
00884 CompAtom *p_i = p.begin();
00885 pExt.resize(numAtoms);
00886 CompAtomExt *pExt_i = pExt.begin();
00887 FullAtom *a_i = atom.begin();
00888 int i; int n = numAtoms;
00889 for ( i=0; i<n; ++i ) {
00890 p_i[i] = a_i[i];
00891 pExt_i[i] = a_i[i];
00892 }
00893
00894 // Measure atom movement to test pairlist validity
00895 doPairlistCheck();
00896
00897 if (flags.doMolly) mollyAverage();
00898
00899 // Must Add Proxy Changes when migration completed!
00900 ProxyListIter pli(proxy);
00901 int *pids = NULL;
00902 int npid;
00903 if (proxySendSpanning == 0) {
00904 npid = proxy.size();
00905 pids = new int[npid];
00906 int *pidi = pids;
00907 int *pide = pids + proxy.size();
00908 int patchNodesLast =
00909 ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00910 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00911 {
00912 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
00913 *(--pide) = pli->node;
00914 } else {
00915 *(pidi++) = pli->node;
00916 }
00917 }
00918 }
00919 else {
00920 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00921 #ifdef USE_NODEPATCHMGR
00922 npid = numNodeChild;
00923 if(numNodeChild>0) {
00924 pids = new int[npid];
00925 memcpy(pids, nodeChildren, sizeof(int)*npid);
00926 }
00927 #else
00928 npid = numChild;
00929 if(numChild>0) {
00930 pids = new int[numChild];
00931 memcpy(pids, children, sizeof(int)*numChild);
00932 }
00933 #endif
00934 #else
00935 npid = nChild;
00936 pids = new int[proxySpanDim];
00937 for (int i=0; i<nChild; i++) pids[i] = child[i];
00938 #endif
00939 }
00940 if (npid) {
00941 #if CMK_PERSISTENT_COMM
00942 if (phsReady == 0)
00943 {
00944 //CmiPrintf("Build on %d phs0:%d\n", CkMyPe(), localphs[0]);
00945 for (int i=0; i<npid; i++) {
00946 localphs[i] = CmiCreatePersistent(pids[i], 30000);
00947 }
00948 nphs = npid;
00949 phsReady = 1;
00950 }
00951 #endif
00952 int seq = flags.sequence;
00953 int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
00954 //begin to prepare proxy msg and send it
00955 int pdMsgPLLen = p.size();
00956 int pdMsgAvgPLLen = 0;
00957 if(flags.doMolly) {
00958 pdMsgAvgPLLen = p_avg.size();
00959 }
00960 int pdMsgPLExtLen = 0;
00961 if(doMigration || isNewProxyAdded) {
00962 pdMsgPLExtLen = pExt.size();
00963 }
00964 ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgPLExtLen, PRIORITY_SIZE) ProxyDataMsg;
00965 SET_PRIORITY(nmsg,seq,priority);
00966 nmsg->patch = patchID;
00967 nmsg->flags = flags;
00968 nmsg->plLen = pdMsgPLLen;
00969 //copying data to the newly created msg
00970 memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
00971 nmsg->avgPlLen = pdMsgAvgPLLen;
00972 if(flags.doMolly) {
00973 memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
00974 }
00975 nmsg->plExtLen = pdMsgPLExtLen;
00976 if(doMigration || isNewProxyAdded){
00977 memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
00978 }
00979
00980 #if NAMD_SeparateWaters != 0
00981 //DMK - Atom Separation (water vs. non-water)
00982 nmsg->numWaterAtoms = numWaterAtoms;
00983 #endif
00984
00985 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
00986 nmsg->isFromImmMsgCall = 0;
00987 #endif
00988
00989 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00990 DebugFileTrace *dft = DebugFileTrace::Object();
00991 dft->openTrace();
00992 dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
00993 for(int i=0; i<npid; i++) {
00994 dft->writeTrace("%d ", pids[i]);
00995 }
00996 dft->writeTrace("\n");
00997 dft->closeTrace();
00998 #endif
00999
01000 if(doMigration) {
01001 ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
01002 }else{
01003 ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
01004 }
01005 #if CMK_PERSISTENT_COMM
01006 CmiUsePersistentHandle(NULL, 0);
01007 #endif
01008 isNewProxyAdded = 0;
01009 }
01010 delete [] pids;
01011 DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
01012
01013 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
01014 positionPtrBegin = p.begin();
01015 positionPtrEnd = p.end();
01016 #endif
01017
01018 if(flags.doMolly) {
01019 avgPositionPtrBegin = p_avg.begin();
01020 avgPositionPtrEnd = p_avg.end();
01021 }
01022 Patch::positionsReady(doMigration);
01023
01024 patchMapRead = 1;
01025
01026 // gzheng
01027 if (useSync) Sync::Object()->PatchReady();
01028 }
|
|
||||||||||||||||
|
Definition at line 1289 of file HomePatch.C. References BigReal, Lattice::c(), SimParameters::fixedAtomsOn, Force, iERROR(), iout, SubmitReduction::item(), iWARN(), Node::molecule, NAMD_bug(), Node::Object(), Lattice::origin(), outer(), SimParameters::pressureProfileSlabs, Molecule::rigid_bond_length(), SimParameters::rigidDie, SimParameters::rigidIter, SimParameters::rigidTol, settle1(), settle1init(), settle1isinitted(), Node::simParameters, simParams, SimParameters::useSettle, SimParameters::watmodel, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz. Referenced by Sequencer::rattle1(). 01291 {
01292 Molecule *mol = Node::Object()->molecule;
01293 SimParameters *simParams = Node::Object()->simParameters;
01294 const int fixedAtomsOn = simParams->fixedAtomsOn;
01295 const int useSettle = simParams->useSettle;
01296 const BigReal dt = timestep / TIMEFACTOR;
01297 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
01298 BigReal tol2 = 2.0 * simParams->rigidTol;
01299 int maxiter = simParams->rigidIter;
01300 int dieOnError = simParams->rigidDie;
01301 int i, iter;
01302 BigReal dsq[10], tmp;
01303 int ial[10], ibl[10];
01304 Vector ref[10]; // reference position
01305 Vector refab[10]; // reference vector
01306 Vector pos[10]; // new position
01307 Vector vel[10]; // new velocity
01308 Vector netdp[10]; // total momentum change from constraint
01309 BigReal rmass[10]; // 1 / mass
01310 int fixed[10]; // is atom fixed?
01311 Tensor wc; // constraint virial
01312 BigReal idz, zmin;
01313 int nslabs;
01314
01315 // Initialize the settle algorithm with water parameters
01316 // settle1() assumes all waters are identical,
01317 // and will generate bad results if they are not.
01318 // XXX this will move to Molecule::build_atom_status when that
01319 // version is debugged
01320 if (!settle1isinitted()) {
01321 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01322 // find a water
01323 if (mol->rigid_bond_length(atom[ig].id) > 0) {
01324 int oatm;
01325 if (simParams->watmodel == WAT_SWM4) {
01326 oatm = ig+3; // skip over Drude and Lonepair
01327 //printf("ig=%d mass_ig=%g oatm=%d mass_oatm=%g\n",
01328 // ig, atom[ig].mass, oatm, atom[oatm].mass);
01329 }
01330 else {
01331 oatm = ig+1;
01332 // Avoid using the Om site to set this by mistake
01333 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
01334 oatm += 1;
01335 }
01336 }
01337
01338 // initialize settle water parameters
01339 settle1init(atom[ig].mass, atom[oatm].mass,
01340 mol->rigid_bond_length(atom[ig].id),
01341 mol->rigid_bond_length(atom[oatm].id));
01342 break; // done with init
01343 }
01344 }
01345 }
01346
01347 if (ppreduction) {
01348 nslabs = simParams->pressureProfileSlabs;
01349 idz = nslabs/lattice.c().z;
01350 zmin = lattice.origin().z - 0.5*lattice.c().z;
01351 }
01352
01353 // Size of a hydrogen group for water
01354 int wathgsize = 3;
01355 int watmodel = simParams->watmodel;
01356 if (watmodel == WAT_TIP4) wathgsize = 4;
01357 else if (watmodel == WAT_SWM4) wathgsize = 5;
01358
01359 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01360 int hgs = atom[ig].hydrogenGroupSize;
01361 if ( hgs == 1 ) continue; // only one atom in group
01362 // cache data in local arrays and integrate positions normally
01363 for ( i = 0; i < hgs; ++i ) {
01364 ref[i] = atom[ig+i].position;
01365 pos[i] = atom[ig+i].position;
01366 vel[i] = atom[ig+i].velocity;
01367 rmass[i] = (atom[ig+1].mass > 0. ? 1. / atom[ig+i].mass : 0.);
01368 //printf("rmass of %i is %f\n", ig+i, rmass[i]);
01369 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
01370 //printf("fixed status of %i is %i\n", i, fixed[i]);
01371 // undo addVelocityToPosition to get proper reference coordinates
01372 if ( fixed[i] ) rmass[i] = 0.; else pos[i] += vel[i] * dt;
01373 }
01374 int icnt = 0;
01375 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) > 0 ) { // for water
01376 if (hgs != wathgsize) {
01377 NAMD_bug("Hydrogen group error caught in rattle1().");
01378 }
01379 // Use SETTLE for water unless some of the water atoms are fixed,
01380 // for speed we test groupFixed rather than the individual atoms
01381 if (useSettle && !atom[ig].groupFixed) {
01382 if (simParams->watmodel == WAT_SWM4) {
01383 // SWM4 ordering: O D LP H1 H2
01384 // do swap(O,LP) and call settle with subarray O H1 H2
01385 // swap back after we return
01386 Vector lp_ref = ref[2];
01387 Vector lp_pos = pos[2];
01388 Vector lp_vel = vel[2];
01389 ref[2] = ref[0];
01390 pos[2] = pos[0];
01391 vel[2] = vel[0];
01392 settle1(ref+2, pos+2, vel+2, invdt);
01393 ref[0] = ref[2];
01394 pos[0] = pos[2];
01395 vel[0] = vel[2];
01396 ref[2] = lp_ref;
01397 pos[2] = lp_pos;
01398 vel[2] = lp_vel;
01399 // determine for LP updated pos and vel
01400 swm4_omrepos(ref, pos, vel, invdt);
01401 }
01402 else {
01403 settle1(ref, pos, vel, invdt);
01404 if (simParams->watmodel == WAT_TIP4) {
01405 tip4_omrepos(ref, pos, vel, invdt);
01406 }
01407 }
01408
01409 // which slab the hydrogen group will belong to
01410 // for pprofile calculations.
01411 int ppoffset, partition;
01412 if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
01413 atom[ig+i].position = pos[i];
01414 } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
01415 atom[ig+i].velocity = vel[i];
01416 } else for ( i = 0; i < wathgsize; ++i ) {
01417 Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
01418 Tensor vir = outer(df, ref[i]);
01419 wc += vir;
01420 f[Results::normal][ig+i] += df;
01421 atom[ig+i].velocity = vel[i];
01422 if (ppreduction) {
01423 // put all the atoms from a water in the same slab. Atom 0
01424 // should be the parent atom.
01425 if (!i) {
01426 BigReal z = pos[i].z;
01427 partition = atom[ig].partition;
01428 int slab = (int)floor((z-zmin)*idz);
01429 if (slab < 0) slab += nslabs;
01430 else if (slab >= nslabs) slab -= nslabs;
01431 ppoffset = 3*(slab + nslabs*partition);
01432 }
01433 ppreduction->item(ppoffset ) += vir.xx;
01434 ppreduction->item(ppoffset+1) += vir.yy;
01435 ppreduction->item(ppoffset+2) += vir.zz;
01436 }
01437 }
01438 continue;
01439 }
01440 if ( !(fixed[1] && fixed[2]) ) {
01441 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01442 }
01443 }
01444 for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
01445 if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) > 0 ) {
01446 if ( !(fixed[0] && fixed[i]) ) {
01447 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01448 }
01449 }
01450 }
01451 if ( icnt == 0 ) continue; // no constraints
01452 for ( i = 0; i < icnt; ++i ) {
01453 refab[i] = ref[ial[i]] - ref[ibl[i]];
01454 }
01455 for ( i = 0; i < hgs; ++i ) {
01456 netdp[i] = 0.;
01457 }
01458 int done;
01459 int consFailure;
01460 for ( iter = 0; iter < maxiter; ++iter ) {
01461 //if (iter > 0) CkPrintf("iteration %d\n", iter);
01462 done = 1;
01463 consFailure = 0;
01464 for ( i = 0; i < icnt; ++i ) {
01465 int a = ial[i]; int b = ibl[i];
01466 Vector pab = pos[a] - pos[b];
01467 BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
01468 BigReal rabsq = dsq[i];
01469 BigReal diffsq = rabsq - pabsq;
01470 if ( fabs(diffsq) > (rabsq * tol2) ) {
01471 Vector &rab = refab[i];
01472 BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
01473 if ( rpab < ( rabsq * 1.0e-6 ) ) {
01474 done = 0;
01475 consFailure = 1;
01476 continue;
01477 }
01478 BigReal rma = rmass[a];
01479 BigReal rmb = rmass[b];
01480 BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
01481 Vector dp = rab * gab;
01482 pos[a] += rma * dp;
01483 pos[b] -= rmb * dp;
01484 if ( invdt != 0. ) {
01485 dp *= invdt;
01486 netdp[a] += dp;
01487 netdp[b] -= dp;
01488 }
01489 done = 0;
01490 }
01491 }
01492 if ( done ) break;
01493 }
01494
01495 if ( consFailure ) {
01496 if ( dieOnError ) {
01497 iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
01498 << (atom[ig].id + 1) << "!\n" << endi;
01499 return -1; // triggers early exit
01500 } else {
01501 iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
01502 << (atom[ig].id + 1) << "!\n" << endi;
01503 }
01504 } else if ( ! done ) {
01505 if ( dieOnError ) {
01506 iout << iERROR << "Exceeded RATTLE iteration limit for atom "
01507 << (atom[ig].id + 1) << "!\n" << endi;
01508 return -1; // triggers early exit
01509 } else {
01510 iout << iWARN << "Exceeded RATTLE iteration limit for atom "
01511 << (atom[ig].id + 1) << "!\n" << endi;
01512 }
01513 }
01514
01515 // store data back to patch
01516 int ppoffset, partition;
01517 if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
01518 atom[ig+i].position = pos[i];
01519 } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
01520 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
01521 } else for ( i = 0; i < hgs; ++i ) {
01522 Force df = netdp[i] * invdt;
01523 Tensor vir = outer(df, ref[i]);
01524 wc += vir;
01525 f[Results::normal][ig+i] += df;
01526 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
01527 if (ppreduction) {
01528 if (!i) {
01529 BigReal z = pos[i].z;
01530 int partition = atom[ig].partition;
01531 int slab = (int)floor((z-zmin)*idz);
01532 if (slab < 0) slab += nslabs;
01533 else if (slab >= nslabs) slab -= nslabs;
01534 ppoffset = 3*(slab + nslabs*partition);
01535 }
01536 ppreduction->item(ppoffset ) += vir.xx;
01537 ppreduction->item(ppoffset+1) += vir.yy;
01538 ppreduction->item(ppoffset+2) += vir.zz;
01539 }
01540 }
01541 }
01542 if ( dt && virial ) *virial += wc;
01543
01544 return 0;
01545 }
|
|
||||||||||||
|
Definition at line 1548 of file HomePatch.C. References BigReal, SimParameters::fixedAtomsOn, iout, iWARN(), Node::molecule, NAMD_bug(), NAMD_die(), Node::Object(), outer(), Molecule::rigid_bond_length(), SimParameters::rigidDie, SimParameters::rigidIter, SimParameters::rigidTol, settle2(), Node::simParameters, simParams, SimParameters::useSettle, SimParameters::watmodel, Vector::x, Vector::y, and Vector::z. Referenced by Sequencer::rattle2(). 01549 {
01550 Molecule *mol = Node::Object()->molecule;
01551 SimParameters *simParams = Node::Object()->simParameters;
01552 const int fixedAtomsOn = simParams->fixedAtomsOn;
01553 const int useSettle = simParams->useSettle;
01554 const BigReal dt = timestep / TIMEFACTOR;
01555 Tensor wc; // constraint virial
01556 BigReal tol = simParams->rigidTol;
01557 int maxiter = simParams->rigidIter;
01558 int dieOnError = simParams->rigidDie;
01559 int i, iter;
01560 BigReal dsqi[10], tmp;
01561 int ial[10], ibl[10];
01562 Vector ref[10]; // reference position
01563 Vector refab[10]; // reference vector
01564 Vector vel[10]; // new velocity
01565 BigReal rmass[10]; // 1 / mass
01566 BigReal redmass[10]; // reduced mass
01567 int fixed[10]; // is atom fixed?
01568
01569 // Size of a hydrogen group for water
01570 int wathgsize = 3;
01571 if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
01572 else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
01573
01574 // CkPrintf("In rattle2!\n");
01575 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01576 // CkPrintf("ig=%d\n",ig);
01577 int hgs = atom[ig].hydrogenGroupSize;
01578 if ( hgs == 1 ) continue; // only one atom in group
01579 // cache data in local arrays and integrate positions normally
01580 for ( i = 0; i < hgs; ++i ) {
01581 ref[i] = atom[ig+i].position;
01582 vel[i] = atom[ig+i].velocity;
01583 rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
01584 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
01585 if ( fixed[i] ) rmass[i] = 0.;
01586 }
01587 int icnt = 0;
01588 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) > 0 ) { // for water
01589 if ( wathgsize != 4 ) {
01590 NAMD_bug("Hydrogen group error caught in rattle2().");
01591 }
01592 // Use SETTLE for water unless some of the water atoms are fixed
01593 if (useSettle && !fixed[0] && !fixed[1] && !fixed[2]) {
01594 settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
01595 for (i=0; i<3; i++) {
01596 atom[ig+i].velocity = vel[i];
01597 }
01598 continue;
01599 }
01600 if ( !(fixed[1] && fixed[2]) ) {
01601 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
01602 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01603 }
01604 }
01605 // CkPrintf("Loop 2\n");
01606 for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
01607 if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) > 0 ) {
01608 if ( !(fixed[0] && fixed[i]) ) {
01609 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
01610 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
01611 ibl[icnt] = i; ++icnt;
01612 }
01613 }
01614 }
01615 if ( icnt == 0 ) continue; // no constraints
01616 // CkPrintf("Loop 3\n");
01617 for ( i = 0; i < icnt; ++i ) {
01618 refab[i] = ref[ial[i]] - ref[ibl[i]];
01619 }
01620 // CkPrintf("Loop 4\n");
01621 int done;
01622 for ( iter = 0; iter < maxiter; ++iter ) {
01623 done = 1;
01624 for ( i = 0; i < icnt; ++i ) {
01625 int a = ial[i]; int b = ibl[i];
01626 Vector vab = vel[a] - vel[b];
01627 Vector &rab = refab[i];
01628 BigReal rabsqi = dsqi[i];
01629 BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
01630 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
01631 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
01632 wc += outer(dp,rab);
01633 vel[a] += rmass[a] * dp;
01634 vel[b] -= rmass[b] * dp;
01635 done = 0;
01636 }
01637 }
01638 if ( done ) break;
01639 //if (done) { if (iter > 0) CkPrintf("iter=%d\n", iter); break; }
01640 }
01641 if ( ! done ) {
01642 if ( dieOnError ) {
01643 NAMD_die("Exceeded maximum number of iterations in rattle2().");
01644 } else {
01645 iout << iWARN <<
01646 "Exceeded maximum number of iterations in rattle2().\n" << endi;
01647 }
01648 }
01649 // store data back to patch
01650 for ( i = 0; i < hgs; ++i ) {
01651 atom[ig+i].velocity = vel[i];
01652 }
01653 }
01654 // CkPrintf("Leaving rattle2!\n");
01655 // check that there isn't a constant needed here!
01656 *virial += wc / ( 0.5 * dt );
01657
01658 }
|
|
|
Definition at line 828 of file HomePatch.C. References ResizeArray< Elem >::begin(), Box< Owner, Data >::close(), DebugM, ResizeArray< Elem >::end(), Results::f, Force, ProxyListElem::forceBox, ProxyCombinedResultMsg::forceList, ProxyListElem::node, ProxyCombinedResultMsg::nodes, Box< Owner, Data >::open(), ResizeArray< Elem >::size(), Vector::x, Vector::y, and Vector::z. 00829 {
00830 DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodes.size()<<")\n");
00831 //CkPrintf("[%d] Homepatch: %d receiveResults from %d nodes\n", CkMyPe(), patchID, n);
00832 for (int i=0; i<msg->nodes.size(); i++) {
00833 int node = msg->nodes[i];
00834 ProxyListElem *pe = proxy.begin();
00835 for ( ; pe->node != node; ++pe );
00836 Results *r = pe->forceBox->open();
00837 if (i==0) {
00838 for ( int k = 0; k < Results::maxNumForces; ++k )
00839 {
00840 Force *f = r->f[k];
00841 register ForceList::iterator f_i, f_e;
00842 f_i = msg->forceList[k].begin();
00843 f_e = msg->forceList[k].end();
00844 //for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
00845
00846 int nf = f_e - f_i;
00847 #ifdef ARCH_POWERPC
00848 #pragma disjoint (*f_i, *f)
00849 #pragma unroll(4)
00850 #endif
00851 for (int count = 0; count < nf; count++) {
00852 f[count].x += f_i[count].x;
00853 f[count].y += f_i[count].y;
00854 f[count].z += f_i[count].z;
00855 }
00856 }
00857 }
00858 pe->forceBox->close(&r);
00859 }
00860
00861 delete msg;
00862 }
|
|
|
Definition at line 809 of file HomePatch.C. References ResizeArray< Elem >::begin(), Box< Owner, Data >::close(), DebugM, ResizeArray< Elem >::end(), Results::f, Force, ProxyListElem::forceBox, ProxyResultMsg::forceList, ProxyListElem::node, ProxyResultMsg::node, and Box< Owner, Data >::open(). 00810 {
00811 DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00812 int n = msg->node;
00813 ProxyListElem *pe = proxy.begin();
00814 for ( ; pe->node != n; ++pe );
00815 Results *r = pe->forceBox->open();
00816 for ( int k = 0; k < Results::maxNumForces; ++k )
00817 {
00818 Force *f = r->f[k];
00819 register ForceList::iterator f_i, f_e;
00820 f_i = msg->forceList[k].begin();
00821 f_e = msg->forceList[k].end();
00822 for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
00823 }
00824 pe->forceBox->close(&r);
00825 delete msg;
00826 }
|
|
|
Definition at line 785 of file HomePatch.C. References ResizeArray< Elem >::begin(), Box< Owner, Data >::close(), DebugM, Results::f, ProxyResultVarsizeMsg::flLen, Force, ProxyResultVarsizeMsg::forceArr, ProxyListElem::forceBox, ProxyResultVarsizeMsg::isZero, ProxyListElem::node, ProxyResultVarsizeMsg::node, and Box< Owner, Data >::open(). Referenced by ProxyMgr::recvResults(). 00785 {
00786 DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00787 int n = msg->node;
00788 ProxyListElem *pe = proxy.begin();
00789 for ( ; pe->node != n; ++pe );
00790 Results *r = pe->forceBox->open();
00791
00792 char *iszeroPtr = msg->isZero;
00793 Force *msgFPtr = msg->forceArr;
00794
00795 for ( int k = 0; k < Results::maxNumForces; ++k )
00796 {
00797 Force *rfPtr = r->f[k];
00798 for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
00799 if((*iszeroPtr)!=1) {
00800 *rfPtr += *msgFPtr;
00801 msgFPtr++;
00802 }
00803 }
00804 }
00805 pe->forceBox->close(&r);
00806 delete msg;
00807 }
|
|
|
Definition at line 605 of file HomePatch.C. Referenced by ProxyMgr::recvNodeAwareSpanningTreeOnHomePatch(). 00605 {}
|
|
||||||||||||
|
Definition at line 611 of file HomePatch.C. References ResizeArray< Elem >::resize(), sendSpanningTree(), and ResizeArray< Elem >::size(). Referenced by ProxyMgr::recvSpanningTreeOnHomePatch(). 00612 {
00613 int i;
00614 nChild=0;
00615 tree.resize(n);
00616 for (i=0; i<n; i++) {
00617 tree[i] = t[i];
00618 }
00619
00620 for (i=1; i<=proxySpanDim; i++) {
00621 if (tree.size() <= i) break;
00622 child[i-1] = tree[i];
00623 nChild++;
00624 }
00625
00626 // send down to children
00627 sendSpanningTree();
00628 }
|
|
|
Definition at line 373 of file HomePatch.C. References ResizeArray< Elem >::add(), ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::checkOut(), DebugM, RegisterProxyMsg::node, and ResizeArray< Elem >::size(). Referenced by ProxyMgr::recvRegisterProxy(). 00373 {
00374 DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
00375 proxy.add(ProxyListElem(msg->node,forceBox.checkOut()));
00376
00377 isNewProxyAdded = 1;
00378
00379 Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
00380 delete msg;
00381 }
|
|
|
Definition at line 1030 of file HomePatch.C. 01031 {
01032 replacementForces = f;
01033 }
|
|
|
Definition at line 1809 of file HomePatch.C. References FullAtomList, and ResizeArray< Elem >::size(). Referenced by Sequencer::algorithm(). 01809 {
01810 FullAtomList tmp_a(&checkpoint_atom); atom = tmp_a;
01811 numAtoms = atom.size();
01812 lattice = checkpoint_lattice;
01813
01814 // DMK - Atom Separation (water vs. non-water)
01815 #if NAMD_SeparateWaters != 0
01816 numWaterAtoms = checkpoint_numWaterAtoms;
01817 #endif
01818 }
|
|
|
Definition at line 282 of file HomePatch.C. References Sequencer::run(). Referenced by Node::run(). 00283 { sequencer->run(); }
|
|
|
Definition at line 1035 of file HomePatch.C. References ResizeArray< Elem >::resize(). Referenced by Sequencer::saveForce(). 01036 {
01037 f_saved[ftag].resize(numAtoms);
01038 for ( int i = 0; i < numAtoms; ++i )
01039 {
01040 f_saved[ftag][i] = f[ftag][i];
01041 }
01042 }
|
|
|
Definition at line 606 of file HomePatch.C. 00606 {}
|
|
|
Definition at line 432 of file HomePatch.C. References ResizeArray< Elem >::add(), ResizeArrayIter< T >::begin(), ResizeArrayIter< T >::end(), NodeIDList, ProxyMgr::Object(), ProxyListIter, ProxyMgr::sendProxies(), and ResizeArray< Elem >::size(). Referenced by ProxyMgr::buildProxySpanningTree2(). 00433 {
00434 ProxyListIter pli(proxy);
00435 NodeIDList list;
00436 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00437 {
00438 list.add(pli->node);
00439 }
00440 ProxyMgr::Object()->sendProxies(patchID, &list[0], list.size());
00441 }
|
|
|
Definition at line 630 of file HomePatch.C. References ProxySpanningTreeMsg::node, ProxyMgr::Object(), ProxySpanningTreeMsg::patch, ProxyMgr::sendSpanningTree(), ResizeArray< Elem >::size(), and ProxySpanningTreeMsg::tree. Referenced by buildSpanningTree(), and recvSpanningTree(). 00631 {
00632 if (tree.size() == 0) return;
00633 ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
00634 msg->patch = patchID;
00635 msg->node = CkMyPe();
00636 msg->tree = tree;
00637 ProxyMgr::Object()->sendSpanningTree(msg);
00638 }
|
|
|
Definition at line 120 of file HomePatch.h. References FullAtomList. Referenced by PatchMgr::fillHomePatchAtomList(). 00120 {
00121 atom = *al;
00122 }
|
|
|
Definition at line 1820 of file HomePatch.C. References LdbCoordinator::Object(), and LdbCoordinator::patchLoad(). Referenced by Sequencer::rebalanceLoad(). 01821 {
01822 LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
01823 }
|
|
|
Definition at line 383 of file HomePatch.C. References ResizeArray< Elem >::begin(), OwnerBox< Owner, Data >::checkIn(), ResizeArray< Elem >::del(), ProxyListElem::forceBox, ProxyListElem::node, and UnregisterProxyMsg::node. Referenced by ProxyMgr::recvUnregisterProxy(). 00383 {
00384 int n = msg->node;
00385 ProxyListElem *pe = proxy.begin();
00386 for ( ; pe->node != n; ++pe );
00387 forceBox.checkIn(pe->forceBox);
00388 proxy.del(pe - proxy.begin());
00389 delete msg;
00390 }
|
|
|
Definition at line 278 of file HomePatch.C. 00279 { sequencer=sequencerPtr; }
|
|
|
Definition at line 40 of file HomePatch.h. |
|
|
Definition at line 38 of file HomePatch.h. |
|
|
Definition at line 39 of file HomePatch.h. |
|
|
Definition at line 154 of file HomePatch.h. Referenced by doAtomMigration(). |
|
|
Definition at line 89 of file HomePatch.h. Referenced by doAtomMigration(), doMarginCheck(), and Sequencer::submitReductions(). |
|
|
Definition at line 156 of file HomePatch.h. Referenced by depositMigration(), and doAtomMigration(). |
|
|
Definition at line 155 of file HomePatch.h. Referenced by depositMigration(). |
1.3.9.1