--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/07/15 17:57:04	614
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/10/28 16:03:37	829
@@ -1,4 +1,4 @@
-#include <cmath>
+#include <math.h>
 #include <iostream>
 using namespace std;
 
@@ -16,12 +16,9 @@ using namespace std;
 #include "mpiSimulation.hpp"
 #endif // is_mpi
 
-
-#define BASE_SEED 123456789
-
 Thermo::Thermo( SimInfo* the_info ) { 
   info = the_info;
-  int baseSeed = BASE_SEED;
+  int baseSeed = the_info->getSeed();
   
   gaussStream = new gaussianSPRNG( baseSeed );
 }
@@ -126,50 +123,79 @@ 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() {
+double Thermo::getVolume() {
 
-  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
-  double u, p, v;
-  double press[3][3];
+  return info->boxVol;
+}
 
-  u = this->getTotalE();
+double Thermo::getPressure() {
 
+  // Relies on the calculation of the full molecular pressure tensor
+  
+  const double p_convert = 1.63882576e8;
+  double press[3][3];
+  double pressure;
+
   this->getPressureTensor(press);
-  p = (press[0][0] + press[1][1] + press[2][2]) / 3.0;
 
-  v = this->getVolume();
+  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
 
-  return (u + (p*v)/e_convert);
+  return pressure;
 }
 
-double Thermo::getVolume() {
+double Thermo::getPressureX() {
 
-  return info->boxVol;
+  // Relies on the calculation of the full molecular pressure tensor
+  
+  const double p_convert = 1.63882576e8;
+  double press[3][3];
+  double pressureX;
+
+  this->getPressureTensor(press);
+
+  pressureX = p_convert * press[0][0];
+
+  return pressureX;
 }
 
-double Thermo::getPressure() {
+double Thermo::getPressureY() {
 
   // Relies on the calculation of the full molecular pressure tensor
   
   const double p_convert = 1.63882576e8;
   double press[3][3];
-  double pressure;
+  double pressureY;
 
   this->getPressureTensor(press);
 
-  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
+  pressureY = p_convert * press[1][1];
 
-  return pressure;
+  return pressureY;
 }
 
+double Thermo::getPressureZ() {
 
+  // Relies on the calculation of the full molecular pressure tensor
+  
+  const double p_convert = 1.63882576e8;
+  double press[3][3];
+  double pressureZ;
+
+  this->getPressureTensor(press);
+
+  pressureZ = p_convert * press[2][2];
+
+  return pressureZ;
+}
+
+
 void Thermo::getPressureTensor(double press[3][3]){
   // returns pressure tensor in units amu*fs^-2*Ang^-1
   // routine derived via viral theorem description in:
@@ -230,7 +256,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];
@@ -251,8 +276,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++){
     
@@ -358,3 +383,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;
+  }
+}