--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/04 19:47:19	458
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/07/10 22:15:53	590
@@ -134,12 +134,99 @@ double Thermo::getPressure(){
   return temperature;
 }
 
-double Thermo::getPressure(){
-  // returns pressure in units amu*fs^-2*Ang^-1
+double Thermo::getEnthalpy() {
+
+  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
+  double u, p, v;
+  double press[3][3];
+
+  u = this->getTotalE();
+
+  this->getPressureTensor(press);
+  p = (press[0][0] + press[1][1] + press[2][2]) / 3.0;
+
+  v = this->getVolume();
+
+  return (u + (p*v)/e_convert);
+}
+
+double Thermo::getVolume() {
+
+  return entry_plug->boxVol;
+}
+
+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);
+
+  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
+
+  return pressure;
+}
+
+
+void Thermo::getPressureTensor(double press[3][3]){
+  // 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
 
-  return 0.0;
+  const double e_convert = 4.184e-4;
+
+  double molmass, volume;
+  double vcom[3];
+  double p_local[9], p_global[9];
+  int i, j, k, l, nMols;
+  Molecule* molecules;
+
+  nMols = entry_plug->n_mol;
+  molecules = entry_plug->molecules;
+  //tau = entry_plug->tau;
+
+  // use velocities of molecular centers of mass and molecular masses:
+  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[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_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+#else
+  for (i=0; i<9; i++) {
+    p_global[i] = p_local[i]; 
+  }
+#endif // is_mpi
+
+  volume = entry_plug->boxVol;
+
+  for(i = 0; i < 3; i++) {
+    for (j = 0; j < 3; j++) {
+      k = 3*i + j;
+      l = 3*j + i;
+      press[i][j] = (p_global[k] - entry_plug->tau[l]*e_convert) / volume;
+    }
+  }
 }
 
 void Thermo::velocitize() {