--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/08 15:20:44	479
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/09 04:06:43	483
@@ -134,16 +134,32 @@ double Thermo::getPressure(){
   return temperature;
 }
 
-double Thermo::getPressure(){
-  // returns pressure in units amu*fs^-2*Ang^-1
+double Thermo::getPressure() {
+  // returns the pressure in units of atm
+  // Relies on the calculation of the full molecular pressure tensor
+  
+  const double p_convert = 1.63882576e8;
+  double press[9];
+  double pressure;
+
+  this->getPressureTensor(press);
+
+  pressure = p_convert * (press[0] + press[4] + press[8]) / 3.0;
+
+  return pressure;
+}
+
+
+void Thermo::getPressureTensor(double press[9]){
+  // returns pressure tensor 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 e_convert = 4.184e-4;
-  const double p_convert = 1.63882576e8;
-  double molmass;
+
+  double molmass, volume;
   double vcom[3];
-  double p_local, p_sum, p_mol, virial;
+  double p_local[9], p_global[9];
   double theBox[3];
   double* tau;
   int i, nMols;
@@ -154,28 +170,42 @@ double Thermo::getPressure(){
   tau = entry_plug->tau;
 
   // use velocities of molecular centers of mass and molecular masses:
-  p_local = 0.0;
+  for (i=0; i < 9; i++) {    
+    p_local[i] = 0.0;
+    p_global[i] = 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;
+
+    p_local[0] += molmass * (vcom[0] * vcom[0]); 
+    p_local[1] += molmass * (vcom[0] * vcom[1]); 
+    p_local[2] += molmass * (vcom[0] * vcom[2]); 
+    p_local[3] += molmass * (vcom[1] * vcom[0]); 
+    p_local[4] += molmass * (vcom[1] * vcom[1]); 
+    p_local[5] += molmass * (vcom[1] * vcom[2]); 
+    p_local[6] += molmass * (vcom[2] * vcom[0]); 
+    p_local[7] += molmass * (vcom[2] * vcom[1]); 
+    p_local[8] += molmass * (vcom[2] * vcom[2]); 
   }
 
   // Get total for entire system from MPI.
 
 #ifdef IS_MPI
-  MPI_Allreduce(&p_local,&p_sum,1,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 #else
-  p_sum = p_local; 
+  for (i=0; i<9; i++) {
+    p_global[i] = p_local[i]; 
+  }
 #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]); 
+  volume = theBox[0] * theBox[1] * theBox[2];
 
-  return p_mol;
+  for(i=0; i<9; i++) {
+    press[i] = (p_global[i] - tau[i]*e_convert) / volume;
+  }
 }
 
 void Thermo::velocitize() {