--- branches/development/src/parallel/ForceMatrixDecomposition.cpp	2011/08/18 20:18:19	1613
+++ branches/development/src/parallel/ForceMatrixDecomposition.cpp	2012/01/06 19:03:05	1668
@@ -36,7 +36,8 @@
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
  * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
 #include "parallel/ForceMatrixDecomposition.hpp"
 #include "math/SquareMatrix3.hpp"
@@ -247,10 +248,9 @@ namespace OpenMD {
       for (int j = 0; j < nLocal_; j++) {
         int jglob = AtomLocalToGlobal[j];
 
-        if (excludes->hasPair(iglob, jglob))           
+        if (excludes->hasPair(iglob, jglob)) 
           excludesForAtom[i].push_back(j);              
         
-        
         if (oneTwo->hasPair(iglob, jglob)) {
           toposForAtom[i].push_back(j);
           topoDist[i].push_back(1);
@@ -836,22 +836,27 @@ namespace OpenMD {
    */
   bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
     int unique_id_1, unique_id_2;
-    
+        
 #ifdef IS_MPI
     // in MPI, we have to look up the unique IDs for each atom
     unique_id_1 = AtomRowToGlobal[atom1];
     unique_id_2 = AtomColToGlobal[atom2];
+#else
+    unique_id_1 = AtomLocalToGlobal[atom1];
+    unique_id_2 = AtomLocalToGlobal[atom2];
+#endif   
 
-    // this situation should only arise in MPI simulations
     if (unique_id_1 == unique_id_2) return true;
-    
+
+#ifdef IS_MPI
     // this prevents us from doing the pair on multiple processors
     if (unique_id_1 < unique_id_2) {
       if ((unique_id_1 + unique_id_2) % 2 == 0) return true;
     } else {
-      if ((unique_id_1 + unique_id_2) % 2 == 1) return true; 
+      if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
     }
 #endif
+    
     return false;
   }
 
@@ -871,7 +876,7 @@ namespace OpenMD {
     
     for (vector<int>::iterator i = excludesForAtom[atom1].begin();
          i != excludesForAtom[atom1].end(); ++i) {
-      if ( (*i) == atom2 )  return true;
+      if ( (*i) == atom2 ) return true;
     }
 
     return false;
@@ -996,8 +1001,8 @@ namespace OpenMD {
   
   void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) {    
 #ifdef IS_MPI
-    pot_row[atom1] += 0.5 *  *(idat.pot);
-    pot_col[atom2] += 0.5 *  *(idat.pot);
+    pot_row[atom1] += RealType(0.5) *  *(idat.pot);
+    pot_col[atom2] += RealType(0.5) *  *(idat.pot);
 
     atomRowData.force[atom1] += *(idat.f1);
     atomColData.force[atom2] -= *(idat.f1);
@@ -1190,17 +1195,24 @@ namespace OpenMD {
                 }
               }
 #else
-              
               for (vector<int>::iterator j1 = cellList_[m1].begin(); 
                    j1 != cellList_[m1].end(); ++j1) {
                 for (vector<int>::iterator j2 = cellList_[m2].begin(); 
                      j2 != cellList_[m2].end(); ++j2) {
-                  
+     
                   // Always do this if we're in different cells or if
-                  // we're in the same cell and the global index of the
-                  // j2 cutoff group is less than the j1 cutoff group
-                  
-                  if (m2 != m1 || (*j2) < (*j1)) {
+                  // we're in the same cell and the global index of
+                  // the j2 cutoff group is greater than or equal to
+                  // the j1 cutoff group.  Note that Rappaport's code
+                  // has a "less than" conditional here, but that
+                  // deals with atom-by-atom computation.  OpenMD
+                  // allows atoms within a single cutoff group to
+                  // interact with each other.
+
+
+
+                  if (m2 != m1 || (*j2) >= (*j1) ) {
+
                     dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)];
                     snap_->wrapVector(dr);
                     cuts = getGroupCutoffs( (*j1), (*j2) );
@@ -1219,7 +1231,7 @@ namespace OpenMD {
       // branch to do all cutoff group pairs
 #ifdef IS_MPI
       for (int j1 = 0; j1 < nGroupsInRow_; j1++) {
-        for (int j2 = 0; j2 < nGroupsInCol_; j2++) {       
+        for (int j2 = 0; j2 < nGroupsInCol_; j2++) {    
           dr = cgColData.position[j2] - cgRowData.position[j1];
           snap_->wrapVector(dr);
           cuts = getGroupCutoffs( j1, j2 );
@@ -1227,18 +1239,20 @@ namespace OpenMD {
             neighborList.push_back(make_pair(j1, j2));
           }
         }
-      }
+      }      
 #else
-      for (int j1 = 0; j1 < nGroups_ - 1; j1++) {
-        for (int j2 = j1 + 1; j2 < nGroups_; j2++) {
+      // include all groups here.
+      for (int j1 = 0; j1 < nGroups_; j1++) {
+        // include self group interactions j2 == j1
+        for (int j2 = j1; j2 < nGroups_; j2++) {
           dr = snap_->cgData.position[j2] - snap_->cgData.position[j1];
           snap_->wrapVector(dr);
           cuts = getGroupCutoffs( j1, j2 );
           if (dr.lengthSquare() < cuts.third) {
             neighborList.push_back(make_pair(j1, j2));
           }
-        }
-      }        
+        }    
+      }
 #endif
     }