| 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 |