--- branches/development/src/flucq/FluctuatingChargeNVT.cpp	2012/06/22 20:01:37	1761
+++ branches/development/src/flucq/FluctuatingChargeNVT.cpp	2013/02/20 15:39:39	1850
@@ -35,7 +35,7 @@
  *                                                                      
  * [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).          
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
  * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
  * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
@@ -51,7 +51,7 @@ namespace OpenMD {
   FluctuatingChargeNVT::FluctuatingChargeNVT(SimInfo* info) : 
     FluctuatingChargePropagator(info), chiTolerance_ (1e-6), 
     maxIterNum_(4), thermo(info),
-    currentSnapshot_(info->getSnapshotManager()->getCurrentSnapshot()) {    
+    snap(info->getSnapshotManager()->getCurrentSnapshot()) {    
   }
 
   void FluctuatingChargeNVT::initialize() {
@@ -68,8 +68,7 @@ namespace OpenMD {
       }
       
       if (!info_->getSimParams()->getUseIntialExtendedSystemState()) {
-        currentSnapshot_->setChiElectronic(0.0);
-        currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
+        snap->setElectronicThermostat(make_pair(0.0, 0.0));
       }
       
       if (!fqParams_->haveTargetTemp()) {
@@ -109,8 +108,9 @@ namespace OpenMD {
     Atom* atom;
     RealType cvel, cpos, cfrc, cmass;
 
-    RealType chi = currentSnapshot_->getChiElectronic();
-    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
+    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
+    RealType chi = thermostat.first;
+    RealType integralOfChidt = thermostat.second;
     RealType instTemp = thermo.getElectronicTemperature();
 
     for (mol = info_->beginMolecule(i); mol != NULL; 
@@ -123,6 +123,7 @@ namespace OpenMD {
         cfrc = atom->getFlucQFrc();
         cmass = atom->getChargeMass();       
 
+        cerr << cpos << "\t" << cvel << "\t" << cfrc << "\t" << instTemp << "\t" <<chi << "\t" << integralOfChidt << "\n";
 	// velocity half step
         cvel += dt2_ * cfrc / cmass - dt2_*chi*cvel;                    
         // position whole step
@@ -137,9 +138,7 @@ namespace OpenMD {
       (tauThermostat_ * tauThermostat_);
 
     integralOfChidt += chi * dt2_;
-    currentSnapshot_->setChiElectronic(chi);
-    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
-       
+    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
   }
 
   void FluctuatingChargeNVT::updateSizes() {
@@ -153,10 +152,11 @@ namespace OpenMD {
     Molecule* mol;
     Atom* atom;
     RealType instTemp;
-    RealType chi = currentSnapshot_->getChiElectronic();
+    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
+    RealType chi = thermostat.first;
     RealType oldChi = chi;
     RealType prevChi;
-    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
+    RealType integralOfChidt = thermostat.second;
     int index;
     RealType cfrc, cvel, cmass;
 
@@ -200,20 +200,19 @@ namespace OpenMD {
         break;
     }
     integralOfChidt += dt2_ * chi;
-    currentSnapshot_->setChiElectronic(chi);
-    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
+    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
   }
   
   void FluctuatingChargeNVT::resetPropagator() {
     if (!hasFlucQ_) return;
-    currentSnapshot_->setChiElectronic(0.0);
-    currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
+    snap->setElectronicThermostat(make_pair(0.0, 0.0));
   }
   
   RealType FluctuatingChargeNVT::calcConservedQuantity() {
     if (!hasFlucQ_) return 0.0;
-    RealType chi = currentSnapshot_->getChiElectronic();
-    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
+    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
+    RealType chi = thermostat.first;
+    RealType integralOfChidt = thermostat.second;
     RealType fkBT = info_->getNFluctuatingCharges() * 
       PhysicalConstants::kB *targetTemp_;