--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/05/25 16:20:37 1568 +++ branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/06/08 16:05:07 1576 @@ -42,6 +42,7 @@ #include "math/SquareMatrix3.hpp" #include "nonbonded/NonBondedInteraction.hpp" #include "brains/SnapshotManager.hpp" +#include "brains/PairList.hpp" using namespace std; namespace OpenMD { @@ -54,20 +55,34 @@ namespace OpenMD { void ForceMatrixDecomposition::distributeInitialData() { snap_ = sman_->getCurrentSnapshot(); storageLayout_ = sman_->getStorageLayout(); + ff_ = info_->getForceField(); nLocal_ = snap_->getNumberOfAtoms(); nGroups_ = snap_->getNumberOfCutoffGroups(); + // gather the information for atomtype IDs (atids): + identsLocal = info_->getIdentArray(); + AtomLocalToGlobal = info_->getGlobalAtomIndices(); + cgLocalToGlobal = info_->getGlobalGroupIndices(); + vector globalGroupMembership = info_->getGlobalGroupMembership(); + vector massFactorsLocal = info_->getMassFactors(); + PairList excludes = info_->getExcludedInteractions(); + PairList oneTwo = info_->getOneTwoInteractions(); + PairList oneThree = info_->getOneThreeInteractions(); + PairList oneFour = info_->getOneFourInteractions(); + #ifdef IS_MPI AtomCommIntRow = new Communicator(nLocal_); AtomCommRealRow = new Communicator(nLocal_); AtomCommVectorRow = new Communicator(nLocal_); AtomCommMatrixRow = new Communicator(nLocal_); + AtomCommPotRow = new Communicator(nLocal_); AtomCommIntColumn = new Communicator(nLocal_); AtomCommRealColumn = new Communicator(nLocal_); AtomCommVectorColumn = new Communicator(nLocal_); AtomCommMatrixColumn = new Communicator(nLocal_); + AtomCommPotColumn = new Communicator(nLocal_); cgCommIntRow = new Communicator(nGroups_); cgCommVectorRow = new Communicator(nGroups_); @@ -88,39 +103,371 @@ namespace OpenMD { cgRowData.setStorageLayout(DataStorage::dslPosition); cgColData.resize(nGroupsInCol_); cgColData.setStorageLayout(DataStorage::dslPosition); - - vector > pot_row(N_INTERACTION_FAMILIES, - vector (nAtomsInRow_, 0.0)); - vector > pot_col(N_INTERACTION_FAMILIES, - vector (nAtomsInCol_, 0.0)); - - - vector pot_local(N_INTERACTION_FAMILIES, 0.0); - - // gather the information for atomtype IDs (atids): - vector identsLocal = info_->getIdentArray(); + identsRow.reserve(nAtomsInRow_); identsCol.reserve(nAtomsInCol_); AtomCommIntRow->gather(identsLocal, identsRow); AtomCommIntColumn->gather(identsLocal, identsCol); - AtomLocalToGlobal = info_->getGlobalAtomIndices(); AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); - cgLocalToGlobal = info_->getGlobalGroupIndices(); cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); - // still need: - // topoDist - // exclude + AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); + AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); + + groupListRow_.clear(); + groupListRow_.reserve(nGroupsInRow_); + for (int i = 0; i < nGroupsInRow_; i++) { + int gid = cgRowToGlobal[i]; + for (int j = 0; j < nAtomsInRow_; j++) { + int aid = AtomRowToGlobal[j]; + if (globalGroupMembership[aid] == gid) + groupListRow_[i].push_back(j); + } + } + + groupListCol_.clear(); + groupListCol_.reserve(nGroupsInCol_); + for (int i = 0; i < nGroupsInCol_; i++) { + int gid = cgColToGlobal[i]; + for (int j = 0; j < nAtomsInCol_; j++) { + int aid = AtomColToGlobal[j]; + if (globalGroupMembership[aid] == gid) + groupListCol_[i].push_back(j); + } + } + + skipsForRowAtom.clear(); + skipsForRowAtom.reserve(nAtomsInRow_); + for (int i = 0; i < nAtomsInRow_; i++) { + int iglob = AtomRowToGlobal[i]; + for (int j = 0; j < nAtomsInCol_; j++) { + int jglob = AtomColToGlobal[j]; + if (excludes.hasPair(iglob, jglob)) + skipsForRowAtom[i].push_back(j); + } + } + + toposForRowAtom.clear(); + toposForRowAtom.reserve(nAtomsInRow_); + for (int i = 0; i < nAtomsInRow_; i++) { + int iglob = AtomRowToGlobal[i]; + int nTopos = 0; + for (int j = 0; j < nAtomsInCol_; j++) { + int jglob = AtomColToGlobal[j]; + if (oneTwo.hasPair(iglob, jglob)) { + toposForRowAtom[i].push_back(j); + topoDistRow[i][nTopos] = 1; + nTopos++; + } + if (oneThree.hasPair(iglob, jglob)) { + toposForRowAtom[i].push_back(j); + topoDistRow[i][nTopos] = 2; + nTopos++; + } + if (oneFour.hasPair(iglob, jglob)) { + toposForRowAtom[i].push_back(j); + topoDistRow[i][nTopos] = 3; + nTopos++; + } + } + } + +#endif + + groupList_.clear(); + groupList_.reserve(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) + groupList_[i].push_back(j); + } + } + + skipsForLocalAtom.clear(); + skipsForLocalAtom.reserve(nLocal_); + + for (int i = 0; i < nLocal_; i++) { + int iglob = AtomLocalToGlobal[i]; + for (int j = 0; j < nLocal_; j++) { + int jglob = AtomLocalToGlobal[j]; + if (excludes.hasPair(iglob, jglob)) + skipsForLocalAtom[i].push_back(j); + } + } + + toposForLocalAtom.clear(); + toposForLocalAtom.reserve(nLocal_); + for (int i = 0; i < nLocal_; i++) { + int iglob = AtomLocalToGlobal[i]; + int nTopos = 0; + for (int j = 0; j < nLocal_; j++) { + int jglob = AtomLocalToGlobal[j]; + if (oneTwo.hasPair(iglob, jglob)) { + toposForLocalAtom[i].push_back(j); + topoDistLocal[i][nTopos] = 1; + nTopos++; + } + if (oneThree.hasPair(iglob, jglob)) { + toposForLocalAtom[i].push_back(j); + topoDistLocal[i][nTopos] = 2; + nTopos++; + } + if (oneFour.hasPair(iglob, jglob)) { + toposForLocalAtom[i].push_back(j); + topoDistLocal[i][nTopos] = 3; + nTopos++; + } + } + } + + } + + void ForceMatrixDecomposition::createGtypeCutoffMap() { + + RealType tol = 1e-6; + RealType rc; + int atid; + set atypes = info_->getSimulatedAtomTypes(); + vector atypeCutoff; + atypeCutoff.reserve( 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++) { + longRangePot_[j] = 0.0; + } + +#ifdef IS_MPI + if (storageLayout_ & DataStorage::dslForce) { + fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero); + fill(atomColData.force.begin(), atomColData.force.end(), V3Zero); + } + + if (storageLayout_ & DataStorage::dslTorque) { + fill(atomRowData.torque.begin(), atomRowData.torque.end(), V3Zero); + fill(atomColData.torque.begin(), atomColData.torque.end(), V3Zero); + } + fill(pot_row.begin(), pot_row.end(), + Vector (0.0)); + fill(pot_col.begin(), pot_col.end(), + Vector (0.0)); + + pot_local = Vector(0.0); + if (storageLayout_ & DataStorage::dslParticlePot) { + fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); + fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0); + } + + if (storageLayout_ & DataStorage::dslDensity) { + fill(atomRowData.density.begin(), atomRowData.density.end(), 0.0); + fill(atomColData.density.begin(), atomColData.density.end(), 0.0); + } + + if (storageLayout_ & DataStorage::dslFunctional) { + fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0); + fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0); + } + + if (storageLayout_ & DataStorage::dslFunctionalDerivative) { + fill(atomRowData.functionalDerivative.begin(), + atomRowData.functionalDerivative.end(), 0.0); + fill(atomColData.functionalDerivative.begin(), + atomColData.functionalDerivative.end(), 0.0); + } + +#else + + if (storageLayout_ & DataStorage::dslParticlePot) { + fill(snap_->atomData.particlePot.begin(), + snap_->atomData.particlePot.end(), 0.0); + } + + if (storageLayout_ & DataStorage::dslDensity) { + fill(snap_->atomData.density.begin(), + snap_->atomData.density.end(), 0.0); + } + if (storageLayout_ & DataStorage::dslFunctional) { + fill(snap_->atomData.functional.begin(), + snap_->atomData.functional.end(), 0.0); + } + if (storageLayout_ & DataStorage::dslFunctionalDerivative) { + fill(snap_->atomData.functionalDerivative.begin(), + snap_->atomData.functionalDerivative.end(), 0.0); + } +#endif + + } + + void ForceMatrixDecomposition::distributeData() { snap_ = sman_->getCurrentSnapshot(); storageLayout_ = sman_->getStorageLayout(); @@ -156,6 +503,9 @@ namespace OpenMD { #endif } + /* collects information obtained during the pre-pair loop onto local + * data structures. + */ void ForceMatrixDecomposition::collectIntermediateData() { snap_ = sman_->getCurrentSnapshot(); storageLayout_ = sman_->getStorageLayout(); @@ -167,14 +517,18 @@ namespace OpenMD { snap_->atomData.density); int n = snap_->atomData.density.size(); - std::vector rho_tmp(n, 0.0); + vector rho_tmp(n, 0.0); AtomCommRealColumn->scatter(atomColData.density, rho_tmp); for (int i = 0; i < n; i++) snap_->atomData.density[i] += rho_tmp[i]; } #endif } - + + /* + * redistributes information obtained during the pre-pair loop out to + * row and column-indexed data structures + */ void ForceMatrixDecomposition::distributeIntermediateData() { snap_ = sman_->getCurrentSnapshot(); storageLayout_ = sman_->getStorageLayout(); @@ -232,18 +586,53 @@ namespace OpenMD { nLocal_ = snap_->getNumberOfAtoms(); - vector > pot_temp(N_INTERACTION_FAMILIES, - vector (nLocal_, 0.0)); + vector pot_temp(nLocal_, + Vector (0.0)); + + // scatter/gather pot_row into the members of my column + + AtomCommPotRow->scatter(pot_row, pot_temp); + + for (int ii = 0; ii < pot_temp.size(); ii++ ) + pot_local += pot_temp[ii]; - for (int i = 0; i < N_INTERACTION_FAMILIES; i++) { - AtomCommRealRow->scatter(pot_row[i], pot_temp[i]); - for (int ii = 0; ii < pot_temp[i].size(); ii++ ) { - pot_local[i] += pot_temp[i][ii]; - } - } + fill(pot_temp.begin(), pot_temp.end(), + Vector (0.0)); + + AtomCommPotColumn->scatter(pot_col, pot_temp); + + for (int ii = 0; ii < pot_temp.size(); ii++ ) + pot_local += pot_temp[ii]; + #endif } + int ForceMatrixDecomposition::getNAtomsInRow() { +#ifdef IS_MPI + return nAtomsInRow_; +#else + return nLocal_; +#endif + } + + /** + * returns the list of atoms belonging to this group. + */ + vector ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ +#ifdef IS_MPI + return groupListRow_[cg1]; +#else + return groupList_[cg1]; +#endif + } + + vector ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ +#ifdef IS_MPI + return groupListCol_[cg2]; +#else + return groupList_[cg2]; +#endif + } Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){ Vector3d d; @@ -285,6 +674,23 @@ namespace OpenMD { snap_->wrapVector(d); return d; } + + RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) { +#ifdef IS_MPI + return massFactorsRow[atom1]; +#else + return massFactorsLocal[atom1]; +#endif + } + + RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { +#ifdef IS_MPI + return massFactorsCol[atom2]; +#else + return massFactorsLocal[atom2]; +#endif + + } Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){ Vector3d d; @@ -299,6 +705,73 @@ namespace OpenMD { return d; } + vector ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) { +#ifdef IS_MPI + return skipsForRowAtom[atom1]; +#else + return skipsForLocalAtom[atom1]; +#endif + } + + /** + * There are a number of reasons to skip a pair or a + * particle. Mostly we do this to exclude atoms who are involved in + * short range interactions (bonds, bends, torsions), but we also + * need to exclude some overcounted interactions that result from + * the parallel decomposition. + */ + bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { + int unique_id_1, unique_id_2; + +#ifdef IS_MPI + // in MPI, we have to look up the unique IDs for each atom + unique_id_1 = AtomRowToGlobal[atom1]; + unique_id_2 = AtomColToGlobal[atom2]; + + // this situation should only arise in MPI simulations + if (unique_id_1 == unique_id_2) return true; + + // this prevents us from doing the pair on multiple processors + if (unique_id_1 < unique_id_2) { + if ((unique_id_1 + unique_id_2) % 2 == 0) return true; + } else { + if ((unique_id_1 + unique_id_2) % 2 == 1) return true; + } +#else + // in the normal loop, the atom numbers are unique + unique_id_1 = atom1; + unique_id_2 = atom2; +#endif + +#ifdef IS_MPI + for (vector::iterator i = skipsForRowAtom[atom1].begin(); + i != skipsForRowAtom[atom1].end(); ++i) { + if ( (*i) == unique_id_2 ) return true; + } +#else + for (vector::iterator i = skipsForLocalAtom[atom1].begin(); + i != skipsForLocalAtom[atom1].end(); ++i) { + if ( (*i) == unique_id_2 ) return true; + } +#endif + } + + int ForceMatrixDecomposition::getTopoDistance(int atom1, int atom2) { + +#ifdef IS_MPI + for (int i = 0; i < toposForRowAtom[atom1].size(); i++) { + if ( toposForRowAtom[atom1][i] == atom2 ) return topoDistRow[atom1][i]; + } +#else + for (int i = 0; i < toposForLocalAtom[atom1].size(); i++) { + if ( toposForLocalAtom[atom1][i] == atom2 ) return topoDistLocal[atom1][i]; + } +#endif + + // zero is default for unconnected (i.e. normal) pair interactions + return 0; + } + void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){ #ifdef IS_MPI atomRowData.force[atom1] += fg; @@ -320,6 +793,11 @@ namespace OpenMD { InteractionData idat; #ifdef IS_MPI + + idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), + ff_->getAtomType(identsCol[atom2]) ); + + if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(atomRowData.aMat[atom1]); idat.A2 = &(atomColData.aMat[atom2]); @@ -338,13 +816,28 @@ namespace OpenMD { if (storageLayout_ & DataStorage::dslDensity) { idat.rho1 = &(atomRowData.density[atom1]); idat.rho2 = &(atomColData.density[atom2]); + } + + if (storageLayout_ & DataStorage::dslFunctional) { + idat.frho1 = &(atomRowData.functional[atom1]); + idat.frho2 = &(atomColData.functional[atom2]); } if (storageLayout_ & DataStorage::dslFunctionalDerivative) { idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]); idat.dfrho2 = &(atomColData.functionalDerivative[atom2]); } + + if (storageLayout_ & DataStorage::dslParticlePot) { + idat.particlePot1 = &(atomRowData.particlePot[atom1]); + idat.particlePot2 = &(atomColData.particlePot[atom2]); + } + #else + + idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), + ff_->getAtomType(identsLocal[atom2]) ); + if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(snap_->atomData.aMat[atom1]); idat.A2 = &(snap_->atomData.aMat[atom2]); @@ -365,18 +858,50 @@ namespace OpenMD { idat.rho2 = &(snap_->atomData.density[atom2]); } + if (storageLayout_ & DataStorage::dslFunctional) { + idat.frho1 = &(snap_->atomData.functional[atom1]); + idat.frho2 = &(snap_->atomData.functional[atom2]); + } + if (storageLayout_ & DataStorage::dslFunctionalDerivative) { idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]); idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]); } + + if (storageLayout_ & DataStorage::dslParticlePot) { + idat.particlePot1 = &(snap_->atomData.particlePot[atom1]); + idat.particlePot2 = &(snap_->atomData.particlePot[atom2]); + } + #endif return idat; + } + + + void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) { +#ifdef IS_MPI + pot_row[atom1] += 0.5 * *(idat.pot); + pot_col[atom2] += 0.5 * *(idat.pot); + + atomRowData.force[atom1] += *(idat.f1); + atomColData.force[atom2] -= *(idat.f1); +#else + longRangePot_ += *(idat.pot); + + snap_->atomData.force[atom1] += *(idat.f1); + snap_->atomData.force[atom2] -= *(idat.f1); +#endif + } + InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ InteractionData idat; #ifdef IS_MPI + idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), + ff_->getAtomType(identsCol[atom2]) ); + if (storageLayout_ & DataStorage::dslElectroFrame) { idat.eFrame1 = &(atomRowData.electroFrame[atom1]); idat.eFrame2 = &(atomColData.electroFrame[atom2]); @@ -385,11 +910,10 @@ namespace OpenMD { idat.t1 = &(atomRowData.torque[atom1]); idat.t2 = &(atomColData.torque[atom2]); } - if (storageLayout_ & DataStorage::dslForce) { - idat.t1 = &(atomRowData.force[atom1]); - idat.t2 = &(atomColData.force[atom2]); - } #else + idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), + ff_->getAtomType(identsLocal[atom2]) ); + if (storageLayout_ & DataStorage::dslElectroFrame) { idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); @@ -398,17 +922,9 @@ namespace OpenMD { idat.t1 = &(snap_->atomData.torque[atom1]); idat.t2 = &(snap_->atomData.torque[atom2]); } - if (storageLayout_ & DataStorage::dslForce) { - idat.t1 = &(snap_->atomData.force[atom1]); - idat.t2 = &(snap_->atomData.force[atom2]); - } -#endif - +#endif } - - - /* * buildNeighborList * @@ -418,6 +934,7 @@ namespace OpenMD { vector > ForceMatrixDecomposition::buildNeighborList() { vector > neighborList; + groupCutoffs cuts; #ifdef IS_MPI cellListRow_.clear(); cellListCol_.clear(); @@ -425,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(); @@ -507,8 +1021,6 @@ namespace OpenMD { } #endif - - for (int m1z = 0; m1z < nCells_.z(); m1z++) { for (int m1y = 0; m1y < nCells_.y(); m1y++) { for (int m1x = 0; m1x < nCells_.x(); m1x++) { @@ -553,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))); } } @@ -572,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))); } }