| 53 |
|
// surrounding cells (not just the 14 upper triangular blocks that |
| 54 |
|
// are used when the processor can see all pairs) |
| 55 |
|
#ifdef IS_MPI |
| 56 |
< |
cellOffsets_.push_back( Vector3i(-1, 0, 0) ); |
| 57 |
< |
cellOffsets_.push_back( Vector3i(-1,-1, 0) ); |
| 58 |
< |
cellOffsets_.push_back( Vector3i( 0,-1, 0) ); |
| 59 |
< |
cellOffsets_.push_back( Vector3i( 1,-1, 0) ); |
| 60 |
< |
cellOffsets_.push_back( Vector3i( 0, 0,-1) ); |
| 61 |
< |
cellOffsets_.push_back( Vector3i(-1, 0, 1) ); |
| 56 |
> |
cellOffsets_.clear(); |
| 57 |
|
cellOffsets_.push_back( Vector3i(-1,-1,-1) ); |
| 58 |
|
cellOffsets_.push_back( Vector3i( 0,-1,-1) ); |
| 59 |
< |
cellOffsets_.push_back( Vector3i( 1,-1,-1) ); |
| 59 |
> |
cellOffsets_.push_back( Vector3i( 1,-1,-1) ); |
| 60 |
> |
cellOffsets_.push_back( Vector3i(-1, 0,-1) ); |
| 61 |
> |
cellOffsets_.push_back( Vector3i( 0, 0,-1) ); |
| 62 |
|
cellOffsets_.push_back( Vector3i( 1, 0,-1) ); |
| 66 |
– |
cellOffsets_.push_back( Vector3i( 1, 1,-1) ); |
| 67 |
– |
cellOffsets_.push_back( Vector3i( 0, 1,-1) ); |
| 63 |
|
cellOffsets_.push_back( Vector3i(-1, 1,-1) ); |
| 64 |
+ |
cellOffsets_.push_back( Vector3i( 0, 1,-1) ); |
| 65 |
+ |
cellOffsets_.push_back( Vector3i( 1, 1,-1) ); |
| 66 |
+ |
cellOffsets_.push_back( Vector3i(-1,-1, 0) ); |
| 67 |
+ |
cellOffsets_.push_back( Vector3i( 0,-1, 0) ); |
| 68 |
+ |
cellOffsets_.push_back( Vector3i( 1,-1, 0) ); |
| 69 |
+ |
cellOffsets_.push_back( Vector3i(-1, 0, 0) ); |
| 70 |
+ |
cellOffsets_.push_back( Vector3i( 0, 0, 0) ); |
| 71 |
+ |
cellOffsets_.push_back( Vector3i( 1, 0, 0) ); |
| 72 |
+ |
cellOffsets_.push_back( Vector3i(-1, 1, 0) ); |
| 73 |
+ |
cellOffsets_.push_back( Vector3i( 0, 1, 0) ); |
| 74 |
+ |
cellOffsets_.push_back( Vector3i( 1, 1, 0) ); |
| 75 |
+ |
cellOffsets_.push_back( Vector3i(-1,-1, 1) ); |
| 76 |
+ |
cellOffsets_.push_back( Vector3i( 0,-1, 1) ); |
| 77 |
+ |
cellOffsets_.push_back( Vector3i( 1,-1, 1) ); |
| 78 |
+ |
cellOffsets_.push_back( Vector3i(-1, 0, 1) ); |
| 79 |
+ |
cellOffsets_.push_back( Vector3i( 0, 0, 1) ); |
| 80 |
+ |
cellOffsets_.push_back( Vector3i( 1, 0, 1) ); |
| 81 |
+ |
cellOffsets_.push_back( Vector3i(-1, 1, 1) ); |
| 82 |
+ |
cellOffsets_.push_back( Vector3i( 0, 1, 1) ); |
| 83 |
+ |
cellOffsets_.push_back( Vector3i( 1, 1, 1) ); |
| 84 |
|
#endif |
| 85 |
|
} |
| 86 |
|
|
| 168 |
|
AtomColToGlobal.resize(nAtomsInCol_); |
| 169 |
|
AtomPlanIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); |
| 170 |
|
AtomPlanIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); |
| 156 |
– |
|
| 157 |
– |
cerr << "Atoms in Local:\n"; |
| 158 |
– |
for (int i = 0; i < AtomLocalToGlobal.size(); i++) { |
| 159 |
– |
cerr << "i =\t" << i << "\t localAt =\t" << AtomLocalToGlobal[i] << "\n"; |
| 160 |
– |
} |
| 161 |
– |
cerr << "Atoms in Row:\n"; |
| 162 |
– |
for (int i = 0; i < AtomRowToGlobal.size(); i++) { |
| 163 |
– |
cerr << "i =\t" << i << "\t rowAt =\t" << AtomRowToGlobal[i] << "\n"; |
| 164 |
– |
} |
| 165 |
– |
cerr << "Atoms in Col:\n"; |
| 166 |
– |
for (int i = 0; i < AtomColToGlobal.size(); i++) { |
| 167 |
– |
cerr << "i =\t" << i << "\t colAt =\t" << AtomColToGlobal[i] << "\n"; |
| 168 |
– |
} |
| 171 |
|
|
| 172 |
|
cgRowToGlobal.resize(nGroupsInRow_); |
| 173 |
|
cgColToGlobal.resize(nGroupsInCol_); |
| 174 |
|
cgPlanIntRow->gather(cgLocalToGlobal, cgRowToGlobal); |
| 175 |
|
cgPlanIntColumn->gather(cgLocalToGlobal, cgColToGlobal); |
| 176 |
|
|
| 175 |
– |
cerr << "Gruops in Local:\n"; |
| 176 |
– |
for (int i = 0; i < cgLocalToGlobal.size(); i++) { |
| 177 |
– |
cerr << "i =\t" << i << "\t localCG =\t" << cgLocalToGlobal[i] << "\n"; |
| 178 |
– |
} |
| 179 |
– |
cerr << "Groups in Row:\n"; |
| 180 |
– |
for (int i = 0; i < cgRowToGlobal.size(); i++) { |
| 181 |
– |
cerr << "i =\t" << i << "\t rowCG =\t" << cgRowToGlobal[i] << "\n"; |
| 182 |
– |
} |
| 183 |
– |
cerr << "Groups in Col:\n"; |
| 184 |
– |
for (int i = 0; i < cgColToGlobal.size(); i++) { |
| 185 |
– |
cerr << "i =\t" << i << "\t colCG =\t" << cgColToGlobal[i] << "\n"; |
| 186 |
– |
} |
| 187 |
– |
|
| 188 |
– |
|
| 177 |
|
massFactorsRow.resize(nAtomsInRow_); |
| 178 |
|
massFactorsCol.resize(nAtomsInCol_); |
| 179 |
|
AtomPlanRealRow->gather(massFactors, massFactorsRow); |
| 233 |
|
} |
| 234 |
|
} |
| 235 |
|
|
| 236 |
< |
#endif |
| 249 |
< |
|
| 250 |
< |
// allocate memory for the parallel objects |
| 251 |
< |
atypesLocal.resize(nLocal_); |
| 252 |
< |
|
| 253 |
< |
for (int i = 0; i < nLocal_; i++) |
| 254 |
< |
atypesLocal[i] = ff_->getAtomType(idents[i]); |
| 255 |
< |
|
| 256 |
< |
groupList_.clear(); |
| 257 |
< |
groupList_.resize(nGroups_); |
| 258 |
< |
for (int i = 0; i < nGroups_; i++) { |
| 259 |
< |
int gid = cgLocalToGlobal[i]; |
| 260 |
< |
for (int j = 0; j < nLocal_; j++) { |
| 261 |
< |
int aid = AtomLocalToGlobal[j]; |
| 262 |
< |
if (globalGroupMembership[aid] == gid) { |
| 263 |
< |
groupList_[i].push_back(j); |
| 264 |
< |
} |
| 265 |
< |
} |
| 266 |
< |
} |
| 267 |
< |
|
| 236 |
> |
#else |
| 237 |
|
excludesForAtom.clear(); |
| 238 |
|
excludesForAtom.resize(nLocal_); |
| 239 |
|
toposForAtom.clear(); |
| 247 |
|
for (int j = 0; j < nLocal_; j++) { |
| 248 |
|
int jglob = AtomLocalToGlobal[j]; |
| 249 |
|
|
| 250 |
< |
if (excludes->hasPair(iglob, jglob)) |
| 250 |
> |
if (excludes->hasPair(iglob, jglob)) |
| 251 |
|
excludesForAtom[i].push_back(j); |
| 252 |
|
|
| 253 |
+ |
|
| 254 |
|
if (oneTwo->hasPair(iglob, jglob)) { |
| 255 |
|
toposForAtom[i].push_back(j); |
| 256 |
|
topoDist[i].push_back(1); |
| 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 |
|
} |
| 566 |
|
|
| 567 |
|
// gather up the cutoff group positions |
| 568 |
|
|
| 579 |
– |
cerr << "before gather\n"; |
| 580 |
– |
for (int i = 0; i < snap_->cgData.position.size(); i++) { |
| 581 |
– |
cerr << "cgpos = " << snap_->cgData.position[i] << "\n"; |
| 582 |
– |
} |
| 583 |
– |
|
| 569 |
|
cgPlanVectorRow->gather(snap_->cgData.position, |
| 570 |
|
cgRowData.position); |
| 571 |
|
|
| 587 |
– |
cerr << "after gather\n"; |
| 588 |
– |
for (int i = 0; i < cgRowData.position.size(); i++) { |
| 589 |
– |
cerr << "cgRpos = " << cgRowData.position[i] << "\n"; |
| 590 |
– |
} |
| 591 |
– |
|
| 572 |
|
cgPlanVectorColumn->gather(snap_->cgData.position, |
| 573 |
|
cgColData.position); |
| 594 |
– |
for (int i = 0; i < cgColData.position.size(); i++) { |
| 595 |
– |
cerr << "cgCpos = " << cgColData.position[i] << "\n"; |
| 596 |
– |
} |
| 574 |
|
|
| 575 |
|
|
| 576 |
|
// if needed, gather the atomic rotation matrices |
| 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(); |
| 709 |
|
|
| 710 |
|
for (int ii = 0; ii < pot_temp.size(); ii++ ) |
| 711 |
|
pairwisePot += pot_temp[ii]; |
| 712 |
+ |
|
| 713 |
+ |
for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
| 714 |
+ |
RealType ploc1 = pairwisePot[ii]; |
| 715 |
+ |
RealType ploc2 = 0.0; |
| 716 |
+ |
MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
| 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 |
|
|
| 736 |
– |
cerr << "pairwisePot = " << pairwisePot << "\n"; |
| 729 |
|
} |
| 730 |
|
|
| 731 |
|
int ForceMatrixDecomposition::getNAtomsInRow() { |
| 760 |
|
|
| 761 |
|
#ifdef IS_MPI |
| 762 |
|
d = cgColData.position[cg2] - cgRowData.position[cg1]; |
| 771 |
– |
cerr << "cg1 = " << cg1 << "\tcg1p = " << cgRowData.position[cg1] << "\n"; |
| 772 |
– |
cerr << "cg2 = " << cg2 << "\tcg2p = " << cgColData.position[cg2] << "\n"; |
| 763 |
|
#else |
| 764 |
|
d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1]; |
| 775 |
– |
cerr << "cg1 = " << cg1 << "\tcg1p = " << snap_->cgData.position[cg1] << "\n"; |
| 776 |
– |
cerr << "cg2 = " << cg2 << "\tcg2p = " << snap_->cgData.position[cg2] << "\n"; |
| 765 |
|
#endif |
| 766 |
|
|
| 767 |
|
snap_->wrapVector(d); |
| 837 |
|
bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { |
| 838 |
|
int unique_id_1, unique_id_2; |
| 839 |
|
|
| 852 |
– |
|
| 853 |
– |
cerr << "sap with atom1, atom2 =\t" << atom1 << "\t" << atom2 << "\n"; |
| 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 |
|
|
| 859 |
– |
cerr << "sap with uid1, uid2 =\t" << unique_id_1 << "\t" << unique_id_2 << "\n"; |
| 845 |
|
// this situation should only arise in MPI simulations |
| 846 |
|
if (unique_id_1 == unique_id_2) return true; |
| 847 |
|
|
| 865 |
|
* field) must still be handled for these pairs. |
| 866 |
|
*/ |
| 867 |
|
bool ForceMatrixDecomposition::excludeAtomPair(int atom1, int atom2) { |
| 868 |
< |
int unique_id_2; |
| 869 |
< |
#ifdef IS_MPI |
| 870 |
< |
// in MPI, we have to look up the unique IDs for the row atom. |
| 886 |
< |
unique_id_2 = AtomColToGlobal[atom2]; |
| 887 |
< |
#else |
| 888 |
< |
// in the normal loop, the atom numbers are unique |
| 889 |
< |
unique_id_2 = atom2; |
| 890 |
< |
#endif |
| 868 |
> |
|
| 869 |
> |
// excludesForAtom was constructed to use row/column indices in the MPI |
| 870 |
> |
// version, and to use local IDs in the non-MPI version: |
| 871 |
|
|
| 872 |
|
for (vector<int>::iterator i = excludesForAtom[atom1].begin(); |
| 873 |
|
i != excludesForAtom[atom1].end(); ++i) { |
| 874 |
< |
if ( (*i) == unique_id_2 ) return true; |
| 874 |
> |
if ( (*i) == atom2 ) return true; |
| 875 |
|
} |
| 876 |
|
|
| 877 |
|
return false; |
| 1111 |
|
// add this cutoff group to the list of groups in this cell; |
| 1112 |
|
cellListCol_[cellIndex].push_back(i); |
| 1113 |
|
} |
| 1114 |
+ |
|
| 1115 |
|
#else |
| 1116 |
|
for (int i = 0; i < nGroups_; i++) { |
| 1117 |
|
rs = snap_->cgData.position[i]; |
| 1137 |
|
// add this cutoff group to the list of groups in this cell; |
| 1138 |
|
cellList_[cellIndex].push_back(i); |
| 1139 |
|
} |
| 1140 |
+ |
|
| 1141 |
|
#endif |
| 1142 |
|
|
| 1143 |
|
for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
| 1150 |
|
os != cellOffsets_.end(); ++os) { |
| 1151 |
|
|
| 1152 |
|
Vector3i m2v = m1v + (*os); |
| 1153 |
< |
|
| 1153 |
> |
|
| 1154 |
> |
|
| 1155 |
|
if (m2v.x() >= nCells_.x()) { |
| 1156 |
|
m2v.x() = 0; |
| 1157 |
|
} else if (m2v.x() < 0) { |
| 1169 |
|
} else if (m2v.z() < 0) { |
| 1170 |
|
m2v.z() = nCells_.z() - 1; |
| 1171 |
|
} |
| 1172 |
< |
|
| 1172 |
> |
|
| 1173 |
|
int m2 = Vlinear (m2v, nCells_); |
| 1174 |
|
|
| 1175 |
|
#ifdef IS_MPI |
| 1178 |
|
for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
| 1179 |
|
j2 != cellListCol_[m2].end(); ++j2) { |
| 1180 |
|
|
| 1181 |
< |
// In parallel, we need to visit *all* pairs of row & |
| 1182 |
< |
// column indicies and will truncate later on. |
| 1181 |
> |
// In parallel, we need to visit *all* pairs of row |
| 1182 |
> |
// & column indicies and will divide labor in the |
| 1183 |
> |
// force evaluation later. |
| 1184 |
|
dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
| 1185 |
|
snap_->wrapVector(dr); |
| 1186 |
|
cuts = getGroupCutoffs( (*j1), (*j2) ); |