| 99 |  | nGroups_ = info_->getNLocalCutoffGroups(); | 
| 100 |  | // gather the information for atomtype IDs (atids): | 
| 101 |  | idents = info_->getIdentArray(); | 
| 102 | + | regions = info_->getRegions(); | 
| 103 |  | AtomLocalToGlobal = info_->getGlobalAtomIndices(); | 
| 104 |  | cgLocalToGlobal = info_->getGlobalGroupIndices(); | 
| 105 |  | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); | 
| 164 |  |  | 
| 165 |  | AtomPlanIntRow->gather(idents, identsRow); | 
| 166 |  | AtomPlanIntColumn->gather(idents, identsCol); | 
| 167 | + |  | 
| 168 | + | regionsRow.resize(nAtomsInRow_); | 
| 169 | + | regionsCol.resize(nAtomsInCol_); | 
| 170 | + |  | 
| 171 | + | AtomPlanIntRow->gather(regions, regionsRow); | 
| 172 | + | AtomPlanIntColumn->gather(regions, regionsCol); | 
| 173 |  |  | 
| 174 |  | // allocate memory for the parallel objects | 
| 175 |  | atypesRow.resize(nAtomsInRow_); | 
| 315 |  |  | 
| 316 |  | void ForceMatrixDecomposition::createGtypeCutoffMap() { | 
| 317 |  |  | 
| 318 | + | GrCut.clear(); | 
| 319 | + | GrCutSq.clear(); | 
| 320 | + | GrlistSq.clear(); | 
| 321 | + |  | 
| 322 |  | RealType tol = 1e-6; | 
| 323 |  | largestRcut_ = 0.0; | 
| 324 |  | int atid; | 
| 430 |  |  | 
| 431 |  | RealType tradRcut = groupMax; | 
| 432 |  |  | 
| 433 | + | GrCut.resize( gTypeCutoffs.size() ); | 
| 434 | + | GrCutSq.resize( gTypeCutoffs.size() ); | 
| 435 | + | GrlistSq.resize( gTypeCutoffs.size() ); | 
| 436 | + |  | 
| 437 | + |  | 
| 438 |  | for (unsigned int i = 0; i < gTypeCutoffs.size();  i++) { | 
| 439 | + | GrCut[i].resize( gTypeCutoffs.size() , 0.0); | 
| 440 | + | GrCutSq[i].resize( gTypeCutoffs.size(), 0.0 ); | 
| 441 | + | GrlistSq[i].resize( gTypeCutoffs.size(), 0.0 ); | 
| 442 | + |  | 
| 443 |  | for (unsigned int j = 0; j < gTypeCutoffs.size();  j++) { | 
| 444 |  | RealType thisRcut; | 
| 445 |  | switch(cutoffPolicy_) { | 
| 462 |  | break; | 
| 463 |  | } | 
| 464 |  |  | 
| 465 | < | pair<int,int> key = make_pair(i,j); | 
| 446 | < | gTypeCutoffMap[key].first = thisRcut; | 
| 465 | > | GrCut[i][j] = thisRcut; | 
| 466 |  | if (thisRcut > largestRcut_) largestRcut_ = thisRcut; | 
| 467 | < | gTypeCutoffMap[key].second = thisRcut*thisRcut; | 
| 468 | < | gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); | 
| 467 | > | GrCutSq[i][j] = thisRcut * thisRcut; | 
| 468 | > | GrlistSq[i][j] = pow(thisRcut + skinThickness_, 2); | 
| 469 | > |  | 
| 470 | > | // pair<int,int> key = make_pair(i,j); | 
| 471 | > | // gTypeCutoffMap[key].first = thisRcut; | 
| 472 | > | // gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2); | 
| 473 |  | // sanity check | 
| 474 |  |  | 
| 475 |  | if (userChoseCutoff_) { | 
| 476 | < | if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { | 
| 476 | > | if (abs(GrCut[i][j] - userCutoff_) > 0.0001) { | 
| 477 |  | sprintf(painCave.errMsg, | 
| 478 |  | "ForceMatrixDecomposition::createGtypeCutoffMap " | 
| 479 |  | "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_); | 
| 486 |  | } | 
| 487 |  | } | 
| 488 |  |  | 
| 489 | < | groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { | 
| 489 | > | void ForceMatrixDecomposition::getGroupCutoffs(int &cg1, int &cg2, RealType &rcut, RealType &rcutsq, RealType &rlistsq) { | 
| 490 |  | int i, j; | 
| 491 |  | #ifdef IS_MPI | 
| 492 |  | i = groupRowToGtype[cg1]; | 
| 495 |  | i = groupToGtype[cg1]; | 
| 496 |  | j = groupToGtype[cg2]; | 
| 497 |  | #endif | 
| 498 | < | return gTypeCutoffMap[make_pair(i,j)]; | 
| 498 | > | rcut = GrCut[i][j]; | 
| 499 | > | rcutsq = GrCutSq[i][j]; | 
| 500 | > | rlistsq = GrlistSq[i][j]; | 
| 501 | > | return; | 
| 502 | > | //return gTypeCutoffMap[make_pair(i,j)]; | 
| 503 |  | } | 
| 504 |  |  | 
| 505 |  | int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { | 
| 1178 |  | idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]); | 
| 1179 |  | idat.atid1 = identsRow[atom1]; | 
| 1180 |  | idat.atid2 = identsCol[atom2]; | 
| 1181 | + |  | 
| 1182 | + | if (regionsRow[atom1] >= 0 && regionsCol[atom2] >= 0) | 
| 1183 | + | idat.sameRegion = (regionsRow[atom1] == regionsCol[atom2]); | 
| 1184 | + |  | 
| 1185 |  | //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | 
| 1186 |  | //                         ff_->getAtomType(identsCol[atom2]) ); | 
| 1187 |  |  | 
| 1241 |  | idat.atid1 = idents[atom1]; | 
| 1242 |  | idat.atid2 = idents[atom2]; | 
| 1243 |  |  | 
| 1244 | + | if (regions[atom1] >= 0 && regions[atom2] >= 0) | 
| 1245 | + | idat.sameRegion = (regions[atom1] == regions[atom2]); | 
| 1246 | + |  | 
| 1247 |  | if (storageLayout_ & DataStorage::dslAmat) { | 
| 1248 |  | idat.A1 = &(snap_->atomData.aMat[atom1]); | 
| 1249 |  | idat.A2 = &(snap_->atomData.aMat[atom2]); | 
| 1361 |  | bool doAllPairs = false; | 
| 1362 |  |  | 
| 1363 |  | RealType rList_ = (largestRcut_ + skinThickness_); | 
| 1364 | + | RealType rcut, rcutsq, rlistsq; | 
| 1365 |  | Snapshot* snap_ = sman_->getCurrentSnapshot(); | 
| 1366 |  | Mat3x3d box; | 
| 1367 |  | Mat3x3d invBox; | 
| 1545 |  | if (usePeriodicBoundaryConditions_) { | 
| 1546 |  | snap_->wrapVector(dr); | 
| 1547 |  | } | 
| 1548 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); | 
| 1549 | < | if (dr.lengthSquare() < cuts.third) { | 
| 1548 | > | getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); | 
| 1549 | > | if (dr.lengthSquare() < rlistsq) { | 
| 1550 |  | neighborList.push_back(make_pair((*j1), (*j2))); | 
| 1551 |  | } | 
| 1552 |  | } | 
| 1572 |  | if (usePeriodicBoundaryConditions_) { | 
| 1573 |  | snap_->wrapVector(dr); | 
| 1574 |  | } | 
| 1575 | < | cuts = getGroupCutoffs( (*j1), (*j2) ); | 
| 1576 | < | if (dr.lengthSquare() < cuts.third) { | 
| 1575 | > | getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq ); | 
| 1576 | > | if (dr.lengthSquare() < rlistsq) { | 
| 1577 |  | neighborList.push_back(make_pair((*j1), (*j2))); | 
| 1578 |  | } | 
| 1579 |  | } | 
| 1593 |  | if (usePeriodicBoundaryConditions_) { | 
| 1594 |  | snap_->wrapVector(dr); | 
| 1595 |  | } | 
| 1596 | < | cuts = getGroupCutoffs( j1, j2 ); | 
| 1597 | < | if (dr.lengthSquare() < cuts.third) { | 
| 1596 | > | getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq); | 
| 1597 | > | if (dr.lengthSquare() < rlistsq) { | 
| 1598 |  | neighborList.push_back(make_pair(j1, j2)); | 
| 1599 |  | } | 
| 1600 |  | } | 
| 1608 |  | if (usePeriodicBoundaryConditions_) { | 
| 1609 |  | snap_->wrapVector(dr); | 
| 1610 |  | } | 
| 1611 | < | cuts = getGroupCutoffs( j1, j2 ); | 
| 1612 | < | if (dr.lengthSquare() < cuts.third) { | 
| 1611 | > | getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq ); | 
| 1612 | > | if (dr.lengthSquare() < rlistsq) { | 
| 1613 |  | neighborList.push_back(make_pair(j1, j2)); | 
| 1614 |  | } | 
| 1615 |  | } |