| 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" |
| 248 |
|
for (int j = 0; j < nLocal_; j++) { |
| 249 |
|
int jglob = AtomLocalToGlobal[j]; |
| 250 |
|
|
| 251 |
< |
if (excludes->hasPair(iglob, jglob)) |
| 251 |
> |
if (excludes->hasPair(iglob, jglob)) |
| 252 |
|
excludesForAtom[i].push_back(j); |
| 253 |
|
|
| 253 |
– |
|
| 254 |
|
if (oneTwo->hasPair(iglob, jglob)) { |
| 255 |
|
toposForAtom[i].push_back(j); |
| 256 |
|
topoDist[i].push_back(1); |
| 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 |
+ |
#else |
| 845 |
+ |
unique_id_1 = AtomLocalToGlobal[atom1]; |
| 846 |
+ |
unique_id_2 = AtomLocalToGlobal[atom2]; |
| 847 |
+ |
#endif |
| 848 |
|
|
| 845 |
– |
// this situation should only arise in MPI simulations |
| 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 |
|
|
| 876 |
|
|
| 877 |
|
for (vector<int>::iterator i = excludesForAtom[atom1].begin(); |
| 878 |
|
i != excludesForAtom[atom1].end(); ++i) { |
| 879 |
< |
if ( (*i) == atom2 ) return true; |
| 879 |
> |
if ( (*i) == atom2 ) return true; |
| 880 |
|
} |
| 881 |
|
|
| 882 |
|
return false; |
| 951 |
|
} |
| 952 |
|
|
| 953 |
|
#else |
| 954 |
+ |
|
| 955 |
+ |
|
| 956 |
+ |
// cerr << "atoms = " << atom1 << " " << atom2 << "\n"; |
| 957 |
+ |
// cerr << "pos1 = " << snap_->atomData.position[atom1] << "\n"; |
| 958 |
+ |
// cerr << "pos2 = " << snap_->atomData.position[atom2] << "\n"; |
| 959 |
|
|
| 960 |
|
idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); |
| 961 |
|
//idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
| 1006 |
|
|
| 1007 |
|
void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) { |
| 1008 |
|
#ifdef IS_MPI |
| 1009 |
< |
pot_row[atom1] += 0.5 * *(idat.pot); |
| 1010 |
< |
pot_col[atom2] += 0.5 * *(idat.pot); |
| 1009 |
> |
pot_row[atom1] += RealType(0.5) * *(idat.pot); |
| 1010 |
> |
pot_col[atom2] += RealType(0.5) * *(idat.pot); |
| 1011 |
|
|
| 1012 |
|
atomRowData.force[atom1] += *(idat.f1); |
| 1013 |
|
atomColData.force[atom2] -= *(idat.f1); |
| 1200 |
|
} |
| 1201 |
|
} |
| 1202 |
|
#else |
| 1193 |
– |
|
| 1203 |
|
for (vector<int>::iterator j1 = cellList_[m1].begin(); |
| 1204 |
|
j1 != cellList_[m1].end(); ++j1) { |
| 1205 |
|
for (vector<int>::iterator j2 = cellList_[m2].begin(); |
| 1206 |
|
j2 != cellList_[m2].end(); ++j2) { |
| 1207 |
< |
|
| 1207 |
> |
|
| 1208 |
|
// Always do this if we're in different cells or if |
| 1209 |
< |
// we're in the same cell and the global index of the |
| 1210 |
< |
// j2 cutoff group is less than the j1 cutoff group |
| 1211 |
< |
|
| 1212 |
< |
if (m2 != m1 || (*j2) < (*j1)) { |
| 1209 |
> |
// we're in the same cell and the global index of |
| 1210 |
> |
// the j2 cutoff group is greater than or equal to |
| 1211 |
> |
// the j1 cutoff group. Note that Rappaport's code |
| 1212 |
> |
// has a "less than" conditional here, but that |
| 1213 |
> |
// deals with atom-by-atom computation. OpenMD |
| 1214 |
> |
// allows atoms within a single cutoff group to |
| 1215 |
> |
// interact with each other. |
| 1216 |
> |
|
| 1217 |
> |
|
| 1218 |
> |
|
| 1219 |
> |
if (m2 != m1 || (*j2) >= (*j1) ) { |
| 1220 |
> |
|
| 1221 |
|
dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
| 1222 |
|
snap_->wrapVector(dr); |
| 1223 |
|
cuts = getGroupCutoffs( (*j1), (*j2) ); |
| 1236 |
|
// branch to do all cutoff group pairs |
| 1237 |
|
#ifdef IS_MPI |
| 1238 |
|
for (int j1 = 0; j1 < nGroupsInRow_; j1++) { |
| 1239 |
< |
for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
| 1239 |
> |
for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
| 1240 |
|
dr = cgColData.position[j2] - cgRowData.position[j1]; |
| 1241 |
|
snap_->wrapVector(dr); |
| 1242 |
|
cuts = getGroupCutoffs( j1, j2 ); |
| 1244 |
|
neighborList.push_back(make_pair(j1, j2)); |
| 1245 |
|
} |
| 1246 |
|
} |
| 1247 |
< |
} |
| 1247 |
> |
} |
| 1248 |
|
#else |
| 1249 |
< |
for (int j1 = 0; j1 < nGroups_ - 1; j1++) { |
| 1250 |
< |
for (int j2 = j1 + 1; j2 < nGroups_; j2++) { |
| 1249 |
> |
// include all groups here. |
| 1250 |
> |
for (int j1 = 0; j1 < nGroups_; j1++) { |
| 1251 |
> |
// include self group interactions j2 == j1 |
| 1252 |
> |
for (int j2 = j1; j2 < nGroups_; j2++) { |
| 1253 |
|
dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; |
| 1254 |
|
snap_->wrapVector(dr); |
| 1255 |
|
cuts = getGroupCutoffs( j1, j2 ); |
| 1256 |
|
if (dr.lengthSquare() < cuts.third) { |
| 1257 |
|
neighborList.push_back(make_pair(j1, j2)); |
| 1258 |
|
} |
| 1259 |
< |
} |
| 1260 |
< |
} |
| 1259 |
> |
} |
| 1260 |
> |
} |
| 1261 |
|
#endif |
| 1262 |
|
} |
| 1263 |
|
|
