--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/05/26 21:56:04 1570 +++ branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/06/03 21:39:49 1575 @@ -55,11 +55,12 @@ 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): - vector identsLocal = info_->getIdentArray(); + identsLocal = info_->getIdentArray(); AtomLocalToGlobal = info_->getGlobalAtomIndices(); cgLocalToGlobal = info_->getGlobalGroupIndices(); vector globalGroupMembership = info_->getGlobalGroupMembership(); @@ -68,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 @@ -76,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_); @@ -101,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_); @@ -147,9 +144,9 @@ namespace OpenMD { skipsForRowAtom.clear(); skipsForRowAtom.reserve(nAtomsInRow_); for (int i = 0; i < nAtomsInRow_; i++) { - int iglob = AtomColToGlobal[i]; + int iglob = AtomRowToGlobal[i]; for (int j = 0; j < nAtomsInCol_; j++) { - int jglob = AtomRowToGlobal[j]; + int jglob = AtomColToGlobal[j]; if (excludes.hasPair(iglob, jglob)) skipsForRowAtom[i].push_back(j); } @@ -158,10 +155,10 @@ namespace OpenMD { toposForRowAtom.clear(); toposForRowAtom.reserve(nAtomsInRow_); for (int i = 0; i < nAtomsInRow_; i++) { - int iglob = AtomColToGlobal[i]; + int iglob = AtomRowToGlobal[i]; int nTopos = 0; for (int j = 0; j < nAtomsInCol_; j++) { - int jglob = AtomRowToGlobal[j]; + int jglob = AtomColToGlobal[j]; if (oneTwo.hasPair(iglob, jglob)) { toposForRowAtom[i].push_back(j); topoDistRow[i][nTopos] = 1; @@ -231,6 +228,77 @@ namespace OpenMD { } } + 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(); @@ -266,6 +334,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(); @@ -277,14 +348,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(); @@ -342,15 +417,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 } @@ -461,11 +545,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; @@ -540,6 +624,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]); @@ -560,12 +649,26 @@ 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]), + ff_->getAtomType(identsLocal[atom2]) ); + if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(snap_->atomData.aMat[atom1]); idat.A2 = &(snap_->atomData.aMat[atom2]); @@ -586,18 +689,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]); @@ -606,11 +741,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]); @@ -619,17 +753,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 * @@ -728,8 +854,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++) {