--- branches/development/src/parallel/ForceDecomposition.cpp	2011/05/29 01:49:17	1572
+++ trunk/src/parallel/ForceDecomposition.cpp	2012/08/22 02:28:28	1782
@@ -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/ForceDecomposition.hpp"
 #include "math/Vector3.hpp"
@@ -47,13 +48,20 @@ namespace OpenMD {
 using namespace std;
 namespace OpenMD {
 
-  ForceDecomposition::ForceDecomposition(SimInfo* info) : info_(info) {
+  ForceDecomposition::ForceDecomposition(SimInfo* info, InteractionManager* iMan) : info_(info), interactionMan_(iMan), needVelocities_(false) {
+
     sman_ = info_->getSnapshotManager();
     storageLayout_ = sman_->getStorageLayout();
     ff_ = info_->getForceField();
+    userChoseCutoff_ = false;
 
     Globals* simParams_ = info_->getSimParams();    
-  
+    if (simParams_->havePrintHeatFlux()) {
+      if (simParams_->getPrintHeatFlux()) {
+        needVelocities_ = true;
+      }
+    }
+
     if (simParams_->haveSkinThickness()) {
       skinThickness_ = simParams_->getSkinThickness();
     } else {      
@@ -86,9 +94,13 @@ namespace OpenMD {
     cellOffsets_.push_back( Vector3i(1, -1,1) );
   }
 
-  SelfData ForceDecomposition::fillSelfData(int atom1) {
-    SelfData sdat;
-    // Still Missing atype, skippedCharge, potVec pot,
+  void ForceDecomposition::fillSelfData(SelfData &sdat, int atom1) {
+
+    sdat.atype = atypesLocal[atom1];
+
+    sdat.pot = &embeddingPot;
+    sdat.excludedPot = &excludedSelfPot;
+
     if (storageLayout_ & DataStorage::dslElectroFrame) {
       sdat.eFrame = &(snap_->atomData.electroFrame[atom1]);
     }
@@ -109,13 +121,29 @@ namespace OpenMD {
       sdat.dfrhodrho = &(snap_->atomData.functionalDerivative[atom1]);
     }
 
-    return sdat;    
+    if (storageLayout_ & DataStorage::dslSkippedCharge) {
+      sdat.skippedCharge = &(snap_->atomData.skippedCharge[atom1]);
+    }
+
+    if (storageLayout_ & DataStorage::dslParticlePot) {
+      sdat.particlePot = &(snap_->atomData.particlePot[atom1]);
+    }
+    
+    if (storageLayout_ & DataStorage::dslFlucQPosition) {
+      sdat.flucQ = &(snap_->atomData.flucQPos[atom1]);
+    }
+    
+    if (storageLayout_ & DataStorage::dslFlucQForce) {
+      sdat.dVdFQ = &(snap_->atomData.flucQFrc[atom1]);
+    }
   }
 
   bool ForceDecomposition::checkNeighborList() {
 
     int nGroups = snap_->cgData.position.size();
-
+    if (needVelocities_) 
+      snap_->cgData.setStorageLayout(DataStorage::dslPosition | DataStorage::dslVelocity);
+    
     // if we have changed the group identities or haven't set up the
     // saved positions we automatically will need a neighbor list update:
 
@@ -137,6 +165,20 @@ namespace OpenMD {
     // a conservative test of list skin crossings
     dispmax = 2.0 * sqrt (3.0 * dispmax * dispmax);
 
-    return (dispmax > skinThickness_);    
+
+    if (dispmax > skinThickness_) 
+      return (dispmax > skinThickness_);   
+
+    return false;
   }
+
+  void ForceDecomposition::addToHeatFlux(Vector3d hf) {
+    Vector3d chf = snap_->getConductiveHeatFlux();
+    chf += hf;
+    snap_->setConductiveHeatFlux(chf);
+  }
+  void ForceDecomposition::setHeatFlux(Vector3d hf) {
+    snap_->setConductiveHeatFlux(hf);
+  }
+
 }