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]); |
1352 |
|
for (int j = 0; j < 3; j++) { |
1353 |
|
scaled[j] -= roundMe(scaled[j]); |
1354 |
|
scaled[j] += 0.5; |
1355 |
+ |
// Handle the special case when an object is exactly on the |
1356 |
+ |
// boundary (a scaled coordinate of 1.0 is the same as |
1357 |
+ |
// scaled coordinate of 0.0) |
1358 |
+ |
if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
1359 |
|
} |
1360 |
|
|
1361 |
|
// find xyz-indices of cell that cutoffGroup is in. |
1380 |
|
for (int j = 0; j < 3; j++) { |
1381 |
|
scaled[j] -= roundMe(scaled[j]); |
1382 |
|
scaled[j] += 0.5; |
1383 |
+ |
// Handle the special case when an object is exactly on the |
1384 |
+ |
// boundary (a scaled coordinate of 1.0 is the same as |
1385 |
+ |
// scaled coordinate of 0.0) |
1386 |
+ |
if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
1387 |
|
} |
1388 |
|
|
1389 |
|
// find xyz-indices of cell that cutoffGroup is in. |
1410 |
|
for (int j = 0; j < 3; j++) { |
1411 |
|
scaled[j] -= roundMe(scaled[j]); |
1412 |
|
scaled[j] += 0.5; |
1413 |
+ |
// Handle the special case when an object is exactly on the |
1414 |
+ |
// boundary (a scaled coordinate of 1.0 is the same as |
1415 |
+ |
// scaled coordinate of 0.0) |
1416 |
+ |
if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
1417 |
|
} |
1418 |
|
|
1419 |
|
// find xyz-indices of cell that cutoffGroup is in. |