--- branches/development/src/parallel/ForceMatrixDecomposition.cpp	2011/11/22 20:38:56	1665
+++ branches/development/src/parallel/ForceMatrixDecomposition.cpp	2012/05/24 14:17:42	1721
@@ -523,6 +523,27 @@ namespace OpenMD {
            atomRowData.skippedCharge.end(), 0.0);
       fill(atomColData.skippedCharge.begin(), 
            atomColData.skippedCharge.end(), 0.0);
+    }
+
+    if (storageLayout_ & DataStorage::dslFlucQForce) {      
+      fill(atomRowData.flucQFrc.begin(), 
+           atomRowData.flucQFrc.end(), 0.0);
+      fill(atomColData.flucQFrc.begin(), 
+           atomColData.flucQFrc.end(), 0.0);
+    }
+
+    if (storageLayout_ & DataStorage::dslElectricField) {    
+      fill(atomRowData.electricField.begin(), 
+           atomRowData.electricField.end(), V3Zero);
+      fill(atomColData.electricField.begin(), 
+           atomColData.electricField.end(), V3Zero);
+    }
+
+    if (storageLayout_ & DataStorage::dslFlucQForce) {    
+      fill(atomRowData.flucQFrc.begin(), atomRowData.flucQFrc.end(),
+           0.0);
+      fill(atomColData.flucQFrc.begin(), atomColData.flucQFrc.end(),
+           0.0);
     }
 
 #endif
@@ -537,19 +558,26 @@ 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);
     }
-    
+
+    if (storageLayout_ & DataStorage::dslElectricField) {      
+      fill(snap_->atomData.electricField.begin(), 
+           snap_->atomData.electricField.end(), V3Zero);
+    }
   }
 
 
@@ -589,6 +617,14 @@ namespace OpenMD {
                                    atomColData.electroFrame);
     }
 
+    // if needed, gather the atomic fluctuating charge values
+    if (storageLayout_ & DataStorage::dslFlucQPosition) {
+      AtomPlanRealRow->gather(snap_->atomData.flucQPos, 
+                              atomRowData.flucQPos);
+      AtomPlanRealColumn->gather(snap_->atomData.flucQPos, 
+                                 atomColData.flucQPos);
+    }
+
 #endif      
   }
   
@@ -611,6 +647,18 @@ namespace OpenMD {
       for (int i = 0; i < n; i++)
         snap_->atomData.density[i] += rho_tmp[i];
     }
+
+    if (storageLayout_ & DataStorage::dslElectricField) {
+      
+      AtomPlanVectorRow->scatter(atomRowData.electricField, 
+                                 snap_->atomData.electricField);
+      
+      int n = snap_->atomData.electricField.size();
+      vector<Vector3d> field_tmp(n, V3Zero);
+      AtomPlanVectorColumn->scatter(atomColData.electricField, field_tmp);
+      for (int i = 0; i < n; i++)
+        snap_->atomData.electricField[i] += field_tmp[i];
+    }
 #endif
   }
 
@@ -690,6 +738,23 @@ namespace OpenMD {
             
     }
     
+    if (storageLayout_ & DataStorage::dslFlucQForce) {
+
+      int nq = snap_->atomData.flucQFrc.size();
+      vector<RealType> fqfrc_tmp(nq, 0.0);
+
+      AtomPlanRealRow->scatter(atomRowData.flucQFrc, fqfrc_tmp);
+      for (int i = 0; i < nq; i++) {
+        snap_->atomData.flucQFrc[i] += fqfrc_tmp[i];
+        fqfrc_tmp[i] = 0.0;
+      }
+      
+      AtomPlanRealColumn->scatter(atomColData.flucQFrc, fqfrc_tmp);
+      for (int i = 0; i < nq; i++) 
+        snap_->atomData.flucQFrc[i] += fqfrc_tmp[i];
+            
+    }
+
     nLocal_ = snap_->getNumberOfAtoms();
 
     vector<potVec> pot_temp(nLocal_, 
@@ -950,8 +1015,18 @@ namespace OpenMD {
       idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
     }
 
+    if (storageLayout_ & DataStorage::dslFlucQPosition) {              
+      idat.flucQ1 = &(atomRowData.flucQPos[atom1]);
+      idat.flucQ2 = &(atomColData.flucQPos[atom2]);
+    }
+
 #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]) );
@@ -995,22 +1070,56 @@ namespace OpenMD {
       idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
       idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
     }
+
+    if (storageLayout_ & DataStorage::dslFlucQPosition) {              
+      idat.flucQ1 = &(snap_->atomData.flucQPos[atom1]);
+      idat.flucQ2 = &(snap_->atomData.flucQPos[atom2]);
+    }
+
 #endif
   }
 
   
   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);
+
+    if (storageLayout_ & DataStorage::dslFlucQForce) {              
+      atomRowData.flucQFrc[atom1] += *(idat.dVdFQ1);
+      atomColData.flucQFrc[atom2] += *(idat.dVdFQ2);
+    }
+
+    if (storageLayout_ & DataStorage::dslElectricField) {              
+      atomRowData.electricField[atom1] += *(idat.eField1);
+      atomColData.electricField[atom2] += *(idat.eField2);
+    }
+
+    // should particle pot be done here also?
 #else
     pairwisePot += *(idat.pot);
 
     snap_->atomData.force[atom1] += *(idat.f1);
     snap_->atomData.force[atom2] -= *(idat.f1);
+
+    if (idat.doParticlePot) {
+      snap_->atomData.particlePot[atom1] += *(idat.vpair) * *(idat.sw);
+      snap_->atomData.particlePot[atom2] -= *(idat.vpair) * *(idat.sw);
+    }
+    
+    if (storageLayout_ & DataStorage::dslFlucQForce) {              
+      snap_->atomData.flucQFrc[atom1] += *(idat.dVdFQ1);
+      snap_->atomData.flucQFrc[atom2] -= *(idat.dVdFQ2);
+    }
+
+    if (storageLayout_ & DataStorage::dslElectricField) {              
+      snap_->atomData.electricField[atom1] += *(idat.eField1);
+      snap_->atomData.electricField[atom2] += *(idat.eField2);
+    }
+
 #endif
     
   }