--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/04 01:57:11	454
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/09 04:06:43	483
@@ -125,54 +125,87 @@ int Thermo::getNDF(){
   return total;
 }
 
-int Thermo::getNDF(){
-  int ndf_local, ndf;
-  
-  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
-    - entry_plug->n_constraints;
+double Thermo::getTemperature(){
 
-#ifdef IS_MPI
-  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
-#else
-  ndf = ndf_local;
-#endif
-
-  ndf = ndf - 3;
-
-  return ndf;
+  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
+  double temperature;
+  
+  temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb );
+  return temperature;
 }
 
-int Thermo::getNDFraw() {
-  int ndfRaw_local, ndfRaw;
-
-  // Raw degrees of freedom that we have to set
-  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;
+double Thermo::getPressure() {
+  // returns the pressure in units of atm
+  // Relies on the calculation of the full molecular pressure tensor
   
-#ifdef IS_MPI
-  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
-#else
-  ndfRaw = ndfRaw_local;
-#endif
+  const double p_convert = 1.63882576e8;
+  double press[9];
+  double pressure;
 
-  return ndfRaw;
-}
+  this->getPressureTensor(press);
 
+  pressure = p_convert * (press[0] + press[4] + press[8]) / 3.0;
 
-double Thermo::getTemperature(){
-
-  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
-  double temperature;
-  
-  temperature = ( 2.0 * this->getKinetic() ) / ( (double)this->getNDF() * kb );
-  return temperature;
+  return pressure;
 }
 
-double Thermo::getPressure(){
-  // returns pressure in units amu*fs^-2*Ang^-1
+
+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
 
-  return 0.0;
+  const double e_convert = 4.184e-4;
+
+  double molmass, volume;
+  double vcom[3];
+  double p_local[9], p_global[9];
+  double theBox[3];
+  double* tau;
+  int i, 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
+
+  entry_plug->getBox(theBox);
+
+  volume = theBox[0] * theBox[1] * theBox[2];
+
+  for(i=0; i<9; i++) {
+    press[i] = (p_global[i] - tau[i]*e_convert) / volume;
+  }
 }
 
 void Thermo::velocitize() {
@@ -186,8 +219,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;
@@ -201,8 +232,8 @@ void Thermo::velocitize() {
   n_oriented    = entry_plug->n_oriented;
   n_constraints = entry_plug->n_constraints;
   
-  kebar = kb * temperature * (double)this->getNDF() / 
-    ( 2.0 * (double)this->getNDFraw() );
+  kebar = kb * temperature * (double)entry_plug->ndf / 
+    ( 2.0 * (double)entry_plug->ndfRaw );
   
   for(vr = 0; vr < n_atoms; vr++){