--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/06/08 16:05:07 1576 +++ branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/06/17 20:16:35 1584 @@ -57,14 +57,16 @@ namespace OpenMD { storageLayout_ = sman_->getStorageLayout(); ff_ = info_->getForceField(); nLocal_ = snap_->getNumberOfAtoms(); - nGroups_ = snap_->getNumberOfCutoffGroups(); + nGroups_ = info_->getNLocalCutoffGroups(); + cerr << "in dId, nGroups = " << nGroups_ << "\n"; // gather the information for atomtype IDs (atids): - identsLocal = info_->getIdentArray(); + idents = info_->getIdentArray(); AtomLocalToGlobal = info_->getGlobalAtomIndices(); cgLocalToGlobal = info_->getGlobalGroupIndices(); vector globalGroupMembership = info_->getGlobalGroupMembership(); - vector massFactorsLocal = info_->getMassFactors(); + massFactors = info_->getMassFactors(); + PairList excludes = info_->getExcludedInteractions(); PairList oneTwo = info_->getOneTwoInteractions(); PairList oneThree = info_->getOneThreeInteractions(); @@ -104,11 +106,11 @@ 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); + AtomCommIntRow->gather(idents, identsRow); + AtomCommIntColumn->gather(idents, identsCol); AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); @@ -116,11 +118,11 @@ namespace OpenMD { cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); - AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); - AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); + AtomCommRealRow->gather(massFactors, massFactorsRow); + AtomCommRealColumn->gather(massFactors, 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 +133,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++) { @@ -141,92 +143,86 @@ namespace OpenMD { } } - skipsForRowAtom.clear(); - skipsForRowAtom.reserve(nAtomsInRow_); + skipsForAtom.clear(); + skipsForAtom.resize(nAtomsInRow_); + toposForAtom.clear(); + toposForAtom.resize(nAtomsInRow_); + topoDist.clear(); + topoDist.resize(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]; + int jglob = AtomColToGlobal[j]; + + if (excludes.hasPair(iglob, jglob)) + skipsForAtom[i].push_back(j); + if (oneTwo.hasPair(iglob, jglob)) { - toposForRowAtom[i].push_back(j); - topoDistRow[i][nTopos] = 1; - nTopos++; + toposForAtom[i].push_back(j); + topoDist[i].push_back(1); + } else { + if (oneThree.hasPair(iglob, jglob)) { + toposForAtom[i].push_back(j); + topoDist[i].push_back(2); + } else { + if (oneFour.hasPair(iglob, jglob)) { + toposForAtom[i].push_back(j); + topoDist[i].push_back(3); + } + } } - 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_); + 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_); + skipsForAtom.clear(); + skipsForAtom.resize(nLocal_); + toposForAtom.clear(); + toposForAtom.resize(nLocal_); + topoDist.clear(); + topoDist.resize(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); - } - } + int jglob = AtomLocalToGlobal[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 (excludes.hasPair(iglob, jglob)) + skipsForAtom[i].push_back(j); + if (oneTwo.hasPair(iglob, jglob)) { - toposForLocalAtom[i].push_back(j); - topoDistLocal[i][nTopos] = 1; - nTopos++; + toposForAtom[i].push_back(j); + topoDist[i].push_back(1); + } else { + if (oneThree.hasPair(iglob, jglob)) { + toposForAtom[i].push_back(j); + topoDist[i].push_back(2); + } else { + if (oneFour.hasPair(iglob, jglob)) { + toposForAtom[i].push_back(j); + topoDist[i].push_back(3); + } + } } - 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++; - } } - } - + } + + createGtypeCutoffMap(); } void ForceMatrixDecomposition::createGtypeCutoffMap() { @@ -236,12 +232,16 @@ namespace OpenMD { 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); + atypeCutoff.resize( atypes.size() ); + + for (set::iterator at = atypes.begin(); + at != atypes.end(); ++at){ atid = (*at)->getIdent(); - atypeCutoff[atid] = rc; + + if (userChoseCutoff_) + atypeCutoff[atid] = userCutoff_; + else + atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at); } vector gTypeCutoffs; @@ -250,6 +250,7 @@ namespace OpenMD { // largest cutoff for any atypes present in this group. #ifdef IS_MPI vector groupCutoffRow(nGroupsInRow_, 0.0); + groupRowToGtype.resize(nGroupsInRow_); for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) { vector atomListRow = getAtomsInGroupRow(cg1); for (vector::iterator ia = atomListRow.begin(); @@ -275,6 +276,7 @@ namespace OpenMD { } vector groupCutoffCol(nGroupsInCol_, 0.0); + groupColToGtype.resize(nGroupsInCol_); for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) { vector atomListCol = getAtomsInGroupColumn(cg2); for (vector::iterator jb = atomListCol.begin(); @@ -298,14 +300,20 @@ namespace OpenMD { } } #else + vector groupCutoff(nGroups_, 0.0); + groupToGtype.resize(nGroups_); + + cerr << "nGroups = " << nGroups_ << "\n"; 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]; + atid = idents[atom1]; if (atypeCutoff[atid] > groupCutoff[cg1]) { groupCutoff[cg1] = atypeCutoff[atid]; } @@ -325,10 +333,10 @@ namespace OpenMD { } #endif + cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n"; // Now we find the maximum group cutoff value present in the simulation - vector::iterator groupMaxLoc = max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); - RealType groupMax = *groupMaxLoc; + RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end()); #ifdef IS_MPI MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX); @@ -337,23 +345,26 @@ namespace OpenMD { RealType tradRcut = groupMax; for (int i = 0; i < gTypeCutoffs.size(); i++) { - for (int j = 0; j < gTypeCutoffs.size(); j++) { - + for (int j = 0; j < gTypeCutoffs.size(); j++) { RealType thisRcut; switch(cutoffPolicy_) { case TRADITIONAL: thisRcut = tradRcut; + break; case MIX: thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]); + break; case MAX: thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]); + break; default: sprintf(painCave.errMsg, "ForceMatrixDecomposition::createGtypeCutoffMap " "hit an unknown cutoff policy!\n"); painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; - simError(); + simError(); + break; } pair key = make_pair(i,j); @@ -371,7 +382,7 @@ namespace OpenMD { if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { sprintf(painCave.errMsg, "ForceMatrixDecomposition::createGtypeCutoffMap " - "user-specified rCut does not match computed group Cutoff\n"); + "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_); painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); @@ -383,26 +394,29 @@ namespace OpenMD { groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { - int i, j; - + int i, j; #ifdef IS_MPI i = groupRowToGtype[cg1]; j = groupColToGtype[cg2]; #else i = groupToGtype[cg1]; j = groupToGtype[cg2]; -#endif - +#endif return gTypeCutoffMap[make_pair(i,j)]; } + int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { + for (int j = 0; j < toposForAtom[atom1].size(); j++) { + if (toposForAtom[atom1][j] == atom2) + return topoDist[atom1][j]; + } + return 0; + } void ForceMatrixDecomposition::zeroWorkArrays() { + pairwisePot = 0.0; + embeddingPot = 0.0; - 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); @@ -418,9 +432,7 @@ namespace OpenMD { Vector (0.0)); fill(pot_col.begin(), pot_col.end(), - Vector (0.0)); - - pot_local = Vector(0.0); + Vector (0.0)); if (storageLayout_ & DataStorage::dslParticlePot) { fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); @@ -594,7 +606,7 @@ namespace OpenMD { AtomCommPotRow->scatter(pot_row, pot_temp); for (int ii = 0; ii < pot_temp.size(); ii++ ) - pot_local += pot_temp[ii]; + pairwisePot += pot_temp[ii]; fill(pot_temp.begin(), pot_temp.end(), Vector (0.0)); @@ -602,9 +614,9 @@ namespace OpenMD { AtomCommPotColumn->scatter(pot_col, pot_temp); for (int ii = 0; ii < pot_temp.size(); ii++ ) - pot_local += pot_temp[ii]; - + pairwisePot += pot_temp[ii]; #endif + } int ForceMatrixDecomposition::getNAtomsInRow() { @@ -679,7 +691,8 @@ namespace OpenMD { #ifdef IS_MPI return massFactorsRow[atom1]; #else - return massFactorsLocal[atom1]; + cerr << "mfs = " << massFactors.size() << " atom1 = " << atom1 << "\n"; + return massFactors[atom1]; #endif } @@ -687,7 +700,7 @@ namespace OpenMD { #ifdef IS_MPI return massFactorsCol[atom2]; #else - return massFactorsLocal[atom2]; + return massFactors[atom2]; #endif } @@ -705,12 +718,8 @@ namespace OpenMD { return d; } - vector ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) { -#ifdef IS_MPI - return skipsForRowAtom[atom1]; -#else - return skipsForLocalAtom[atom1]; -#endif + vector ForceMatrixDecomposition::getSkipsForAtom(int atom1) { + return skipsForAtom[atom1]; } /** @@ -743,35 +752,15 @@ namespace OpenMD { unique_id_2 = atom2; #endif -#ifdef IS_MPI - for (vector::iterator i = skipsForRowAtom[atom1].begin(); - i != skipsForRowAtom[atom1].end(); ++i) { + for (vector::iterator i = skipsForAtom[atom1].begin(); + i != skipsForAtom[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; + return false; } + void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){ #ifdef IS_MPI atomRowData.force[atom1] += fg; @@ -789,14 +778,12 @@ namespace OpenMD { } // filling interaction blocks with pointers - InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) { - InteractionData idat; - + void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, + int atom1, int atom2) { #ifdef IS_MPI idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), ff_->getAtomType(identsCol[atom2]) ); - if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(atomRowData.aMat[atom1]); @@ -835,8 +822,8 @@ namespace OpenMD { #else - idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), - ff_->getAtomType(identsLocal[atom2]) ); + idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), + ff_->getAtomType(idents[atom2]) ); if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(snap_->atomData.aMat[atom1]); @@ -853,7 +840,7 @@ namespace OpenMD { idat.t2 = &(snap_->atomData.torque[atom2]); } - if (storageLayout_ & DataStorage::dslDensity) { + if (storageLayout_ & DataStorage::dslDensity) { idat.rho1 = &(snap_->atomData.density[atom1]); idat.rho2 = &(snap_->atomData.density[atom2]); } @@ -874,11 +861,10 @@ namespace OpenMD { } #endif - return idat; } - void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) { + 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); @@ -886,8 +872,8 @@ namespace OpenMD { atomRowData.force[atom1] += *(idat.f1); atomColData.force[atom2] -= *(idat.f1); #else - longRangePot_ += *(idat.pot); - + pairwisePot += *(idat.pot); + snap_->atomData.force[atom1] += *(idat.f1); snap_->atomData.force[atom2] -= *(idat.f1); #endif @@ -895,9 +881,8 @@ namespace OpenMD { } - InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ - - InteractionData idat; + void ForceMatrixDecomposition::fillSkipData(InteractionData &idat, + int atom1, int atom2) { #ifdef IS_MPI idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), ff_->getAtomType(identsCol[atom2]) ); @@ -906,25 +891,49 @@ namespace OpenMD { idat.eFrame1 = &(atomRowData.electroFrame[atom1]); idat.eFrame2 = &(atomColData.electroFrame[atom2]); } + if (storageLayout_ & DataStorage::dslTorque) { idat.t1 = &(atomRowData.torque[atom1]); idat.t2 = &(atomColData.torque[atom2]); } + + if (storageLayout_ & DataStorage::dslSkippedCharge) { + idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); + idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); + } #else - idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), - ff_->getAtomType(identsLocal[atom2]) ); + idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), + ff_->getAtomType(idents[atom2]) ); if (storageLayout_ & DataStorage::dslElectroFrame) { idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); } + if (storageLayout_ & DataStorage::dslTorque) { idat.t1 = &(snap_->atomData.torque[atom1]); idat.t2 = &(snap_->atomData.torque[atom2]); } + + if (storageLayout_ & DataStorage::dslSkippedCharge) { + idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); + idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); + } #endif } + + void ForceMatrixDecomposition::unpackSkipData(InteractionData &idat, int atom1, int atom2) { +#ifdef IS_MPI + pot_row[atom1] += 0.5 * *(idat.pot); + pot_col[atom2] += 0.5 * *(idat.pot); +#else + pairwisePot += *(idat.pot); +#endif + + } + + /* * buildNeighborList * @@ -958,16 +967,28 @@ namespace OpenMD { Vector3d rs, scaled, dr; Vector3i whichCell; int cellIndex; + int nCtot = nCells_.x() * nCells_.y() * nCells_.z(); #ifdef IS_MPI + cellListRow_.resize(nCtot); + cellListCol_.resize(nCtot); +#else + cellList_.resize(nCtot); +#endif + +#ifdef IS_MPI for (int i = 0; i < nGroupsInRow_; i++) { rs = cgRowData.position[i]; + // scaled positions relative to the box vectors scaled = invHmat * rs; + // wrap the vector back into the unit box by subtracting integer box // numbers - for (int j = 0; j < 3; j++) + for (int j = 0; j < 3; j++) { scaled[j] -= roundMe(scaled[j]); + scaled[j] += 0.5; + } // find xyz-indices of cell that cutoffGroup is in. whichCell.x() = nCells_.x() * scaled.x(); @@ -976,18 +997,23 @@ namespace OpenMD { // find single index of this cell: cellIndex = Vlinear(whichCell, nCells_); + // add this cutoff group to the list of groups in this cell; cellListRow_[cellIndex].push_back(i); } for (int i = 0; i < nGroupsInCol_; i++) { rs = cgColData.position[i]; + // scaled positions relative to the box vectors scaled = invHmat * rs; + // wrap the vector back into the unit box by subtracting integer box // numbers - for (int j = 0; j < 3; j++) + for (int j = 0; j < 3; j++) { scaled[j] -= roundMe(scaled[j]); + scaled[j] += 0.5; + } // find xyz-indices of cell that cutoffGroup is in. whichCell.x() = nCells_.x() * scaled.x(); @@ -996,18 +1022,23 @@ namespace OpenMD { // find single index of this cell: cellIndex = Vlinear(whichCell, nCells_); + // add this cutoff group to the list of groups in this cell; cellListCol_[cellIndex].push_back(i); } #else for (int i = 0; i < nGroups_; i++) { rs = snap_->cgData.position[i]; + // scaled positions relative to the box vectors scaled = invHmat * rs; + // wrap the vector back into the unit box by subtracting integer box // numbers - for (int j = 0; j < 3; j++) + for (int j = 0; j < 3; j++) { scaled[j] -= roundMe(scaled[j]); + scaled[j] += 0.5; + } // find xyz-indices of cell that cutoffGroup is in. whichCell.x() = nCells_.x() * scaled.x(); @@ -1015,7 +1046,8 @@ namespace OpenMD { whichCell.z() = nCells_.z() * scaled.z(); // find single index of this cell: - cellIndex = Vlinear(whichCell, nCells_); + cellIndex = Vlinear(whichCell, nCells_); + // add this cutoff group to the list of groups in this cell; cellList_[cellIndex].push_back(i); } @@ -1073,11 +1105,12 @@ namespace OpenMD { } } #else + for (vector::iterator j1 = cellList_[m1].begin(); j1 != cellList_[m1].end(); ++j1) { for (vector::iterator j2 = cellList_[m2].begin(); j2 != cellList_[m2].end(); ++j2) { - + // Always do this if we're in different cells or if // we're in the same cell and the global index of the // j2 cutoff group is less than the j1 cutoff group @@ -1097,13 +1130,13 @@ namespace OpenMD { } } } - + // save the local cutoff group positions for the check that is // done on each loop: saved_CG_positions_.clear(); for (int i = 0; i < nGroups_; i++) saved_CG_positions_.push_back(snap_->cgData.position[i]); - + return neighborList; } } //end namespace OpenMD