--- trunk/OOPSE/libmdtools/NPTi.cpp	2003/09/19 14:22:29	769
+++ 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,25 +397,15 @@ 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 Three_NkBT;
   double Energy;
   double thermostat_kinetic;
   double thermostat_potential;
   double barostat_kinetic;
   double barostat_potential;
   double tb2;
-  double eta2;  
+  double eta2; 
 
-  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
-
   Energy = tStats->getTotalE();
 
   thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / 
@@ -416,7 +414,7 @@ template<typename T> double NPTi<T>::getConservedQuant
   thermostat_potential = fkBT* integralOfChidt / eConvert;
 
 
-  barostat_kinetic = fkBT * tauBarostat * tauBarostat * eta * eta / 
+  barostat_kinetic = 3.0 * NkBT * tauBarostat * tauBarostat * eta * eta / 
     (2.0 * eConvert);
   
   barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /