--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/03/21 17:42:12	378
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/03/25 22:54:16	401
@@ -10,6 +10,8 @@ using namespace std;
 #include "Thermo.hpp"
 #include "SRI.hpp"
 #include "Integrator.hpp"
+#define __C
+//#include "mpiSimulation.hpp"
 
 #define BASE_SEED 123456789
 
@@ -77,28 +79,26 @@ double Thermo::getPotential(){
 
 double Thermo::getPotential(){
   
+  double potential_local;
   double potential;
-  double potential_global;
   int el, nSRI;
   SRI** sris;
 
   sris = entry_plug->sr_interactions;
   nSRI = entry_plug->n_SRI;
 
-  potential = 0.0;
-  potential_global = 0.0;
-  potential += entry_plug->lrPot;
+  potential_local = 0.0;
+  potential_local += entry_plug->lrPot;
 
-  for( el=0; el<nSRI; el++ ){
-    
-    potential += sris[el]->get_potential();
+  for( el=0; el<nSRI; el++ ){    
+    potential_local += sris[el]->get_potential();
   }
 
   // Get total potential for entire system from MPI.
 #ifdef IS_MPI
-  MPI::COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
-  potential = potential_global;
- 
+  MPI::COMM_WORLD.Allreduce(&potential_local,&potential,1,MPI_DOUBLE,MPI_SUM);
+#else
+  potential = potential_local; 
 #endif // is_mpi
 
   return potential;
@@ -116,10 +116,19 @@ double Thermo::getTemperature(){
 
   const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
   double temperature;
+  int ndf_local, ndf;
   
-  int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
-    - entry_plug->n_constraints - 3;
+  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
+    - entry_plug->n_constraints;
 
+#ifdef IS_MPI
+  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
+#else
+  ndf = ndf_local;
+#endif
+
+  ndf = ndf - 3;
+  
   temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
   return temperature;
 }
@@ -139,8 +148,7 @@ void Thermo::velocitize() {
   double vx, vy, vz;
   double jx, jy, jz;
   int i, vr, vd; // velocity randomizer loop counters
-  double vdrift[3];
-  double mtot = 0.0;
+  double *vdrift;
   double vbar;
   const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
   double av2;
@@ -185,35 +193,20 @@ void Thermo::velocitize() {
     atoms[vr]->set_vy( vy );
     atoms[vr]->set_vz( vz );
   }
+
+  // Get the Center of Mass drift velocity.
+
+  vdrift = getCOMVel();
   
   //  Corrects for the center of mass drift.
   // sums all the momentum and divides by total mass.
-  
-  mtot = 0.0;
-  vdrift[0] = 0.0;
-  vdrift[1] = 0.0;
-  vdrift[2] = 0.0;
-  for(vd = 0; vd < n_atoms; vd++){
-    
-    vdrift[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
-    vdrift[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
-    vdrift[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
-    
-    mtot += atoms[vd]->getMass();
-  }
-  
-  for (vd = 0; vd < 3; vd++) {
-    vdrift[vd] = vdrift[vd] / mtot;
-  }
-  
 
   for(vd = 0; vd < n_atoms; vd++){
     
     vx = atoms[vd]->get_vx();
     vy = atoms[vd]->get_vy();
     vz = atoms[vd]->get_vz();
-    
-    
+        
     vx -= vdrift[0];
     vy -= vdrift[1];
     vz -= vdrift[2];
@@ -246,3 +239,50 @@ void Thermo::velocitize() {
     }	
   }
 }
+
+double* Thermo::getCOMVel(){
+
+  double mtot, mtot_local;
+  double* vdrift;
+  double vdrift_local[3];
+  int vd, n_atoms;
+  Atom** atoms;
+
+  vdrift = new double[3];
+  // 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 = entry_plug->n_atoms;  
+  atoms   = entry_plug->atoms;
+
+  mtot_local = 0.0;
+  vdrift_local[0] = 0.0;
+  vdrift_local[1] = 0.0;
+  vdrift_local[2] = 0.0;
+  
+  for(vd = 0; vd < n_atoms; vd++){
+    
+    vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
+    vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
+    vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
+    
+    mtot_local += atoms[vd]->getMass();
+  }
+
+#ifdef IS_MPI
+  MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM);
+  MPI::COMM_WORLD.Allreduce(&vdrift_local,&vdrift,3,MPI_DOUBLE,MPI_SUM);
+#else
+  mtot = mtot_local;
+  for(vd = 0; vd < 3; vd++) {
+    vdrift[vd] = vdrift_local[vd];
+  }
+#endif
+    
+  for (vd = 0; vd < 3; vd++) {
+    vdrift[vd] = vdrift[vd] / mtot;
+  }
+  
+  return vdrift;
+}
+