--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/09/09 20:35:25	755
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2004/04/19 22:13:01	1125
@@ -1,4 +1,4 @@
-#include <cmath>
+#include <math.h>
 #include <iostream>
 using namespace std;
 
@@ -33,46 +33,43 @@ double Thermo::getKinetic(){
   double kinetic;
   double amass;
   double aVel[3], aJ[3], I[3][3];
-  int j, kl;
+  int i, j, k, kl;
 
-  DirectionalAtom *dAtom;
-
-  int n_atoms;
   double kinetic_global;
-  Atom** atoms;
-
+  vector<StuntDouble *> integrableObjects = info->integrableObjects;
   
-  n_atoms = info->n_atoms;
-  atoms = info->atoms;
-
   kinetic = 0.0;
   kinetic_global = 0.0;
-  for( kl=0; kl < n_atoms; kl++ ){
-    
-    atoms[kl]->getVel(aVel);
-    amass = atoms[kl]->getMass();
-    
-    for (j=0; j < 3; j++) 
-      kinetic += amass * aVel[j] * aVel[j];
 
-    if( atoms[kl]->isDirectional() ){
-            
-      dAtom = (DirectionalAtom *)atoms[kl];
+  for (kl=0; kl<integrableObjects.size(); kl++) {
+    integrableObjects[kl]->getVel(aVel);
+    amass = integrableObjects[kl]->getMass();
 
-      dAtom->getJ( aJ );
-      dAtom->getI( I );
-      
-      for (j=0; j<3; j++) 
-        kinetic += aJ[j]*aJ[j] / I[j][j];
-      
-    }
+   for(j=0; j<3; j++) 
+      kinetic += amass*aVel[j]*aVel[j];
+
+   if (integrableObjects[kl]->isDirectional()){
+ 
+      integrableObjects[kl]->getJ( aJ );
+      integrableObjects[kl]->getI( I );
+
+      if (integrableObjects[kl]->isLinear()) {
+        i = integrableObjects[kl]->linearAxis();
+        j = (i+1)%3;
+        k = (i+2)%3;
+        kinetic += aJ[j]*aJ[j]/I[j][j] + aJ[k]*aJ[k]/I[k][k];
+      } else {
+        for (j=0; j<3; j++) 
+          kinetic += aJ[j]*aJ[j] / I[j][j];
+      }
+   }
   }
 #ifdef IS_MPI
   MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
 		MPI_SUM, MPI_COMM_WORLD);
   kinetic = kinetic_global;
 #endif //is_mpi
-
+  
   kinetic = kinetic * 0.5 / e_convert;
 
   return kinetic;
@@ -104,12 +101,6 @@ double Thermo::getPotential(){
   potential = potential_local; 
 #endif // is_mpi
 
-#ifdef IS_MPI
-  /*
-  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
-  */
-#endif
-
   return potential;
 }
 
@@ -123,27 +114,11 @@ double Thermo::getTemperature(){
 
 double Thermo::getTemperature(){
 
-  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
+  const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K)
   double temperature;
-  
+
   temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
   return temperature;
-}
-
-double Thermo::getEnthalpy() {
-
-  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
-  double u, p, v;
-  double press[3][3];
-
-  u = this->getTotalE();
-
-  this->getPressureTensor(press);
-  p = (press[0][0] + press[1][1] + press[2][2]) / 3.0;
-
-  v = this->getVolume();
-
-  return (u + (p*v)/e_convert);
 }
 
 double Thermo::getVolume() {
@@ -272,7 +247,6 @@ void Thermo::velocitize() {
 
 void Thermo::velocitize() {
   
-  double x,y;
   double aVel[3], aJ[3], I[3][3];
   int i, j, vr, vd; // velocity randomizer loop counters
   double vdrift[3];
@@ -293,8 +267,8 @@ void Thermo::velocitize() {
   n_oriented    = info->n_oriented;
   n_constraints = info->n_constraints;
   
-  kebar = kb * temperature * (double)info->ndf / 
-    ( 2.0 * (double)info->ndfRaw );
+  kebar = kb * temperature * (double)info->ndfRaw / 
+    ( 2.0 * (double)info->ndf );
   
   for(vr = 0; vr < n_atoms; vr++){
     
@@ -302,9 +276,7 @@ void Thermo::velocitize() {
 
     av2 = 2.0 * kebar / atoms[vr]->getMass();
     vbar = sqrt( av2 );
- 
-//     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
-    
+
     // picks random velocities from a gaussian distribution
     // centered on vbar
 
@@ -400,3 +372,47 @@ void Thermo::getCOMVel(double vdrift[3]){
   
 }
 
+void Thermo::getCOM(double COM[3]){
+
+  double mtot, mtot_local;
+  double aPos[3], amass;
+  double COM_local[3];
+  int i, n_atoms, j;
+  Atom** atoms;
+
+  // We are very careless here with the distinction between n_atoms and n_local
+  // We should really fix this before someone pokes an eye out.
+
+  n_atoms = info->n_atoms;  
+  atoms   = info->atoms;
+
+  mtot_local = 0.0;
+  COM_local[0] = 0.0;
+  COM_local[1] = 0.0;
+  COM_local[2] = 0.0;
+  
+  for(i = 0; i < n_atoms; i++){
+    
+    amass = atoms[i]->getMass();
+    atoms[i]->getPos( aPos );
+
+    for(j = 0; j < 3; j++) 
+      COM_local[j] += aPos[j] * amass;
+    
+    mtot_local += amass;
+  }
+
+#ifdef IS_MPI
+  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+#else
+  mtot = mtot_local;
+  for(i = 0; i < 3; i++) {
+    COM[i] = COM_local[i];
+  }
+#endif
+    
+  for (i = 0; i < 3; i++) {
+    COM[i] = COM[i] / mtot;
+  }
+}