--- branches/development/src/parallel/ForceMatrixDecomposition.cpp	2012/05/24 20:59:54	1723
+++ branches/development/src/parallel/ForceMatrixDecomposition.cpp	2012/06/18 18:23:20	1756
@@ -461,7 +461,6 @@ namespace OpenMD {
     }
   }
 
-
   groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
     int i, j;   
 #ifdef IS_MPI
@@ -852,13 +851,6 @@ namespace OpenMD {
       pairwisePot[ii] = ploc2;
     }
 
-    for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
-      RealType ploc1 = embeddingPot[ii];
-      RealType ploc2 = 0.0;
-      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
-      embeddingPot[ii] = ploc2;
-    }
-    
     // Here be dragons.
     MPI::Intracomm col = colComm.getComm();
 
@@ -868,8 +860,29 @@ namespace OpenMD {
 
 
 #endif
+
+  }
+
+  /** 
+   * Collects information obtained during the post-pair (and embedding
+   * functional) loops onto local data structures.
+   */
+  void ForceMatrixDecomposition::collectSelfData() {
+    snap_ = sman_->getCurrentSnapshot();
+    storageLayout_ = sman_->getStorageLayout();
 
+#ifdef IS_MPI
+    for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
+      RealType ploc1 = embeddingPot[ii];
+      RealType ploc2 = 0.0;
+      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
+      embeddingPot[ii] = ploc2;
+    }    
+#endif
+    
   }
+
+
 
   int ForceMatrixDecomposition::getNAtomsInRow() {   
 #ifdef IS_MPI
@@ -993,16 +1006,20 @@ namespace OpenMD {
    * We need to exclude some overcounted interactions that result from
    * the parallel decomposition.
    */
-  bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
-    int unique_id_1, unique_id_2;
+  bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) {
+    int unique_id_1, unique_id_2, group1, group2;
         
 #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];
+    group1 = cgRowToGlobal[cg1];
+    group2 = cgColToGlobal[cg2];
 #else
     unique_id_1 = AtomLocalToGlobal[atom1];
     unique_id_2 = AtomLocalToGlobal[atom2];
+    group1 = cgLocalToGlobal[cg1];
+    group2 = cgLocalToGlobal[cg2];
 #endif   
 
     if (unique_id_1 == unique_id_2) return true;
@@ -1014,6 +1031,12 @@ namespace OpenMD {
     } else {
       if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
     }
+#endif    
+
+#ifndef IS_MPI
+    if (group1 == group2) {
+      if (unique_id_1 < unique_id_2) return true;
+    }
 #endif
     
     return false;
@@ -1116,14 +1139,7 @@ namespace OpenMD {
 
 #else
     
-
-    // cerr << "atoms = " << atom1 << " " << atom2 << "\n";
-    // cerr << "pos1 = " << snap_->atomData.position[atom1] << "\n";
-    // cerr << "pos2 = " << snap_->atomData.position[atom2] << "\n";
-
     idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
-    //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), 
-    //                         ff_->getAtomType(idents[atom2]) );
 
     if (storageLayout_ & DataStorage::dslAmat) {
       idat.A1 = &(snap_->atomData.aMat[atom1]);
@@ -1183,8 +1199,8 @@ namespace OpenMD {
     atomColData.force[atom2] -= *(idat.f1);
 
     if (storageLayout_ & DataStorage::dslFlucQForce) {              
-      atomRowData.flucQFrc[atom1] += *(idat.dVdFQ1);
-      atomColData.flucQFrc[atom2] += *(idat.dVdFQ2);
+      atomRowData.flucQFrc[atom1] -= *(idat.dVdFQ1);
+      atomColData.flucQFrc[atom2] -= *(idat.dVdFQ2);
     }
 
     if (storageLayout_ & DataStorage::dslElectricField) {              
@@ -1208,7 +1224,7 @@ namespace OpenMD {
     }
     
     if (storageLayout_ & DataStorage::dslFlucQForce) {              
-      snap_->atomData.flucQFrc[atom1] += *(idat.dVdFQ1);
+      snap_->atomData.flucQFrc[atom1] -= *(idat.dVdFQ1);
       snap_->atomData.flucQFrc[atom2] -= *(idat.dVdFQ2);
     }