--- branches/development/src/flucq/FluctuatingChargeNVT.cpp	2012/06/05 17:50:53	1735
+++ trunk/src/flucq/FluctuatingChargeNVT.cpp	2012/08/22 02:28:28	1782
@@ -45,84 +45,59 @@
 #include "utils/simError.h"
 #include "utils/PhysicalConstants.hpp"
 
+
 namespace OpenMD {
 
   FluctuatingChargeNVT::FluctuatingChargeNVT(SimInfo* info) : 
-    FluctuatingChargePropagator(info), chiTolerance_ (1e-6), maxIterNum_(4),
-    thermo(info), 
-    currentSnapshot_(info->getSnapshotManager()->getCurrentSnapshot()) {
-
-
-    if (info_->usesFluctuatingCharges()) {
-      if (info_->getNFluctuatingCharges() > 0) {
-
-        hasFlucQ_ = true;
-        Globals* simParams = info_->getSimParams();
-        FluctuatingChargeParameters* fqParams = simParams->getFluctuatingChargeParameters();
-
-        if (simParams->haveDt()) {
-          dt_ = simParams->getDt();
-          dt2_ = dt_ * 0.5;
-        } else {
-          sprintf(painCave.errMsg,
-                  "FluctuatingChargeNVT Error: dt is not set\n");
-          painCave.isFatal = 1;
-          simError();
-        }
-        
-        if (!simParams->getUseIntialExtendedSystemState()) {
-          currentSnapshot_->setChiElectronic(0.0);
-          currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
-        }
-        
-        if (!fqParams->haveTargetTemp()) {
-          sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT "
-                  "propagator without a flucQ.targetTemp!\n");
-          painCave.isFatal = 1;
-          painCave.severity = OPENMD_ERROR;
-          simError();
-        } else {
-          targetTemp_ = fqParams->getTargetTemp();
-        }
-        
-        // We must set tauThermostat.
-        
-        if (!fqParams->haveTauThermostat()) {
-          sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n"
-                  "\tpropagator, you must set flucQ.tauThermostat .\n");
-          
-          painCave.severity = OPENMD_ERROR;
-          painCave.isFatal = 1;
-          simError();
-        } else {
-          tauThermostat_ = fqParams->getTauThermostat();
-        }
-        updateSizes();
-      }
-    }       
+    FluctuatingChargePropagator(info), chiTolerance_ (1e-6), 
+    maxIterNum_(4), thermo(info),
+    snap(info->getSnapshotManager()->getCurrentSnapshot()) {    
   }
 
   void FluctuatingChargeNVT::initialize() {
-
-    if (!hasFlucQ_) return;
-
-    SimInfo::MoleculeIterator i;
-    Molecule::FluctuatingChargeIterator  j;
-    Molecule* mol;
-    Atom* atom;
-    
-    for (mol = info_->beginMolecule(i); mol != NULL; 
-         mol = info_->nextMolecule(i)) {
-      for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
-           atom = mol->nextFluctuatingCharge(j)) {
-        atom->setFlucQPos(0.0);
-        atom->setFlucQVel(0.0);
+    FluctuatingChargePropagator::initialize();  
+    if (hasFlucQ_) {
+      if (info_->getSimParams()->haveDt()) {
+        dt_ = info_->getSimParams()->getDt();
+        dt2_ = dt_ * 0.5;
+      } else {
+        sprintf(painCave.errMsg,
+                "FluctuatingChargeNVT Error: dt is not set\n");
+        painCave.isFatal = 1;
+        simError();
       }
+      
+      if (!info_->getSimParams()->getUseIntialExtendedSystemState()) {
+        snap->setElectronicThermostat(make_pair(0.0, 0.0));
+      }
+      
+      if (!fqParams_->haveTargetTemp()) {
+        sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT "
+                "propagator without a flucQ.targetTemp!\n");
+        painCave.isFatal = 1;
+        painCave.severity = OPENMD_ERROR;
+        simError();
+      } else {
+        targetTemp_ = fqParams_->getTargetTemp();
+      }
+      
+      // We must set tauThermostat.
+      
+      if (!fqParams_->haveTauThermostat()) {
+        sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n"
+                "\tpropagator, you must set flucQ.tauThermostat .\n");
+        
+        painCave.severity = OPENMD_ERROR;
+        painCave.isFatal = 1;
+        simError();
+      } else {
+        tauThermostat_ = fqParams_->getTauThermostat();
+      }
+      updateSizes(); 
     }
-    
-    cerr << "Yeah, you should probably implement this\n";
   }
 
+
   void FluctuatingChargeNVT::moveA() {
 
     if (!hasFlucQ_) return;
@@ -133,11 +108,11 @@ 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; 
          mol = info_->nextMolecule(i)) {
       for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
@@ -146,20 +121,14 @@ namespace OpenMD {
         cvel = atom->getFlucQVel();
         cpos = atom->getFlucQPos();
         cfrc = atom->getFlucQFrc();
-        cmass = atom->getChargeMass();
-        
-	cerr << atom->getType() << "\n";
-	cerr << "Before\n";
-	cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\tcfrc: " << cfrc << "\tcmass: " << cmass << "\n";
-        
+        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
         cpos += dt_ * cvel;
 	
-	cerr << "After\n";
-	cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\n\n";        
-
         atom->setFlucQVel(cvel);
         atom->setFlucQPos(cpos);
       }
@@ -169,14 +138,10 @@ namespace OpenMD {
       (tauThermostat_ * tauThermostat_);
 
     integralOfChidt += chi * dt2_;
-    cerr << "Move A instTemp: " << instTemp << "\n";
-    currentSnapshot_->setChiElectronic(chi);
-    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
-       
+    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
   }
 
   void FluctuatingChargeNVT::updateSizes() {
-    if (!hasFlucQ_) return;
     oldVel_.resize(info_->getNFluctuatingCharges());
   }
 
@@ -187,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;
 
@@ -210,7 +176,6 @@ namespace OpenMD {
     for(int k = 0; k < maxIterNum_; k++) {
       index = 0;
       instTemp = thermo.getElectronicTemperature();
-      cerr << "MoveB instTemp: " << instTemp << "\n";
       // evolve chi another half step using the temperature at t + dt/2
       prevChi = chi;
       chi = oldChi + dt2_ * (instTemp / targetTemp_ - 1.0) / 
@@ -227,7 +192,6 @@ namespace OpenMD {
           
           // velocity half step
           cvel = oldVel_[index] + dt2_ * cfrc / cmass - dt2_*chi*oldVel_[index];
-          
           atom->setFlucQVel(cvel);      
           ++index;          
         }
@@ -236,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_;