--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/05/27 16:45:44 1571 +++ branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/06/08 16:05:07 1576 @@ -69,7 +69,6 @@ namespace OpenMD { PairList oneTwo = info_->getOneTwoInteractions(); PairList oneThree = info_->getOneThreeInteractions(); PairList oneFour = info_->getOneFourInteractions(); - vector pot_local(N_INTERACTION_FAMILIES, 0.0); #ifdef IS_MPI @@ -77,11 +76,13 @@ namespace OpenMD { 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_); @@ -102,12 +103,7 @@ 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)); - + identsRow.reserve(nAtomsInRow_); identsCol.reserve(nAtomsInCol_); @@ -229,9 +225,249 @@ namespace OpenMD { 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(); @@ -267,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(); @@ -278,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(); @@ -343,15 +586,24 @@ 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 } @@ -462,11 +714,11 @@ namespace OpenMD { } /** - * 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. + * 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; @@ -545,6 +797,7 @@ namespace OpenMD { 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]); @@ -565,11 +818,21 @@ namespace OpenMD { 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]), @@ -595,14 +858,43 @@ 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; @@ -618,10 +910,6 @@ 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]) ); @@ -634,10 +922,6 @@ 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 } @@ -650,6 +934,7 @@ namespace OpenMD { vector > ForceMatrixDecomposition::buildNeighborList() { vector > neighborList; + groupCutoffs cuts; #ifdef IS_MPI cellListRow_.clear(); cellListCol_.clear(); @@ -657,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(); @@ -739,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++) { @@ -785,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))); } } @@ -804,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))); } }