nano_mpi/src/dynamics_nose_hoover.F90

! MTK (Martina Tuckerman and Klein) Nosie' Hoover chains
! See Molecular Physics, 1996, Vol. 87, No. 5, 1117-1157 for details of
! the algorithm.
!
! author: Charles F. Vardeman II
! date:   6-12-02
!
! Interfaces:
! Called from dynamics module:
! public: nose_hoover_dynamics():  controls nose-hoover chains dynamics
! for the length of the simulation time:
!     No Arguments:

! private:
! nose_hoovera: updates nose_hoover chains and then advances particle positions and velocities
!     by velocity verlet.
! nose_hooverb: updates velocities, second half of velocity verlet. and thermostat velocities and
!     positions.
! init_hoover: initalize hoover structures
! apply_NHC_par: performs details of nose-hoover chains.


module dynamics_nose_hoover
  use definitions,            ONLY : DP,NDIM
  use constants
  use simulation
  use parameter
  use force_module,           ONLY : calc_forces
  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
#ifdef MPI
  use mpi_module,             ONLY : mpi_allreduce,mpi_comm_world,mpi_err,&
       mpi_double_precision, mpi_sum
#endif

  implicit none
  PRIVATE


  public :: nose_hoover_dynamics

  type (sim_status) :: nose_status


contains


  subroutine nose_hoover_dynamics()
    use velocity, ONLY : calc_temp


    logical                          :: update_nlist
    logical                          :: do_pot
    logical                          :: not_done
    logical                          :: nmflag

    integer                          :: step, i

    ! initialize status (ke,pe,time, etc.)
    call init_status(nose_status)
    call init_dynamics_loop( nose_status )

    !     Start the main simulation loop.

    do_pot = .true.
    not_done = .true.
    nmflag  =  .false.
    step = 0

    dynamics: do
       if (.not.not_done) exit dynamics

       step = step + 1


       if (checktemp.and.(mod(step,check_temp_steps).eq.0)) then
          call calc_temp(nose_status%temp)
          if (dabs(nose_status%temp-target_temp).gt.therm_variance) then
             call scale_velocities(target_temp)
          end if
       end if

       call nose_hoovera()

       call check(update_nlist)


       if (do_pot) then
          call calc_forces(update_nlist,nmflag,pe = nose_status%pot_e)
       else
          call calc_forces(update_nlist,nmflag)
       endif


       call nose_hooverb()

       call evolve_time(step,nose_status,do_pot,not_done)


#ifdef MPI
       call mpi_barrier(mpi_comm_world,mpi_err)
#endif
    end do dynamics

  end subroutine nose_hoover_dynamics


! nose_hoovera(): advances velocities from T to T + DT/2 and positions from T to T + DT
  subroutine nose_hoovera()
    use velocity, ONLY : remove_cm_momentum,remove_angular_momentum


    integer                 ::  i, dim
    real( kind = DP)        :: ke, v2, instatemp, chi
    real( kind = DP)        :: kebar

    real( kind = DP ), dimension(ndim) :: v_dim_i

    real(kind=DP) :: ke_local

    !     units for time are femtosec  (10^-15 sec)
    !     units for distance are angstroms (10^-10 m)
    !     units for velocity are angstroms femtosec^-1
    !     units for mass are a.m.u. (1.661 * 10^-27 kg)
    !     units for force are kcal mol^-1 angstrom^-1
    !
    !     converter will put the final terms into angstroms.
    !       or angstrom/femtosecond.

    !  converter = 1.0d0/2.390664d3


    ! call nose_hoover chains
    call apply_NHC_par()

    ! advance velocities from t to t + dt/2
    do i = 1, nlocal
       v(1:ndim,i) = v(1:ndim,i) + &
            (dt2*f(1:ndim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
    end do

    ! advance positions from t to t+dt
    do i = 1, nlocal
       q(1:ndim,i) = q(1:ndim,i) + dt*v(1:ndim,i)
    end do


  end subroutine nose_hoovera


! advances velocities another half step from t+dt/2 to t + dt
  subroutine nose_hooverb( )
    use velocity, ONLY : remove_cm_momentum,remove_angular_momentum
    integer i

    !     units for time are femtosec  (10^-15 sec)
    !     units for distance are angstroms (10^-10 m)
    !     units for velocity are angstroms femtosec^-1
    !     units for mass are a.m.u. (1.661 * 10^-27 kg)
    !     units for force are kcal mol^-1 angstrom^-1

    do i = 1, nlocal
       v(1:ndim,i) = v(1:ndim,i) + &
            (dt2*f(1:ndim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
    end do

    ! apply nose-hoover thermostat
    call apply_NHC_par()

    if (sim_type  == "cluster") then
       call remove_cm_momentum()
       call remove_angular_momentum()
    endif
  end subroutine nose_hooverb


! does the Nose-Hoover part of the integration from t = 0 to t=DT/2
      ! 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()
    use velocity, ONLY : calc_temp
    real(kind = DP) :: system_ke
    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

    integer ::  i

    call calc_temp(system_temp,system_ke)
    
    G_NkT = natoms*ndim*K_BOLTZ*target_temp
    ! and aren't the system_temp and system_ke similar to chris's kt and ke?

    ! Q_mass is a function of the omega parameter
    Q_mass = G_NkT / (omega * omega)

    ! advance the zeta term to zeta(t + dt)
    zeta = zeta + dt*((system_temp*2 - G_NkT)/Q_mass)
    zetaScale = zeta * dt

    ! 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
    enddo

 end subroutine apply_NHC_par


!!$    ! other code borrowed from Chris Fennell's code of
!!$    ! G.J. Martyna et al. Mol. Phys., 1996, Vol. 87, No. 5, 1117-1157
!!$
!!$  subroutine apply_NHC_par()
!!$    use velocity, ONLY : calc_temp
!!$    real(kind = DP) :: scale
!!$    real(kind = DP) :: system_ke
!!$    real(kind = DP) :: system_temp
!!$    real(kind = DP) :: G_NkT
!!$    real(kind = DP) :: omega2
!!$    real(kind = DP) :: value1, value2, value3, value4
!!$    real(kind = DP) :: wdti2(n_ysval), wdti4(n_ysval), wdti8(n_ysval)
!!$    real(kind = DP) :: G_val(N_nos), vval(N_nos), xival(N_nos)
!!$    real(kind = DP) :: Q_mass
!!$
!!$    integer ::  i, j
!!$
!!$    scale = 1._DP
!!$
!!$    call calc_temp(system_temp,system_ke)
!!$    
!!$      G_NkT = natoms*ndim*K_BOLTZ*my_temp
!!$     
!!$    ! we need to get the w numbers for these values
!!$    ! but i think that chris has calculated the values 
!!$       ! because they are values of constants
!!$    value1 =  0.67560359598d0 * dt
!!$    value2 =  -0.85120719196d0 * dt
!!$    value3 =  0.207245385897d0 * dt
!!$    value4 =  -0.328981543589d0 * dt
!!$    if (n_ysval.eq.3) then
!!$       wdti2(1) = value1
!!$       wdti2(2) = value2
!!$       wdti2(3) = value1
!!$       wdti4(1) = 0.5_DP * wdti2(1)
!!$       wdti4(2) = 0.5_DP * wdti2(2)
!!$       wdti4(3) = wdti4(1)
!!$       wdti8(1) = 0.5_DP * wdti4(1)
!!$       wdti8(2) = 0.5_DP * wdti4(2)
!!$       wdti8(3) = wdti8(1)
!!$    else
!!$       write(*,*) 'Error: n_syval must be 3'
!!$       stop
!!$    endif
!!$
!!$    omega2 = omega * omega
!!$
!!$    ! set the Q_mass values
!!$    Q_mass(1) = G_NkT/omega2
!!$    do i = 2, N_nos-1
!!$       Q_mass(i) = system_temp/omega2
!!$    enddo
!!$
!!$    qmass(N_nos) = 2 * omega
!!$
!!$    ! zero the intial xivals and vvals
!!$   do i = 1, N_nos
!!$       xival(i) = 0._DP
!!$       vval(i) = 0._DP
!!$    enddo
!!$
!!$    ! Start the Nose_Hoover chain stepping
!!$    ! update the forces
!!$    G_val(1) = (system_ke - G_NkT)/Q_mass(1)
!!$
!!$    ! start the multiple time stepping
!!$    do i = 1, n_cval
!!$       do j = 1, n_ysval
!!$        
!!$          ! update the thermostat velocities
!!$          vval(N_nos) = vval(N_nos) + G_val(N_nos)*wdti4(j)/i
!!$          do mcount = 1, N_nos-1
!!$             aa = exp((-wdti8(j)/i) * vval((N_nos+1)-mcount))
!!$             vval(N_nos-mcount) = vval(N_nos-mcount)*aa*aa &
!!$                  + wdti4(j)*G_val(N_nos)*wdti4(j)/(i*i)
!!$          enddo
!!$        
!!$          ! update particle velocities
!!$          aa = exp((-wdti2(j)/i)*vval(1))
!!$          scale = scale*aa
!!$
!!$          ! update the forces
!!$          G_val(1) = (scale*scale*ke -G_NkT)/Q_mass(1)
!!$
!!$          ! update the thermostat positions
!!$          do mcount = 1, N_nos-1
!!$             xival(mcount) = xival(mcount) + vval(mcount)*(wdti2(j)/i)
!!$          enddo
!!$
!!$          ! update the thermostat velocities
!!$          do mcount = 1, N_nos-1
!!$             aa = exp((-wdti8(j)/i)*vval(mcount+1))
!!$             vval(mcount) = vval(mcount)*aa*aa &
!!$                  + (wdti4(j)/i)*G_val(mcount)*aa
!!$             G_val(mcount+1) = (Q_mass(mcount)*vval(mcount)*vval(mcount) &
!!$                  - system_temp)/Q_mass(mcount+1)
!!$          enddo
!!$
!!$          vval(N_nos) = vval(N_nos) + G_val(N_nos)*(wdti4(j)/i)
!!$       enddo
!!$    enddo
!!$
!!$    ! scale both linear velocities and angular momenta
!!$    do i = 1, nlocal
!!$       v(1:ndim,i) = v(1:ndim,i)*scale
!!$       !something here for what we use for angular momenta
!!$    enddo
!!$
!!$  end subroutine apply_NHC_par

        
end module dynamics_nose_hoover
Revision:	8
Committed:	Wed Jun 19 18:48:10 2002 UTC (22 years ago) by pconfort
File size:	9537 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	pconfort	8	! MTK (Martina Tuckerman and Klein) Nosie' Hoover chains
2			! See Molecular Physics, 1996, Vol. 87, No. 5, 1117-1157 for details of
3			! the algorithm.
4			!
5			! author: Charles F. Vardeman II
6			! date: 6-12-02
7			!
8			! Interfaces:
9			! Called from dynamics module:
10			! public: nose_hoover_dynamics(): controls nose-hoover chains dynamics
11			! for the length of the simulation time:
12			! No Arguments:
13
14			! private:
15			! nose_hoovera: updates nose_hoover chains and then advances particle positions and velocities
16			! by velocity verlet.
17			! nose_hooverb: updates velocities, second half of velocity verlet. and thermostat velocities and
18			! positions.
19			! init_hoover: initalize hoover structures
20			! apply_NHC_par: performs details of nose-hoover chains.
21
22
23	chuckv	4	module dynamics_nose_hoover
24			use definitions, ONLY : DP,NDIM
25			use constants
26			use simulation
27			use parameter
28			use force_module, ONLY : calc_forces
29			use force_utilities, ONLY : check
30			use velocity, ONLY : scale_velocities
31			use dynamics_utilities, ONLY : evolve_time, sim_status, init_dynamics_loop, &
32			DT2,DTSQRT,init_status
33			#ifdef MPI
34			use mpi_module, ONLY : mpi_allreduce,mpi_comm_world,mpi_err,&
35			mpi_double_precision, mpi_sum
36			#endif
37
38			implicit none
39			PRIVATE
40
41
42
43			public :: nose_hoover_dynamics
44	pconfort	8
45	chuckv	4	type (sim_status) :: nose_status
46
47
48			contains
49
50
51	pconfort	8	subroutine nose_hoover_dynamics()
52			use velocity, ONLY : calc_temp
53	chuckv	4
54
55
56	pconfort	8	logical :: update_nlist
57			logical :: do_pot
58			logical :: not_done
59			logical :: nmflag
60	chuckv	4
61	pconfort	8	integer :: step, i
62	chuckv	4
63	pconfort	8	! initialize status (ke,pe,time, etc.)
64			call init_status(nose_status)
65			call init_dynamics_loop( nose_status )
66	chuckv	4
67	pconfort	8	! Start the main simulation loop.
68	chuckv	4
69	pconfort	8	do_pot = .true.
70			not_done = .true.
71			nmflag = .false.
72			step = 0
73	chuckv	4
74	pconfort	8	dynamics: do
75			if (.not.not_done) exit dynamics
76	chuckv	4
77	pconfort	8	step = step + 1
78	chuckv	4
79	pconfort	8
80			if (checktemp.and.(mod(step,check_temp_steps).eq.0)) then
81	chuckv	4	call calc_temp(nose_status%temp)
82	pconfort	8	if (dabs(nose_status%temp-target_temp).gt.therm_variance) then
83			call scale_velocities(target_temp)
84			end if
85			end if
86	chuckv	4
87	pconfort	8	call nose_hoovera()
88	chuckv	4
89	pconfort	8	call check(update_nlist)
90	chuckv	4
91
92	pconfort	8	if (do_pot) then
93			call calc_forces(update_nlist,nmflag,pe = nose_status%pot_e)
94			else
95			call calc_forces(update_nlist,nmflag)
96			endif
97	chuckv	4
98	pconfort	8
99			call nose_hooverb()
100
101			call evolve_time(step,nose_status,do_pot,not_done)
102
103
104	chuckv	4	#ifdef MPI
105	pconfort	8	call mpi_barrier(mpi_comm_world,mpi_err)
106	chuckv	4	#endif
107	pconfort	8	end do dynamics
108	chuckv	4
109	pconfort	8	end subroutine nose_hoover_dynamics
110	chuckv	4
111
112	pconfort	8	! nose_hoovera(): advances velocities from T to T + DT/2 and positions from T to T + DT
113			subroutine nose_hoovera()
114			use velocity, ONLY : remove_cm_momentum,remove_angular_momentum
115	chuckv	4
116
117
118
119	pconfort	8	integer :: i, dim
120			real( kind = DP) :: ke, v2, instatemp, chi
121			real( kind = DP) :: kebar
122	chuckv	4
123	pconfort	8	real( kind = DP ), dimension(ndim) :: v_dim_i
124	chuckv	4
125	pconfort	8	real(kind=DP) :: ke_local
126	chuckv	4
127	pconfort	8	! units for time are femtosec (10^-15 sec)
128			! units for distance are angstroms (10^-10 m)
129			! units for velocity are angstroms femtosec^-1
130			! units for mass are a.m.u. (1.661 * 10^-27 kg)
131			! units for force are kcal mol^-1 angstrom^-1
132			!
133			! converter will put the final terms into angstroms.
134			! or angstrom/femtosecond.
135	chuckv	4
136	pconfort	8	! converter = 1.0d0/2.390664d3
137	chuckv	4
138
139	pconfort	8	! call nose_hoover chains
140			call apply_NHC_par()
141	chuckv	4
142	pconfort	8	! advance velocities from t to t + dt/2
143			do i = 1, nlocal
144			v(1:ndim,i) = v(1:ndim,i) + &
145			(dt2f(1:ndim,i)/mass(i))KCALMOL_TO_AMUA2FS2
146			end do
147	chuckv	4
148	pconfort	8	! advance positions from t to t+dt
149			do i = 1, nlocal
150			q(1:ndim,i) = q(1:ndim,i) + dt*v(1:ndim,i)
151			end do
152	chuckv	4
153
154	pconfort	8	end subroutine nose_hoovera
155	chuckv	4
156
157	pconfort	8	! advances velocities another half step from t+dt/2 to t + dt
158			subroutine nose_hooverb( )
159			use velocity, ONLY : remove_cm_momentum,remove_angular_momentum
160			integer i
161	chuckv	4
162	pconfort	8	! units for time are femtosec (10^-15 sec)
163			! units for distance are angstroms (10^-10 m)
164			! units for velocity are angstroms femtosec^-1
165			! units for mass are a.m.u. (1.661 * 10^-27 kg)
166			! units for force are kcal mol^-1 angstrom^-1
167	chuckv	4
168	pconfort	8	do i = 1, nlocal
169			v(1:ndim,i) = v(1:ndim,i) + &
170			(dt2f(1:ndim,i)/mass(i))KCALMOL_TO_AMUA2FS2
171			end do
172	chuckv	4
173	pconfort	8	! apply nose-hoover thermostat
174			call apply_NHC_par()
175	chuckv	4
176	pconfort	8	if (sim_type == "cluster") then
177			call remove_cm_momentum()
178			call remove_angular_momentum()
179			endif
180			end subroutine nose_hooverb
181	chuckv	4
182
183	pconfort	8	! does the Nose-Hoover part of the integration from t = 0 to t=DT/2
184			! borrowed from Chris Fennel's code of thermostating
185			! by way of Hoover, Phys.Rev.A 1985, Vol.31 (3) 1695-1697
186	chuckv	4
187	pconfort	8	subroutine apply_NHC_par()
188			use velocity, ONLY : calc_temp
189			real(kind = DP) :: system_ke
190			real(kind = DP) :: system_temp
191			real(kind = DP) :: G_NkT
192			real(kind = DP) :: Q_mass
193			real(kind = DP) :: zeta
194			real(kind = DP) :: zetaScale
195			real(kind = DP),dimension(ndim) :: vi
196	chuckv	4
197	pconfort	8	integer :: i
198	chuckv	4
199	pconfort	8	call calc_temp(system_temp,system_ke)
200
201			G_NkT = natomsndimK_BOLTZ*target_temp
202			! and aren't the system_temp and system_ke similar to chris's kt and ke?
203
204			! Q_mass is a function of the omega parameter
205			Q_mass = G_NkT / (omega * omega)
206
207			! advance the zeta term to zeta(t + dt)
208			zeta = zeta + dt((system_temp2 - G_NkT)/Q_mass)
209			zetaScale = zeta * dt
210
211			! perform thermostating on velocities and angular momenta
212			do i = 1, nlocal
213			vi(1:ndim) = v(1:ndim,i)*zetaScale
214			! similar thing here for our angular momenta
215
216			v(1:ndim,i) = v(1:ndim,i) - vi(1:ndim)
217			! similar thing here for our angular momenta
218			enddo
219
220			end subroutine apply_NHC_par
221
222
223
224			!!$ ! other code borrowed from Chris Fennell's code of
225			!!$ ! G.J. Martyna et al. Mol. Phys., 1996, Vol. 87, No. 5, 1117-1157
226			!!$
227			!!$ subroutine apply_NHC_par()
228			!!$ use velocity, ONLY : calc_temp
229			!!$ real(kind = DP) :: scale
230			!!$ real(kind = DP) :: system_ke
231			!!$ real(kind = DP) :: system_temp
232			!!$ real(kind = DP) :: G_NkT
233			!!$ real(kind = DP) :: omega2
234			!!$ real(kind = DP) :: value1, value2, value3, value4
235			!!$ real(kind = DP) :: wdti2(n_ysval), wdti4(n_ysval), wdti8(n_ysval)
236			!!$ real(kind = DP) :: G_val(N_nos), vval(N_nos), xival(N_nos)
237			!!$ real(kind = DP) :: Q_mass
238			!!$
239			!!$ integer :: i, j
240			!!$
241			!!$ scale = 1._DP
242			!!$
243			!!$ call calc_temp(system_temp,system_ke)
244			!!$
245			!!$ G_NkT = natomsndimK_BOLTZ*my_temp
246			!!$
247			!!$ ! we need to get the w numbers for these values
248			!!$ ! but i think that chris has calculated the values
249			!!$ ! because they are values of constants
250			!!$ value1 = 0.67560359598d0 * dt
251			!!$ value2 = -0.85120719196d0 * dt
252			!!$ value3 = 0.207245385897d0 * dt
253			!!$ value4 = -0.328981543589d0 * dt
254			!!$ if (n_ysval.eq.3) then
255			!!$ wdti2(1) = value1
256			!!$ wdti2(2) = value2
257			!!$ wdti2(3) = value1
258			!!$ wdti4(1) = 0.5_DP * wdti2(1)
259			!!$ wdti4(2) = 0.5_DP * wdti2(2)
260			!!$ wdti4(3) = wdti4(1)
261			!!$ wdti8(1) = 0.5_DP * wdti4(1)
262			!!$ wdti8(2) = 0.5_DP * wdti4(2)
263			!!$ wdti8(3) = wdti8(1)
264			!!$ else
265			!!$ write(,) 'Error: n_syval must be 3'
266			!!$ stop
267			!!$ endif
268			!!$
269			!!$ omega2 = omega * omega
270			!!$
271			!!$ ! set the Q_mass values
272			!!$ Q_mass(1) = G_NkT/omega2
273			!!$ do i = 2, N_nos-1
274			!!$ Q_mass(i) = system_temp/omega2
275			!!$ enddo
276			!!$
277			!!$ qmass(N_nos) = 2 * omega
278			!!$
279			!!$ ! zero the intial xivals and vvals
280			!!$ do i = 1, N_nos
281			!!$ xival(i) = 0._DP
282			!!$ vval(i) = 0._DP
283			!!$ enddo
284			!!$
285			!!$ ! Start the Nose_Hoover chain stepping
286			!!$ ! update the forces
287			!!$ G_val(1) = (system_ke - G_NkT)/Q_mass(1)
288			!!$
289			!!$ ! start the multiple time stepping
290			!!$ do i = 1, n_cval
291			!!$ do j = 1, n_ysval
292			!!$
293			!!$ ! update the thermostat velocities
294			!!$ vval(N_nos) = vval(N_nos) + G_val(N_nos)*wdti4(j)/i
295			!!$ do mcount = 1, N_nos-1
296			!!$ aa = exp((-wdti8(j)/i) * vval((N_nos+1)-mcount))
297			!!$ vval(N_nos-mcount) = vval(N_nos-mcount)aaaa &
298			!!$ + wdti4(j)G_val(N_nos)wdti4(j)/(i*i)
299			!!$ enddo
300			!!$
301			!!$ ! update particle velocities
302			!!$ aa = exp((-wdti2(j)/i)*vval(1))
303			!!$ scale = scale*aa
304			!!$
305			!!$ ! update the forces
306			!!$ G_val(1) = (scalescaleke -G_NkT)/Q_mass(1)
307			!!$
308			!!$ ! update the thermostat positions
309			!!$ do mcount = 1, N_nos-1
310			!!$ xival(mcount) = xival(mcount) + vval(mcount)*(wdti2(j)/i)
311			!!$ enddo
312			!!$
313			!!$ ! update the thermostat velocities
314			!!$ do mcount = 1, N_nos-1
315			!!$ aa = exp((-wdti8(j)/i)*vval(mcount+1))
316			!!$ vval(mcount) = vval(mcount)aaaa &
317			!!$ + (wdti4(j)/i)G_val(mcount)aa
318			!!$ G_val(mcount+1) = (Q_mass(mcount)vval(mcount)vval(mcount) &
319			!!$ - system_temp)/Q_mass(mcount+1)
320			!!$ enddo
321			!!$
322			!!$ vval(N_nos) = vval(N_nos) + G_val(N_nos)*(wdti4(j)/i)
323			!!$ enddo
324			!!$ enddo
325			!!$
326			!!$ ! scale both linear velocities and angular momenta
327			!!$ do i = 1, nlocal
328			!!$ v(1:ndim,i) = v(1:ndim,i)*scale
329			!!$ !something here for what we use for angular momenta
330			!!$ enddo
331			!!$
332			!!$ end subroutine apply_NHC_par
333
334
335
336	chuckv	4	end module dynamics_nose_hoover