--- trunk/nano_mpi/src/dynamics_nose_hoover.F90	2002/06/19 18:48:10	8
+++ trunk/nano_mpi/src/dynamics_nose_hoover.F90	2002/07/01 15:27:33	9
@@ -29,7 +29,7 @@ module dynamics_nose_hoover
   use force_utilities,        ONLY : check
   use velocity,               ONLY : scale_velocities
   use dynamics_utilities,     ONLY : evolve_time, sim_status, init_dynamics_loop, &
-       DT2,DTSQRT,init_status
+       DT2,DTSQRT,DTSQ2,init_status
 #ifdef MPI
   use mpi_module,             ONLY : mpi_allreduce,mpi_comm_world,mpi_err,&
        mpi_double_precision, mpi_sum
@@ -37,14 +37,14 @@ module dynamics_nose_hoover
 
   implicit none
   PRIVATE
+  SAVE
 
 
-
   public :: nose_hoover_dynamics
 
   type (sim_status) :: nose_status
 
-
+ 
 contains
 
 
@@ -58,6 +58,8 @@ contains
     logical                          :: not_done
     logical                          :: nmflag
 
+
+
     integer                          :: step, i
 
     ! initialize status (ke,pe,time, etc.)
@@ -137,7 +139,8 @@ contains
 
 
     ! call nose_hoover chains
-    call apply_NHC_par()
+    call apply_NHC_par(dt2)
+    
 
     ! advance velocities from t to t + dt/2
     do i = 1, nlocal
@@ -155,7 +158,7 @@ contains
 
 
 ! advances velocities another half step from t+dt/2 to t + dt
-  subroutine nose_hooverb( )
+  subroutine nose_hooverb()
     use velocity, ONLY : remove_cm_momentum,remove_angular_momentum
     integer i
 
@@ -171,12 +174,12 @@ contains
     end do
 
     ! apply nose-hoover thermostat
-    call apply_NHC_par()
+    call apply_NHC_par(dt2)
 
-    if (sim_type  == "cluster") then
-       call remove_cm_momentum()
-       call remove_angular_momentum()
-    endif
+!    if (sim_type  == "cluster") then
+!       call remove_cm_momentum()
+!       call remove_angular_momentum()
+!    endif
   end subroutine nose_hooverb
 
 
@@ -184,37 +187,50 @@ contains
       ! borrowed from Chris Fennel's code of thermostating
       ! by way of Hoover, Phys.Rev.A 1985, Vol.31 (3) 1695-1697
 
-  subroutine apply_NHC_par()
+  subroutine apply_NHC_par(this_time)
     use velocity, ONLY : calc_temp
-    real(kind = DP) :: system_ke
+    real(kind = DP) :: system_ke, ke_temp
     real(kind = DP) :: system_temp
     real(kind = DP) :: G_NkT
     real(kind = DP) :: Q_mass
-    real(kind = DP) :: zeta
     real(kind = DP) :: zetaScale
     real(kind = DP),dimension(ndim) :: vi
+    real(kind = DP),intent(out) :: this_time
 
     integer ::  i
-
-    call calc_temp(system_temp,system_ke)
+    integer ::  j
     
-    G_NkT = natoms*ndim*K_BOLTZ*target_temp
-    ! and aren't the system_temp and system_ke similar to chris's kt and ke?
+    call calc_temp(system_temp, system_ke)
+
+    G_NkT = natoms*ndim * 8.31451d-7 * target_temp
+    ke_temp = system_ke*KCALMOL_TO_AMUA2FS2
 
+!    write(*,*) 'System_temp: ', system_temp
+!    write(*,*) 'ke_temp: ',ke_temp
+
     ! Q_mass is a function of the omega parameter
-    Q_mass = G_NkT / (omega * omega)
+    ! Q_mass = G_NkT / (omeg * omeg)
+     Q_mass = 50000._DP
 
     ! advance the zeta term to zeta(t + dt)
-    zeta = zeta + dt*((system_temp*2 - G_NkT)/Q_mass)
-    zetaScale = zeta * dt
+    nzeta = nzeta + this_time*((ke_temp*2 - G_NkT)/Q_mass)
 
+!    write(*,*) 'zeta: ',zeta
+!    write(*,*) 'dt: ',dt
+!    write(*,*) 'ke_temp: ',ke_temp
+!    write(*,*) 'G_NkT: ',G_NkT
+!    write(*,*) 'Q_mass: ',Q_mass
+
+    zetaScale = nzeta * this_time
+
+!    write(*,*) 'zetaScale: ',zetaScale
+
     ! perform thermostating on velocities and angular momenta
     do i = 1, nlocal
-       vi(1:ndim) = v(1:ndim,i)*zetaScale
-       ! similar thing here for our angular momenta
-    
-       v(1:ndim,i) = v(1:ndim,i) - vi(1:ndim)
-       ! similar thing here for our angular momenta
+       do j = 1, ndim
+          vi(j) = v(j,i)*zetaScale          
+          v(j,i) = v(j,i) - vi(j)
+       enddo
     enddo
 
  end subroutine apply_NHC_par