--- trunk/OOPSE/libmdtools/NPTi.cpp	2003/09/17 14:22:15	768
+++ trunk/OOPSE/libmdtools/NPTi.cpp	2003/09/19 14:22:29	769
@@ -389,40 +389,48 @@ template<typename T> double NPTi<T>::getConservedQuant
 template<typename T> double NPTi<T>::getConservedQuantity(void){
 
   double conservedQuantity;
+  double LkBT;
+  double fkBT;
+  double f1kBT;
+  double f2kBT;
+  double NkBT;
+  double Energy;
+  double thermostat_kinetic;
+  double thermostat_potential;
+  double barostat_kinetic;
+  double barostat_potential;
   double tb2;
   double eta2;  
-  double E_NPT;
-  double U;
-  double TS;
-  double PV;
-  double extra;
+
+  LkBT = (double)(info->getNDF() + 4) * kB * targetTemp;  // 3N + 1
+  fkBT = (double)(info->getNDF()    ) * kB * targetTemp;  // 3N - 3
+  f1kBT = (double)(info->getNDF()+ 1) * kB * targetTemp;  // 3N - 3 + 1
+  NkBT = (double)(info->getNDF() + 3) * kB * targetTemp;  // 3N
+  f2kBT = (double)(info->getNDF()+ 2) * kB * targetTemp;  // 3N - 3 + 1
 
-  U = tStats->getTotalE();
+  Energy = tStats->getTotalE();
 
-  TS = fkBT * 
-    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
-
-  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
-
-  tb2 = tauBarostat * tauBarostat;
-  eta2 = eta * eta;
+  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / 
+    (2.0 * eConvert);
 
+  thermostat_potential = fkBT* integralOfChidt / eConvert;
 
-  extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert;
 
+  barostat_kinetic = fkBT * tauBarostat * tauBarostat * eta * eta / 
+    (2.0 * eConvert);
+  
+  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / 
+    eConvert;
+
+  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
+    barostat_kinetic + barostat_potential;
+  
   cout.width(8);
   cout.precision(8);
 
-  
-//   cout << info->getTime() << "\t" 
-//        << chi << "\t"
-//        << eta << "\t"
-//        << U << "\t" 
-//        << TS << "\t"
-//        << PV << "\t"
-//        << extra << "\t"
-//        << U+TS+PV+extra << endl;
+  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << 
+      "\t" << thermostat_potential << "\t" << barostat_kinetic << 
+      "\t" << barostat_potential << "\t" << conservedQuantity << endl;
 
-  conservedQuantity = U+TS+PV+extra;
   return conservedQuantity; 
 }