--- trunk/OOPSE/libmdtools/NPTi.cpp	2003/09/17 14:22:15	768
+++ trunk/OOPSE/libmdtools/NPTi.cpp	2003/09/19 14:55:41	770
@@ -215,7 +215,6 @@ template<typename T> void NPTi<T>::moveB( void ){
   
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
-
 
   // Set things up for the iteration:
 
@@ -370,7 +369,7 @@ template<typename T> int NPTi<T>::readyCheck() {
     simError();
   } 
 
-    if (!have_eta_tolerance) {
+  if (!have_eta_tolerance) {
     sprintf( painCave.errMsg,
              "NPTi warning: setting eta tolerance to 1e-6\n");
     etaTolerance = 1e-6;
@@ -378,9 +377,18 @@ template<typename T> int NPTi<T>::readyCheck() {
     painCave.isFatal = 0;
     simError();
   } 
-  // We need NkBT a lot, so just set it here:
-
+  
+  
+  // We need NkBT a lot, so just set it here: This is the RAW number
+  // of particles, so no subtraction or addition of constraints or
+  // orientational degrees of freedom:
+  
   NkBT = (double)Nparticles * kB * targetTemp;
+  
+  // fkBT is used because the thermostat operates on more degrees of freedom
+  // than the barostat (when there are particles with orientational degrees
+  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
+  
   fkBT = (double)info->ndf * kB * targetTemp;
 
   return 1;
@@ -389,40 +397,38 @@ template<typename T> double NPTi<T>::getConservedQuant
 template<typename T> double NPTi<T>::getConservedQuantity(void){
 
   double conservedQuantity;
+  double Three_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;
+  double eta2; 
 
-  U = tStats->getTotalE();
+  Energy = tStats->getTotalE();
 
-  TS = fkBT * 
-    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
+  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / 
+    (2.0 * eConvert);
 
-  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
+  thermostat_potential = fkBT* integralOfChidt / eConvert;
 
-  tb2 = tauBarostat * tauBarostat;
-  eta2 = eta * eta;
 
+  barostat_kinetic = 3.0 * NkBT * tauBarostat * tauBarostat * eta * eta / 
+    (2.0 * eConvert);
+  
+  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / 
+    eConvert;
 
-  extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / 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; 
 }