--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/06/03 21:39:49 1575 +++ branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/06/08 20:26:56 1577 @@ -57,8 +57,8 @@ namespace OpenMD { storageLayout_ = sman_->getStorageLayout(); ff_ = info_->getForceField(); nLocal_ = snap_->getNumberOfAtoms(); - nGroups_ = snap_->getNumberOfCutoffGroups(); + nGroups_ = info_->getNLocalCutoffGroups(); // gather the information for atomtype IDs (atids): identsLocal = info_->getIdentArray(); AtomLocalToGlobal = info_->getGlobalAtomIndices(); @@ -104,8 +104,8 @@ namespace OpenMD { cgColData.resize(nGroupsInCol_); cgColData.setStorageLayout(DataStorage::dslPosition); - identsRow.reserve(nAtomsInRow_); - identsCol.reserve(nAtomsInCol_); + identsRow.resize(nAtomsInRow_); + identsCol.resize(nAtomsInCol_); AtomCommIntRow->gather(identsLocal, identsRow); AtomCommIntColumn->gather(identsLocal, identsCol); @@ -120,7 +120,7 @@ namespace OpenMD { AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); groupListRow_.clear(); - groupListRow_.reserve(nGroupsInRow_); + groupListRow_.resize(nGroupsInRow_); for (int i = 0; i < nGroupsInRow_; i++) { int gid = cgRowToGlobal[i]; for (int j = 0; j < nAtomsInRow_; j++) { @@ -131,7 +131,7 @@ namespace OpenMD { } groupListCol_.clear(); - groupListCol_.reserve(nGroupsInCol_); + groupListCol_.resize(nGroupsInCol_); for (int i = 0; i < nGroupsInCol_; i++) { int gid = cgColToGlobal[i]; for (int j = 0; j < nAtomsInCol_; j++) { @@ -142,7 +142,7 @@ namespace OpenMD { } skipsForRowAtom.clear(); - skipsForRowAtom.reserve(nAtomsInRow_); + skipsForRowAtom.resize(nAtomsInRow_); for (int i = 0; i < nAtomsInRow_; i++) { int iglob = AtomRowToGlobal[i]; for (int j = 0; j < nAtomsInCol_; j++) { @@ -153,7 +153,7 @@ namespace OpenMD { } toposForRowAtom.clear(); - toposForRowAtom.reserve(nAtomsInRow_); + toposForRowAtom.resize(nAtomsInRow_); for (int i = 0; i < nAtomsInRow_; i++) { int iglob = AtomRowToGlobal[i]; int nTopos = 0; @@ -178,20 +178,21 @@ namespace OpenMD { } #endif - groupList_.clear(); - groupList_.reserve(nGroups_); + groupList_.resize(nGroups_); for (int i = 0; i < nGroups_; i++) { int gid = cgLocalToGlobal[i]; for (int j = 0; j < nLocal_; j++) { int aid = AtomLocalToGlobal[j]; - if (globalGroupMembership[aid] == gid) + if (globalGroupMembership[aid] == gid) { groupList_[i].push_back(j); + + } } } skipsForLocalAtom.clear(); - skipsForLocalAtom.reserve(nLocal_); + skipsForLocalAtom.resize(nLocal_); for (int i = 0; i < nLocal_; i++) { int iglob = AtomLocalToGlobal[i]; @@ -201,9 +202,8 @@ namespace OpenMD { skipsForLocalAtom[i].push_back(j); } } - toposForLocalAtom.clear(); - toposForLocalAtom.reserve(nLocal_); + toposForLocalAtom.resize(nLocal_); for (int i = 0; i < nLocal_; i++) { int iglob = AtomLocalToGlobal[i]; int nTopos = 0; @@ -225,9 +225,178 @@ namespace OpenMD { nTopos++; } } - } + } + } + void ForceMatrixDecomposition::createGtypeCutoffMap() { + + RealType tol = 1e-6; + RealType rc; + int atid; + set atypes = info_->getSimulatedAtomTypes(); + vector atypeCutoff; + atypeCutoff.resize( atypes.size() ); + + for (set::iterator at = atypes.begin(); at != atypes.end(); ++at){ + rc = interactionMan_->getSuggestedCutoffRadius(*at); + atid = (*at)->getIdent(); + atypeCutoff[atid] = rc; + } + + vector gTypeCutoffs; + + // first we do a single loop over the cutoff groups to find the + // largest cutoff for any atypes present in this group. +#ifdef IS_MPI + vector groupCutoffRow(nGroupsInRow_, 0.0); + for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) { + vector atomListRow = getAtomsInGroupRow(cg1); + for (vector::iterator ia = atomListRow.begin(); + ia != atomListRow.end(); ++ia) { + int atom1 = (*ia); + atid = identsRow[atom1]; + if (atypeCutoff[atid] > groupCutoffRow[cg1]) { + groupCutoffRow[cg1] = atypeCutoff[atid]; + } + } + + bool gTypeFound = false; + for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { + if (abs(groupCutoffRow[cg1] - gTypeCutoffs[gt]) < tol) { + groupRowToGtype[cg1] = gt; + gTypeFound = true; + } + } + if (!gTypeFound) { + gTypeCutoffs.push_back( groupCutoffRow[cg1] ); + groupRowToGtype[cg1] = gTypeCutoffs.size() - 1; + } + + } + vector groupCutoffCol(nGroupsInCol_, 0.0); + for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) { + vector atomListCol = getAtomsInGroupColumn(cg2); + for (vector::iterator jb = atomListCol.begin(); + jb != atomListCol.end(); ++jb) { + int atom2 = (*jb); + atid = identsCol[atom2]; + if (atypeCutoff[atid] > groupCutoffCol[cg2]) { + groupCutoffCol[cg2] = atypeCutoff[atid]; + } + } + bool gTypeFound = false; + for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { + if (abs(groupCutoffCol[cg2] - gTypeCutoffs[gt]) < tol) { + groupColToGtype[cg2] = gt; + gTypeFound = true; + } + } + if (!gTypeFound) { + gTypeCutoffs.push_back( groupCutoffCol[cg2] ); + groupColToGtype[cg2] = gTypeCutoffs.size() - 1; + } + } +#else + vector groupCutoff(nGroups_, 0.0); + for (int cg1 = 0; cg1 < nGroups_; cg1++) { + groupCutoff[cg1] = 0.0; + vector atomList = getAtomsInGroupRow(cg1); + for (vector::iterator ia = atomList.begin(); + ia != atomList.end(); ++ia) { + int atom1 = (*ia); + atid = identsLocal[atom1]; + if (atypeCutoff[atid] > groupCutoff[cg1]) { + groupCutoff[cg1] = atypeCutoff[atid]; + } + } + + bool gTypeFound = false; + for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { + if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) { + groupToGtype[cg1] = gt; + gTypeFound = true; + } + } + if (!gTypeFound) { + gTypeCutoffs.push_back( groupCutoff[cg1] ); + groupToGtype[cg1] = gTypeCutoffs.size() - 1; + } + } +#endif + + // Now we find the maximum group cutoff value present in the simulation + + vector::iterator groupMaxLoc = max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); + RealType groupMax = *groupMaxLoc; + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX); +#endif + + RealType tradRcut = groupMax; + + for (int i = 0; i < gTypeCutoffs.size(); i++) { + for (int j = 0; j < gTypeCutoffs.size(); j++) { + + RealType thisRcut; + switch(cutoffPolicy_) { + case TRADITIONAL: + thisRcut = tradRcut; + case MIX: + thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]); + case MAX: + thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]); + default: + sprintf(painCave.errMsg, + "ForceMatrixDecomposition::createGtypeCutoffMap " + "hit an unknown cutoff policy!\n"); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + } + + pair key = make_pair(i,j); + gTypeCutoffMap[key].first = thisRcut; + + if (thisRcut > largestRcut_) largestRcut_ = thisRcut; + + gTypeCutoffMap[key].second = thisRcut*thisRcut; + + gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); + + // sanity check + + if (userChoseCutoff_) { + if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { + sprintf(painCave.errMsg, + "ForceMatrixDecomposition::createGtypeCutoffMap " + "user-specified rCut does not match computed group Cutoff\n"); + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + } + } + } + } + } + + + groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { + int i, j; + +#ifdef IS_MPI + i = groupRowToGtype[cg1]; + j = groupColToGtype[cg2]; +#else + i = groupToGtype[cg1]; + j = groupToGtype[cg2]; +#endif + + return gTypeCutoffMap[make_pair(i,j)]; + } + + void ForceMatrixDecomposition::zeroWorkArrays() { for (int j = 0; j < N_INTERACTION_FAMILIES; j++) { @@ -765,6 +934,7 @@ namespace OpenMD { vector > ForceMatrixDecomposition::buildNeighborList() { vector > neighborList; + groupCutoffs cuts; #ifdef IS_MPI cellListRow_.clear(); cellListCol_.clear(); @@ -772,10 +942,7 @@ namespace OpenMD { cellList_.clear(); #endif - // dangerous to not do error checking. - RealType rCut_; - - RealType rList_ = (rCut_ + skinThickness_); + RealType rList_ = (largestRcut_ + skinThickness_); RealType rl2 = rList_ * rList_; Snapshot* snap_ = sman_->getCurrentSnapshot(); Mat3x3d Hmat = snap_->getHmat(); @@ -898,7 +1065,8 @@ namespace OpenMD { if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; snap_->wrapVector(dr); - if (dr.lengthSquare() < rl2) { + cuts = getGroupCutoffs( (*j1), (*j2) ); + if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair((*j1), (*j2))); } } @@ -917,7 +1085,8 @@ namespace OpenMD { if (m2 != m1 || (*j2) < (*j1)) { dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; snap_->wrapVector(dr); - if (dr.lengthSquare() < rl2) { + cuts = getGroupCutoffs( (*j1), (*j2) ); + if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair((*j1), (*j2))); } }