--- trunk/src/parallel/ForceMatrixDecomposition.cpp	2013/06/19 17:19:07	1893
+++ trunk/src/parallel/ForceMatrixDecomposition.cpp	2013/08/19 13:12:00	1929
@@ -99,6 +99,7 @@ namespace OpenMD {
     nGroups_ = info_->getNLocalCutoffGroups();
     // gather the information for atomtype IDs (atids):
     idents = info_->getIdentArray();
+    regions = info_->getRegions();
     AtomLocalToGlobal = info_->getGlobalAtomIndices();
     cgLocalToGlobal = info_->getGlobalGroupIndices();
     vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
@@ -163,6 +164,12 @@ namespace OpenMD {
     
     AtomPlanIntRow->gather(idents, identsRow);
     AtomPlanIntColumn->gather(idents, identsCol);
+
+    regionsRow.resize(nAtomsInRow_);
+    regionsCol.resize(nAtomsInCol_);
+    
+    AtomPlanIntRow->gather(regions, regionsRow);
+    AtomPlanIntColumn->gather(regions, regionsCol);
     
     // allocate memory for the parallel objects
     atypesRow.resize(nAtomsInRow_);
@@ -308,6 +315,10 @@ namespace OpenMD {
    
   void ForceMatrixDecomposition::createGtypeCutoffMap() {
     
+    GrCut.clear();
+    GrCutSq.clear();
+    GrlistSq.clear();
+
     RealType tol = 1e-6;
     largestRcut_ = 0.0;
     int atid;
@@ -419,7 +430,16 @@ namespace OpenMD {
     
     RealType tradRcut = groupMax;
 
+    GrCut.resize( gTypeCutoffs.size() );
+    GrCutSq.resize( gTypeCutoffs.size() );
+    GrlistSq.resize( gTypeCutoffs.size() );
+
+
     for (unsigned int i = 0; i < gTypeCutoffs.size();  i++) {
+      GrCut[i].resize( gTypeCutoffs.size() , 0.0);
+      GrCutSq[i].resize( gTypeCutoffs.size(), 0.0 );
+      GrlistSq[i].resize( gTypeCutoffs.size(), 0.0 );
+
       for (unsigned int j = 0; j < gTypeCutoffs.size();  j++) {       
         RealType thisRcut;
         switch(cutoffPolicy_) {
@@ -442,15 +462,18 @@ namespace OpenMD {
           break;
         }
 
-        pair<int,int> key = make_pair(i,j);
-        gTypeCutoffMap[key].first = thisRcut;
+        GrCut[i][j] = thisRcut;
         if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
-        gTypeCutoffMap[key].second = thisRcut*thisRcut;
-        gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
+        GrCutSq[i][j] = thisRcut * thisRcut;
+        GrlistSq[i][j] = pow(thisRcut + skinThickness_, 2);
+
+        // pair<int,int> key = make_pair(i,j);
+        // gTypeCutoffMap[key].first = thisRcut;
+        // gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
         // sanity check
         
         if (userChoseCutoff_) {
-          if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) {
+          if (abs(GrCut[i][j] - userCutoff_) > 0.0001) {
             sprintf(painCave.errMsg,
                     "ForceMatrixDecomposition::createGtypeCutoffMap " 
                     "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_);
@@ -463,7 +486,7 @@ namespace OpenMD {
     }
   }
 
-  groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
+  void ForceMatrixDecomposition::getGroupCutoffs(int &cg1, int &cg2, RealType &rcut, RealType &rcutsq, RealType &rlistsq) {
     int i, j;   
 #ifdef IS_MPI
     i = groupRowToGtype[cg1];
@@ -472,7 +495,11 @@ namespace OpenMD {
     i = groupToGtype[cg1];
     j = groupToGtype[cg2];
 #endif    
-    return gTypeCutoffMap[make_pair(i,j)];
+    rcut = GrCut[i][j];
+    rcutsq = GrCutSq[i][j];
+    rlistsq = GrlistSq[i][j];
+    return;
+    //return gTypeCutoffMap[make_pair(i,j)];
   }
 
   int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) {
@@ -1149,6 +1176,12 @@ namespace OpenMD {
    
 #ifdef IS_MPI
     idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]);
+    idat.atid1 = identsRow[atom1];
+    idat.atid2 = identsCol[atom2];
+
+    if (regionsRow[atom1] >= 0 && regionsCol[atom2] >= 0) 
+      idat.sameRegion = (regionsRow[atom1] == regionsCol[atom2]);
+       
     //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), 
     //                         ff_->getAtomType(identsCol[atom2]) );
     
@@ -1205,7 +1238,12 @@ namespace OpenMD {
 #else
     
     idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
+    idat.atid1 = idents[atom1];
+    idat.atid2 = idents[atom2];
 
+    if (regions[atom1] >= 0 && regions[atom2] >= 0) 
+      idat.sameRegion = (regions[atom1] == regions[atom2]);
+
     if (storageLayout_ & DataStorage::dslAmat) {
       idat.A1 = &(snap_->atomData.aMat[atom1]);
       idat.A2 = &(snap_->atomData.aMat[atom2]);
@@ -1316,13 +1354,14 @@ namespace OpenMD {
    * first element of pair is row-indexed CutoffGroup
    * second element of pair is column-indexed CutoffGroup
    */
-  vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
-      
-    vector<pair<int, int> > neighborList;
+  void ForceMatrixDecomposition::buildNeighborList(vector<pair<int,int> >& neighborList) {
+    
+    neighborList.clear();
     groupCutoffs cuts;
     bool doAllPairs = false;
 
     RealType rList_ = (largestRcut_ + skinThickness_);
+    RealType rcut, rcutsq, rlistsq;
     Snapshot* snap_ = sman_->getCurrentSnapshot();
     Mat3x3d box;
     Mat3x3d invBox;
@@ -1506,8 +1545,8 @@ namespace OpenMD {
                   if (usePeriodicBoundaryConditions_) {
                     snap_->wrapVector(dr);
                   }
-                  cuts = getGroupCutoffs( (*j1), (*j2) );
-                  if (dr.lengthSquare() < cuts.third) {
+                  getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq );
+                  if (dr.lengthSquare() < rlistsq) {
                     neighborList.push_back(make_pair((*j1), (*j2)));
                   }                  
                 }
@@ -1533,8 +1572,8 @@ namespace OpenMD {
                     if (usePeriodicBoundaryConditions_) {
                       snap_->wrapVector(dr);
                     }
-                    cuts = getGroupCutoffs( (*j1), (*j2) );
-                    if (dr.lengthSquare() < cuts.third) {
+                    getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq );
+                    if (dr.lengthSquare() < rlistsq) {
                       neighborList.push_back(make_pair((*j1), (*j2)));
                     }
                   }
@@ -1554,8 +1593,8 @@ namespace OpenMD {
           if (usePeriodicBoundaryConditions_) {
             snap_->wrapVector(dr);
           }
-          cuts = getGroupCutoffs( j1, j2 );
-          if (dr.lengthSquare() < cuts.third) {
+          getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq);
+          if (dr.lengthSquare() < rlistsq) {
             neighborList.push_back(make_pair(j1, j2));
           }
         }
@@ -1569,8 +1608,8 @@ namespace OpenMD {
           if (usePeriodicBoundaryConditions_) {
             snap_->wrapVector(dr);
           }
-          cuts = getGroupCutoffs( j1, j2 );
-          if (dr.lengthSquare() < cuts.third) {
+          getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq );
+          if (dr.lengthSquare() < rlistsq) {
             neighborList.push_back(make_pair(j1, j2));
           }
         }    
@@ -1583,7 +1622,5 @@ namespace OpenMD {
     saved_CG_positions_.clear();
     for (int i = 0; i < nGroups_; i++)
       saved_CG_positions_.push_back(snap_->cgData.position[i]);
-    
-    return neighborList;
   }
 } //end namespace OpenMD