--- branches/new_design/OOPSE-3.0/src/integrators/NPTxyz.cpp	2004/12/07 14:43:58	1866
+++ branches/new_design/OOPSE-3.0/src/integrators/NPTxyz.cpp	2004/12/07 23:08:14	1867
@@ -22,44 +22,54 @@ double NPTxyz::calcConservedQuantity(){
     
 double NPTxyz::calcConservedQuantity(){
 
-  double conservedQuantity;
-  double totalEnergy;
-  double thermostat_kinetic;
-  double thermostat_potential;
-  double barostat_kinetic;
-  double barostat_potential;
-  double trEta;
+    // 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
+    // orientational degrees of freedom:
+    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;
 
-  totalEnergy = thermo.getTotalE();
+    // fkBT is used because the thermostat operates on more degrees of freedom
+    // than the barostat (when there are particles with orientational degrees
+    // of freedom).  
+    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;        
 
-  thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
+    double conservedQuantity;
+    double totalEnergy;
+    double thermostat_kinetic;
+    double thermostat_potential;
+    double barostat_kinetic;
+    double barostat_potential;
+    double trEta;
 
-  thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
+    totalEnergy = thermo.getTotalE();
 
+    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
+
+    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
+
     SquareMatrix<double, 3> tmp = eta.transpose() * eta;
     trEta = tmp.trace();
 
-  barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
+    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
 
-  barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
+    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
 
-  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
-    barostat_kinetic + barostat_potential;
+    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
+        barostat_kinetic + barostat_potential;
 
 
-  return conservedQuantity;
+    return conservedQuantity;
 
 }
 
     
 void NPTxyz::scaleSimBox(){
 
-  int i,j,k;
-  Mat3x3d scaleMat;
-  double eta2ij, scaleFactor;
-  double bigScale, smallScale, offDiagMax;
-  Mat3x3d hm;
-  Mat3x3d hmnew;
+    int i,j,k;
+    Mat3x3d scaleMat;
+    double eta2ij, scaleFactor;
+    double bigScale, smallScale, offDiagMax;
+    Mat3x3d hm;
+    Mat3x3d hmnew;
 
 
 
@@ -68,47 +78,54 @@ void NPTxyz::scaleSimBox(){
   // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
   //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
 
-  bigScale = 1.0;
-  smallScale = 1.0;
-  offDiagMax = 0.0;
+    bigScale = 1.0;
+    smallScale = 1.0;
+    offDiagMax = 0.0;
 
-  for(i=0; i<3; i++){
-    for(j=0; j<3; j++){
-      scaleMat(i, j) = 0.0;
-      if(i==j) scaleMat(i, j) = 1.0;
+    for(i=0; i<3; i++){
+        for(j=0; j<3; j++){
+            scaleMat(i, j) = 0.0;
+            if(i==j) {
+                scaleMat(i, j) = 1.0;
+            }
+        }
     }
-  }
 
-  for(i=0;i<3;i++){
+    for(i=0;i<3;i++){
 
     // calculate the scaleFactors
 
-    scaleFactor = exp(dt*eta(i, i));
+        scaleFactor = exp(dt*eta(i, i));
 
-    scaleMat(i, i) = scaleFactor;
+        scaleMat(i, i) = scaleFactor;
 
-    if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
-    if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
-  }
+        if (scaleMat(i, i) > bigScale) {
+            bigScale = scaleMat(i, i);
+        }
+        
+        if (scaleMat(i, i) < smallScale) {
+            smallScale = scaleMat(i, i);
+        }
+    }
 
-  if ((bigScale > 1.1) || (smallScale < 0.9)) {
-    sprintf( painCave.errMsg,
-             "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
-             " Check your tauBarostat, as it is probably too small!\n\n"
-             " scaleMat = [%lf\t%lf\t%lf]\n"
-             "            [%lf\t%lf\t%lf]\n"
-             "            [%lf\t%lf\t%lf]\n",
-             scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
-             scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
-             scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
-    painCave.isFatal = 1;
-    simError();
-  } else {
+    if ((bigScale > 1.1) || (smallScale < 0.9)) {
+        sprintf( painCave.errMsg,
+            "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
+            " Check your tauBarostat, as it is probably too small!\n\n"
+            " scaleMat = [%lf\t%lf\t%lf]\n"
+            "            [%lf\t%lf\t%lf]\n"
+            "            [%lf\t%lf\t%lf]\n",
+        scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
+        scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
+        scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
+        painCave.isFatal = 1;
+        simError();
+    } else {
 
         Mat3x3d hmat = currentSnapshot_->getHmat();
         hmat = hmat *scaleMat;
         currentSnapshot_->setHmat(hmat);
-  }
+    }
 }
 
 void NPTxyz::loadEta() {