| 36 |
|
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
| 37 |
|
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
| 38 |
|
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
| 39 |
< |
* [4] Vardeman & Gezelter, in progress (2009). |
| 39 |
> |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
| 40 |
> |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 42 |
|
#include "parallel/ForceMatrixDecomposition.hpp" |
| 43 |
|
#include "math/SquareMatrix3.hpp" |
| 230 |
|
topoDist[i].push_back(3); |
| 231 |
|
} |
| 232 |
|
} |
| 232 |
– |
} |
| 233 |
– |
} |
| 234 |
– |
} |
| 235 |
– |
|
| 236 |
– |
#endif |
| 237 |
– |
|
| 238 |
– |
// allocate memory for the parallel objects |
| 239 |
– |
atypesLocal.resize(nLocal_); |
| 240 |
– |
|
| 241 |
– |
for (int i = 0; i < nLocal_; i++) |
| 242 |
– |
atypesLocal[i] = ff_->getAtomType(idents[i]); |
| 243 |
– |
|
| 244 |
– |
groupList_.clear(); |
| 245 |
– |
groupList_.resize(nGroups_); |
| 246 |
– |
for (int i = 0; i < nGroups_; i++) { |
| 247 |
– |
int gid = cgLocalToGlobal[i]; |
| 248 |
– |
for (int j = 0; j < nLocal_; j++) { |
| 249 |
– |
int aid = AtomLocalToGlobal[j]; |
| 250 |
– |
if (globalGroupMembership[aid] == gid) { |
| 251 |
– |
groupList_[i].push_back(j); |
| 233 |
|
} |
| 234 |
|
} |
| 235 |
|
} |
| 236 |
|
|
| 237 |
+ |
#else |
| 238 |
|
excludesForAtom.clear(); |
| 239 |
|
excludesForAtom.resize(nLocal_); |
| 240 |
|
toposForAtom.clear(); |
| 267 |
|
} |
| 268 |
|
} |
| 269 |
|
} |
| 270 |
< |
|
| 270 |
> |
#endif |
| 271 |
> |
|
| 272 |
> |
// allocate memory for the parallel objects |
| 273 |
> |
atypesLocal.resize(nLocal_); |
| 274 |
> |
|
| 275 |
> |
for (int i = 0; i < nLocal_; i++) |
| 276 |
> |
atypesLocal[i] = ff_->getAtomType(idents[i]); |
| 277 |
> |
|
| 278 |
> |
groupList_.clear(); |
| 279 |
> |
groupList_.resize(nGroups_); |
| 280 |
> |
for (int i = 0; i < nGroups_; i++) { |
| 281 |
> |
int gid = cgLocalToGlobal[i]; |
| 282 |
> |
for (int j = 0; j < nLocal_; j++) { |
| 283 |
> |
int aid = AtomLocalToGlobal[j]; |
| 284 |
> |
if (globalGroupMembership[aid] == gid) { |
| 285 |
> |
groupList_[i].push_back(j); |
| 286 |
> |
} |
| 287 |
> |
} |
| 288 |
> |
} |
| 289 |
> |
|
| 290 |
> |
|
| 291 |
|
createGtypeCutoffMap(); |
| 292 |
|
|
| 293 |
|
} |
| 685 |
|
} |
| 686 |
|
|
| 687 |
|
AtomPlanRealColumn->scatter(atomColData.skippedCharge, skch_tmp); |
| 688 |
< |
for (int i = 0; i < ns; i++) |
| 688 |
> |
for (int i = 0; i < ns; i++) |
| 689 |
|
snap_->atomData.skippedCharge[i] += skch_tmp[i]; |
| 690 |
+ |
|
| 691 |
|
} |
| 692 |
|
|
| 693 |
|
nLocal_ = snap_->getNumberOfAtoms(); |
| 717 |
|
pairwisePot[ii] = ploc2; |
| 718 |
|
} |
| 719 |
|
|
| 720 |
+ |
for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
| 721 |
+ |
RealType ploc1 = embeddingPot[ii]; |
| 722 |
+ |
RealType ploc2 = 0.0; |
| 723 |
+ |
MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
| 724 |
+ |
embeddingPot[ii] = ploc2; |
| 725 |
+ |
} |
| 726 |
+ |
|
| 727 |
|
#endif |
| 728 |
|
|
| 729 |
|
} |
| 836 |
|
*/ |
| 837 |
|
bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { |
| 838 |
|
int unique_id_1, unique_id_2; |
| 839 |
< |
|
| 839 |
> |
|
| 840 |
|
#ifdef IS_MPI |
| 841 |
|
// in MPI, we have to look up the unique IDs for each atom |
| 842 |
|
unique_id_1 = AtomRowToGlobal[atom1]; |
| 843 |
|
unique_id_2 = AtomColToGlobal[atom2]; |
| 844 |
< |
|
| 845 |
< |
// this situation should only arise in MPI simulations |
| 844 |
> |
#else |
| 845 |
> |
unique_id_1 = AtomLocalToGlobal[atom1]; |
| 846 |
> |
unique_id_2 = AtomLocalToGlobal[atom2]; |
| 847 |
> |
#endif |
| 848 |
> |
|
| 849 |
|
if (unique_id_1 == unique_id_2) return true; |
| 850 |
< |
|
| 850 |
> |
|
| 851 |
> |
#ifdef IS_MPI |
| 852 |
|
// this prevents us from doing the pair on multiple processors |
| 853 |
|
if (unique_id_1 < unique_id_2) { |
| 854 |
|
if ((unique_id_1 + unique_id_2) % 2 == 0) return true; |
| 855 |
|
} else { |
| 856 |
< |
if ((unique_id_1 + unique_id_2) % 2 == 1) return true; |
| 856 |
> |
if ((unique_id_1 + unique_id_2) % 2 == 1) return true; |
| 857 |
|
} |
| 858 |
|
#endif |
| 859 |
+ |
|
| 860 |
|
return false; |
| 861 |
|
} |
| 862 |
|
|
| 870 |
|
* field) must still be handled for these pairs. |
| 871 |
|
*/ |
| 872 |
|
bool ForceMatrixDecomposition::excludeAtomPair(int atom1, int atom2) { |
| 873 |
< |
int unique_id_2; |
| 874 |
< |
#ifdef IS_MPI |
| 875 |
< |
// in MPI, we have to look up the unique IDs for the row atom. |
| 861 |
< |
unique_id_2 = AtomColToGlobal[atom2]; |
| 862 |
< |
#else |
| 863 |
< |
// in the normal loop, the atom numbers are unique |
| 864 |
< |
unique_id_2 = atom2; |
| 865 |
< |
#endif |
| 873 |
> |
|
| 874 |
> |
// excludesForAtom was constructed to use row/column indices in the MPI |
| 875 |
> |
// version, and to use local IDs in the non-MPI version: |
| 876 |
|
|
| 877 |
|
for (vector<int>::iterator i = excludesForAtom[atom1].begin(); |
| 878 |
|
i != excludesForAtom[atom1].end(); ++i) { |
| 879 |
< |
if ( (*i) == unique_id_2 ) return true; |
| 879 |
> |
if ( (*i) == atom2 ) return true; |
| 880 |
|
} |
| 881 |
|
|
| 882 |
|
return false; |
| 1001 |
|
|
| 1002 |
|
void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) { |
| 1003 |
|
#ifdef IS_MPI |
| 1004 |
< |
pot_row[atom1] += 0.5 * *(idat.pot); |
| 1005 |
< |
pot_col[atom2] += 0.5 * *(idat.pot); |
| 1004 |
> |
pot_row[atom1] += RealType(0.5) * *(idat.pot); |
| 1005 |
> |
pot_col[atom2] += RealType(0.5) * *(idat.pot); |
| 1006 |
|
|
| 1007 |
|
atomRowData.force[atom1] += *(idat.f1); |
| 1008 |
|
atomColData.force[atom2] -= *(idat.f1); |
| 1195 |
|
} |
| 1196 |
|
} |
| 1197 |
|
#else |
| 1188 |
– |
|
| 1198 |
|
for (vector<int>::iterator j1 = cellList_[m1].begin(); |
| 1199 |
|
j1 != cellList_[m1].end(); ++j1) { |
| 1200 |
|
for (vector<int>::iterator j2 = cellList_[m2].begin(); |
| 1201 |
|
j2 != cellList_[m2].end(); ++j2) { |
| 1202 |
< |
|
| 1202 |
> |
|
| 1203 |
|
// Always do this if we're in different cells or if |
| 1204 |
< |
// we're in the same cell and the global index of the |
| 1205 |
< |
// j2 cutoff group is less than the j1 cutoff group |
| 1206 |
< |
|
| 1207 |
< |
if (m2 != m1 || (*j2) < (*j1)) { |
| 1204 |
> |
// we're in the same cell and the global index of |
| 1205 |
> |
// the j2 cutoff group is greater than or equal to |
| 1206 |
> |
// the j1 cutoff group. Note that Rappaport's code |
| 1207 |
> |
// has a "less than" conditional here, but that |
| 1208 |
> |
// deals with atom-by-atom computation. OpenMD |
| 1209 |
> |
// allows atoms within a single cutoff group to |
| 1210 |
> |
// interact with each other. |
| 1211 |
> |
|
| 1212 |
> |
|
| 1213 |
> |
|
| 1214 |
> |
if (m2 != m1 || (*j2) >= (*j1) ) { |
| 1215 |
> |
|
| 1216 |
|
dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
| 1217 |
|
snap_->wrapVector(dr); |
| 1218 |
|
cuts = getGroupCutoffs( (*j1), (*j2) ); |
| 1231 |
|
// branch to do all cutoff group pairs |
| 1232 |
|
#ifdef IS_MPI |
| 1233 |
|
for (int j1 = 0; j1 < nGroupsInRow_; j1++) { |
| 1234 |
< |
for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
| 1234 |
> |
for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
| 1235 |
|
dr = cgColData.position[j2] - cgRowData.position[j1]; |
| 1236 |
|
snap_->wrapVector(dr); |
| 1237 |
|
cuts = getGroupCutoffs( j1, j2 ); |
| 1239 |
|
neighborList.push_back(make_pair(j1, j2)); |
| 1240 |
|
} |
| 1241 |
|
} |
| 1242 |
< |
} |
| 1242 |
> |
} |
| 1243 |
|
#else |
| 1244 |
< |
for (int j1 = 0; j1 < nGroups_ - 1; j1++) { |
| 1245 |
< |
for (int j2 = j1 + 1; j2 < nGroups_; j2++) { |
| 1244 |
> |
// include all groups here. |
| 1245 |
> |
for (int j1 = 0; j1 < nGroups_; j1++) { |
| 1246 |
> |
// include self group interactions j2 == j1 |
| 1247 |
> |
for (int j2 = j1; j2 < nGroups_; j2++) { |
| 1248 |
|
dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; |
| 1249 |
|
snap_->wrapVector(dr); |
| 1250 |
|
cuts = getGroupCutoffs( j1, j2 ); |
| 1251 |
|
if (dr.lengthSquare() < cuts.third) { |
| 1252 |
|
neighborList.push_back(make_pair(j1, j2)); |
| 1253 |
|
} |
| 1254 |
< |
} |
| 1255 |
< |
} |
| 1254 |
> |
} |
| 1255 |
> |
} |
| 1256 |
|
#endif |
| 1257 |
|
} |
| 1258 |
|
|
