--- trunk/OOPSE/libmdtools/NPTf.cpp	2003/09/04 21:48:35	746
+++ trunk/OOPSE/libmdtools/NPTf.cpp	2003/09/15 16:52:02	763
@@ -25,6 +25,7 @@ template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
 {
   int i, j;
   chi = 0.0;
+  integralOfChidt = 0.0;
 
   for(i = 0; i < 3; i++) 
     for (j = 0; j < 3; j++) 
@@ -72,7 +73,7 @@ template<typename T> void NPTf<T>::moveA() {
           (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
         
         vScale[i][j] = eta[i][j] + chi;
-        
+          
       } else {
         
         eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
@@ -189,7 +190,6 @@ template<typename T> void NPTf<T>::moveA() {
       if (i != j)
         if (fabs(scaleMat[i][j]) > offDiagMax) 
           offDiagMax = fabs(scaleMat[i][j]);
-      
     }
 
     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
@@ -378,3 +378,28 @@ template<typename T> int NPTf<T>::readyCheck() {
 
   return 1;
 }
+
+template<typename T> double NPTf<T>::getConservedQuantity(void){
+
+  double conservedQuantity;
+  double tb2;
+  double eta2[3][3];  
+  double trEta;
+
+  //HNVE
+  conservedQuantity = tStats->getTotalE();
+
+  //HNVT
+  conservedQuantity += (info->getNDF() * kB * targetTemp * 
+    (integralOfChidt + tauThermostat * tauThermostat * chi * chi /2)) / eConvert;
+
+  //HNPT
+  tb2 = tauBarostat *tauBarostat;
+
+  trEta = info->matTrace3(eta);
+  
+  conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
+			3*NkBT/2 * tb2 * trEta * trEta) / eConvert;
+  
+  return conservedQuantity; 
+}