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,DTSQ2,init_status
#ifdef MPI
  use mpi_module,             ONLY : mpi_allreduce,mpi_comm_world,mpi_err,&
       mpi_double_precision, mpi_sum
#endif

  implicit none
  PRIVATE
  SAVE


  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(dt2)
    

    ! 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(dt2)

!    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(this_time)
    use velocity, ONLY : calc_temp
    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) :: zetaScale
    real(kind = DP),dimension(ndim) :: vi
    real(kind = DP),intent(out) :: this_time

    integer ::  i
    integer ::  j
    
    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 / (omeg * omeg)
     Q_mass = 50000._DP

    ! advance the zeta term to zeta(t + 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
       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


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