--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/03/25 22:54:16	401
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/03 20:21:54	447
@@ -4,15 +4,19 @@ using namespace std;
 
 #ifdef IS_MPI
 #include <mpi.h>
-#include <mpi++.h>
 #endif //is_mpi
 
 #include "Thermo.hpp"
 #include "SRI.hpp"
 #include "Integrator.hpp"
+#include "simError.h"
+
+#ifdef IS_MPI
 #define __C
-//#include "mpiSimulation.hpp"
+#include "mpiSimulation.hpp"
+#endif // is_mpi
 
+
 #define BASE_SEED 123456789
 
 Thermo::Thermo( SimInfo* the_entry_plug ) { 
@@ -68,7 +72,8 @@ double Thermo::getKinetic(){
     }
   }
 #ifdef IS_MPI
-  MPI::COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
+  MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
+		MPI_SUM, MPI_COMM_WORLD);
   kinetic = kinetic_global;
 #endif //is_mpi
 
@@ -82,25 +87,33 @@ double Thermo::getPotential(){
   double potential_local;
   double potential;
   int el, nSRI;
-  SRI** sris;
+  Molecule* molecules;
 
-  sris = entry_plug->sr_interactions;
+  molecules = entry_plug->molecules;
   nSRI = entry_plug->n_SRI;
 
   potential_local = 0.0;
+  potential = 0.0;
   potential_local += entry_plug->lrPot;
 
-  for( el=0; el<nSRI; el++ ){    
-    potential_local += sris[el]->get_potential();
+  for( el=0; el<entry_plug->n_mol; el++ ){    
+    potential_local += molecules[el].getPotential();
   }
 
   // Get total potential for entire system from MPI.
 #ifdef IS_MPI
-  MPI::COMM_WORLD.Allreduce(&potential_local,&potential,1,MPI_DOUBLE,MPI_SUM);
+  MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE,
+		MPI_SUM, MPI_COMM_WORLD);
 #else
   potential = potential_local; 
 #endif // is_mpi
 
+#ifdef IS_MPI
+  /*
+  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
+  */
+#endif
+
   return potential;
 }
 
@@ -122,7 +135,7 @@ double Thermo::getTemperature(){
     - entry_plug->n_constraints;
 
 #ifdef IS_MPI
-  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
+  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
 #else
   ndf = ndf_local;
 #endif
@@ -134,10 +147,9 @@ double Thermo::getPressure(){
 }
 
 double Thermo::getPressure(){
-
-//  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
-// const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
-//  const double conv_A_m = 1.0E-10; //convert A -> m
+  // returns pressure in units amu*fs^-2*Ang^-1
+  // routine derived via viral theorem description in:
+  // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322
 
   return 0.0;
 }
@@ -148,13 +160,13 @@ void Thermo::velocitize() {
   double vx, vy, vz;
   double jx, jy, jz;
   int i, vr, vd; // velocity randomizer loop counters
-  double *vdrift;
+  double vdrift[3];
   double vbar;
   const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
   double av2;
   double kebar;
-  int ndf; // number of degrees of freedom
-  int ndfRaw; // the raw number of degrees of freedom
+  int ndf, ndf_local; // number of degrees of freedom
+  int ndfRaw, ndfRaw_local; // the raw number of degrees of freedom
   int n_atoms;
   Atom** atoms;
   DirectionalAtom* dAtom;
@@ -168,9 +180,22 @@ void Thermo::velocitize() {
   n_oriented    = entry_plug->n_oriented;
   n_constraints = entry_plug->n_constraints;
   
+  // Raw degrees of freedom that we have to set
+  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;
 
-  ndfRaw = 3 * n_atoms + 3 * n_oriented;
-  ndf = ndfRaw - n_constraints - 3;
+  // Degrees of freedom that can contain kinetic energy
+  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
+    - entry_plug->n_constraints;
+  
+#ifdef IS_MPI
+  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+#else
+  ndfRaw = ndfRaw_local;
+  ndf = ndf_local;
+#endif
+  ndf = ndf - 3;
+
   kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
   
   for(vr = 0; vr < n_atoms; vr++){
@@ -179,7 +204,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
@@ -196,7 +221,7 @@ void Thermo::velocitize() {
 
   // Get the Center of Mass drift velocity.
 
-  vdrift = getCOMVel();
+  getCOMVel(vdrift);
   
   //  Corrects for the center of mass drift.
   // sums all the momentum and divides by total mass.
@@ -240,15 +265,13 @@ double* Thermo::getCOMVel(){
   }
 }
 
-double* Thermo::getCOMVel(){
+void Thermo::getCOMVel(double vdrift[3]){
 
   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.
 
@@ -270,8 +293,8 @@ double* Thermo::getCOMVel(){
   }
 
 #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);
+  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
 #else
   mtot = mtot_local;
   for(vd = 0; vd < 3; vd++) {
@@ -283,6 +306,5 @@ double* Thermo::getCOMVel(){
     vdrift[vd] = vdrift[vd] / mtot;
   }
   
-  return vdrift;
 }