| 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) { |
| 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 |
|
|
| 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]); |
| 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; |
| 1531 |
|
if (usePeriodicBoundaryConditions_) { |
| 1532 |
|
snap_->wrapVector(dr); |
| 1533 |
|
} |
| 1534 |
< |
cuts = getGroupCutoffs( (*j1), (*j2) ); |
| 1535 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 1534 |
> |
getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); |
| 1535 |
> |
if (dr.lengthSquare() < rlistsq) { |
| 1536 |
|
neighborList.push_back(make_pair((*j1), (*j2))); |
| 1537 |
|
} |
| 1538 |
|
} |
| 1558 |
|
if (usePeriodicBoundaryConditions_) { |
| 1559 |
|
snap_->wrapVector(dr); |
| 1560 |
|
} |
| 1561 |
< |
cuts = getGroupCutoffs( (*j1), (*j2) ); |
| 1562 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 1561 |
> |
getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); |
| 1562 |
> |
if (dr.lengthSquare() < rlistsq) { |
| 1563 |
|
neighborList.push_back(make_pair((*j1), (*j2))); |
| 1564 |
|
} |
| 1565 |
|
} |
| 1579 |
|
if (usePeriodicBoundaryConditions_) { |
| 1580 |
|
snap_->wrapVector(dr); |
| 1581 |
|
} |
| 1582 |
< |
cuts = getGroupCutoffs( j1, j2 ); |
| 1583 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 1582 |
> |
getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq); |
| 1583 |
> |
if (dr.lengthSquare() < rlistsq) { |
| 1584 |
|
neighborList.push_back(make_pair(j1, j2)); |
| 1585 |
|
} |
| 1586 |
|
} |
| 1594 |
|
if (usePeriodicBoundaryConditions_) { |
| 1595 |
|
snap_->wrapVector(dr); |
| 1596 |
|
} |
| 1597 |
< |
cuts = getGroupCutoffs( j1, j2 ); |
| 1598 |
< |
if (dr.lengthSquare() < cuts.third) { |
| 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]); |
| 1586 |
– |
|
| 1587 |
– |
return neighborList; |
| 1611 |
|
} |
| 1612 |
|
} //end namespace OpenMD |