--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2012/06/18 18:23:20 1756 +++ trunk/src/parallel/ForceMatrixDecomposition.cpp 2013/06/16 15:15:42 1879 @@ -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). */ @@ -175,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_); @@ -307,7 +310,6 @@ namespace OpenMD { RealType tol = 1e-6; largestRcut_ = 0.0; - RealType rc; int atid; set atypes = info_->getSimulatedAtomTypes(); @@ -392,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; @@ -417,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: @@ -474,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]; } @@ -484,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) { @@ -500,6 +504,12 @@ namespace OpenMD { Vector (0.0)); 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) { @@ -549,13 +559,6 @@ namespace OpenMD { atomColData.electricField.end(), V3Zero); } - if (storageLayout_ & DataStorage::dslFlucQForce) { - fill(atomRowData.flucQFrc.begin(), atomRowData.flucQFrc.end(), - 0.0); - fill(atomColData.flucQFrc.begin(), atomColData.flucQFrc.end(), - 0.0); - } - #endif // even in parallel, we need to zero out the local arrays: @@ -629,15 +632,22 @@ 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, @@ -669,6 +679,8 @@ namespace OpenMD { 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, @@ -676,7 +688,8 @@ namespace OpenMD { int n = snap_->atomData.electricField.size(); vector field_tmp(n, V3Zero); - AtomPlanVectorColumn->scatter(atomColData.electricField, field_tmp); + AtomPlanVectorColumn->scatter(atomColData.electricField, + field_tmp); for (int i = 0; i < n; i++) snap_->atomData.electricField[i] += field_tmp[i]; } @@ -776,17 +789,40 @@ namespace OpenMD { } + 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 @@ -804,12 +840,18 @@ namespace OpenMD { 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 @@ -851,6 +893,13 @@ namespace OpenMD { pairwisePot[ii] = ploc2; } + for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { + RealType ploc1 = excludedPot[ii]; + RealType ploc2 = 0.0; + MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); + excludedPot[ii] = ploc2; + } + // Here be dragons. MPI::Intracomm col = colComm.getComm(); @@ -877,6 +926,12 @@ namespace OpenMD { 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 @@ -920,7 +975,9 @@ namespace OpenMD { d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1]; #endif - snap_->wrapVector(d); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(d); + } return d; } @@ -950,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; } @@ -963,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; } @@ -993,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; } @@ -1007,19 +1067,19 @@ namespace OpenMD { * the parallel decomposition. */ bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) { - int unique_id_1, unique_id_2, group1, group2; + 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]; + // group1 = cgRowToGlobal[cg1]; + // group2 = cgColToGlobal[cg2]; #else unique_id_1 = AtomLocalToGlobal[atom1]; unique_id_2 = AtomLocalToGlobal[atom2]; - group1 = cgLocalToGlobal[cg1]; - group2 = cgLocalToGlobal[cg2]; + int group1 = cgLocalToGlobal[cg1]; + int group2 = cgLocalToGlobal[cg2]; #endif if (unique_id_1 == unique_id_2) return true; @@ -1097,16 +1157,21 @@ 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) { idat.rho1 = &(atomRowData.density[atom1]); idat.rho2 = &(atomColData.density[atom2]); @@ -1146,16 +1211,21 @@ namespace OpenMD { 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]); @@ -1194,6 +1264,8 @@ 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); @@ -1210,6 +1282,7 @@ namespace OpenMD { #else pairwisePot += *(idat.pot); + excludedPot += *(idat.excludedPot); snap_->atomData.force[atom1] += *(idat.f1); snap_->atomData.force[atom2] -= *(idat.f1); @@ -1249,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. @@ -1318,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. @@ -1338,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. @@ -1409,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))); @@ -1431,12 +1527,12 @@ 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); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } cuts = getGroupCutoffs( (*j1), (*j2) ); if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair((*j1), (*j2))); @@ -1455,7 +1551,9 @@ 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); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } cuts = getGroupCutoffs( j1, j2 ); if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair(j1, j2)); @@ -1468,7 +1566,9 @@ 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); + if (usePeriodicBoundaryConditions_) { + snap_->wrapVector(dr); + } cuts = getGroupCutoffs( j1, j2 ); if (dr.lengthSquare() < cuts.third) { neighborList.push_back(make_pair(j1, j2));