--- trunk/src/integrators/NPTi.cpp	2009/11/25 20:02:06	1390
+++ trunk/src/integrators/NPTi.cpp	2013/06/16 15:15:42	1879
@@ -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 "NPTi.hpp"
@@ -77,7 +78,7 @@ namespace OpenMD {
   }
 
   void NPTi::calcVelScale() {
-    vScale = chi + eta;
+    vScale = thermostat.first + eta;
   }
 
   void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
@@ -111,9 +112,9 @@ namespace OpenMD {
       painCave.isFatal = 1;
       simError();
     } else {
-      Mat3x3d hmat = currentSnapshot_->getHmat();
+      Mat3x3d hmat = snap->getHmat();
       hmat *= scaleFactor;
-      currentSnapshot_->setHmat(hmat);
+      snap->setHmat(hmat);
     }
 
   }
@@ -125,8 +126,7 @@ namespace OpenMD {
 
   RealType NPTi::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
     // of integrableObjects, so no subtraction or addition of constraints or
@@ -145,11 +145,12 @@ namespace OpenMD {
     RealType barostat_kinetic;
     RealType barostat_potential;
 
-    Energy =thermo.getTotalE();
+    Energy =thermo.getTotalEnergy();
 
-    thermostat_kinetic = fkBT* tt2 * chi * chi / (2.0 * PhysicalConstants::energyConvert);
+    thermostat_kinetic = fkBT* tt2 * thermostat.first * 
+      thermostat.first / (2.0 * PhysicalConstants::energyConvert);
 
-    thermostat_potential = fkBT* integralOfChidt / PhysicalConstants::energyConvert;
+    thermostat_potential = fkBT* thermostat.second / PhysicalConstants::energyConvert;
 
 
     barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * PhysicalConstants::energyConvert);
@@ -164,7 +165,7 @@ namespace OpenMD {
   }
 
   void NPTi::loadEta() {
-    Mat3x3d etaMat = currentSnapshot_->getEta();
+    Mat3x3d etaMat = snap->getBarostat();
     eta = etaMat(0,0);
     //if (fabs(etaMat(1,1) - eta) >= OpenMD::epsilon || fabs(etaMat(1,1) - eta) >= OpenMD::epsilon || !etaMat.isDiagonal()) {
     //    sprintf( painCave.errMsg,
@@ -179,6 +180,6 @@ namespace OpenMD {
     etaMat(0, 0) = eta;
     etaMat(1, 1) = eta;
     etaMat(2, 2) = eta;
-    currentSnapshot_->setEta(etaMat);
+    snap->setBarostat(etaMat);
   }
 }