--- branches/development/src/parallel/ForceMatrixDecomposition.cpp	2011/08/30 15:45:35	1616
+++ branches/development/src/parallel/ForceMatrixDecomposition.cpp	2012/04/27 20:44:16	1706
@@ -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"
@@ -536,19 +537,21 @@ namespace OpenMD {
       fill(snap_->atomData.density.begin(), 
            snap_->atomData.density.end(), 0.0);
     }
+
     if (storageLayout_ & DataStorage::dslFunctional) {
       fill(snap_->atomData.functional.begin(), 
            snap_->atomData.functional.end(), 0.0);
     }
+
     if (storageLayout_ & DataStorage::dslFunctionalDerivative) {      
       fill(snap_->atomData.functionalDerivative.begin(), 
            snap_->atomData.functionalDerivative.end(), 0.0);
     }
+
     if (storageLayout_ & DataStorage::dslSkippedCharge) {      
       fill(snap_->atomData.skippedCharge.begin(), 
            snap_->atomData.skippedCharge.end(), 0.0);
     }
-    
   }
 
 
@@ -950,6 +953,11 @@ 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]), 
@@ -1000,8 +1008,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);