--- trunk/OOPSE/libmdtools/ExtendedSystem.cpp	2003/04/03 23:49:20	453
+++ trunk/OOPSE/libmdtools/ExtendedSystem.cpp	2003/04/04 01:57:11	454
@@ -1,4 +1,9 @@
 #include <math.h>
+#include "Atom.hpp"
+#include "Molecule.hpp"
+#include "SimInfo.hpp"
+#include "Thermo.hpp"
+#include "ExtendedSystem.hpp"
 
 ExtendedSystem::ExtendedSystem( SimInfo &info ) {
 
@@ -9,6 +14,7 @@ ExtendedSystem::ExtendedSystem( SimInfo &info ) {
   atoms = info.atoms;
   nMols = info.n_mol;
   molecules = info.molecules;
+  zeta = 0;
 
 }
 
@@ -16,48 +22,60 @@ void ExtendedSystem::nose_hoover_nvt( double ke, doubl
 }
 
 
-void ExtendedSystem::nose_hoover_nvt( double ke, double dt, double temp ){
+void ExtendedSystem::NoseHooverNVT( double dt ){
 
   // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
   
-  int i, j, degrees_freedom;
-  double ke, dt, temp, kB;
-  double keconverter, NkBT, zetaScale, ke_temp;
-  double vxi, vyi, vzi, jxi, jyi, jzi;
+  int i;
+  double kB, keconverter, NkBT, zetaScale, ke_temp;
+  double vx, vy, vz, jx, jy, jz;
   
-  degrees_freedom = 6*nmol; // number of degrees of freedom for the system
   kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K
   keconverter = 4.184e-4; // to convert ke from kcal/mol to amu*Ang^2*fs^-2 / K
   
-  ke_temp = ke * keconverter;
-  NkBT = degrees_freedom*kB*temp;
+  ke_temp = getKinetic() * keconverter;
+  NkBT = (double)getNDF() * kB * targetTemp;
 
   // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin &
   // qmass is set in the parameter file
-  zeta = zeta + dt*((ke_temp*2 - NkBT)/qmass);
+  zeta += dt*((ke_temp*2 - NkBT)/qmass);
   zetaScale = zeta * dt;
 
   // perform thermostat scaling on linear velocities and angular momentum
-  for(i = 0, i < nmol; i++ ) {
-    vxi = vx(i)*zetaScale;
-    vyi = vy(i)*zetaScale;
-    vzi = vz(i)*zetaScale;
-    jxi = jx(i)*zetaScale;
-    jyi = jy(i)*zetaScale;
-    jzi = jz(i)*zetaScale;
-
-    vx(i) = vx(i) - vxi;
-    vy(i) = vy(i) - vyi;
-    vz(i) = vz(i) - vzi;
-    jx(i) = jx(i) - jxi;
-    jy(i) = jy(i) - jyi;
-    jz(i) = jz(i) - jzi;
+  
+  for(i = 0; i < n_atoms; i++){
+    
+    vx = atoms[i]->get_vx();
+    vy = atoms[i]->get_vy();
+    vz = atoms[i]->get_vz();
+    
+    atoms[i]->set_vx(vx - zetaScale * vx);
+    atoms[i]->set_vy(vy - zetaScale * vy);
+    atoms[i]->set_vz(vz - zetaScale * vz);
   }
+  if( n_oriented ){
+    
+    for( i=0; i < n_atoms; i++ ){
+      
+      if( atoms[i]->isDirectional() ){
+	
+	dAtom = (DirectionalAtom *)atoms[i];
+        
+        jx = dAtom->getJx();
+        jy = dAtom->getJy();
+        jz = dAtom->getJz();
+	
+	dAtom->setJx( jx - zetaScale * jx);
+	dAtom->setJy( jy - zetaScale * jy);
+	dAtom->setJz( jz - zetaScale * jz);
+      }
+    }	
+  }
 }
 
 
-void ExtendedSystem::nose_hoover_anderson_npt(double pressure, double ke, double dt, 
-                                   double temp ) {
+void ExtendedSystem::NoseHooverAndersonNPT(double pressure, double ke, 
+                                           double dt, double temp ) {
 
   // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 
   // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500
@@ -82,7 +100,7 @@ void ExtendedSystem::nose_hoover_anderson_npt(double p
   // propogate the strain rate
 
   epsilon_dot +=  dt * ( (p_mol - pressure_ext)*volume 
-                         / (tau_relax*tau_relax * kB * temp) );
+                         / (tau_relax*tau_relax * kB * targetTemp) );
 
   // determine the change in cell volume
   scale = pow( (1.0 + dt * 3.0 * epsilon_dot), (1.0 / 3.0));
@@ -101,7 +119,7 @@ void ExtendedSystem::nose_hoover_anderson_npt(double p
   boxy = boxy_old*scale;
   boxz = boxz_old*scale;
 
-  epsilonScale = epsilon_dot * dt;
+  epsilonScale = epsilonDot * dt;
 
   // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin 
   // qmass is set in the parameter file
@@ -109,32 +127,37 @@ void ExtendedSystem::nose_hoover_anderson_npt(double p
   zetaScale = zeta * dt;
   
   // apply barostating and thermostating to velocities and angular momenta
-  
-  for (i=0; i < nmol; i++) {
-
-    vxi = vx(i)*epsilonScale;
-    vyi = vy(i)*epsilonScale;
-    vzi = vz(i)*epsilonScale;
-    vxi = vxi + vx(i)*zetaScale;
-    vyi = vyi + vy(i)*zetaScale;
-    vzi = vzi + vz(i)*zetaScale;
-    jxi = jx(i)*zetaScale;
-    jyi = jy(i)*zetaScale;
-    jzi = jz(i)*zetaScale;
-
-    vx(i) = vx(i) - vxi;
-    vy(i) = vy(i) - vyi;
-    vz(i) = vz(i) - vzi;
-    jx(i) = jx(i) - jxi;
-    jy(i) = jy(i) - jyi;
-    jz(i) = jz(i) - jzi;
+  for(i = 0; i < n_atoms; i++){
+    
+    vx = atoms[i]->get_vx();
+    vy = atoms[i]->get_vy();
+    vz = atoms[i]->get_vz();
+    
+    atoms[i]->set_vx(vx * (1.0 - zetaScale * epsilonScale));
+    atoms[i]->set_vy(vy * (1.0 - zetaScale * epsilonScale));
+    atoms[i]->set_vz(vz * (1.0 - zetaScale * epsilonScale));
   }
-  
-
-
+  if( n_oriented ){
+    
+    for( i=0; i < n_atoms; i++ ){
+      
+      if( atoms[i]->isDirectional() ){
+	
+	dAtom = (DirectionalAtom *)atoms[i];
+        
+        jx = dAtom->getJx();
+        jy = dAtom->getJy();
+        jz = dAtom->getJz();
+	
+	dAtom->setJx( jx * (1.0 - zetaScale));
+	dAtom->setJy( jy * (1.0 - zetaScale));
+	dAtom->setJz( jz * (1.0 - zetaScale));
+      }
+    }	
+  }
 }
 
-void ExtendedSystem::affine_transform( double scale ){
+void ExtendedSystem::AffineTransform( double scale ){
 
   int i;
   double boxx_old, boxy_old, boxz_old, percentScale;