--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2011/08/18 20:18:19 1613 +++ trunk/src/parallel/ForceMatrixDecomposition.cpp 2013/06/16 15:15:42 1879 @@ -35,8 +35,9 @@ * * [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). - * [4] Vardeman & Gezelter, in progress (2009). + * [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). */ #include "parallel/ForceMatrixDecomposition.hpp" #include "math/SquareMatrix3.hpp" @@ -94,7 +95,7 @@ namespace OpenMD { storageLayout_ = sman_->getStorageLayout(); ff_ = info_->getForceField(); nLocal_ = snap_->getNumberOfAtoms(); - + nGroups_ = info_->getNLocalCutoffGroups(); // gather the information for atomtype IDs (atids): idents = info_->getIdentArray(); @@ -108,7 +109,13 @@ 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(); @@ -144,8 +151,13 @@ 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_); @@ -163,6 +175,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_); @@ -247,10 +262,9 @@ namespace OpenMD { for (int j = 0; j < nLocal_; j++) { int jglob = AtomLocalToGlobal[j]; - if (excludes->hasPair(iglob, jglob)) + if (excludes->hasPair(iglob, jglob)) excludesForAtom[i].push_back(j); - if (oneTwo->hasPair(iglob, jglob)) { toposForAtom[i].push_back(j); topoDist[i].push_back(1); @@ -296,7 +310,6 @@ namespace OpenMD { RealType tol = 1e-6; largestRcut_ = 0.0; - RealType rc; int atid; set atypes = info_->getSimulatedAtomTypes(); @@ -381,7 +394,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; @@ -406,8 +419,8 @@ namespace OpenMD { RealType tradRcut = groupMax; - for (int i = 0; i < gTypeCutoffs.size(); i++) { - for (int j = 0; j < gTypeCutoffs.size(); j++) { + for (unsigned int i = 0; i < gTypeCutoffs.size(); i++) { + for (unsigned int j = 0; j < gTypeCutoffs.size(); j++) { RealType thisRcut; switch(cutoffPolicy_) { case TRADITIONAL: @@ -450,7 +463,6 @@ namespace OpenMD { } } - groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { int i, j; #ifdef IS_MPI @@ -464,7 +476,7 @@ namespace OpenMD { } 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]; } @@ -474,6 +486,8 @@ namespace OpenMD { void ForceMatrixDecomposition::zeroWorkArrays() { pairwisePot = 0.0; embeddingPot = 0.0; + excludedPot = 0.0; + excludedSelfPot = 0.0; #ifdef IS_MPI if (storageLayout_ & DataStorage::dslForce) { @@ -492,6 +506,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 +545,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: @@ -537,19 +571,26 @@ namespace OpenMD { 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); } + if (storageLayout_ & DataStorage::dslSkippedCharge) { 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); + } } @@ -572,6 +613,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) { @@ -580,15 +632,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 } @@ -611,6 +678,21 @@ namespace OpenMD { 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 } @@ -690,26 +772,120 @@ 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; @@ -718,16 +894,51 @@ namespace OpenMD { } 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; + excludedPot[ii] = ploc2; } + // Here be dragons. + MPI::Intracomm col = colComm.getComm(); + + col.Allreduce(MPI::IN_PLACE, + &snap_->frameData.conductiveHeatFlux[0], 3, + MPI::REALTYPE, MPI::SUM); + + #endif } + /** + * 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::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); + embeddingPot[ii] = ploc2; + } + for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { + RealType ploc1 = excludedSelfPot[ii]; + RealType ploc2 = 0.0; + MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); + excludedSelfPot[ii] = ploc2; + } +#endif + + } + + + int ForceMatrixDecomposition::getNAtomsInRow() { #ifdef IS_MPI return nAtomsInRow_; @@ -764,11 +975,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; @@ -778,8 +1007,9 @@ namespace OpenMD { #else d = snap_->cgData.position[cg1] - snap_->atomData.position[atom1]; #endif - - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } @@ -791,8 +1021,9 @@ namespace OpenMD { #else d = snap_->cgData.position[cg2] - snap_->atomData.position[atom2]; #endif - - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } @@ -821,8 +1052,9 @@ namespace OpenMD { #else d = snap_->atomData.position[atom2] - snap_->atomData.position[atom1]; #endif - - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } @@ -834,24 +1066,39 @@ 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 - // this situation should only arise in MPI simulations if (unique_id_1 == unique_id_2) return true; - + +#ifdef IS_MPI // 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; + 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; } @@ -871,7 +1118,7 @@ namespace OpenMD { for (vector::iterator i = excludesForAtom[atom1].begin(); i != excludesForAtom[atom1].end(); ++i) { - if ( (*i) == atom2 ) return true; + if ( (*i) == atom2 ) return true; } return false; @@ -910,14 +1157,19 @@ namespace OpenMD { 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) { @@ -945,27 +1197,35 @@ namespace OpenMD { idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); } -#else + if (storageLayout_ & DataStorage::dslFlucQPosition) { + idat.flucQ1 = &(atomRowData.flucQPos[atom1]); + idat.flucQ2 = &(atomColData.flucQPos[atom2]); + } +#else + 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]); @@ -990,22 +1250,62 @@ 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 } 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); + 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 } @@ -1022,58 +1322,71 @@ namespace OpenMD { groupCutoffs cuts; bool doAllPairs = false; + RealType rList_ = (largestRcut_ + skinThickness_); + 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. @@ -1091,13 +1404,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. @@ -1111,19 +1428,23 @@ 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. @@ -1182,7 +1503,9 @@ namespace OpenMD { // & column indicies and will divide labor in the // force evaluation later. dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; - snap_->wrapVector(dr); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } cuts = getGroupCutoffs( (*j1), (*j2) ); if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair((*j1), (*j2))); @@ -1190,19 +1513,26 @@ 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 - - if (m2 != m1 || (*j2) < (*j1)) { + // we're in the same cell and the global index of + // the j2 cutoff group is greater than or equal to + // the j1 cutoff group. Note that Rappaport's code + // has a "less than" conditional here, but that + // deals with atom-by-atom computation. 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); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } cuts = getGroupCutoffs( (*j1), (*j2) ); if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair((*j1), (*j2))); @@ -1219,26 +1549,32 @@ namespace OpenMD { // branch to do all cutoff group pairs #ifdef IS_MPI for (int j1 = 0; j1 < nGroupsInRow_; j1++) { - for (int j2 = 0; j2 < nGroupsInCol_; j2++) { + for (int j2 = 0; j2 < nGroupsInCol_; j2++) { dr = cgColData.position[j2] - cgRowData.position[j1]; - snap_->wrapVector(dr); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } cuts = getGroupCutoffs( j1, j2 ); if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair(j1, j2)); } } - } + } #else - for (int j1 = 0; j1 < nGroups_ - 1; j1++) { - for (int j2 = j1 + 1; j2 < nGroups_; j2++) { + // include all groups here. + for (int j1 = 0; j1 < nGroups_; j1++) { + // 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); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } cuts = getGroupCutoffs( j1, j2 ); if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair(j1, j2)); } - } - } + } + } #endif }