--- trunk/src/integrators/NgammaT.cpp	2010/05/10 17:28:26	1442
+++ trunk/src/integrators/NgammaT.cpp	2015/03/05 16:30:23	2069
@@ -35,8 +35,9 @@
  *                                                                      
  * [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).                        
+ * [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).
  */
  
 #include "brains/SimInfo.hpp"
@@ -62,7 +63,7 @@ namespace OpenMD {
 
   }
   void NgammaT::evolveEtaA() {
-    Mat3x3d hmat = currentSnapshot_->getHmat();
+    Mat3x3d hmat = snap->getHmat();
     RealType hz = hmat(2, 2);
     RealType Axy = hmat(0,0) * hmat(1, 1);
     RealType sx = -hz * (press(0, 0) - targetPressure/PhysicalConstants::pressureConvert);
@@ -74,7 +75,7 @@ namespace OpenMD {
   }
 
   void NgammaT::evolveEtaB() {
-    Mat3x3d hmat = currentSnapshot_->getHmat();
+    Mat3x3d hmat = snap->getHmat();
     RealType hz = hmat(2, 2);
     RealType Axy = hmat(0,0) * hmat(1, 1);
     prevEta = eta;
@@ -92,7 +93,7 @@ namespace OpenMD {
 	vScale(i, j) = eta(i, j);
 
 	if (i == j) {
-	  vScale(i, j) += chi;
+	  vScale(i, j) += thermostat.first;
 	}
       }
     }
@@ -119,9 +120,9 @@ namespace OpenMD {
     scaleMat(0, 0) = exp(dt*eta(0, 0));
     scaleMat(1, 1) = exp(dt*eta(1, 1));    
     scaleMat(2, 2) = exp(dt*eta(2, 2));
-    Mat3x3d hmat = currentSnapshot_->getHmat();
+    Mat3x3d hmat = snap->getHmat();
     hmat = hmat *scaleMat;
-    currentSnapshot_->setHmat(hmat);
+    snap->setHmat(hmat);
 
   }
 
@@ -141,8 +142,7 @@ namespace OpenMD {
 
   RealType NgammaT::calcConservedQuantity(){
 
-    chi= currentSnapshot_->getChi();
-    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
+    thermostat = snap->getThermostat();
     loadEta();
     
     // We need NkBT a lot, so just set it here: This is the RAW number
@@ -156,11 +156,12 @@ namespace OpenMD {
     fkBT = info_->getNdf()*PhysicalConstants::kB *targetTemp;    
     
 
-    RealType totalEnergy = thermo.getTotalE();
+    RealType totalEnergy = thermo.getTotalEnergy();
 
-    RealType thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * PhysicalConstants::energyConvert);
+    RealType thermostat_kinetic = fkBT * tt2 * thermostat.first *
+      thermostat.first /(2.0 * PhysicalConstants::energyConvert);
 
-    RealType thermostat_potential = fkBT* integralOfChidt / PhysicalConstants::energyConvert;
+    RealType thermostat_potential = fkBT* thermostat.second / PhysicalConstants::energyConvert;
 
     SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
     RealType trEta = tmp.trace();
@@ -169,19 +170,19 @@ namespace OpenMD {
 
     RealType barostat_potential = (targetPressure * thermo.getVolume() / PhysicalConstants::pressureConvert) /PhysicalConstants::energyConvert;
 
-    Mat3x3d hmat = currentSnapshot_->getHmat();
-    RealType hz = hmat(2, 2);
+    Mat3x3d hmat = snap->getHmat();
     RealType area = hmat(0,0) * hmat(1, 1);
 
-    RealType conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
-      barostat_kinetic + barostat_potential - surfaceTension * area/ PhysicalConstants::energyConvert;
+    RealType conservedQuantity = totalEnergy + thermostat_kinetic
+      + thermostat_potential + barostat_kinetic + barostat_potential
+      - surfaceTension * area/ PhysicalConstants::energyConvert;
 
     return conservedQuantity;
 
   }
 
   void NgammaT::loadEta() {
-    eta= currentSnapshot_->getEta();
+    eta= snap->getBarostat();
 
     //if (!eta.isDiagonal()) {
     //    sprintf( painCave.errMsg,
@@ -192,7 +193,7 @@ namespace OpenMD {
   }
 
   void NgammaT::saveEta() {
-    currentSnapshot_->setEta(eta);
+    snap->setBarostat(eta);
   }
 
 }