--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2012/04/27 20:44:16 1706 +++ trunk/src/parallel/ForceMatrixDecomposition.cpp 2014/02/28 13:25:13 1971 @@ -35,7 +35,7 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ @@ -95,10 +95,11 @@ namespace OpenMD { storageLayout_ = sman_->getStorageLayout(); ff_ = info_->getForceField(); nLocal_ = snap_->getNumberOfAtoms(); - + nGroups_ = info_->getNLocalCutoffGroups(); // gather the information for atomtype IDs (atids): idents = info_->getIdentArray(); + regions = info_->getRegions(); AtomLocalToGlobal = info_->getGlobalAtomIndices(); cgLocalToGlobal = info_->getGlobalGroupIndices(); vector globalGroupMembership = info_->getGlobalGroupMembership(); @@ -109,11 +110,17 @@ namespace OpenMD { PairList* oneTwo = info_->getOneTwoInteractions(); PairList* oneThree = info_->getOneThreeInteractions(); PairList* oneFour = info_->getOneFourInteractions(); - + + if (needVelocities_) + snap_->cgData.setStorageLayout(DataStorage::dslPosition | + DataStorage::dslVelocity); + else + snap_->cgData.setStorageLayout(DataStorage::dslPosition); + #ifdef IS_MPI - MPI::Intracomm row = rowComm.getComm(); - MPI::Intracomm col = colComm.getComm(); + MPI_Comm row = rowComm.getComm(); + MPI_Comm col = colComm.getComm(); AtomPlanIntRow = new Plan(row, nLocal_); AtomPlanRealRow = new Plan(row, nLocal_); @@ -145,14 +152,25 @@ namespace OpenMD { cgRowData.resize(nGroupsInRow_); cgRowData.setStorageLayout(DataStorage::dslPosition); cgColData.resize(nGroupsInCol_); - cgColData.setStorageLayout(DataStorage::dslPosition); - + if (needVelocities_) + // we only need column velocities if we need them. + cgColData.setStorageLayout(DataStorage::dslPosition | + DataStorage::dslVelocity); + else + cgColData.setStorageLayout(DataStorage::dslPosition); + identsRow.resize(nAtomsInRow_); identsCol.resize(nAtomsInCol_); AtomPlanIntRow->gather(idents, identsRow); AtomPlanIntColumn->gather(idents, identsCol); + + regionsRow.resize(nAtomsInRow_); + regionsCol.resize(nAtomsInCol_); + AtomPlanIntRow->gather(regions, regionsRow); + AtomPlanIntColumn->gather(regions, regionsCol); + // allocate memory for the parallel objects atypesRow.resize(nAtomsInRow_); atypesCol.resize(nAtomsInCol_); @@ -165,6 +183,9 @@ namespace OpenMD { pot_row.resize(nAtomsInRow_); pot_col.resize(nAtomsInCol_); + expot_row.resize(nAtomsInRow_); + expot_col.resize(nAtomsInCol_); + AtomRowToGlobal.resize(nAtomsInRow_); AtomColToGlobal.resize(nAtomsInCol_); AtomPlanIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); @@ -294,9 +315,12 @@ namespace OpenMD { void ForceMatrixDecomposition::createGtypeCutoffMap() { + GrCut.clear(); + GrCutSq.clear(); + GrlistSq.clear(); + RealType tol = 1e-6; largestRcut_ = 0.0; - RealType rc; int atid; set atypes = info_->getSimulatedAtomTypes(); @@ -381,7 +405,7 @@ namespace OpenMD { } bool gTypeFound = false; - for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { + for (unsigned int gt = 0; gt < gTypeCutoffs.size(); gt++) { if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) { groupToGtype[cg1] = gt; gTypeFound = true; @@ -400,14 +424,23 @@ namespace OpenMD { gTypeCutoffs.end()); #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, - MPI::MAX); + MPI_Allreduce(MPI_IN_PLACE, &groupMax, 1, MPI_REALTYPE, + MPI_MAX, MPI_COMM_WORLD); #endif RealType tradRcut = groupMax; - for (int i = 0; i < gTypeCutoffs.size(); i++) { - for (int j = 0; j < gTypeCutoffs.size(); j++) { + GrCut.resize( gTypeCutoffs.size() ); + GrCutSq.resize( gTypeCutoffs.size() ); + GrlistSq.resize( gTypeCutoffs.size() ); + + + for (unsigned int i = 0; i < gTypeCutoffs.size(); i++) { + GrCut[i].resize( gTypeCutoffs.size() , 0.0); + GrCutSq[i].resize( gTypeCutoffs.size(), 0.0 ); + GrlistSq[i].resize( gTypeCutoffs.size(), 0.0 ); + + for (unsigned int j = 0; j < gTypeCutoffs.size(); j++) { RealType thisRcut; switch(cutoffPolicy_) { case TRADITIONAL: @@ -429,15 +462,18 @@ namespace OpenMD { break; } - pair key = make_pair(i,j); - gTypeCutoffMap[key].first = thisRcut; + GrCut[i][j] = thisRcut; if (thisRcut > largestRcut_) largestRcut_ = thisRcut; - gTypeCutoffMap[key].second = thisRcut*thisRcut; - gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); + GrCutSq[i][j] = thisRcut * thisRcut; + GrlistSq[i][j] = pow(thisRcut + skinThickness_, 2); + + // pair key = make_pair(i,j); + // gTypeCutoffMap[key].first = thisRcut; + // gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); // sanity check if (userChoseCutoff_) { - if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { + if (abs(GrCut[i][j] - userCutoff_) > 0.0001) { sprintf(painCave.errMsg, "ForceMatrixDecomposition::createGtypeCutoffMap " "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_); @@ -450,8 +486,7 @@ namespace OpenMD { } } - - groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { + void ForceMatrixDecomposition::getGroupCutoffs(int &cg1, int &cg2, RealType &rcut, RealType &rcutsq, RealType &rlistsq) { int i, j; #ifdef IS_MPI i = groupRowToGtype[cg1]; @@ -460,20 +495,26 @@ namespace OpenMD { i = groupToGtype[cg1]; j = groupToGtype[cg2]; #endif - return gTypeCutoffMap[make_pair(i,j)]; + rcut = GrCut[i][j]; + rcutsq = GrCutSq[i][j]; + rlistsq = GrlistSq[i][j]; + return; + //return gTypeCutoffMap[make_pair(i,j)]; } int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { - for (int j = 0; j < toposForAtom[atom1].size(); j++) { + for (unsigned 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; + excludedPot = 0.0; + excludedSelfPot = 0.0; #ifdef IS_MPI if (storageLayout_ & DataStorage::dslForce) { @@ -492,6 +533,12 @@ namespace OpenMD { fill(pot_col.begin(), pot_col.end(), Vector (0.0)); + fill(expot_row.begin(), expot_row.end(), + Vector (0.0)); + + fill(expot_col.begin(), expot_col.end(), + Vector (0.0)); + if (storageLayout_ & DataStorage::dslParticlePot) { fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); @@ -525,6 +572,20 @@ namespace OpenMD { atomColData.skippedCharge.end(), 0.0); } + if (storageLayout_ & DataStorage::dslFlucQForce) { + fill(atomRowData.flucQFrc.begin(), + atomRowData.flucQFrc.end(), 0.0); + fill(atomColData.flucQFrc.begin(), + atomColData.flucQFrc.end(), 0.0); + } + + if (storageLayout_ & DataStorage::dslElectricField) { + fill(atomRowData.electricField.begin(), + atomRowData.electricField.end(), V3Zero); + fill(atomColData.electricField.begin(), + atomColData.electricField.end(), V3Zero); + } + #endif // even in parallel, we need to zero out the local arrays: @@ -552,6 +613,11 @@ namespace OpenMD { fill(snap_->atomData.skippedCharge.begin(), snap_->atomData.skippedCharge.end(), 0.0); } + + if (storageLayout_ & DataStorage::dslElectricField) { + fill(snap_->atomData.electricField.begin(), + snap_->atomData.electricField.end(), V3Zero); + } } @@ -574,6 +640,17 @@ namespace OpenMD { cgPlanVectorColumn->gather(snap_->cgData.position, cgColData.position); + + + if (needVelocities_) { + // gather up the atomic velocities + AtomPlanVectorColumn->gather(snap_->atomData.velocity, + atomColData.velocity); + + cgPlanVectorColumn->gather(snap_->cgData.velocity, + cgColData.velocity); + } + // if needed, gather the atomic rotation matrices if (storageLayout_ & DataStorage::dslAmat) { @@ -582,15 +659,30 @@ namespace OpenMD { AtomPlanMatrixColumn->gather(snap_->atomData.aMat, atomColData.aMat); } - - // if needed, gather the atomic eletrostatic frames - if (storageLayout_ & DataStorage::dslElectroFrame) { - AtomPlanMatrixRow->gather(snap_->atomData.electroFrame, - atomRowData.electroFrame); - AtomPlanMatrixColumn->gather(snap_->atomData.electroFrame, - atomColData.electroFrame); + + // if needed, gather the atomic eletrostatic information + if (storageLayout_ & DataStorage::dslDipole) { + AtomPlanVectorRow->gather(snap_->atomData.dipole, + atomRowData.dipole); + AtomPlanVectorColumn->gather(snap_->atomData.dipole, + atomColData.dipole); } + if (storageLayout_ & DataStorage::dslQuadrupole) { + AtomPlanMatrixRow->gather(snap_->atomData.quadrupole, + atomRowData.quadrupole); + AtomPlanMatrixColumn->gather(snap_->atomData.quadrupole, + atomColData.quadrupole); + } + + // if needed, gather the atomic fluctuating charge values + if (storageLayout_ & DataStorage::dslFlucQPosition) { + AtomPlanRealRow->gather(snap_->atomData.flucQPos, + atomRowData.flucQPos); + AtomPlanRealColumn->gather(snap_->atomData.flucQPos, + atomColData.flucQPos); + } + #endif } @@ -612,6 +704,21 @@ namespace OpenMD { AtomPlanRealColumn->scatter(atomColData.density, rho_tmp); for (int i = 0; i < n; i++) snap_->atomData.density[i] += rho_tmp[i]; + } + + // this isn't necessary if we don't have polarizable atoms, but + // we'll leave it here for now. + if (storageLayout_ & DataStorage::dslElectricField) { + + AtomPlanVectorRow->scatter(atomRowData.electricField, + snap_->atomData.electricField); + + int n = snap_->atomData.electricField.size(); + vector field_tmp(n, V3Zero); + AtomPlanVectorColumn->scatter(atomColData.electricField, + field_tmp); + for (int i = 0; i < n; i++) + snap_->atomData.electricField[i] += field_tmp[i]; } #endif } @@ -692,46 +799,175 @@ namespace OpenMD { } + if (storageLayout_ & DataStorage::dslFlucQForce) { + + int nq = snap_->atomData.flucQFrc.size(); + vector fqfrc_tmp(nq, 0.0); + + AtomPlanRealRow->scatter(atomRowData.flucQFrc, fqfrc_tmp); + for (int i = 0; i < nq; i++) { + snap_->atomData.flucQFrc[i] += fqfrc_tmp[i]; + fqfrc_tmp[i] = 0.0; + } + + AtomPlanRealColumn->scatter(atomColData.flucQFrc, fqfrc_tmp); + for (int i = 0; i < nq; i++) + snap_->atomData.flucQFrc[i] += fqfrc_tmp[i]; + + } + + if (storageLayout_ & DataStorage::dslElectricField) { + + int nef = snap_->atomData.electricField.size(); + vector efield_tmp(nef, V3Zero); + + AtomPlanVectorRow->scatter(atomRowData.electricField, efield_tmp); + for (int i = 0; i < nef; i++) { + snap_->atomData.electricField[i] += efield_tmp[i]; + efield_tmp[i] = 0.0; + } + + AtomPlanVectorColumn->scatter(atomColData.electricField, efield_tmp); + for (int i = 0; i < nef; i++) + snap_->atomData.electricField[i] += efield_tmp[i]; + } + + nLocal_ = snap_->getNumberOfAtoms(); vector pot_temp(nLocal_, Vector (0.0)); + vector expot_temp(nLocal_, + Vector (0.0)); // scatter/gather pot_row into the members of my column AtomPlanPotRow->scatter(pot_row, pot_temp); + AtomPlanPotRow->scatter(expot_row, expot_temp); - for (int ii = 0; ii < pot_temp.size(); ii++ ) + for (int ii = 0; ii < pot_temp.size(); ii++ ) pairwisePot += pot_temp[ii]; - + + for (int ii = 0; ii < expot_temp.size(); ii++ ) + excludedPot += expot_temp[ii]; + + if (storageLayout_ & DataStorage::dslParticlePot) { + // This is the pairwise contribution to the particle pot. The + // embedding contribution is added in each of the low level + // non-bonded routines. In single processor, this is done in + // unpackInteractionData, not in collectData. + for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { + for (int i = 0; i < nLocal_; i++) { + // factor of two is because the total potential terms are divided + // by 2 in parallel due to row/ column scatter + snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii); + } + } + } + fill(pot_temp.begin(), pot_temp.end(), Vector (0.0)); + fill(expot_temp.begin(), expot_temp.end(), + Vector (0.0)); AtomPlanPotColumn->scatter(pot_col, pot_temp); + AtomPlanPotColumn->scatter(expot_col, expot_temp); for (int ii = 0; ii < pot_temp.size(); ii++ ) pairwisePot += pot_temp[ii]; + + for (int ii = 0; ii < expot_temp.size(); ii++ ) + excludedPot += expot_temp[ii]; + + if (storageLayout_ & DataStorage::dslParticlePot) { + // This is the pairwise contribution to the particle pot. The + // embedding contribution is added in each of the low level + // non-bonded routines. In single processor, this is done in + // unpackInteractionData, not in collectData. + for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { + for (int i = 0; i < nLocal_; i++) { + // factor of two is because the total potential terms are divided + // by 2 in parallel due to row/ column scatter + snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii); + } + } + } + if (storageLayout_ & DataStorage::dslParticlePot) { + int npp = snap_->atomData.particlePot.size(); + vector ppot_temp(npp, 0.0); + + // This is the direct or embedding contribution to the particle + // pot. + + AtomPlanRealRow->scatter(atomRowData.particlePot, ppot_temp); + for (int i = 0; i < npp; i++) { + snap_->atomData.particlePot[i] += ppot_temp[i]; + } + + fill(ppot_temp.begin(), ppot_temp.end(), 0.0); + + AtomPlanRealColumn->scatter(atomColData.particlePot, ppot_temp); + for (int i = 0; i < npp; i++) { + snap_->atomData.particlePot[i] += ppot_temp[i]; + } + } + for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { RealType ploc1 = pairwisePot[ii]; RealType ploc2 = 0.0; - MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); pairwisePot[ii] = ploc2; } for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { - RealType ploc1 = embeddingPot[ii]; + RealType ploc1 = excludedPot[ii]; RealType ploc2 = 0.0; - MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); - embeddingPot[ii] = ploc2; + MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + excludedPot[ii] = ploc2; } + // Here be dragons. + MPI_Comm col = colComm.getComm(); + + MPI_Allreduce(MPI_IN_PLACE, + &snap_->frameData.conductiveHeatFlux[0], 3, + MPI_REALTYPE, MPI_SUM, col); + + #endif } - int ForceMatrixDecomposition::getNAtomsInRow() { + /** + * Collects information obtained during the post-pair (and embedding + * functional) loops onto local data structures. + */ + void ForceMatrixDecomposition::collectSelfData() { + snap_ = sman_->getCurrentSnapshot(); + storageLayout_ = sman_->getStorageLayout(); + #ifdef IS_MPI + for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { + RealType ploc1 = embeddingPot[ii]; + RealType ploc2 = 0.0; + MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + embeddingPot[ii] = ploc2; + } + for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { + RealType ploc1 = excludedSelfPot[ii]; + RealType ploc2 = 0.0; + MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + excludedSelfPot[ii] = ploc2; + } +#endif + + } + + + + int& ForceMatrixDecomposition::getNAtomsInRow() { +#ifdef IS_MPI return nAtomsInRow_; #else return nLocal_; @@ -741,7 +977,7 @@ namespace OpenMD { /** * returns the list of atoms belonging to this group. */ - vector ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ + vector& ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ #ifdef IS_MPI return groupListRow_[cg1]; #else @@ -749,7 +985,7 @@ namespace OpenMD { #endif } - vector ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ + vector& ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ #ifdef IS_MPI return groupListCol_[cg2]; #else @@ -766,11 +1002,29 @@ namespace OpenMD { d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1]; #endif - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } + Vector3d& ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){ +#ifdef IS_MPI + return cgColData.velocity[cg2]; +#else + return snap_->cgData.velocity[cg2]; +#endif + } + Vector3d& ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ +#ifdef IS_MPI + return atomColData.velocity[atom2]; +#else + return snap_->atomData.velocity[atom2]; +#endif + } + + Vector3d ForceMatrixDecomposition::getAtomToGroupVectorRow(int atom1, int cg1){ Vector3d d; @@ -780,8 +1034,9 @@ namespace OpenMD { #else d = snap_->cgData.position[cg1] - snap_->atomData.position[atom1]; #endif - - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } @@ -793,12 +1048,13 @@ namespace OpenMD { #else d = snap_->cgData.position[cg2] - snap_->atomData.position[atom2]; #endif - - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } - RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) { + RealType& ForceMatrixDecomposition::getMassFactorRow(int atom1) { #ifdef IS_MPI return massFactorsRow[atom1]; #else @@ -806,7 +1062,7 @@ namespace OpenMD { #endif } - RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { + RealType& ForceMatrixDecomposition::getMassFactorColumn(int atom2) { #ifdef IS_MPI return massFactorsCol[atom2]; #else @@ -823,12 +1079,13 @@ namespace OpenMD { #else d = snap_->atomData.position[atom2] - snap_->atomData.position[atom1]; #endif - - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } - vector ForceMatrixDecomposition::getExcludesForAtom(int atom1) { + vector& ForceMatrixDecomposition::getExcludesForAtom(int atom1) { return excludesForAtom[atom1]; } @@ -836,16 +1093,20 @@ namespace OpenMD { * We need to exclude some overcounted interactions that result from * the parallel decomposition. */ - bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { + bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) { 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]; + // group1 = cgRowToGlobal[cg1]; + // group2 = cgColToGlobal[cg2]; #else unique_id_1 = AtomLocalToGlobal[atom1]; unique_id_2 = AtomLocalToGlobal[atom2]; + int group1 = cgLocalToGlobal[cg1]; + int group2 = cgLocalToGlobal[cg2]; #endif if (unique_id_1 == unique_id_2) return true; @@ -857,6 +1118,12 @@ namespace OpenMD { } else { if ((unique_id_1 + unique_id_2) % 2 == 1) return true; } +#endif + +#ifndef IS_MPI + if (group1 == group2) { + if (unique_id_1 < unique_id_2) return true; + } #endif return false; @@ -908,25 +1175,36 @@ namespace OpenMD { idat.excluded = excludeAtomPair(atom1, atom2); #ifdef IS_MPI - idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]); - //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), - // ff_->getAtomType(identsCol[atom2]) ); - + //idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]); + idat.atid1 = identsRow[atom1]; + idat.atid2 = identsCol[atom2]; + + if (regionsRow[atom1] >= 0 && regionsCol[atom2] >= 0) { + idat.sameRegion = (regionsRow[atom1] == regionsCol[atom2]); + } else { + idat.sameRegion = false; + } + if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(atomRowData.aMat[atom1]); idat.A2 = &(atomColData.aMat[atom2]); } - if (storageLayout_ & DataStorage::dslElectroFrame) { - 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::dslDipole) { + idat.dipole1 = &(atomRowData.dipole[atom1]); + idat.dipole2 = &(atomColData.dipole[atom2]); + } + + if (storageLayout_ & DataStorage::dslQuadrupole) { + idat.quadrupole1 = &(atomRowData.quadrupole[atom1]); + idat.quadrupole2 = &(atomColData.quadrupole[atom2]); + } + if (storageLayout_ & DataStorage::dslDensity) { idat.rho1 = &(atomRowData.density[atom1]); idat.rho2 = &(atomColData.density[atom2]); @@ -952,32 +1230,43 @@ namespace OpenMD { idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); } + if (storageLayout_ & DataStorage::dslFlucQPosition) { + idat.flucQ1 = &(atomRowData.flucQPos[atom1]); + idat.flucQ2 = &(atomColData.flucQPos[atom2]); + } + #else + //idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); + idat.atid1 = idents[atom1]; + idat.atid2 = idents[atom2]; - // cerr << "atoms = " << atom1 << " " << atom2 << "\n"; - // cerr << "pos1 = " << snap_->atomData.position[atom1] << "\n"; - // cerr << "pos2 = " << snap_->atomData.position[atom2] << "\n"; + if (regions[atom1] >= 0 && regions[atom2] >= 0) { + idat.sameRegion = (regions[atom1] == regions[atom2]); + } else { + idat.sameRegion = false; + } - idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); - //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), - // ff_->getAtomType(idents[atom2]) ); - if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(snap_->atomData.aMat[atom1]); idat.A2 = &(snap_->atomData.aMat[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::dslDipole) { + idat.dipole1 = &(snap_->atomData.dipole[atom1]); + idat.dipole2 = &(snap_->atomData.dipole[atom2]); + } + + if (storageLayout_ & DataStorage::dslQuadrupole) { + idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]); + idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]); + } + if (storageLayout_ & DataStorage::dslDensity) { idat.rho1 = &(snap_->atomData.density[atom1]); idat.rho2 = &(snap_->atomData.density[atom2]); @@ -1002,6 +1291,12 @@ namespace OpenMD { idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); } + + if (storageLayout_ & DataStorage::dslFlucQPosition) { + idat.flucQ1 = &(snap_->atomData.flucQPos[atom1]); + idat.flucQ2 = &(snap_->atomData.flucQPos[atom2]); + } + #endif } @@ -1010,14 +1305,48 @@ namespace OpenMD { #ifdef IS_MPI pot_row[atom1] += RealType(0.5) * *(idat.pot); pot_col[atom2] += RealType(0.5) * *(idat.pot); + expot_row[atom1] += RealType(0.5) * *(idat.excludedPot); + expot_col[atom2] += RealType(0.5) * *(idat.excludedPot); atomRowData.force[atom1] += *(idat.f1); atomColData.force[atom2] -= *(idat.f1); + + if (storageLayout_ & DataStorage::dslFlucQForce) { + atomRowData.flucQFrc[atom1] -= *(idat.dVdFQ1); + atomColData.flucQFrc[atom2] -= *(idat.dVdFQ2); + } + + if (storageLayout_ & DataStorage::dslElectricField) { + atomRowData.electricField[atom1] += *(idat.eField1); + atomColData.electricField[atom2] += *(idat.eField2); + } + #else pairwisePot += *(idat.pot); + excludedPot += *(idat.excludedPot); snap_->atomData.force[atom1] += *(idat.f1); snap_->atomData.force[atom2] -= *(idat.f1); + + if (idat.doParticlePot) { + // This is the pairwise contribution to the particle pot. The + // embedding contribution is added in each of the low level + // non-bonded routines. In parallel, this calculation is done + // in collectData, not in unpackInteractionData. + snap_->atomData.particlePot[atom1] += *(idat.vpair) * *(idat.sw); + snap_->atomData.particlePot[atom2] += *(idat.vpair) * *(idat.sw); + } + + if (storageLayout_ & DataStorage::dslFlucQForce) { + snap_->atomData.flucQFrc[atom1] -= *(idat.dVdFQ1); + snap_->atomData.flucQFrc[atom2] -= *(idat.dVdFQ2); + } + + if (storageLayout_ & DataStorage::dslElectricField) { + snap_->atomData.electricField[atom1] += *(idat.eField1); + snap_->atomData.electricField[atom2] += *(idat.eField2); + } + #endif } @@ -1028,64 +1357,78 @@ namespace OpenMD { * first element of pair is row-indexed CutoffGroup * second element of pair is column-indexed CutoffGroup */ - vector > ForceMatrixDecomposition::buildNeighborList() { - - vector > neighborList; + void ForceMatrixDecomposition::buildNeighborList(vector >& neighborList) { + + neighborList.clear(); groupCutoffs cuts; bool doAllPairs = false; + RealType rList_ = (largestRcut_ + skinThickness_); + RealType rcut, rcutsq, rlistsq; + Snapshot* snap_ = sman_->getCurrentSnapshot(); + Mat3x3d box; + Mat3x3d invBox; + + Vector3d rs, scaled, dr; + Vector3i whichCell; + int cellIndex; + #ifdef IS_MPI cellListRow_.clear(); cellListCol_.clear(); #else cellList_.clear(); #endif - - RealType rList_ = (largestRcut_ + skinThickness_); - RealType rl2 = rList_ * rList_; - Snapshot* snap_ = sman_->getCurrentSnapshot(); - Mat3x3d Hmat = snap_->getHmat(); - Vector3d Hx = Hmat.getColumn(0); - Vector3d Hy = Hmat.getColumn(1); - Vector3d Hz = Hmat.getColumn(2); - - nCells_.x() = (int) ( Hx.length() )/ rList_; - nCells_.y() = (int) ( Hy.length() )/ rList_; - nCells_.z() = (int) ( Hz.length() )/ rList_; - + + if (!usePeriodicBoundaryConditions_) { + box = snap_->getBoundingBox(); + invBox = snap_->getInvBoundingBox(); + } else { + box = snap_->getHmat(); + invBox = snap_->getInvHmat(); + } + + Vector3d boxX = box.getColumn(0); + Vector3d boxY = box.getColumn(1); + Vector3d boxZ = box.getColumn(2); + + nCells_.x() = int( boxX.length() / rList_ ); + nCells_.y() = int( boxY.length() / rList_ ); + nCells_.z() = int( boxZ.length() / rList_ ); + // handle small boxes where the cell offsets can end up repeating cells if (nCells_.x() < 3) doAllPairs = true; if (nCells_.y() < 3) doAllPairs = true; if (nCells_.z() < 3) doAllPairs = true; - - Mat3x3d invHmat = snap_->getInvHmat(); - 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 - + if (!doAllPairs) { #ifdef IS_MPI - + for (int i = 0; i < nGroupsInRow_; i++) { rs = cgRowData.position[i]; // scaled positions relative to the box vectors - scaled = invHmat * rs; + scaled = invBox * rs; // wrap the vector back into the unit box by subtracting integer box // numbers for (int j = 0; j < 3; j++) { scaled[j] -= roundMe(scaled[j]); scaled[j] += 0.5; + // Handle the special case when an object is exactly on the + // boundary (a scaled coordinate of 1.0 is the same as + // scaled coordinate of 0.0) + if (scaled[j] >= 1.0) scaled[j] -= 1.0; } // find xyz-indices of cell that cutoffGroup is in. @@ -1103,13 +1446,17 @@ namespace OpenMD { rs = cgColData.position[i]; // scaled positions relative to the box vectors - scaled = invHmat * rs; + scaled = invBox * rs; // wrap the vector back into the unit box by subtracting integer box // numbers for (int j = 0; j < 3; j++) { scaled[j] -= roundMe(scaled[j]); scaled[j] += 0.5; + // Handle the special case when an object is exactly on the + // boundary (a scaled coordinate of 1.0 is the same as + // scaled coordinate of 0.0) + if (scaled[j] >= 1.0) scaled[j] -= 1.0; } // find xyz-indices of cell that cutoffGroup is in. @@ -1123,25 +1470,29 @@ namespace OpenMD { // 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; + scaled = invBox * rs; // wrap the vector back into the unit box by subtracting integer box // numbers for (int j = 0; j < 3; j++) { scaled[j] -= roundMe(scaled[j]); scaled[j] += 0.5; + // Handle the special case when an object is exactly on the + // boundary (a scaled coordinate of 1.0 is the same as + // scaled coordinate of 0.0) + if (scaled[j] >= 1.0) scaled[j] -= 1.0; } // find xyz-indices of cell that cutoffGroup is in. - whichCell.x() = nCells_.x() * scaled.x(); - whichCell.y() = nCells_.y() * scaled.y(); - whichCell.z() = nCells_.z() * scaled.z(); + whichCell.x() = int(nCells_.x() * scaled.x()); + whichCell.y() = int(nCells_.y() * scaled.y()); + whichCell.z() = int(nCells_.z() * scaled.z()); // find single index of this cell: cellIndex = Vlinear(whichCell, nCells_); @@ -1194,9 +1545,11 @@ namespace OpenMD { // & column indicies and will divide labor in the // force evaluation later. dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; - snap_->wrapVector(dr); - cuts = getGroupCutoffs( (*j1), (*j2) ); - if (dr.lengthSquare() < cuts.third) { + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } + getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair((*j1), (*j2))); } } @@ -1216,14 +1569,14 @@ namespace OpenMD { // allows atoms within a single cutoff group to // interact with each other. - - if (m2 != m1 || (*j2) >= (*j1) ) { dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; - snap_->wrapVector(dr); - cuts = getGroupCutoffs( (*j1), (*j2) ); - if (dr.lengthSquare() < cuts.third) { + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } + getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair((*j1), (*j2))); } } @@ -1240,9 +1593,11 @@ namespace OpenMD { for (int j1 = 0; j1 < nGroupsInRow_; j1++) { for (int j2 = 0; j2 < nGroupsInCol_; j2++) { dr = cgColData.position[j2] - cgRowData.position[j1]; - snap_->wrapVector(dr); - cuts = getGroupCutoffs( j1, j2 ); - if (dr.lengthSquare() < cuts.third) { + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } + getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair(j1, j2)); } } @@ -1253,9 +1608,11 @@ namespace OpenMD { // include self group interactions j2 == j1 for (int j2 = j1; j2 < nGroups_; j2++) { dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; - snap_->wrapVector(dr); - cuts = getGroupCutoffs( j1, j2 ); - if (dr.lengthSquare() < cuts.third) { + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } + getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq ); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair(j1, j2)); } } @@ -1268,7 +1625,5 @@ namespace OpenMD { 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