--- trunk/src/parallel/ForceMatrixDecomposition.cpp 2013/06/19 17:19:07 1893 +++ trunk/src/parallel/ForceMatrixDecomposition.cpp 2014/02/26 14:14:50 1969 @@ -99,6 +99,7 @@ namespace OpenMD { 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(); @@ -118,8 +119,8 @@ namespace OpenMD { #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_); @@ -163,6 +164,12 @@ namespace OpenMD { 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_); @@ -308,6 +315,10 @@ namespace OpenMD { void ForceMatrixDecomposition::createGtypeCutoffMap() { + GrCut.clear(); + GrCutSq.clear(); + GrlistSq.clear(); + RealType tol = 1e-6; largestRcut_ = 0.0; int atid; @@ -413,13 +424,22 @@ namespace OpenMD { gTypeCutoffs.end()); #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, - MPI::MAX); + MPI_Allreduce(&groupMax, &groupMax, 1, MPI_REALTYPE, + MPI_MAX, MPI_COMM_WORLD); #endif RealType tradRcut = groupMax; + 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_) { @@ -442,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_); @@ -463,7 +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]; @@ -472,7 +495,11 @@ 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) { @@ -889,23 +916,23 @@ namespace OpenMD { 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 = excludedPot[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); excludedPot[ii] = ploc2; } // Here be dragons. - MPI::Intracomm col = colComm.getComm(); + MPI_Comm col = colComm.getComm(); - col.Allreduce(MPI::IN_PLACE, + MPI_Allreduce(MPI_IN_PLACE, &snap_->frameData.conductiveHeatFlux[0], 3, - MPI::REALTYPE, MPI::SUM); + MPI_REALTYPE, MPI_SUM, col); #endif @@ -924,13 +951,13 @@ namespace OpenMD { 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); + 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::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); excludedSelfPot[ii] = ploc2; } #endif @@ -1148,10 +1175,16 @@ 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]); @@ -1204,8 +1237,16 @@ namespace OpenMD { #else - idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); + //idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); + idat.atid1 = idents[atom1]; + idat.atid2 = idents[atom2]; + if (regions[atom1] >= 0 && regions[atom2] >= 0) { + idat.sameRegion = (regions[atom1] == regions[atom2]); + } else { + idat.sameRegion = false; + } + if (storageLayout_ & DataStorage::dslAmat) { idat.A1 = &(snap_->atomData.aMat[atom1]); idat.A2 = &(snap_->atomData.aMat[atom2]); @@ -1316,13 +1357,14 @@ 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; @@ -1350,9 +1392,9 @@ namespace OpenMD { 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_; + 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 @@ -1448,9 +1490,9 @@ namespace OpenMD { } // 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_); @@ -1506,8 +1548,8 @@ namespace OpenMD { if (usePeriodicBoundaryConditions_) { snap_->wrapVector(dr); } - cuts = getGroupCutoffs( (*j1), (*j2) ); - if (dr.lengthSquare() < cuts.third) { + getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair((*j1), (*j2))); } } @@ -1533,8 +1575,8 @@ namespace OpenMD { if (usePeriodicBoundaryConditions_) { snap_->wrapVector(dr); } - cuts = getGroupCutoffs( (*j1), (*j2) ); - if (dr.lengthSquare() < cuts.third) { + getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair((*j1), (*j2))); } } @@ -1554,8 +1596,8 @@ namespace OpenMD { if (usePeriodicBoundaryConditions_) { snap_->wrapVector(dr); } - cuts = getGroupCutoffs( j1, j2 ); - if (dr.lengthSquare() < cuts.third) { + getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair(j1, j2)); } } @@ -1569,8 +1611,8 @@ namespace OpenMD { if (usePeriodicBoundaryConditions_) { snap_->wrapVector(dr); } - cuts = getGroupCutoffs( j1, j2 ); - if (dr.lengthSquare() < cuts.third) { + getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq ); + if (dr.lengthSquare() < rlistsq) { neighborList.push_back(make_pair(j1, j2)); } } @@ -1583,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