--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/03 13:43:02	444
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/08 15:20:44	479
@@ -4,7 +4,6 @@ using namespace std;
 
 #ifdef IS_MPI
 #include <mpi.h>
-#include <mpi++.h>
 #endif //is_mpi
 
 #include "Thermo.hpp"
@@ -73,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
 
@@ -100,17 +100,10 @@ double Thermo::getPotential(){
     potential_local += molecules[el].getPotential();
   }
 
-#ifdef IS_MPI
-  /*
-  std::cerr << "node " << worldRank << ": before LONG RANGE pot = " << entry_plug->lrPot 
-	    << "; pot_local = " << potential_local 
-	    << "; pot = " << potential << "\n";
-  */
-#endif
-
   // 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
@@ -136,30 +129,53 @@ double Thermo::getTemperature(){
 
   const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
   double temperature;
-  int ndf_local, ndf;
   
-  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 );
+  temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb );
   return temperature;
 }
 
 double Thermo::getPressure(){
+  // 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
 
-//  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
+  const double e_convert = 4.184e-4;
+  const double p_convert = 1.63882576e8;
+  double molmass;
+  double vcom[3];
+  double p_local, p_sum, p_mol, virial;
+  double theBox[3];
+  double* tau;
+  int i, nMols;
+  Molecule* molecules;
 
-  return 0.0;
+  nMols = entry_plug->n_mol;
+  molecules = entry_plug->molecules;
+  tau = entry_plug->tau;
+
+  // use velocities of molecular centers of mass and molecular masses:
+  p_local = 0.0;
+
+  for (i=0; i < nMols; i++) {
+    molmass = molecules[i].getCOMvel(vcom);
+    p_local += (vcom[0]*vcom[0] + vcom[1]*vcom[1] + vcom[2]*vcom[2]) * molmass;
+  }
+
+  // Get total for entire system from MPI.
+
+#ifdef IS_MPI
+  MPI_Allreduce(&p_local,&p_sum,1,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+#else
+  p_sum = p_local; 
+#endif // is_mpi
+
+  virial = tau[0] + tau[4] + tau[8];
+  entry_plug->getBox(theBox);
+
+  p_mol = p_convert * (p_sum - virial*e_convert) / 
+    (3.0 * theBox[0] * theBox[1]* theBox[2]); 
+
+  return p_mol;
 }
 
 void Thermo::velocitize() {
@@ -173,8 +189,6 @@ void Thermo::velocitize() {
   const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
   double av2;
   double kebar;
-  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;
@@ -188,24 +202,9 @@ 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;
-
-  // Degrees of freedom that can contain kinetic energy
-  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
-    - entry_plug->n_constraints;
+  kebar = kb * temperature * (double)entry_plug->ndf / 
+    ( 2.0 * (double)entry_plug->ndfRaw );
   
-#ifdef IS_MPI
-  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
-  MPI::COMM_WORLD.Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM);
-#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++){
     
     // uses equipartition theory to solve for vbar in angstrom/fs
@@ -261,7 +260,7 @@ void Thermo::velocitize() {
 
 	vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
 	jy = vbar * gaussStream->getGaussian();
-
+        
 	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
 	jz = vbar * gaussStream->getGaussian();
 	
@@ -301,8 +300,8 @@ void Thermo::getCOMVel(double vdrift[3]){
   }
 
 #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++) {