| 310 |
|
|
| 311 |
|
RealType tol = 1e-6; |
| 312 |
|
largestRcut_ = 0.0; |
| 313 |
– |
RealType rc; |
| 313 |
|
int atid; |
| 314 |
|
set<AtomType*> atypes = info_->getSimulatedAtomTypes(); |
| 315 |
|
|
| 394 |
|
} |
| 395 |
|
|
| 396 |
|
bool gTypeFound = false; |
| 397 |
< |
for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { |
| 397 |
> |
for (unsigned int gt = 0; gt < gTypeCutoffs.size(); gt++) { |
| 398 |
|
if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) { |
| 399 |
|
groupToGtype[cg1] = gt; |
| 400 |
|
gTypeFound = true; |
| 419 |
|
|
| 420 |
|
RealType tradRcut = groupMax; |
| 421 |
|
|
| 422 |
< |
for (int i = 0; i < gTypeCutoffs.size(); i++) { |
| 423 |
< |
for (int j = 0; j < gTypeCutoffs.size(); j++) { |
| 422 |
> |
for (unsigned int i = 0; i < gTypeCutoffs.size(); i++) { |
| 423 |
> |
for (unsigned int j = 0; j < gTypeCutoffs.size(); j++) { |
| 424 |
|
RealType thisRcut; |
| 425 |
|
switch(cutoffPolicy_) { |
| 426 |
|
case TRADITIONAL: |
| 476 |
|
} |
| 477 |
|
|
| 478 |
|
int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { |
| 479 |
< |
for (int j = 0; j < toposForAtom[atom1].size(); j++) { |
| 479 |
> |
for (unsigned int j = 0; j < toposForAtom[atom1].size(); j++) { |
| 480 |
|
if (toposForAtom[atom1][j] == atom2) |
| 481 |
|
return topoDist[atom1][j]; |
| 482 |
|
} |
| 487 |
|
pairwisePot = 0.0; |
| 488 |
|
embeddingPot = 0.0; |
| 489 |
|
excludedPot = 0.0; |
| 490 |
+ |
excludedSelfPot = 0.0; |
| 491 |
|
|
| 492 |
|
#ifdef IS_MPI |
| 493 |
|
if (storageLayout_ & DataStorage::dslForce) { |
| 559 |
|
atomColData.electricField.end(), V3Zero); |
| 560 |
|
} |
| 561 |
|
|
| 562 |
– |
if (storageLayout_ & DataStorage::dslFlucQForce) { |
| 563 |
– |
fill(atomRowData.flucQFrc.begin(), atomRowData.flucQFrc.end(), |
| 564 |
– |
0.0); |
| 565 |
– |
fill(atomColData.flucQFrc.begin(), atomColData.flucQFrc.end(), |
| 566 |
– |
0.0); |
| 567 |
– |
} |
| 568 |
– |
|
| 562 |
|
#endif |
| 563 |
|
// even in parallel, we need to zero out the local arrays: |
| 564 |
|
|
| 632 |
|
AtomPlanMatrixColumn->gather(snap_->atomData.aMat, |
| 633 |
|
atomColData.aMat); |
| 634 |
|
} |
| 635 |
< |
|
| 636 |
< |
// if needed, gather the atomic eletrostatic frames |
| 637 |
< |
if (storageLayout_ & DataStorage::dslElectroFrame) { |
| 638 |
< |
AtomPlanMatrixRow->gather(snap_->atomData.electroFrame, |
| 639 |
< |
atomRowData.electroFrame); |
| 640 |
< |
AtomPlanMatrixColumn->gather(snap_->atomData.electroFrame, |
| 641 |
< |
atomColData.electroFrame); |
| 635 |
> |
|
| 636 |
> |
// if needed, gather the atomic eletrostatic information |
| 637 |
> |
if (storageLayout_ & DataStorage::dslDipole) { |
| 638 |
> |
AtomPlanVectorRow->gather(snap_->atomData.dipole, |
| 639 |
> |
atomRowData.dipole); |
| 640 |
> |
AtomPlanVectorColumn->gather(snap_->atomData.dipole, |
| 641 |
> |
atomColData.dipole); |
| 642 |
|
} |
| 643 |
|
|
| 644 |
+ |
if (storageLayout_ & DataStorage::dslQuadrupole) { |
| 645 |
+ |
AtomPlanMatrixRow->gather(snap_->atomData.quadrupole, |
| 646 |
+ |
atomRowData.quadrupole); |
| 647 |
+ |
AtomPlanMatrixColumn->gather(snap_->atomData.quadrupole, |
| 648 |
+ |
atomColData.quadrupole); |
| 649 |
+ |
} |
| 650 |
+ |
|
| 651 |
|
// if needed, gather the atomic fluctuating charge values |
| 652 |
|
if (storageLayout_ & DataStorage::dslFlucQPosition) { |
| 653 |
|
AtomPlanRealRow->gather(snap_->atomData.flucQPos, |
| 686 |
|
|
| 687 |
|
int n = snap_->atomData.electricField.size(); |
| 688 |
|
vector<Vector3d> field_tmp(n, V3Zero); |
| 689 |
< |
AtomPlanVectorColumn->scatter(atomColData.electricField, field_tmp); |
| 689 |
> |
AtomPlanVectorColumn->scatter(atomColData.electricField, |
| 690 |
> |
field_tmp); |
| 691 |
|
for (int i = 0; i < n; i++) |
| 692 |
|
snap_->atomData.electricField[i] += field_tmp[i]; |
| 693 |
|
} |
| 907 |
|
RealType ploc2 = 0.0; |
| 908 |
|
MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
| 909 |
|
embeddingPot[ii] = ploc2; |
| 910 |
+ |
} |
| 911 |
+ |
for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
| 912 |
+ |
RealType ploc1 = excludedSelfPot[ii]; |
| 913 |
+ |
RealType ploc2 = 0.0; |
| 914 |
+ |
MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
| 915 |
+ |
excludedSelfPot[ii] = ploc2; |
| 916 |
|
} |
| 917 |
|
#endif |
| 918 |
|
|
| 1043 |
|
* the parallel decomposition. |
| 1044 |
|
*/ |
| 1045 |
|
bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) { |
| 1046 |
< |
int unique_id_1, unique_id_2, group1, group2; |
| 1046 |
> |
int unique_id_1, unique_id_2; |
| 1047 |
|
|
| 1048 |
|
#ifdef IS_MPI |
| 1049 |
|
// in MPI, we have to look up the unique IDs for each atom |
| 1050 |
|
unique_id_1 = AtomRowToGlobal[atom1]; |
| 1051 |
|
unique_id_2 = AtomColToGlobal[atom2]; |
| 1052 |
< |
group1 = cgRowToGlobal[cg1]; |
| 1053 |
< |
group2 = cgColToGlobal[cg2]; |
| 1052 |
> |
// group1 = cgRowToGlobal[cg1]; |
| 1053 |
> |
// group2 = cgColToGlobal[cg2]; |
| 1054 |
|
#else |
| 1055 |
|
unique_id_1 = AtomLocalToGlobal[atom1]; |
| 1056 |
|
unique_id_2 = AtomLocalToGlobal[atom2]; |
| 1057 |
< |
group1 = cgLocalToGlobal[cg1]; |
| 1058 |
< |
group2 = cgLocalToGlobal[cg2]; |
| 1057 |
> |
int group1 = cgLocalToGlobal[cg1]; |
| 1058 |
> |
int group2 = cgLocalToGlobal[cg2]; |
| 1059 |
|
#endif |
| 1060 |
|
|
| 1061 |
|
if (unique_id_1 == unique_id_2) return true; |
| 1133 |
|
idat.A2 = &(atomColData.aMat[atom2]); |
| 1134 |
|
} |
| 1135 |
|
|
| 1129 |
– |
if (storageLayout_ & DataStorage::dslElectroFrame) { |
| 1130 |
– |
idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
| 1131 |
– |
idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
| 1132 |
– |
} |
| 1133 |
– |
|
| 1136 |
|
if (storageLayout_ & DataStorage::dslTorque) { |
| 1137 |
|
idat.t1 = &(atomRowData.torque[atom1]); |
| 1138 |
|
idat.t2 = &(atomColData.torque[atom2]); |
| 1139 |
|
} |
| 1140 |
|
|
| 1141 |
+ |
if (storageLayout_ & DataStorage::dslDipole) { |
| 1142 |
+ |
idat.dipole1 = &(atomRowData.dipole[atom1]); |
| 1143 |
+ |
idat.dipole2 = &(atomColData.dipole[atom2]); |
| 1144 |
+ |
} |
| 1145 |
+ |
|
| 1146 |
+ |
if (storageLayout_ & DataStorage::dslQuadrupole) { |
| 1147 |
+ |
idat.quadrupole1 = &(atomRowData.quadrupole[atom1]); |
| 1148 |
+ |
idat.quadrupole2 = &(atomColData.quadrupole[atom2]); |
| 1149 |
+ |
} |
| 1150 |
+ |
|
| 1151 |
|
if (storageLayout_ & DataStorage::dslDensity) { |
| 1152 |
|
idat.rho1 = &(atomRowData.density[atom1]); |
| 1153 |
|
idat.rho2 = &(atomColData.density[atom2]); |
| 1187 |
|
idat.A2 = &(snap_->atomData.aMat[atom2]); |
| 1188 |
|
} |
| 1189 |
|
|
| 1190 |
< |
if (storageLayout_ & DataStorage::dslElectroFrame) { |
| 1179 |
< |
idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
| 1180 |
< |
idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
| 1181 |
< |
} |
| 1190 |
> |
RealType ct = dot(idat.A1->getColumn(2), idat.A2->getColumn(2)); |
| 1191 |
|
|
| 1192 |
|
if (storageLayout_ & DataStorage::dslTorque) { |
| 1193 |
|
idat.t1 = &(snap_->atomData.torque[atom1]); |
| 1194 |
|
idat.t2 = &(snap_->atomData.torque[atom2]); |
| 1195 |
|
} |
| 1196 |
|
|
| 1197 |
+ |
if (storageLayout_ & DataStorage::dslDipole) { |
| 1198 |
+ |
idat.dipole1 = &(snap_->atomData.dipole[atom1]); |
| 1199 |
+ |
idat.dipole2 = &(snap_->atomData.dipole[atom2]); |
| 1200 |
+ |
} |
| 1201 |
+ |
|
| 1202 |
+ |
if (storageLayout_ & DataStorage::dslQuadrupole) { |
| 1203 |
+ |
idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]); |
| 1204 |
+ |
idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]); |
| 1205 |
+ |
} |
| 1206 |
+ |
|
| 1207 |
|
if (storageLayout_ & DataStorage::dslDensity) { |
| 1208 |
|
idat.rho1 = &(snap_->atomData.density[atom1]); |
| 1209 |
|
idat.rho2 = &(snap_->atomData.density[atom2]); |
| 1308 |
|
#endif |
| 1309 |
|
|
| 1310 |
|
RealType rList_ = (largestRcut_ + skinThickness_); |
| 1292 |
– |
RealType rl2 = rList_ * rList_; |
| 1311 |
|
Snapshot* snap_ = sman_->getCurrentSnapshot(); |
| 1312 |
|
Mat3x3d Hmat = snap_->getHmat(); |
| 1313 |
|
Vector3d Hx = Hmat.getColumn(0); |
| 1351 |
|
for (int j = 0; j < 3; j++) { |
| 1352 |
|
scaled[j] -= roundMe(scaled[j]); |
| 1353 |
|
scaled[j] += 0.5; |
| 1354 |
+ |
// Handle the special case when an object is exactly on the |
| 1355 |
+ |
// boundary (a scaled coordinate of 1.0 is the same as |
| 1356 |
+ |
// scaled coordinate of 0.0) |
| 1357 |
+ |
if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
| 1358 |
|
} |
| 1359 |
|
|
| 1360 |
|
// find xyz-indices of cell that cutoffGroup is in. |
| 1379 |
|
for (int j = 0; j < 3; j++) { |
| 1380 |
|
scaled[j] -= roundMe(scaled[j]); |
| 1381 |
|
scaled[j] += 0.5; |
| 1382 |
+ |
// Handle the special case when an object is exactly on the |
| 1383 |
+ |
// boundary (a scaled coordinate of 1.0 is the same as |
| 1384 |
+ |
// scaled coordinate of 0.0) |
| 1385 |
+ |
if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
| 1386 |
|
} |
| 1387 |
|
|
| 1388 |
|
// find xyz-indices of cell that cutoffGroup is in. |
| 1409 |
|
for (int j = 0; j < 3; j++) { |
| 1410 |
|
scaled[j] -= roundMe(scaled[j]); |
| 1411 |
|
scaled[j] += 0.5; |
| 1412 |
+ |
// Handle the special case when an object is exactly on the |
| 1413 |
+ |
// boundary (a scaled coordinate of 1.0 is the same as |
| 1414 |
+ |
// scaled coordinate of 0.0) |
| 1415 |
+ |
if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
| 1416 |
|
} |
| 1417 |
|
|
| 1418 |
|
// find xyz-indices of cell that cutoffGroup is in. |