--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/07/15 17:10:50	611
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2004/06/07 14:09:02	1251
@@ -1,4 +1,4 @@
-#include <cmath>
+#include <math.h>
 #include <iostream>
 using namespace std;
 
@@ -10,18 +10,20 @@ using namespace std;
 #include "SRI.hpp"
 #include "Integrator.hpp"
 #include "simError.h"
+#include "MatVec3.h"
 
 #ifdef IS_MPI
 #define __C
 #include "mpiSimulation.hpp"
 #endif // is_mpi
 
+inline double roundMe( double x ){
+	  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
+}
 
-#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 );
 }
@@ -36,46 +38,43 @@ double Thermo::getKinetic(){
   double kinetic;
   double amass;
   double aVel[3], aJ[3], I[3][3];
-  int j, kl;
+  int i, j, k, kl;
 
-  DirectionalAtom *dAtom;
-
-  int n_atoms;
   double kinetic_global;
-  Atom** atoms;
-
+  vector<StuntDouble *> integrableObjects = info->integrableObjects;
   
-  n_atoms = entry_plug->n_atoms;
-  atoms = entry_plug->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];
 
-    if( atoms[kl]->isDirectional() ){
-            
-      dAtom = (DirectionalAtom *)atoms[kl];
+  for (kl=0; kl<integrableObjects.size(); kl++) {
+    integrableObjects[kl]->getVel(aVel);
+    amass = integrableObjects[kl]->getMass();
 
-      dAtom->getJ( aJ );
-      dAtom->getI( I );
-      
-      for (j=0; j<3; j++) 
-        kinetic += aJ[j]*aJ[j] / I[j][j];
-      
-    }
+   for(j=0; j<3; j++) 
+      kinetic += amass*aVel[j]*aVel[j];
+
+   if (integrableObjects[kl]->isDirectional()){
+ 
+      integrableObjects[kl]->getJ( aJ );
+      integrableObjects[kl]->getI( I );
+
+      if (integrableObjects[kl]->isLinear()) {
+        i = integrableObjects[kl]->linearAxis();
+        j = (i+1)%3;
+        k = (i+2)%3;
+        kinetic += aJ[j]*aJ[j]/I[j][j] + aJ[k]*aJ[k]/I[k][k];
+      } else {
+        for (j=0; j<3; j++) 
+          kinetic += aJ[j]*aJ[j] / I[j][j];
+      }
+   }
   }
 #ifdef IS_MPI
   MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,
 		MPI_SUM, MPI_COMM_WORLD);
   kinetic = kinetic_global;
 #endif //is_mpi
-
+  
   kinetic = kinetic * 0.5 / e_convert;
 
   return kinetic;
@@ -88,14 +87,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();
   }
 
@@ -107,12 +106,6 @@ double Thermo::getPotential(){
   potential = potential_local; 
 #endif // is_mpi
 
-#ifdef IS_MPI
-  /*
-  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
-  */
-#endif
-
   return potential;
 }
 
@@ -126,47 +119,76 @@ 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;
-  
-  temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb );
+
+  temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb );
   return temperature;
 }
 
-double Thermo::getEnthalpy() {
+double Thermo::getVolume() {
 
-  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
-  double u, p, v;
+  return info->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;
 
-  u = this->getTotalE();
+  this->getPressureTensor(press);
 
+  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
+
+  return pressure;
+}
+
+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);
-  p = (press[0][0] + press[1][1] + press[2][2]) / 3.0;
 
-  v = this->getVolume();
+  pressureX = p_convert * press[0][0];
 
-  return (u + (p*v)/e_convert);
+  return pressureX;
 }
 
-double Thermo::getVolume() {
+double Thermo::getPressureY() {
 
-  return entry_plug->boxVol;
+  // 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::getPressure() {
+double Thermo::getPressureZ() {
 
   // Relies on the calculation of the full molecular pressure tensor
   
   const double p_convert = 1.63882576e8;
   double press[3][3];
-  double pressure;
+  double pressureZ;
 
   this->getPressureTensor(press);
 
-  pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
+  pressureZ = p_convert * press[2][2];
 
-  return pressure;
+  return pressureZ;
 }
 
 
@@ -180,22 +202,23 @@ void Thermo::getPressureTensor(double press[3][3]){
   double molmass, volume;
   double vcom[3];
   double p_local[9], p_global[9];
-  int i, j, k, l, nMols;
-  Molecule* molecules;
+  int i, j, k;
 
-  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);
+  // use velocities of integrableObjects and their masses:  
 
+  for (i=0; i < info->integrableObjects.size(); i++) {
+
+    molmass = info->integrableObjects[i]->getMass();
+    
+    info->integrableObjects[i]->getVel(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]); 
@@ -205,10 +228,11 @@ void Thermo::getPressureTensor(double press[3][3]){
     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
@@ -222,54 +246,72 @@ void Thermo::getPressureTensor(double press[3][3]){
   for(i = 0; i < 3; i++) {
     for (j = 0; j < 3; j++) {
       k = 3*i + j;
-      press[i][j] = (p_global[k] + entry_plug->tau[k]*e_convert) / volume;
+      press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume;
     }
   }
 }
 
 void Thermo::velocitize() {
   
-  double x,y;
   double aVel[3], aJ[3], I[3][3];
-  int i, j, vr, vd; // velocity randomizer loop counters
+  int i, j, l, m, n, 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 n_atoms;
-  Atom** atoms;
-  DirectionalAtom* dAtom;
   double temperature;
-  int n_oriented;
-  int n_constraints;
+  int nobj;
 
-  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;
+  nobj = info->integrableObjects.size();
   
-  kebar = kb * temperature * (double)entry_plug->ndf / 
-    ( 2.0 * (double)entry_plug->ndfRaw );
+  temperature   = info->target_temp;
   
-  for(vr = 0; vr < n_atoms; vr++){
+  kebar = kb * temperature * (double)info->ndfRaw / 
+    ( 2.0 * (double)info->ndf );
+  
+  for(vr = 0; vr < nobj; vr++){
     
     // uses equipartition theory to solve for vbar in angstrom/fs
 
-    av2 = 2.0 * kebar / atoms[vr]->getMass();
+    av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass();
     vbar = sqrt( av2 );
- 
-//     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
-    
+
     // picks random velocities from a gaussian distribution
     // centered on vbar
 
     for (j=0; j<3; j++) 
       aVel[j] = vbar * gaussStream->getGaussian();
     
-    atoms[vr]->setVel( aVel );
+    info->integrableObjects[vr]->setVel( aVel );
+    
+    if(info->integrableObjects[vr]->isDirectional()){
 
+      info->integrableObjects[vr]->getI( I );
+
+      if (info->integrableObjects[vr]->isLinear()) {
+
+        l= info->integrableObjects[vr]->linearAxis();
+        m = (l+1)%3;
+        n = (l+2)%3;
+
+        aJ[l] = 0.0;
+        vbar = sqrt( 2.0 * kebar * I[m][m] );
+        aJ[m] = vbar * gaussStream->getGaussian();
+        vbar = sqrt( 2.0 * kebar * I[n][n] );
+        aJ[n] = vbar * gaussStream->getGaussian();
+        
+      } else {
+        for (j = 0 ; j < 3; j++) {
+          vbar = sqrt( 2.0 * kebar * I[j][j] );
+          aJ[j] = vbar * gaussStream->getGaussian();
+        }	
+      } // else isLinear
+
+      info->integrableObjects[vr]->setJ( aJ ); 
+      
+    }//isDirectional 
+
   }
 
   // Get the Center of Mass drift velocity.
@@ -279,36 +321,16 @@ void Thermo::velocitize() {
   //  Corrects for the center of mass drift.
   // sums all the momentum and divides by total mass.
 
-  for(vd = 0; vd < n_atoms; vd++){
+  for(vd = 0; vd < nobj; vd++){
     
-    atoms[vd]->getVel(aVel);
+    info->integrableObjects[vd]->getVel(aVel);
     
     for (j=0; j < 3; j++) 
       aVel[j] -= vdrift[j];
         
-    atoms[vd]->setVel( aVel );
+    info->integrableObjects[vd]->setVel( aVel );
   }
-  if( n_oriented ){
-  
-    for( i=0; i<n_atoms; i++ ){
-      
-      if( atoms[i]->isDirectional() ){
-	
-	dAtom = (DirectionalAtom *)atoms[i];
-        dAtom->getI( I );
-        
-        for (j = 0 ; j < 3; j++) {
 
-          vbar = sqrt( 2.0 * kebar * I[j][j] );
-          aJ[j] = vbar * gaussStream->getGaussian();
-
-        }	
-
-        dAtom->setJ( aJ );
-
-      }
-    }	
-  }
 }
 
 void Thermo::getCOMVel(double vdrift[3]){
@@ -316,24 +338,20 @@ void Thermo::getCOMVel(double vdrift[3]){
   double mtot, mtot_local;
   double aVel[3], amass;
   double vdrift_local[3];
-  int vd, n_atoms, j;
-  Atom** atoms;
+  int vd, j;
+  int nobj;
 
-  // 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.
+  nobj   = info->integrableObjects.size();
 
-  n_atoms = entry_plug->n_atoms;  
-  atoms   = entry_plug->atoms;
-
   mtot_local = 0.0;
   vdrift_local[0] = 0.0;
   vdrift_local[1] = 0.0;
   vdrift_local[2] = 0.0;
   
-  for(vd = 0; vd < n_atoms; vd++){
+  for(vd = 0; vd < nobj; vd++){
     
-    amass = atoms[vd]->getMass();
-    atoms[vd]->getVel( aVel );
+    amass = info->integrableObjects[vd]->getMass();
+    info->integrableObjects[vd]->getVel( aVel );
 
     for(j = 0; j < 3; j++) 
       vdrift_local[j] += aVel[j] * amass;
@@ -357,3 +375,66 @@ 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, j;
+  int nobj;
+
+  mtot_local = 0.0;
+  COM_local[0] = 0.0;
+  COM_local[1] = 0.0;
+  COM_local[2] = 0.0;
+
+  nobj = info->integrableObjects.size();
+  for(i = 0; i < nobj; i++){
+    
+    amass = info->integrableObjects[i]->getMass();
+    info->integrableObjects[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;
+  }
+}
+
+void Thermo::removeCOMdrift() {
+  double vdrift[3], aVel[3];
+  int vd, j, nobj;
+
+  nobj = info->integrableObjects.size();
+
+  // Get the Center of Mass drift velocity.
+
+  getCOMVel(vdrift);
+  
+  //  Corrects for the center of mass drift.
+  // sums all the momentum and divides by total mass.
+
+  for(vd = 0; vd < nobj; vd++){
+    
+    info->integrableObjects[vd]->getVel(aVel);
+    
+    for (j=0; j < 3; j++) 
+      aVel[j] -= vdrift[j];
+        
+    info->integrableObjects[vd]->setVel( aVel );
+  }
+}