--- trunk/OOPSE/libmdtools/Thermo.cpp	2004/04/19 22:13:01	1125
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2004/04/22 21:33:55	1131
@@ -10,6 +10,7 @@ using namespace std;
 #include "SRI.hpp"
 #include "Integrator.hpp"
 #include "simError.h"
+#include "MatVec3.h"
 
 #ifdef IS_MPI
 #define __C
\ No newline at end of file
@@ -195,7 +196,7 @@ void Thermo::getPressureTensor(double press[3][3]){
   const double e_convert = 4.184e-4;
 
   double molmass, volume;
-  double vcom[3];
+  double vcom[3], pcom[3], fcom[3], scaled[3];
   double p_local[9], p_global[9];
   int i, j, k, nMols;
   Molecule* molecules;
\ No newline at end of file
@@ -210,18 +211,33 @@ void Thermo::getPressureTensor(double press[3][3]){
     p_global[i] = 0.0;
   }
 
-  for (i=0; i < nMols; i++) {
-    molmass = molecules[i].getCOMvel(vcom);
+  for (i=0; i < info->integrableObjects.size(); i++) {
 
-    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]); 
+    molmass = info->integrableObjects[i]->getMass();
+    
+    info->integrableObjects[i]->getVel(vcom);
+    info->integrableObjects[i]->getPos(pcom);
+    info->integrableObjects[i]->getFrc(fcom);
+
+    matVecMul3(info->HmatInv, pcom, scaled);
+  
+    for(j=0; j<3; j++)
+      scaled[j] -= roundMe(scaled[j]);
+
+    // calc the wrapped real coordinates from the wrapped scaled coordinates
+  
+    matVecMul3(info->Hmat, scaled, pcom);
+    
+    p_local[0] += molmass * (vcom[0] * vcom[0]) + fcom[0]*pcom[0]*eConvert; 
+    p_local[1] += molmass * (vcom[0] * vcom[1]) + fcom[0]*pcom[1]*eConvert; 
+    p_local[2] += molmass * (vcom[0] * vcom[2]) + fcom[0]*pcom[2]*eConvert; 
+    p_local[3] += molmass * (vcom[1] * vcom[0]) + fcom[1]*pcom[0]*eConvert; 
+    p_local[4] += molmass * (vcom[1] * vcom[1]) + fcom[1]*pcom[1]*eConvert; 
+    p_local[5] += molmass * (vcom[1] * vcom[2]) + fcom[1]*pcom[2]*eConvert; 
+    p_local[6] += molmass * (vcom[2] * vcom[0]) + fcom[2]*pcom[0]*eConvert; 
+    p_local[7] += molmass * (vcom[2] * vcom[1]) + fcom[2]*pcom[1]*eConvert; 
+    p_local[8] += molmass * (vcom[2] * vcom[2]) + fcom[2]*pcom[2]*eConvert; 
+    
   }
 
   // Get total for entire system from MPI.
\ No newline at end of file
@@ -239,7 +255,7 @@ 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] + info->tau[k]*e_convert) / volume;
+      press[i][j] = p_global[k] /  volume;
 
     }
   }
\ No newline at end of file
@@ -248,33 +264,27 @@ void Thermo::velocitize() {
 void Thermo::velocitize() {
   
   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;
+  double temperature;
+  int nobj;
 
-  atoms         = info->atoms;
-  n_atoms       = info->n_atoms;
+  nobj = info->integrableObjects.size();
+  
   temperature   = info->target_temp;
-  n_oriented    = info->n_oriented;
-  n_constraints = info->n_constraints;
   
   kebar = kb * temperature * (double)info->ndfRaw / 
     ( 2.0 * (double)info->ndf );
   
-  for(vr = 0; vr < n_atoms; vr++){
+  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 );
 
     // picks random velocities from a gaussian distribution
\ No newline at end of file
@@ -283,8 +293,35 @@ void Thermo::velocitize() {
     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.
\ No newline at end of file
@@ -294,36 +331,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]){
\ No newline at end of file
@@ -331,24 +348,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 = info->n_atoms;  
-  atoms   = info->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;
\ No newline at end of file
@@ -377,24 +390,19 @@ 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;
+  int i, 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.
-
-  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++){
+
+  nobj = info->integrableObjects.size();
+  for(i = 0; i < nobj; i++){
     
-    amass = atoms[i]->getMass();
-    atoms[i]->getPos( aPos );
+    amass = info->integrableObjects[i]->getMass();
+    info->integrableObjects[i]->getPos( aPos );
 
     for(j = 0; j < 3; j++) 
       COM_local[j] += aPos[j] * amass;
\ No newline at end of file
@@ -416,3 +424,27 @@ void Thermo::getCOM(double COM[3]){
     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 );
+  }
+}
\ No newline at end of file