| 35 |
|
* |
| 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). |
| 38 |
> |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
| 39 |
|
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
| 40 |
|
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 308 |
|
|
| 309 |
|
void ForceMatrixDecomposition::createGtypeCutoffMap() { |
| 310 |
|
|
| 311 |
+ |
GrCut.clear(); |
| 312 |
+ |
GrCutSq.clear(); |
| 313 |
+ |
GrlistSq.clear(); |
| 314 |
+ |
|
| 315 |
|
RealType tol = 1e-6; |
| 316 |
|
largestRcut_ = 0.0; |
| 317 |
|
int atid; |
| 422 |
|
#endif |
| 423 |
|
|
| 424 |
|
RealType tradRcut = groupMax; |
| 425 |
+ |
|
| 426 |
+ |
GrCut.resize( gTypeCutoffs.size() ); |
| 427 |
+ |
GrCutSq.resize( gTypeCutoffs.size() ); |
| 428 |
+ |
GrlistSq.resize( gTypeCutoffs.size() ); |
| 429 |
|
|
| 430 |
+ |
|
| 431 |
|
for (unsigned int i = 0; i < gTypeCutoffs.size(); i++) { |
| 432 |
+ |
GrCut[i].resize( gTypeCutoffs.size() , 0.0); |
| 433 |
+ |
GrCutSq[i].resize( gTypeCutoffs.size(), 0.0 ); |
| 434 |
+ |
GrlistSq[i].resize( gTypeCutoffs.size(), 0.0 ); |
| 435 |
+ |
|
| 436 |
|
for (unsigned int j = 0; j < gTypeCutoffs.size(); j++) { |
| 437 |
|
RealType thisRcut; |
| 438 |
|
switch(cutoffPolicy_) { |
| 455 |
|
break; |
| 456 |
|
} |
| 457 |
|
|
| 458 |
< |
pair<int,int> key = make_pair(i,j); |
| 446 |
< |
gTypeCutoffMap[key].first = thisRcut; |
| 458 |
> |
GrCut[i][j] = thisRcut; |
| 459 |
|
if (thisRcut > largestRcut_) largestRcut_ = thisRcut; |
| 460 |
< |
gTypeCutoffMap[key].second = thisRcut*thisRcut; |
| 461 |
< |
gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); |
| 460 |
> |
GrCutSq[i][j] = thisRcut * thisRcut; |
| 461 |
> |
GrlistSq[i][j] = pow(thisRcut + skinThickness_, 2); |
| 462 |
> |
|
| 463 |
> |
// pair<int,int> key = make_pair(i,j); |
| 464 |
> |
// gTypeCutoffMap[key].first = thisRcut; |
| 465 |
> |
// gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); |
| 466 |
|
// sanity check |
| 467 |
|
|
| 468 |
|
if (userChoseCutoff_) { |
| 469 |
< |
if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { |
| 469 |
> |
if (abs(GrCut[i][j] - userCutoff_) > 0.0001) { |
| 470 |
|
sprintf(painCave.errMsg, |
| 471 |
|
"ForceMatrixDecomposition::createGtypeCutoffMap " |
| 472 |
|
"user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_); |
| 479 |
|
} |
| 480 |
|
} |
| 481 |
|
|
| 482 |
< |
groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { |
| 482 |
> |
void ForceMatrixDecomposition::getGroupCutoffs(int &cg1, int &cg2, RealType &rcut, RealType &rcutsq, RealType &rlistsq) { |
| 483 |
|
int i, j; |
| 484 |
|
#ifdef IS_MPI |
| 485 |
|
i = groupRowToGtype[cg1]; |
| 488 |
|
i = groupToGtype[cg1]; |
| 489 |
|
j = groupToGtype[cg2]; |
| 490 |
|
#endif |
| 491 |
< |
return gTypeCutoffMap[make_pair(i,j)]; |
| 491 |
> |
rcut = GrCut[i][j]; |
| 492 |
> |
rcutsq = GrCutSq[i][j]; |
| 493 |
> |
rlistsq = GrlistSq[i][j]; |
| 494 |
> |
return; |
| 495 |
> |
//return gTypeCutoffMap[make_pair(i,j)]; |
| 496 |
|
} |
| 497 |
|
|
| 498 |
|
int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { |
| 499 |
|
for (unsigned int j = 0; j < toposForAtom[atom1].size(); j++) { |
| 500 |
|
if (toposForAtom[atom1][j] == atom2) |
| 501 |
|
return topoDist[atom1][j]; |
| 502 |
< |
} |
| 502 |
> |
} |
| 503 |
|
return 0; |
| 504 |
|
} |
| 505 |
|
|
| 579 |
|
atomColData.electricField.end(), V3Zero); |
| 580 |
|
} |
| 581 |
|
|
| 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 |
– |
|
| 582 |
|
#endif |
| 583 |
|
// even in parallel, we need to zero out the local arrays: |
| 584 |
|
|
| 652 |
|
AtomPlanMatrixColumn->gather(snap_->atomData.aMat, |
| 653 |
|
atomColData.aMat); |
| 654 |
|
} |
| 655 |
< |
|
| 656 |
< |
// if needed, gather the atomic eletrostatic frames |
| 657 |
< |
if (storageLayout_ & DataStorage::dslElectroFrame) { |
| 658 |
< |
AtomPlanMatrixRow->gather(snap_->atomData.electroFrame, |
| 659 |
< |
atomRowData.electroFrame); |
| 660 |
< |
AtomPlanMatrixColumn->gather(snap_->atomData.electroFrame, |
| 661 |
< |
atomColData.electroFrame); |
| 655 |
> |
|
| 656 |
> |
// if needed, gather the atomic eletrostatic information |
| 657 |
> |
if (storageLayout_ & DataStorage::dslDipole) { |
| 658 |
> |
AtomPlanVectorRow->gather(snap_->atomData.dipole, |
| 659 |
> |
atomRowData.dipole); |
| 660 |
> |
AtomPlanVectorColumn->gather(snap_->atomData.dipole, |
| 661 |
> |
atomColData.dipole); |
| 662 |
|
} |
| 663 |
|
|
| 664 |
+ |
if (storageLayout_ & DataStorage::dslQuadrupole) { |
| 665 |
+ |
AtomPlanMatrixRow->gather(snap_->atomData.quadrupole, |
| 666 |
+ |
atomRowData.quadrupole); |
| 667 |
+ |
AtomPlanMatrixColumn->gather(snap_->atomData.quadrupole, |
| 668 |
+ |
atomColData.quadrupole); |
| 669 |
+ |
} |
| 670 |
+ |
|
| 671 |
|
// if needed, gather the atomic fluctuating charge values |
| 672 |
|
if (storageLayout_ & DataStorage::dslFlucQPosition) { |
| 673 |
|
AtomPlanRealRow->gather(snap_->atomData.flucQPos, |
| 699 |
|
snap_->atomData.density[i] += rho_tmp[i]; |
| 700 |
|
} |
| 701 |
|
|
| 702 |
+ |
// this isn't necessary if we don't have polarizable atoms, but |
| 703 |
+ |
// we'll leave it here for now. |
| 704 |
|
if (storageLayout_ & DataStorage::dslElectricField) { |
| 705 |
|
|
| 706 |
|
AtomPlanVectorRow->scatter(atomRowData.electricField, |
| 708 |
|
|
| 709 |
|
int n = snap_->atomData.electricField.size(); |
| 710 |
|
vector<Vector3d> field_tmp(n, V3Zero); |
| 711 |
< |
AtomPlanVectorColumn->scatter(atomColData.electricField, field_tmp); |
| 711 |
> |
AtomPlanVectorColumn->scatter(atomColData.electricField, |
| 712 |
> |
field_tmp); |
| 713 |
|
for (int i = 0; i < n; i++) |
| 714 |
|
snap_->atomData.electricField[i] += field_tmp[i]; |
| 715 |
|
} |
| 807 |
|
for (int i = 0; i < nq; i++) |
| 808 |
|
snap_->atomData.flucQFrc[i] += fqfrc_tmp[i]; |
| 809 |
|
|
| 810 |
+ |
} |
| 811 |
+ |
|
| 812 |
+ |
if (storageLayout_ & DataStorage::dslElectricField) { |
| 813 |
+ |
|
| 814 |
+ |
int nef = snap_->atomData.electricField.size(); |
| 815 |
+ |
vector<Vector3d> efield_tmp(nef, V3Zero); |
| 816 |
+ |
|
| 817 |
+ |
AtomPlanVectorRow->scatter(atomRowData.electricField, efield_tmp); |
| 818 |
+ |
for (int i = 0; i < nef; i++) { |
| 819 |
+ |
snap_->atomData.electricField[i] += efield_tmp[i]; |
| 820 |
+ |
efield_tmp[i] = 0.0; |
| 821 |
+ |
} |
| 822 |
+ |
|
| 823 |
+ |
AtomPlanVectorColumn->scatter(atomColData.electricField, efield_tmp); |
| 824 |
+ |
for (int i = 0; i < nef; i++) |
| 825 |
+ |
snap_->atomData.electricField[i] += efield_tmp[i]; |
| 826 |
|
} |
| 827 |
|
|
| 828 |
+ |
|
| 829 |
|
nLocal_ = snap_->getNumberOfAtoms(); |
| 830 |
|
|
| 831 |
|
vector<potVec> pot_temp(nLocal_, |
| 959 |
|
|
| 960 |
|
|
| 961 |
|
|
| 962 |
< |
int ForceMatrixDecomposition::getNAtomsInRow() { |
| 962 |
> |
int& ForceMatrixDecomposition::getNAtomsInRow() { |
| 963 |
|
#ifdef IS_MPI |
| 964 |
|
return nAtomsInRow_; |
| 965 |
|
#else |
| 970 |
|
/** |
| 971 |
|
* returns the list of atoms belonging to this group. |
| 972 |
|
*/ |
| 973 |
< |
vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ |
| 973 |
> |
vector<int>& ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ |
| 974 |
|
#ifdef IS_MPI |
| 975 |
|
return groupListRow_[cg1]; |
| 976 |
|
#else |
| 978 |
|
#endif |
| 979 |
|
} |
| 980 |
|
|
| 981 |
< |
vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ |
| 981 |
> |
vector<int>& ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ |
| 982 |
|
#ifdef IS_MPI |
| 983 |
|
return groupListCol_[cg2]; |
| 984 |
|
#else |
| 995 |
|
d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1]; |
| 996 |
|
#endif |
| 997 |
|
|
| 998 |
< |
snap_->wrapVector(d); |
| 998 |
> |
if (usePeriodicBoundaryConditions_) { |
| 999 |
> |
snap_->wrapVector(d); |
| 1000 |
> |
} |
| 1001 |
|
return d; |
| 1002 |
|
} |
| 1003 |
|
|
| 1004 |
< |
Vector3d ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){ |
| 1004 |
> |
Vector3d& ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){ |
| 1005 |
|
#ifdef IS_MPI |
| 1006 |
|
return cgColData.velocity[cg2]; |
| 1007 |
|
#else |
| 1009 |
|
#endif |
| 1010 |
|
} |
| 1011 |
|
|
| 1012 |
< |
Vector3d ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ |
| 1012 |
> |
Vector3d& ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ |
| 1013 |
|
#ifdef IS_MPI |
| 1014 |
|
return atomColData.velocity[atom2]; |
| 1015 |
|
#else |
| 1027 |
|
#else |
| 1028 |
|
d = snap_->cgData.position[cg1] - snap_->atomData.position[atom1]; |
| 1029 |
|
#endif |
| 1030 |
< |
|
| 1031 |
< |
snap_->wrapVector(d); |
| 1030 |
> |
if (usePeriodicBoundaryConditions_) { |
| 1031 |
> |
snap_->wrapVector(d); |
| 1032 |
> |
} |
| 1033 |
|
return d; |
| 1034 |
|
} |
| 1035 |
|
|
| 1041 |
|
#else |
| 1042 |
|
d = snap_->cgData.position[cg2] - snap_->atomData.position[atom2]; |
| 1043 |
|
#endif |
| 1044 |
< |
|
| 1045 |
< |
snap_->wrapVector(d); |
| 1044 |
> |
if (usePeriodicBoundaryConditions_) { |
| 1045 |
> |
snap_->wrapVector(d); |
| 1046 |
> |
} |
| 1047 |
|
return d; |
| 1048 |
|
} |
| 1049 |
|
|
| 1050 |
< |
RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) { |
| 1050 |
> |
RealType& ForceMatrixDecomposition::getMassFactorRow(int atom1) { |
| 1051 |
|
#ifdef IS_MPI |
| 1052 |
|
return massFactorsRow[atom1]; |
| 1053 |
|
#else |
| 1055 |
|
#endif |
| 1056 |
|
} |
| 1057 |
|
|
| 1058 |
< |
RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { |
| 1058 |
> |
RealType& ForceMatrixDecomposition::getMassFactorColumn(int atom2) { |
| 1059 |
|
#ifdef IS_MPI |
| 1060 |
|
return massFactorsCol[atom2]; |
| 1061 |
|
#else |
| 1072 |
|
#else |
| 1073 |
|
d = snap_->atomData.position[atom2] - snap_->atomData.position[atom1]; |
| 1074 |
|
#endif |
| 1075 |
< |
|
| 1076 |
< |
snap_->wrapVector(d); |
| 1075 |
> |
if (usePeriodicBoundaryConditions_) { |
| 1076 |
> |
snap_->wrapVector(d); |
| 1077 |
> |
} |
| 1078 |
|
return d; |
| 1079 |
|
} |
| 1080 |
|
|
| 1081 |
< |
vector<int> ForceMatrixDecomposition::getExcludesForAtom(int atom1) { |
| 1081 |
> |
vector<int>& ForceMatrixDecomposition::getExcludesForAtom(int atom1) { |
| 1082 |
|
return excludesForAtom[atom1]; |
| 1083 |
|
} |
| 1084 |
|
|
| 1169 |
|
|
| 1170 |
|
#ifdef IS_MPI |
| 1171 |
|
idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]); |
| 1172 |
+ |
idat.atid1 = identsRow[atom1]; |
| 1173 |
+ |
idat.atid2 = identsCol[atom2]; |
| 1174 |
|
//idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
| 1175 |
|
// ff_->getAtomType(identsCol[atom2]) ); |
| 1176 |
|
|
| 1179 |
|
idat.A2 = &(atomColData.aMat[atom2]); |
| 1180 |
|
} |
| 1181 |
|
|
| 1135 |
– |
if (storageLayout_ & DataStorage::dslElectroFrame) { |
| 1136 |
– |
idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
| 1137 |
– |
idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
| 1138 |
– |
} |
| 1139 |
– |
|
| 1182 |
|
if (storageLayout_ & DataStorage::dslTorque) { |
| 1183 |
|
idat.t1 = &(atomRowData.torque[atom1]); |
| 1184 |
|
idat.t2 = &(atomColData.torque[atom2]); |
| 1185 |
|
} |
| 1186 |
|
|
| 1187 |
+ |
if (storageLayout_ & DataStorage::dslDipole) { |
| 1188 |
+ |
idat.dipole1 = &(atomRowData.dipole[atom1]); |
| 1189 |
+ |
idat.dipole2 = &(atomColData.dipole[atom2]); |
| 1190 |
+ |
} |
| 1191 |
+ |
|
| 1192 |
+ |
if (storageLayout_ & DataStorage::dslQuadrupole) { |
| 1193 |
+ |
idat.quadrupole1 = &(atomRowData.quadrupole[atom1]); |
| 1194 |
+ |
idat.quadrupole2 = &(atomColData.quadrupole[atom2]); |
| 1195 |
+ |
} |
| 1196 |
+ |
|
| 1197 |
|
if (storageLayout_ & DataStorage::dslDensity) { |
| 1198 |
|
idat.rho1 = &(atomRowData.density[atom1]); |
| 1199 |
|
idat.rho2 = &(atomColData.density[atom2]); |
| 1227 |
|
#else |
| 1228 |
|
|
| 1229 |
|
idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); |
| 1230 |
+ |
idat.atid1 = idents[atom1]; |
| 1231 |
+ |
idat.atid2 = idents[atom2]; |
| 1232 |
|
|
| 1233 |
|
if (storageLayout_ & DataStorage::dslAmat) { |
| 1234 |
|
idat.A1 = &(snap_->atomData.aMat[atom1]); |
| 1235 |
|
idat.A2 = &(snap_->atomData.aMat[atom2]); |
| 1236 |
|
} |
| 1237 |
|
|
| 1184 |
– |
if (storageLayout_ & DataStorage::dslElectroFrame) { |
| 1185 |
– |
idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
| 1186 |
– |
idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
| 1187 |
– |
} |
| 1188 |
– |
|
| 1238 |
|
if (storageLayout_ & DataStorage::dslTorque) { |
| 1239 |
|
idat.t1 = &(snap_->atomData.torque[atom1]); |
| 1240 |
|
idat.t2 = &(snap_->atomData.torque[atom2]); |
| 1241 |
|
} |
| 1242 |
|
|
| 1243 |
+ |
if (storageLayout_ & DataStorage::dslDipole) { |
| 1244 |
+ |
idat.dipole1 = &(snap_->atomData.dipole[atom1]); |
| 1245 |
+ |
idat.dipole2 = &(snap_->atomData.dipole[atom2]); |
| 1246 |
+ |
} |
| 1247 |
+ |
|
| 1248 |
+ |
if (storageLayout_ & DataStorage::dslQuadrupole) { |
| 1249 |
+ |
idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]); |
| 1250 |
+ |
idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]); |
| 1251 |
+ |
} |
| 1252 |
+ |
|
| 1253 |
|
if (storageLayout_ & DataStorage::dslDensity) { |
| 1254 |
|
idat.rho1 = &(snap_->atomData.density[atom1]); |
| 1255 |
|
idat.rho2 = &(snap_->atomData.density[atom2]); |
| 1340 |
|
* first element of pair is row-indexed CutoffGroup |
| 1341 |
|
* second element of pair is column-indexed CutoffGroup |
| 1342 |
|
*/ |
| 1343 |
< |
vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() { |
| 1344 |
< |
|
| 1345 |
< |
vector<pair<int, int> > neighborList; |
| 1343 |
> |
void ForceMatrixDecomposition::buildNeighborList(vector<pair<int,int> >& neighborList) { |
| 1344 |
> |
|
| 1345 |
> |
neighborList.clear(); |
| 1346 |
|
groupCutoffs cuts; |
| 1347 |
|
bool doAllPairs = false; |
| 1348 |
|
|
| 1349 |
+ |
RealType rList_ = (largestRcut_ + skinThickness_); |
| 1350 |
+ |
RealType rcut, rcutsq, rlistsq; |
| 1351 |
+ |
Snapshot* snap_ = sman_->getCurrentSnapshot(); |
| 1352 |
+ |
Mat3x3d box; |
| 1353 |
+ |
Mat3x3d invBox; |
| 1354 |
+ |
|
| 1355 |
+ |
Vector3d rs, scaled, dr; |
| 1356 |
+ |
Vector3i whichCell; |
| 1357 |
+ |
int cellIndex; |
| 1358 |
+ |
|
| 1359 |
|
#ifdef IS_MPI |
| 1360 |
|
cellListRow_.clear(); |
| 1361 |
|
cellListCol_.clear(); |
| 1362 |
|
#else |
| 1363 |
|
cellList_.clear(); |
| 1364 |
|
#endif |
| 1365 |
< |
|
| 1366 |
< |
RealType rList_ = (largestRcut_ + skinThickness_); |
| 1367 |
< |
Snapshot* snap_ = sman_->getCurrentSnapshot(); |
| 1368 |
< |
Mat3x3d Hmat = snap_->getHmat(); |
| 1369 |
< |
Vector3d Hx = Hmat.getColumn(0); |
| 1370 |
< |
Vector3d Hy = Hmat.getColumn(1); |
| 1371 |
< |
Vector3d Hz = Hmat.getColumn(2); |
| 1372 |
< |
|
| 1373 |
< |
nCells_.x() = (int) ( Hx.length() )/ rList_; |
| 1374 |
< |
nCells_.y() = (int) ( Hy.length() )/ rList_; |
| 1375 |
< |
nCells_.z() = (int) ( Hz.length() )/ rList_; |
| 1376 |
< |
|
| 1365 |
> |
|
| 1366 |
> |
if (!usePeriodicBoundaryConditions_) { |
| 1367 |
> |
box = snap_->getBoundingBox(); |
| 1368 |
> |
invBox = snap_->getInvBoundingBox(); |
| 1369 |
> |
} else { |
| 1370 |
> |
box = snap_->getHmat(); |
| 1371 |
> |
invBox = snap_->getInvHmat(); |
| 1372 |
> |
} |
| 1373 |
> |
|
| 1374 |
> |
Vector3d boxX = box.getColumn(0); |
| 1375 |
> |
Vector3d boxY = box.getColumn(1); |
| 1376 |
> |
Vector3d boxZ = box.getColumn(2); |
| 1377 |
> |
|
| 1378 |
> |
nCells_.x() = (int) ( boxX.length() )/ rList_; |
| 1379 |
> |
nCells_.y() = (int) ( boxY.length() )/ rList_; |
| 1380 |
> |
nCells_.z() = (int) ( boxZ.length() )/ rList_; |
| 1381 |
> |
|
| 1382 |
|
// handle small boxes where the cell offsets can end up repeating cells |
| 1383 |
|
|
| 1384 |
|
if (nCells_.x() < 3) doAllPairs = true; |
| 1385 |
|
if (nCells_.y() < 3) doAllPairs = true; |
| 1386 |
|
if (nCells_.z() < 3) doAllPairs = true; |
| 1387 |
< |
|
| 1314 |
< |
Mat3x3d invHmat = snap_->getInvHmat(); |
| 1315 |
< |
Vector3d rs, scaled, dr; |
| 1316 |
< |
Vector3i whichCell; |
| 1317 |
< |
int cellIndex; |
| 1387 |
> |
|
| 1388 |
|
int nCtot = nCells_.x() * nCells_.y() * nCells_.z(); |
| 1389 |
< |
|
| 1389 |
> |
|
| 1390 |
|
#ifdef IS_MPI |
| 1391 |
|
cellListRow_.resize(nCtot); |
| 1392 |
|
cellListCol_.resize(nCtot); |
| 1393 |
|
#else |
| 1394 |
|
cellList_.resize(nCtot); |
| 1395 |
|
#endif |
| 1396 |
< |
|
| 1396 |
> |
|
| 1397 |
|
if (!doAllPairs) { |
| 1398 |
|
#ifdef IS_MPI |
| 1399 |
< |
|
| 1399 |
> |
|
| 1400 |
|
for (int i = 0; i < nGroupsInRow_; i++) { |
| 1401 |
|
rs = cgRowData.position[i]; |
| 1402 |
|
|
| 1403 |
|
// scaled positions relative to the box vectors |
| 1404 |
< |
scaled = invHmat * rs; |
| 1404 |
> |
scaled = invBox * rs; |
| 1405 |
|
|
| 1406 |
|
// wrap the vector back into the unit box by subtracting integer box |
| 1407 |
|
// numbers |
| 1429 |
|
rs = cgColData.position[i]; |
| 1430 |
|
|
| 1431 |
|
// scaled positions relative to the box vectors |
| 1432 |
< |
scaled = invHmat * rs; |
| 1432 |
> |
scaled = invBox * rs; |
| 1433 |
|
|
| 1434 |
|
// wrap the vector back into the unit box by subtracting integer box |
| 1435 |
|
// numbers |
| 1453 |
|
// add this cutoff group to the list of groups in this cell; |
| 1454 |
|
cellListCol_[cellIndex].push_back(i); |
| 1455 |
|
} |
| 1456 |
< |
|
| 1456 |
> |
|
| 1457 |
|
#else |
| 1458 |
|
for (int i = 0; i < nGroups_; i++) { |
| 1459 |
|
rs = snap_->cgData.position[i]; |
| 1460 |
|
|
| 1461 |
|
// scaled positions relative to the box vectors |
| 1462 |
< |
scaled = invHmat * rs; |
| 1462 |
> |
scaled = invBox * rs; |
| 1463 |
|
|
| 1464 |
|
// wrap the vector back into the unit box by subtracting integer box |
| 1465 |
|
// numbers |
| 1528 |
|
// & column indicies and will divide labor in the |
| 1529 |
|
// force evaluation later. |
| 1530 |
|
dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
| 1531 |
< |
snap_->wrapVector(dr); |
| 1532 |
< |
cuts = getGroupCutoffs( (*j1), (*j2) ); |
| 1533 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 1531 |
> |
if (usePeriodicBoundaryConditions_) { |
| 1532 |
> |
snap_->wrapVector(dr); |
| 1533 |
> |
} |
| 1534 |
> |
getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); |
| 1535 |
> |
if (dr.lengthSquare() < rlistsq) { |
| 1536 |
|
neighborList.push_back(make_pair((*j1), (*j2))); |
| 1537 |
|
} |
| 1538 |
|
} |
| 1552 |
|
// allows atoms within a single cutoff group to |
| 1553 |
|
// interact with each other. |
| 1554 |
|
|
| 1483 |
– |
|
| 1484 |
– |
|
| 1555 |
|
if (m2 != m1 || (*j2) >= (*j1) ) { |
| 1556 |
|
|
| 1557 |
|
dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
| 1558 |
< |
snap_->wrapVector(dr); |
| 1559 |
< |
cuts = getGroupCutoffs( (*j1), (*j2) ); |
| 1560 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 1558 |
> |
if (usePeriodicBoundaryConditions_) { |
| 1559 |
> |
snap_->wrapVector(dr); |
| 1560 |
> |
} |
| 1561 |
> |
getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); |
| 1562 |
> |
if (dr.lengthSquare() < rlistsq) { |
| 1563 |
|
neighborList.push_back(make_pair((*j1), (*j2))); |
| 1564 |
|
} |
| 1565 |
|
} |
| 1576 |
|
for (int j1 = 0; j1 < nGroupsInRow_; j1++) { |
| 1577 |
|
for (int j2 = 0; j2 < nGroupsInCol_; j2++) { |
| 1578 |
|
dr = cgColData.position[j2] - cgRowData.position[j1]; |
| 1579 |
< |
snap_->wrapVector(dr); |
| 1580 |
< |
cuts = getGroupCutoffs( j1, j2 ); |
| 1581 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 1579 |
> |
if (usePeriodicBoundaryConditions_) { |
| 1580 |
> |
snap_->wrapVector(dr); |
| 1581 |
> |
} |
| 1582 |
> |
getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq); |
| 1583 |
> |
if (dr.lengthSquare() < rlistsq) { |
| 1584 |
|
neighborList.push_back(make_pair(j1, j2)); |
| 1585 |
|
} |
| 1586 |
|
} |
| 1591 |
|
// include self group interactions j2 == j1 |
| 1592 |
|
for (int j2 = j1; j2 < nGroups_; j2++) { |
| 1593 |
|
dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; |
| 1594 |
< |
snap_->wrapVector(dr); |
| 1595 |
< |
cuts = getGroupCutoffs( j1, j2 ); |
| 1596 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 1594 |
> |
if (usePeriodicBoundaryConditions_) { |
| 1595 |
> |
snap_->wrapVector(dr); |
| 1596 |
> |
} |
| 1597 |
> |
getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq ); |
| 1598 |
> |
if (dr.lengthSquare() < rlistsq) { |
| 1599 |
|
neighborList.push_back(make_pair(j1, j2)); |
| 1600 |
|
} |
| 1601 |
|
} |
| 1608 |
|
saved_CG_positions_.clear(); |
| 1609 |
|
for (int i = 0; i < nGroups_; i++) |
| 1610 |
|
saved_CG_positions_.push_back(snap_->cgData.position[i]); |
| 1535 |
– |
|
| 1536 |
– |
return neighborList; |
| 1611 |
|
} |
| 1612 |
|
} //end namespace OpenMD |