--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/07/10 22:15:53	590
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/07/15 14:45:09	608
@@ -33,10 +33,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;
 
@@ -51,24 +51,23 @@ double Thermo::getKinetic(){
   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
@@ -232,9 +231,8 @@ void Thermo::velocitize() {
 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.
@@ -268,13 +266,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.
@@ -286,17 +282,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 ){
   
@@ -305,19 +296,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 );
+
       }
     }	
   }
@@ -326,8 +315,9 @@ 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
@@ -343,11 +333,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