#include <HomePatch.h>
Inheritance diagram for HomePatch:
|
|
Definition at line 347 of file HomePatch.C. 00348 {
00349 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00350 ptnTree.resize(0);
00351 delete [] children;
00352 #else
00353 delete [] child;
00354 #endif
00355 }
|
|
||||||||||||||||||||
|
Definition at line 1145 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, SimParameters::watmodel, Vector::x, Vector::y, and Vector::z. Referenced by Sequencer::addForceToMomentum(). 01147 {
01148 SimParameters *simParams = Node::Object()->simParameters;
01149 const BigReal dt = timestep / TIMEFACTOR;
01150 ForceList *f_use;
01151
01152 // Apply corrections for massless particles
01153 // For now this is specific to TIP4P water
01154 if (simParams->watmodel == WAT_TIP4) {
01155
01156 ForceList *f_mod = ( useSaved ? f_saved : f );
01157 for (int i=0; i<numAtoms; i++) {
01158 if (atom[i].mass < 0.01) {
01159 redistrib_tip4p_forces( f_mod[ftag][i-3], f_mod[ftag][i-2], f_mod[ftag][i-1], f_mod[ftag][i], i-3 , virial);
01160 }
01161 }
01162
01163 f_use = f_mod;
01164
01165 } else {
01166 f_use = ( useSaved ? f_saved : f );
01167 }
01168
01169 if ( simParams->fixedAtomsOn ) {
01170 for ( int i = 0; i < numAtoms; ++i ) {
01171 if ( atom[i].atomFixed || atom[i].mass <= 0. ) atom[i].velocity = 0;
01172 else atom[i].velocity += f_use[ftag][i] * ( dt * namd_reciprocal (atom[i].mass) );
01173 }
01174 } else {
01175 FullAtom *atom_arr = atom.begin();
01176 const Force *force_arr = f_use[ftag].const_begin();
01177 #ifdef ARCH_POWERPC
01178 #pragma disjoint (*force_arr, *atom_arr)
01179 #endif
01180 for ( int i = 0; i < numAtoms; ++i ) {
01181 if (atom[i].mass == 0.) continue;
01182 BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.);
01183 //printf("Taking reciprocal of mass %f\n", atom[i].mass);
01184 atom_arr[i].velocity.x += force_arr[i].x * recip_val;
01185 atom_arr[i].velocity.y += force_arr[i].y * recip_val;
01186 atom_arr[i].velocity.z += force_arr[i].z * recip_val;
01187 }
01188 }
01189 }
|
|
|
Definition at line 1191 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(). 01192 {
01193 SimParameters *simParams = Node::Object()->simParameters;
01194 const BigReal dt = timestep / TIMEFACTOR;
01195 if ( simParams->fixedAtomsOn ) {
01196 for ( int i = 0; i < numAtoms; ++i ) {
01197 if ( ! atom[i].atomFixed ) atom[i].position += atom[i].velocity * dt;
01198 }
01199 } else {
01200 FullAtom *atom_arr = atom.begin();
01201 for ( int i = 0; i < numAtoms; ++i ) {
01202 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
01203 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
01204 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
01205 }
01206 }
01207 }
|
|
|
Implements Patch. Definition at line 358 of file HomePatch.C. References Sequencer::awaken(), DebugM, ExtForce::force, ExtForce::replace, and Sequencer::thread. 00359 {
00360 if ( ! --boxesOpen )
00361 {
00362 if ( replacementForces ) {
00363 for ( int i = 0; i < numAtoms; ++i ) {
00364 if ( replacementForces[i].replace ) {
00365 for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
00366 f[Results::normal][i] = replacementForces[i].force;
00367 }
00368 }
00369 replacementForces = 0;
00370 }
00371 DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
00372 << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
00373 // only awaken suspended threads. Then say it is suspended.
00374 sequencer->awaken();
00375 return;
00376 }
00377 else
00378 {
00379 DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
00380 }
00381 }
|
|
|
Definition at line 617 of file HomePatch.C. References ResizeArrayIter< T >::begin(), compDistance(), ResizeArrayIter< T >::end(), ResizeArray< Elem >::find(), 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(). 00618 {
00619 nChild = 0;
00620 int psize = proxy.size();
00621 if (psize == 0) return;
00622 NodeIDList oldtree = tree;
00623 int oldsize = oldtree.size();
00624 tree.resize(psize + 1);
00625 tree.setall(-1);
00626 tree[0] = CkMyPe();
00627 int s=1, e=psize+1;
00628 ProxyListIter pli(proxy);
00629 int patchNodesLast =
00630 ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00631 int nNonPatch = 0;
00632 #if 0
00633 int* pelists = new int[psize];
00634 for (int i=0; i<psize; i++) pelists[i] = -1;
00635 for ( pli = pli.begin(); pli != pli.end(); ++pli ) {
00636 int idx = rand()%psize;
00637 while (pelists[idx] != -1) { idx++; if (idx == psize) idx = 0; }
00638 pelists[idx] = pli->node;
00639 }
00640 for ( int i=0; i<psize; i++)
00641 {
00642 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pelists[i]) ) {
00643 tree[--e] = pelists[i];
00644 } else {
00645 tree[s++] = pelists[i];
00646 nNonPatch++;
00647 }
00648 }
00649 delete [] pelists;
00650 #else
00651 // try to put it to the same old tree
00652 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00653 {
00654 int oldindex = oldtree.find(pli->node);
00655 if (oldindex != -1 && oldindex < psize) {
00656 tree[oldindex] = pli->node;
00657 }
00658 }
00659 s=1; e=psize;
00660 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00661 {
00662 if (tree.find(pli->node) != -1) continue; // already assigned
00663 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
00664 while (tree[e] != -1) { e--; if (e==-1) e = psize; }
00665 tree[e] = pli->node;
00666 } else {
00667 while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
00668 tree[s] = pli->node;
00669 nNonPatch++;
00670 }
00671 }
00672 #if 1
00673 if (oldsize==0 && nNonPatch) {
00674 // first time, sort by distance
00675 qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
00676 }
00677 #endif
00678
00679 //CkPrintf("home: %d:(%d) %d %d %d %d %d\n", patchID, tree.size(),tree[0],tree[1],tree[2],tree[3],tree[4]);
00680
00681 #if USE_TOPOMAP && USE_SPANNING_TREE
00682
00683 //Right now only works for spanning trees with two levels
00684 int *treecopy = new int [psize];
00685 int subroots[proxySpanDim];
00686 int subsizes[proxySpanDim];
00687 int subtrees[proxySpanDim][proxySpanDim];
00688 int idxes[proxySpanDim];
00689 int i = 0;
00690
00691 for(i=0;i<proxySpanDim;i++){
00692 subsizes[i] = 0;
00693 idxes[i] = i;
00694 }
00695
00696 for(i=0;i<psize;i++){
00697 treecopy[i] = tree[i+1];
00698 }
00699
00700 TopoManager tmgr;
00701 tmgr.sortRanksByHops(treecopy,nNonPatch);
00702 tmgr.sortRanksByHops(treecopy+nNonPatch,
00703 psize-nNonPatch);
00704
00705 /* build tree and subtrees */
00706 nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
00707 delete [] treecopy;
00708
00709 for(int i=1;i<psize+1;i++){
00710 int isSubroot=0;
00711 for(int j=0;j<nChild;j++)
00712 if(tree[i]==subroots[j]){
00713 isSubroot=1;
00714 break;
00715 }
00716 if(isSubroot) continue;
00717
00718 int bAdded = 0;
00719 tmgr.sortIndexByHops(tree[i], subroots,
00720 idxes, proxySpanDim);
00721 for(int j=0;j<proxySpanDim;j++){
00722 if(subsizes[idxes[j]]<proxySpanDim){
00723 subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
00724 bAdded = 1;
00725 break;
00726 }
00727 }
00728 if( psize > proxySpanDim && ! bAdded ) {
00729 NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
00730 }
00731 }
00732
00733 #else /* USE_TOPOMAP && USE_SPANNING_TREE */
00734
00735 for (int i=1; i<=proxySpanDim; i++) {
00736 if (tree.size() <= i) break;
00737 child[i-1] = tree[i];
00738 nChild++;
00739 }
00740 #endif
00741 #endif
00742
00743 #if 0
00744 // for debugging
00745 CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
00746 for (int i=0; i<psize+1; i++) {
00747 CkPrintf("%d ", tree[i]);
00748 }
00749 CkPrintf("\n");
00750 #endif
00751 // send to children nodes
00752 sendSpanningTree();
00753 }
|
|
|
Definition at line 1686 of file HomePatch.C. References FullAtomList. Referenced by Sequencer::algorithm(). 01686 {
01687 FullAtomList tmp_a(&atom); checkpoint_atom = tmp_a;
01688 checkpoint_lattice = lattice;
01689
01690 // DMK - Atom Separation (water vs. non-water)
01691 #if NAMD_SeparateWaters != 0
01692 checkpoint_numWaterAtoms = numWaterAtoms;
01693 #endif
01694 }
|
|
|
Definition at line 2034 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(). 02035 {
02036
02037 if (!inMigration) { // We have to buffer changes due to migration
02038 // until our patch is in migration mode
02039 msgbuf[numMlBuf++] = msg;
02040 return;
02041 }
02042
02043
02044 // DMK - Atom Separation (water vs. non-water)
02045 #if NAMD_SeparateWaters != 0
02046
02047
02048 // Merge the incoming list of atoms with the current list of
02049 // atoms. Note that mergeSeparatedAtomList() will apply any
02050 // required transformations to the incoming atoms as it is
02051 // separating them.
02052 mergeAtomList(msg->migrationList);
02053
02054
02055 #else
02056
02057 // Merge the incoming list of atoms with the current list of
02058 // atoms. Apply transformations to the incoming atoms as they are
02059 // added to this patch's list.
02060 {
02061 MigrationList &migrationList = msg->migrationList;
02062 MigrationList::iterator mi;
02063 Transform mother_transform;
02064 for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
02065 DebugM(1,"Migrating atom " << mi->id << " to patch "
02066 << patchID << " with position " << mi->position << "\n");
02067 if ( mi->hydrogenGroupSize ) {
02068 mi->position = lattice.nearest(mi->position,center,&(mi->transform));
02069 mother_transform = mi->transform;
02070 } else {
02071 mi->position = lattice.reverse_transform(mi->position,mi->transform);
02072 mi->position = lattice.apply_transform(mi->position,mother_transform);
02073 mi->transform = mother_transform;
02074 }
02075 atom.add(*mi);
02076 }
02077 }
02078
02079
02080 #endif // if (NAMD_SeparateWaters != 0)
02081
02082
02083 numAtoms = atom.size();
02084 delete msg;
02085
02086 DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
02087 if (!--patchMigrationCounter) {
02088 DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
02089 allMigrationIn = true;
02090 patchMigrationCounter = numNeighbors;
02091 if (migrationSuspended) {
02092 DebugM(3,"patch "<<patchID<<" is being awakened\n");
02093 migrationSuspended = false;
02094 sequencer->awaken();
02095 return;
02096 }
02097 else {
02098 DebugM(3,"patch "<<patchID<<" was not suspended\n");
02099 }
02100 }
02101 }
|
|
|
Definition at line 1909 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(). 01910 {
01911 int i;
01912
01913 for (i=0; i<numNeighbors; i++) {
01914 realInfo[i].mList.resize(0);
01915 }
01916
01917 // Purge the AtomMap
01918 AtomMap::Object()->unregisterIDs(patchID,p.begin(),p.end());
01919
01920 // Determine atoms that need to migrate
01921
01922 BigReal minx = min.x;
01923 BigReal miny = min.y;
01924 BigReal minz = min.z;
01925 BigReal maxx = max.x;
01926 BigReal maxy = max.y;
01927 BigReal maxz = max.z;
01928
01929 int xdev, ydev, zdev;
01930 int delnum = 0;
01931
01932 FullAtomList::iterator atom_i = atom.begin();
01933 FullAtomList::iterator atom_e = atom.end();
01934
01935 // DMK - Atom Separation (water vs. non-water)
01936 #if NAMD_SeparateWaters != 0
01937 FullAtomList::iterator atom_first = atom_i;
01938 int numLostWaterAtoms = 0;
01939 #endif
01940
01941 while ( atom_i != atom_e ) {
01942 if ( atom_i->hydrogenGroupSize ) {
01943
01944 ScaledPosition s = lattice.scale(atom_i->position);
01945
01946 // check if atom is within bounds
01947 if (s.x < minx) xdev = 0;
01948 else if (maxx <= s.x) xdev = 2;
01949 else xdev = 1;
01950
01951 if (s.y < miny) ydev = 0;
01952 else if (maxy <= s.y) ydev = 2;
01953 else ydev = 1;
01954
01955 if (s.z < minz) zdev = 0;
01956 else if (maxz <= s.z) zdev = 2;
01957 else zdev = 1;
01958
01959 }
01960
01961 if (mInfo[xdev][ydev][zdev]) { // process atom for migration
01962 // Don't migrate if destination is myself
01963
01964 // See if we have a migration list already
01965 MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
01966 DebugM(3,"Migrating atom " << atom_i->id << " from patch "
01967 << patchID << " with position " << atom_i->position << "\n");
01968 mCur.add(*atom_i);
01969
01970 ++delnum;
01971
01972
01973 // DMK - Atom Separation (water vs. non-water)
01974 #if NAMD_SeparateWaters != 0
01975 // Check to see if this atom is part of a water molecule. If
01976 // so, numWaterAtoms needs to be adjusted to refect the lost of
01977 // this atom.
01978 // NOTE: The atom separation code assumes that if the oxygen
01979 // atom of the hydrogen group making up the water molecule is
01980 // migrated to another HomePatch, the hydrogens will also
01981 // move!!!
01982 int atomIndex = atom_i - atom_first;
01983 if (atomIndex < numWaterAtoms)
01984 numLostWaterAtoms++;
01985 #endif
01986
01987
01988 } else {
01989
01990 if ( delnum ) { *(atom_i-delnum) = *atom_i; }
01991
01992 }
01993
01994 ++atom_i;
01995 }
01996
01997 // DMK - Atom Separation (water vs. non-water)
01998 #if NAMD_SeparateWaters != 0
01999 numWaterAtoms -= numLostWaterAtoms;
02000 #endif
02001
02002
02003 int delpos = numAtoms - delnum;
02004 DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
02005 atom.del(delpos,delnum);
02006
02007 numAtoms = atom.size();
02008
02009 PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
02010
02011 // signal depositMigration() that we are inMigration mode
02012 inMigration = true;
02013
02014 // Drain the migration message buffer
02015 for (i=0; i<numMlBuf; i++) {
02016 DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
02017 depositMigration(msgbuf[i]);
02018 }
02019 numMlBuf = 0;
02020
02021 if (!allMigrationIn) {
02022 DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
02023 migrationSuspended = true;
02024 sequencer->suspend();
02025 migrationSuspended = false;
02026 }
02027 allMigrationIn = false;
02028
02029 inMigration = false;
02030 marginViolations = 0;
02031 }
|
|
|
Definition at line 1788 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(). 01789 {
01790 if ( ! flags.doNonbonded ) return;
01791
01792 SimParameters *simParams = Node::Object()->simParameters;
01793 BigReal hgcut = 0.5 * simParams->hgroupCutoff; hgcut *= hgcut;
01794 BigReal maxrad2 = 0.;
01795
01796 FullAtomList::iterator p_i = atom.begin();
01797 FullAtomList::iterator p_e = atom.end();
01798
01799 while ( p_i != p_e ) {
01800 int hgs = p_i->hydrogenGroupSize;
01801 p_i->nonbondedGroupIsAtom = 0;
01802 BigReal x = p_i->position.x;
01803 BigReal y = p_i->position.y;
01804 BigReal z = p_i->position.z;
01805 ++p_i;
01806 int oversize = 0;
01807 // limit spatial extent
01808 for ( int i = 1; i < hgs; ++i ) {
01809 p_i->nonbondedGroupIsAtom = 0;
01810 BigReal dx = p_i->position.x - x;
01811 BigReal dy = p_i->position.y - y;
01812 BigReal dz = p_i->position.z - z;
01813 BigReal r2 = dx * dx + dy * dy + dz * dz;
01814 ++p_i;
01815 if ( r2 > hgcut ) oversize = 1;
01816 else if ( r2 > maxrad2 ) maxrad2 = r2;
01817 }
01818 // also limit to at most 4 atoms per group
01819 if ( oversize || hgs > 4 ) {
01820 p_i -= hgs;
01821 for ( int i = 0; i < hgs; ++i ) {
01822 p_i->nonbondedGroupIsAtom = 1;
01823 ++p_i;
01824 }
01825 }
01826 }
01827
01828 flags.maxGroupRadius = sqrt(maxrad2);
01829
01830 }
|
|
|
Definition at line 1832 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(). 01833 {
01834 SimParameters *simParams = Node::Object()->simParameters;
01835
01836 BigReal sysdima = lattice.a_r().unit() * lattice.a();
01837 BigReal sysdimb = lattice.b_r().unit() * lattice.b();
01838 BigReal sysdimc = lattice.c_r().unit() * lattice.c();
01839
01840 BigReal minSize = simParams->patchDimension - simParams->margin;
01841
01842 if ( ( (max.x - min.x)*aAway*sysdima < minSize*0.9999 ) ||
01843 ( (max.y - min.y)*bAway*sysdimb < minSize*0.9999 ) ||
01844 ( (max.z - min.z)*cAway*sysdimc < minSize*0.9999 ) ) {
01845
01846 NAMD_die("Periodic cell has become too small for original patch grid!\n"
01847 "Possible solutions are to restart from a recent checkpoint,\n"
01848 "increase margin, or disable useFlexibleCell for liquid simulation.");
01849 }
01850
01851 BigReal cutoff = simParams->cutoff;
01852
01853 BigReal margina = 0.5 * ( (max.x - min.x) * aAway - cutoff / sysdima );
01854 BigReal marginb = 0.5 * ( (max.y - min.y) * bAway - cutoff / sysdimb );
01855 BigReal marginc = 0.5 * ( (max.z - min.z) * cAway - cutoff / sysdimc );
01856
01857 if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
01858 NAMD_die("Periodic cell has become too small for original patch grid!\n"
01859 "There are probably many margin violations already reported.\n"
01860 "Possible solutions are to restart from a recent checkpoint,\n"
01861 "increase margin, or disable useFlexibleCell for liquid simulation.");
01862 }
01863
01864 BigReal minx = min.x - margina;
01865 BigReal miny = min.y - marginb;
01866 BigReal minz = min.z - marginc;
01867 BigReal maxx = max.x + margina;
01868 BigReal maxy = max.y + marginb;
01869 BigReal maxz = max.z + marginc;
01870
01871 int xdev, ydev, zdev;
01872 int problemCount = 0;
01873
01874 FullAtomList::iterator p_i = atom.begin();
01875 FullAtomList::iterator p_e = atom.end();
01876 for ( ; p_i != p_e; ++p_i ) {
01877
01878 ScaledPosition s = lattice.scale(p_i->position);
01879
01880 // check if atom is within bounds
01881 if (s.x < minx) xdev = 0;
01882 else if (maxx <= s.x) xdev = 2;
01883 else xdev = 1;
01884
01885 if (s.y < miny) ydev = 0;
01886 else if (maxy <= s.y) ydev = 2;
01887 else ydev = 1;
01888
01889 if (s.z < minz) zdev = 0;
01890 else if (maxz <= s.z) zdev = 2;
01891 else zdev = 1;
01892
01893 if (mInfo[xdev][ydev][zdev]) { // somewhere else to be
01894 ++problemCount;
01895 }
01896
01897 }
01898
01899 marginViolations = problemCount;
01900 // if ( problemCount ) {
01901 // iout << iERROR <<
01902 // "Found " << problemCount << " margin violations!\n" << endi;
01903 // }
01904
01905 }
|
|
|
Definition at line 1713 of file HomePatch.C. References ResizeArray< Elem >::begin(), BigReal, 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(). 01714 {
01715 SimParameters *simParams = Node::Object()->simParameters;
01716
01717 if ( ! flags.usePairlists ) {
01718 flags.pairlistTolerance = 0.;
01719 flags.maxAtomMovement = 99999.;
01720 return;
01721 }
01722
01723 int i; int n = numAtoms;
01724 CompAtom *p_i = p.begin();
01725
01726 if ( flags.savePairlists ) {
01727 flags.pairlistTolerance = doPairlistCheck_newTolerance;
01728 flags.maxAtomMovement = 0.;
01729 doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
01730 doPairlistCheck_lattice = lattice;
01731 doPairlistCheck_positions.resize(numAtoms);
01732 CompAtom *psave_i = doPairlistCheck_positions.begin();
01733 for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
01734 return;
01735 }
01736
01737 Lattice &lattice_old = doPairlistCheck_lattice;
01738 Position center_cur = lattice.unscale(center);
01739 Position center_old = lattice_old.unscale(center);
01740 Vector center_delta = center_cur - center_old;
01741
01742 // find max deviation to corner (any neighbor shares a corner)
01743 BigReal max_cd = 0.;
01744 for ( i=0; i<2; ++i ) {
01745 for ( int j=0; j<2; ++j ) {
01746 for ( int k=0; k<2; ++k ) {
01747 ScaledPosition corner( i ? min.x : max.x ,
01748 j ? min.y : max.y ,
01749 k ? min.z : max.z );
01750 Vector corner_delta =
01751 lattice.unscale(corner) - lattice_old.unscale(corner);
01752 corner_delta -= center_delta;
01753 BigReal cd = corner_delta.length2();
01754 if ( cd > max_cd ) max_cd = cd;
01755 }
01756 }
01757 }
01758 max_cd = sqrt(max_cd);
01759
01760 // find max deviation of atoms relative to center
01761 BigReal max_pd = 0.;
01762 CompAtom *p_old_i = doPairlistCheck_positions.begin();
01763 for ( i=0; i<n; ++i ) {
01764 Vector p_delta = p_i[i].position - p_old_i[i].position;
01765 p_delta -= center_delta;
01766 BigReal pd = p_delta.length2();
01767 if ( pd > max_pd ) max_pd = pd;
01768 }
01769 max_pd = sqrt(max_pd);
01770
01771 BigReal max_tol = max_pd + max_cd;
01772
01773 flags.maxAtomMovement = max_tol;
01774
01775 // if ( max_tol > flags.pairlistTolerance ) iout << "tolerance " << max_tol << " > " << flags.pairlistTolerance << "\n" << endi;
01776
01777 if ( max_tol > ( (1. - simParams->pairlistTrigger) *
01778 doPairlistCheck_newTolerance ) ) {
01779 doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
01780 }
01781
01782 if ( max_tol > doPairlistCheck_newTolerance ) {
01783 doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
01784 }
01785
01786 }
|
|
|
Definition at line 120 of file HomePatch.h. References FullAtomList. Referenced by ComputeHomePatch::doWork(), and dumpbench(). 00120 { return (atom); }
|
|
|
Definition at line 1549 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(). 01550 {
01551 Molecule *mol = Node::Object()->molecule;
01552 SimParameters *simParams = Node::Object()->simParameters;
01553 BigReal tol = simParams->mollyTol;
01554 int maxiter = simParams->mollyIter;
01555 int i, iter;
01556 HGArrayBigReal dsq;
01557 BigReal tmp;
01558 HGArrayInt ial, ibl;
01559 HGArrayVector ref; // reference position
01560 HGArrayVector refab; // reference vector
01561 HGArrayBigReal rmass; // 1 / mass
01562 BigReal *lambda; // Lagrange multipliers
01563 CompAtom *avg; // averaged position
01564 int numLambdas = 0;
01565 HGArrayInt fixed; // is atom fixed?
01566
01567 // iout<<iINFO << "mollyAverage: "<<std::endl<<endi;
01568 p_avg.resize(numAtoms);
01569 for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
01570
01571 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01572 int hgs = atom[ig].hydrogenGroupSize;
01573 if ( hgs == 1 ) continue; // only one atom in group
01574 for ( i = 0; i < hgs; ++i ) {
01575 ref[i] = atom[ig+i].position;
01576 rmass[i] = 1. / atom[ig+i].mass;
01577 fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
01578 if ( fixed[i] ) rmass[i] = 0.;
01579 }
01580 avg = &(p_avg[ig]);
01581 int icnt = 0;
01582
01583 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) ) { // for water
01584 if ( hgs != 3 ) {
01585 NAMD_die("Hydrogen group error caught in mollyAverage(). It's a bug!\n");
01586 }
01587 if ( !(fixed[1] && fixed[2]) ) {
01588 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01589 }
01590 }
01591 for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
01592 if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) ) {
01593 if ( !(fixed[0] && fixed[i]) ) {
01594 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01595 }
01596 }
01597 }
01598 if ( icnt == 0 ) continue; // no constraints
01599 numLambdas += icnt;
01600 molly_lambda.resize(numLambdas);
01601 lambda = &(molly_lambda[numLambdas - icnt]);
01602 for ( i = 0; i < icnt; ++i ) {
01603 refab[i] = ref[ial[i]] - ref[ibl[i]];
01604 }
01605 // iout<<iINFO<<"hgs="<<hgs<<" m="<<icnt<<std::endl<<endi;
01606 iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
01607 if ( iter == maxiter ) {
01608 iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
01609 }
01610 }
01611
01612 // for ( i=0; i<numAtoms; ++i ) {
01613 // if ( ( p_avg[i].position - p[i].position ).length2() > 1.0 ) {
01614 // iout << iERROR << "MOLLY moved atom " << (p[i].id + 1) << " from "
01615 // << p[i].position << " to " << p_avg[i].position << "\n" << endi;
01616 // }
01617 // }
01618
01619 }
|
|
|
Definition at line 1623 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(). 01624 {
01625 Molecule *mol = Node::Object()->molecule;
01626 SimParameters *simParams = Node::Object()->simParameters;
01627 Tensor wc; // constraint virial
01628 int i;
01629 HGArrayInt ial, ibl;
01630 HGArrayVector ref; // reference position
01631 CompAtom *avg; // averaged position
01632 HGArrayVector refab; // reference vector
01633 HGArrayVector force; // new force
01634 HGArrayBigReal rmass; // 1 / mass
01635 BigReal *lambda; // Lagrange multipliers
01636 int numLambdas = 0;
01637 HGArrayInt fixed; // is atom fixed?
01638
01639 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01640 int hgs = atom[ig].hydrogenGroupSize;
01641 if (hgs == 1 ) continue; // only one atom in group
01642 for ( i = 0; i < hgs; ++i ) {
01643 ref[i] = atom[ig+i].position;
01644 force[i] = f[Results::slow][ig+i];
01645 rmass[i] = 1. / atom[ig+i].mass;
01646 fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
01647 if ( fixed[i] ) rmass[i] = 0.;
01648 }
01649 int icnt = 0;
01650 // c-ji I'm only going to mollify water for now
01651 if ( ( mol->rigid_bond_length(atom[ig].id) ) ) { // for water
01652 if ( hgs != 3 ) {
01653 NAMD_die("Hydrogen group error caught in mollyMollify(). It's a bug!\n");
01654 }
01655 if ( !(fixed[1] && fixed[2]) ) {
01656 ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01657 }
01658 }
01659 for ( i = 1; i < hgs; ++i ) { // normal bonds to mother atom
01660 if ( ( mol->rigid_bond_length(atom[ig+i].id) ) ) {
01661 if ( !(fixed[0] && fixed[i]) ) {
01662 ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01663 }
01664 }
01665 }
01666
01667 if ( icnt == 0 ) continue; // no constraints
01668 lambda = &(molly_lambda[numLambdas]);
01669 numLambdas += icnt;
01670 for ( i = 0; i < icnt; ++i ) {
01671 refab[i] = ref[ial[i]] - ref[ibl[i]];
01672 }
01673 avg = &(p_avg[ig]);
01674 mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
01675 // store data back to patch
01676 for ( i = 0; i < hgs; ++i ) {
01677 wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
01678 f[Results::slow][ig+i] = force[i];
01679 }
01680 }
01681 // check that there isn't a constant needed here!
01682 *virial += wc;
01683 p_avg.resize(0);
01684 }
|
|
|
Reimplemented from Patch. Definition at line 835 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(). 00836 {
00837 flags.sequence++;
00838
00839 if (!patchMapRead) { readPatchMap(); }
00840
00841 if (numNeighbors) {
00842 if (doMigration) {
00843 doAtomMigration();
00844 } else {
00845 doMarginCheck();
00846 }
00847 }
00848 doMigration = (doMigration && numNeighbors) || ! patchMapRead;
00849
00850 // Workaround for oversize groups
00851 doGroupSizeCheck();
00852
00853 // Copy information needed by computes and proxys to Patch::p.
00854 p.resize(numAtoms);
00855 CompAtom *p_i = p.begin();
00856 #ifdef MEM_OPT_VERSION
00857 pExt.resize(numAtoms);
00858 CompAtomExt *pExt_i = pExt.begin();
00859 #endif
00860 FullAtom *a_i = atom.begin();
00861 int i; int n = numAtoms;
00862 for ( i=0; i<n; ++i ) {
00863 p_i[i] = a_i[i];
00864 #ifdef MEM_OPT_VERSION
00865 pExt_i[i] = a_i[i];
00866 #endif
00867 }
00868
00869 // Measure atom movement to test pairlist validity
00870 doPairlistCheck();
00871
00872 if (flags.doMolly) mollyAverage();
00873
00874 // Must Add Proxy Changes when migration completed!
00875 ProxyListIter pli(proxy);
00876 int *pids = NULL;
00877 int npid;
00878 if (proxySendSpanning == 0) {
00879 npid = proxy.size();
00880 pids = new int[npid];
00881 int *pidi = pids;
00882 int *pide = pids + proxy.size();
00883 int patchNodesLast =
00884 ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00885 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00886 {
00887 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
00888 *(--pide) = pli->node;
00889 } else {
00890 *(pidi++) = pli->node;
00891 }
00892 }
00893 }
00894 else {
00895 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00896 npid = numChild;
00897 if(numChild>0) {
00898 pids = new int[numChild];
00899 memcpy(pids, children, sizeof(int)*numChild);
00900 }
00901 #else
00902 npid = nChild;
00903 pids = new int[proxySpanDim];
00904 for (int i=0; i<nChild; i++) pids[i] = child[i];
00905 #endif
00906 }
00907 if (npid) {
00908 #if CMK_PERSISTENT_COMM
00909 if (phsReady == 0)
00910 {
00911 //CmiPrintf("Build on %d phs0:%d\n", CkMyPe(), localphs[0]);
00912 for (int i=0; i<npid; i++) {
00913 localphs[i] = CmiCreatePersistent(pids[i], 30000);
00914 }
00915 nphs = npid;
00916 phsReady = 1;
00917 }
00918 #endif
00919 int seq = flags.sequence;
00920 int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
00921 //begin to prepare proxy msg and send it
00922 int pdMsgPLLen = p.size();
00923 int pdMsgAvgPLLen = 0;
00924 if(flags.doMolly) {
00925 pdMsgAvgPLLen = p_avg.size();
00926 }
00927 int pdMsgPLExtLen = 0;
00928 #ifdef MEM_OPT_VERSION
00929 if(doMigration || isNewProxyAdded) {
00930 pdMsgPLExtLen = pExt.size();
00931 }
00932 #endif
00933 ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgPLExtLen, PRIORITY_SIZE) ProxyDataMsg;
00934 SET_PRIORITY(nmsg,seq,priority);
00935 nmsg->patch = patchID;
00936 nmsg->flags = flags;
00937 nmsg->plLen = pdMsgPLLen;
00938 //copying data to the newly created msg
00939 memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
00940 nmsg->avgPlLen = pdMsgAvgPLLen;
00941 if(flags.doMolly) {
00942 memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
00943 }
00944 nmsg->plExtLen = pdMsgPLExtLen;
00945 #ifdef MEM_OPT_VERSION
00946 if(doMigration || isNewProxyAdded){
00947 memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
00948 }
00949 #endif
00950
00951 #if NAMD_SeparateWaters != 0
00952 //DMK - Atom Separation (water vs. non-water)
00953 nmsg->numWaterAtoms = numWaterAtoms;
00954 #endif
00955
00956
00957 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00958 DebugFileTrace *dft = DebugFileTrace::Object();
00959 dft->openTrace();
00960 dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
00961 for(int i=0; i<npid; i++) {
00962 dft->writeTrace("%d ", pids[i]);
00963 }
00964 dft->writeTrace("\n");
00965 dft->closeTrace();
00966 #endif
00967
00968 if(doMigration) {
00969 ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
00970 }else{
00971 ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
00972 }
00973 #if CMK_PERSISTENT_COMM
00974 CmiUsePersistentHandle(NULL, 0);
00975 #endif
00976 #ifdef MEM_OPT_VERSION
00977 isNewProxyAdded = 0;
00978 #endif
00979 }
00980 delete [] pids;
00981 DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
00982
00983 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
00984 positionPtrBegin = p.begin();
00985 positionPtrEnd = p.end();
00986 #endif
00987
00988 if(flags.doMolly) {
00989 avgPositionPtrBegin = p_avg.begin();
00990 avgPositionPtrEnd = p_avg.end();
00991 }
00992 Patch::positionsReady(doMigration);
00993
00994 patchMapRead = 1;
00995
00996 // gzheng
00997 if (useSync) Sync::Object()->PatchReady();
00998 }
|
|
||||||||||||||||
|
Definition at line 1210 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(). 01212 {
01213 Molecule *mol = Node::Object()->molecule;
01214 SimParameters *simParams = Node::Object()->simParameters;
01215 const int fixedAtomsOn = simParams->fixedAtomsOn;
01216 const int useSettle = simParams->useSettle;
01217 const BigReal dt = timestep / TIMEFACTOR;
01218 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt; // precalc 1/dt
01219 BigReal tol2 = 2.0 * simParams->rigidTol;
01220 int maxiter = simParams->rigidIter;
01221 int dieOnError = simParams->rigidDie;
01222 int i, iter;
01223 BigReal dsq[10], tmp;
01224 int ial[10], ibl[10];
01225 Vector ref[10]; // reference position
01226 Vector refab[10]; // reference vector
01227 Vector pos[10]; // new position
01228 Vector vel[10]; // new velocity
01229 Vector netdp[10]; // total momentum change from constraint
01230 BigReal rmass[10]; // 1 / mass
01231 int fixed[10]; // is atom fixed?
01232 Tensor wc; // constraint virial
01233 BigReal idz, zmin;
01234 int nslabs;
01235
01236 // Initialize the settle algorithm with water parameters
01237 // settle1() assumes all waters are identical,
01238 // and will generate bad results if they are not.
01239 // XXX this will move to Molecule::build_atom_status when that
01240 // version is debugged
01241 if (!settle1isinitted()) {
01242 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01243 // find a water
01244 if (mol->rigid_bond_length(atom[ig].id) > 0) {
01245 int oatm = ig+1;
01246
01247 // Avoid using the Om site to set this by mistake
01248 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
01249 oatm += 1;
01250 }
01251
01252 // initialize settle water parameters
01253 settle1init(atom[ig].mass, atom[oatm].mass,
01254 mol->rigid_bond_length(atom[ig].id),
01255 mol->rigid_bond_length(atom[oatm].id));
01256 break; // done with init
01257 }
01258 }
01259 }
01260
01261 if (ppreduction) {
01262 nslabs = simParams->pressureProfileSlabs;
01263 idz = nslabs/lattice.c().z;
01264 zmin = lattice.origin().z - 0.5*lattice.c().z;
01265 }
01266
01267 // Size of a hydrogen group for water
01268 int wathgsize = 3;
01269 int watmodel = simParams->watmodel;
01270 if (watmodel == WAT_TIP4) wathgsize = 4;
01271
01272 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01273 int hgs = atom[ig].hydrogenGroupSize;
01274 if ( hgs == 1 ) continue; // only one atom in group
01275 // cache data in local arrays and integrate positions normally
01276 for ( i = 0; i < hgs; ++i ) {
01277 ref[i] = atom[ig+i].position;
01278 pos[i] = atom[ig+i].position;
01279 vel[i] = atom[ig+i].velocity;
01280 rmass[i] = (atom[ig+1].mass > 0. ? 1. / atom[ig+i].mass : 0.);
01281 //printf("rmass of %i is %f\n", ig+i, rmass[i]);
01282 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
01283 // undo addVelocityToPosition to get proper reference coordinates
01284 if ( fixed[i] ) rmass[i] = 0.; else pos[i] += vel[i] * dt;
01285 }
01286 int icnt = 0;
01287 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) > 0 ) { // for water
01288 if (hgs != wathgsize) {
01289 NAMD_bug("Hydrogen group error caught in rattle1().");
01290 }
01291 // Use SETTLE for water unless some of the water atoms are fixed,
01292 // for speed we test groupFixed rather than the individual atoms
01293 if (useSettle && !atom[ig].groupFixed) {
01294 settle1(ref, pos, vel, invdt);
01295 if (simParams->watmodel == WAT_TIP4) tip4_omrepos(ref, pos, vel, invdt);
01296
01297 // which slab the hydrogen group will belong to
01298 // for pprofile calculations.
01299 int ppoffset, partition;
01300 if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
01301 atom[ig+i].position = pos[i];
01302 } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
01303 atom[ig+i].velocity = vel[i];
01304 } |