--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/03 20:21:54	447
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/10/03 22:11:53	799
@@ -16,12 +16,9 @@ using namespace std;
 #include "mpiSimulation.hpp"
 #endif // is_mpi
 
-
-#define BASE_SEED 123456789
-
-Thermo::Thermo( SimInfo* the_entry_plug ) { 
-  entry_plug = the_entry_plug;
-  int baseSeed = BASE_SEED;
+Thermo::Thermo( SimInfo* the_info ) { 
+  info = the_info;
+  int baseSeed = the_info->getSeed();
   
   gaussStream = new gaussianSPRNG( baseSeed );
 }
@@ -33,10 +30,10 @@ double Thermo::getKinetic(){
 double Thermo::getKinetic(){
 
   const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
-  double vx2, vy2, vz2;
-  double kinetic, v_sqr;
-  int kl;
-  double jx2, jy2, jz2; // the square of the angular momentums
+  double kinetic;
+  double amass;
+  double aVel[3], aJ[3], I[3][3];
+  int j, kl;
 
   DirectionalAtom *dAtom;
 
@@ -45,30 +42,29 @@ double Thermo::getKinetic(){
   Atom** atoms;
 
   
-  n_atoms = entry_plug->n_atoms;
-  atoms = entry_plug->atoms;
+  n_atoms = info->n_atoms;
+  atoms = info->atoms;
 
   kinetic = 0.0;
   kinetic_global = 0.0;
   for( kl=0; kl < n_atoms; kl++ ){
+    
+    atoms[kl]->getVel(aVel);
+    amass = atoms[kl]->getMass();
+    
+    for (j=0; j < 3; j++) 
+      kinetic += amass * aVel[j] * aVel[j];
 
-    vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
-    vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy();
-    vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz();
-
-    v_sqr = vx2 + vy2 + vz2;
-    kinetic += atoms[kl]->getMass() * v_sqr;
-
     if( atoms[kl]->isDirectional() ){
             
       dAtom = (DirectionalAtom *)atoms[kl];
+
+      dAtom->getJ( aJ );
+      dAtom->getI( I );
       
-      jx2 = dAtom->getJx() * dAtom->getJx();	
-      jy2 = dAtom->getJy() * dAtom->getJy();
-      jz2 = dAtom->getJz() * dAtom->getJz();
+      for (j=0; j<3; j++) 
+        kinetic += aJ[j]*aJ[j] / I[j][j];
       
-      kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy()) 
-	+ (jz2 / dAtom->getIzz());
     }
   }
 #ifdef IS_MPI
@@ -89,14 +85,14 @@ double Thermo::getPotential(){
   int el, nSRI;
   Molecule* molecules;
 
-  molecules = entry_plug->molecules;
-  nSRI = entry_plug->n_SRI;
+  molecules = info->molecules;
+  nSRI = info->n_SRI;
 
   potential_local = 0.0;
   potential = 0.0;
-  potential_local += entry_plug->lrPot;
+  potential_local += info->lrPot;
 
-  for( el=0; el<entry_plug->n_mol; el++ ){    
+  for( el=0; el<info->n_mol; el++ ){    
     potential_local += molecules[el].getPotential();
   }
 
@@ -127,46 +123,146 @@ double Thermo::getTemperature(){
 
 double Thermo::getTemperature(){
 
-  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
+  const double kb = 1.9872156E-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;
+  temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
+  return temperature;
+}
 
-#ifdef IS_MPI
-  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
-#else
-  ndf = ndf_local;
-#endif
+double Thermo::getVolume() {
 
-  ndf = ndf - 3;
+  return info->boxVol;
+}
+
+double Thermo::getPressure() {
+
+  // Relies on the calculation of the full molecular pressure tensor
   
-  temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
-  return temperature;
+  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;
 }
 
-double Thermo::getPressure(){
-  // returns pressure in units amu*fs^-2*Ang^-1
+double Thermo::getPressureX() {
+
+  // Relies on the calculation of the full molecular pressure tensor
+  
+  const double p_convert = 1.63882576e8;
+  double press[3][3];
+  double pressureX;
+
+  this->getPressureTensor(press);
+
+  pressureX = p_convert * press[0][0];
+
+  return pressureX;
+}
+
+double Thermo::getPressureY() {
+
+  // Relies on the calculation of the full molecular pressure tensor
+  
+  const double p_convert = 1.63882576e8;
+  double press[3][3];
+  double pressureY;
+
+  this->getPressureTensor(press);
+
+  pressureY = p_convert * press[1][1];
+
+  return pressureY;
+}
+
+double Thermo::getPressureZ() {
+
+  // Relies on the calculation of the full molecular pressure tensor
+  
+  const double p_convert = 1.63882576e8;
+  double press[3][3];
+  double pressureZ;
+
+  this->getPressureTensor(press);
+
+  pressureZ = p_convert * press[2][2];
+
+  return pressureZ;
+}
+
+
+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, nMols;
+  Molecule* molecules;
+
+  nMols = info->n_mol;
+  molecules = info->molecules;
+  //tau = info->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 = this->getVolume();
+
+  for(i = 0; i < 3; i++) {
+    for (j = 0; j < 3; j++) {
+      k = 3*i + j;
+      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
+
+    }
+  }
 }
 
 void Thermo::velocitize() {
   
-  double x,y;
-  double vx, vy, vz;
-  double jx, jy, jz;
-  int i, vr, vd; // velocity randomizer loop counters
+  double aVel[3], aJ[3], I[3][3];
+  int i, j, vr, vd; // velocity randomizer loop counters
   double vdrift[3];
   double vbar;
   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;
@@ -174,30 +270,15 @@ void Thermo::velocitize() {
   int n_oriented;
   int n_constraints;
 
-  atoms         = entry_plug->atoms;
-  n_atoms       = entry_plug->n_atoms;
-  temperature   = entry_plug->target_temp;
-  n_oriented    = entry_plug->n_oriented;
-  n_constraints = entry_plug->n_constraints;
+  atoms         = info->atoms;
+  n_atoms       = info->n_atoms;
+  temperature   = info->target_temp;
+  n_oriented    = info->n_oriented;
+  n_constraints = info->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)info->ndfRaw / 
+    ( 2.0 * (double)info->ndf );
   
-#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++){
     
     // uses equipartition theory to solve for vbar in angstrom/fs
@@ -210,13 +291,11 @@ void Thermo::velocitize() {
     // picks random velocities from a gaussian distribution
     // centered on vbar
 
-    vx = vbar * gaussStream->getGaussian();
-    vy = vbar * gaussStream->getGaussian();
-    vz = vbar * gaussStream->getGaussian();
+    for (j=0; j<3; j++) 
+      aVel[j] = vbar * gaussStream->getGaussian();
+    
+    atoms[vr]->setVel( aVel );
 
-    atoms[vr]->set_vx( vx ); 
-    atoms[vr]->set_vy( vy );
-    atoms[vr]->set_vz( vz );
   }
 
   // Get the Center of Mass drift velocity.
@@ -228,17 +307,12 @@ void Thermo::velocitize() {
 
   for(vd = 0; vd < n_atoms; vd++){
     
-    vx = atoms[vd]->get_vx();
-    vy = atoms[vd]->get_vy();
-    vz = atoms[vd]->get_vz();
-        
-    vx -= vdrift[0];
-    vy -= vdrift[1];
-    vz -= vdrift[2];
+    atoms[vd]->getVel(aVel);
     
-    atoms[vd]->set_vx(vx);
-    atoms[vd]->set_vy(vy);
-    atoms[vd]->set_vz(vz);
+    for (j=0; j < 3; j++) 
+      aVel[j] -= vdrift[j];
+        
+    atoms[vd]->setVel( aVel );
   }
   if( n_oriented ){
   
@@ -247,19 +321,17 @@ void Thermo::velocitize() {
       if( atoms[i]->isDirectional() ){
 	
 	dAtom = (DirectionalAtom *)atoms[i];
+        dAtom->getI( I );
+        
+        for (j = 0 ; j < 3; j++) {
 
-	vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
-	jx = vbar * gaussStream->getGaussian();
+          vbar = sqrt( 2.0 * kebar * I[j][j] );
+          aJ[j] = vbar * gaussStream->getGaussian();
 
-	vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
-	jy = vbar * gaussStream->getGaussian();
+        }	
 
-	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
-	jz = vbar * gaussStream->getGaussian();
-	
-	dAtom->setJx( jx );
-	dAtom->setJy( jy );
-	dAtom->setJz( jz );
+        dAtom->setJ( aJ );
+
       }
     }	
   }
@@ -268,15 +340,16 @@ void Thermo::getCOMVel(double vdrift[3]){
 void Thermo::getCOMVel(double vdrift[3]){
 
   double mtot, mtot_local;
+  double aVel[3], amass;
   double vdrift_local[3];
-  int vd, n_atoms;
+  int vd, n_atoms, j;
   Atom** atoms;
 
   // 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.
 
-  n_atoms = entry_plug->n_atoms;  
-  atoms   = entry_plug->atoms;
+  n_atoms = info->n_atoms;  
+  atoms   = info->atoms;
 
   mtot_local = 0.0;
   vdrift_local[0] = 0.0;
@@ -285,11 +358,13 @@ void Thermo::getCOMVel(double vdrift[3]){
   
   for(vd = 0; vd < n_atoms; vd++){
     
-    vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
-    vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
-    vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
+    amass = atoms[vd]->getMass();
+    atoms[vd]->getVel( aVel );
+
+    for(j = 0; j < 3; j++) 
+      vdrift_local[j] += aVel[j] * amass;
     
-    mtot_local += atoms[vd]->getMass();
+    mtot_local += amass;
   }
 
 #ifdef IS_MPI
@@ -308,3 +383,47 @@ void Thermo::getCOMVel(double vdrift[3]){
   
 }
 
+void Thermo::getCOM(double COM[3]){
+
+  double mtot, mtot_local;
+  double aPos[3], amass;
+  double COM_local[3];
+  int i, n_atoms, j;
+  Atom** atoms;
+
+  // 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.
+
+  n_atoms = info->n_atoms;  
+  atoms   = info->atoms;
+
+  mtot_local = 0.0;
+  COM_local[0] = 0.0;
+  COM_local[1] = 0.0;
+  COM_local[2] = 0.0;
+  
+  for(i = 0; i < n_atoms; i++){
+    
+    amass = atoms[i]->getMass();
+    atoms[i]->getPos( aPos );
+
+    for(j = 0; j < 3; j++) 
+      COM_local[j] += aPos[j] * amass;
+    
+    mtot_local += amass;
+  }
+
+#ifdef IS_MPI
+  MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+#else
+  mtot = mtot_local;
+  for(i = 0; i < 3; i++) {
+    COM[i] = COM_local[i];
+  }
+#endif
+    
+  for (i = 0; i < 3; i++) {
+    COM[i] = COM[i] / mtot;
+  }
+}