--- branches/development/src/parallel/ForceMatrixDecomposition.cpp 2012/06/22 20:01:37 1761 +++ trunk/src/parallel/ForceMatrixDecomposition.cpp 2013/07/01 21:09:37 1895 @@ -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). */ @@ -310,7 +310,6 @@ namespace OpenMD { RealType tol = 1e-6; largestRcut_ = 0.0; - RealType rc; int atid; set atypes = info_->getSimulatedAtomTypes(); @@ -395,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; @@ -420,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: @@ -477,10 +476,10 @@ 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]; - } + } return 0; } @@ -558,13 +557,6 @@ namespace OpenMD { atomRowData.electricField.end(), V3Zero); fill(atomColData.electricField.begin(), 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 @@ -640,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, @@ -680,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, @@ -687,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]; } @@ -785,8 +787,25 @@ namespace OpenMD { 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_, @@ -920,7 +939,7 @@ namespace OpenMD { - int ForceMatrixDecomposition::getNAtomsInRow() { + int& ForceMatrixDecomposition::getNAtomsInRow() { #ifdef IS_MPI return nAtomsInRow_; #else @@ -931,7 +950,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 @@ -939,7 +958,7 @@ namespace OpenMD { #endif } - vector ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ + vector& ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ #ifdef IS_MPI return groupListCol_[cg2]; #else @@ -956,11 +975,13 @@ 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){ + Vector3d& ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){ #ifdef IS_MPI return cgColData.velocity[cg2]; #else @@ -968,7 +989,7 @@ namespace OpenMD { #endif } - Vector3d ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ + Vector3d& ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ #ifdef IS_MPI return atomColData.velocity[atom2]; #else @@ -986,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; } @@ -999,12 +1021,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 @@ -1012,7 +1035,7 @@ namespace OpenMD { #endif } - RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { + RealType& ForceMatrixDecomposition::getMassFactorColumn(int atom2) { #ifdef IS_MPI return massFactorsCol[atom2]; #else @@ -1029,12 +1052,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]; } @@ -1043,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; @@ -1125,6 +1149,8 @@ namespace OpenMD { #ifdef IS_MPI idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]); + idat.atid1 = identsRow[atom1]; + idat.atid2 = identsCol[atom2]; //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), // ff_->getAtomType(identsCol[atom2]) ); @@ -1133,16 +1159,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]); @@ -1176,22 +1207,29 @@ namespace OpenMD { #else idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); + idat.atid1 = idents[atom1]; + idat.atid2 = 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]); @@ -1282,64 +1320,77 @@ 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_); + 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. @@ -1357,13 +1408,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. @@ -1377,19 +1432,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. @@ -1448,7 +1507,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))); @@ -1470,12 +1531,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))); @@ -1494,7 +1555,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)); @@ -1507,7 +1570,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)); @@ -1522,7 +1587,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