ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/nano_mpi/src/dynamics_nose_hoover.F90
(Generate patch)

Comparing trunk/nano_mpi/src/dynamics_nose_hoover.F90 (file contents):
Revision 7 by chuckv, Mon Jun 10 17:18:36 2002 UTC vs.
Revision 8 by pconfort, Wed Jun 19 18:48:10 2002 UTC

# Line 1 | Line 1
1 + ! 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   module dynamics_nose_hoover
24    use definitions,            ONLY : DP,NDIM
25    use constants
# Line 19 | Line 41 | module dynamics_nose_hoover
41  
42  
43    public :: nose_hoover_dynamics
44 <
44 >
45    type (sim_status) :: nose_status
46  
47  
48   contains
49  
50  
51 < subroutine nose_hoover_dynamics()
52 <  use velocity, ONLY : calc_temp
51 >  subroutine nose_hoover_dynamics()
52 >    use velocity, ONLY : calc_temp
53  
54  
55  
56 <  logical                          :: update_nlist
57 <  logical                          :: do_pot
58 <  logical                          :: not_done
59 <  logical                          :: nmflag
56 >    logical                          :: update_nlist
57 >    logical                          :: do_pot
58 >    logical                          :: not_done
59 >    logical                          :: nmflag
60  
61 <  integer                          :: step, i
40 <
41 <  call init_status(nose_status)
42 <  call init_dynamics_loop( nose_status )
61 >    integer                          :: step, i
62  
63 <  !     Start the main simulation loop.
63 >    ! initialize status (ke,pe,time, etc.)
64 >    call init_status(nose_status)
65 >    call init_dynamics_loop( nose_status )
66  
67 <  do_pot = .true.
47 <  not_done = .true.
48 <  nmflag  =  .false.
49 <  step = 0
67 >    !     Start the main simulation loop.
68  
69 <  dynamics: do  
70 <     if (.not.not_done) exit dynamics
69 >    do_pot = .true.
70 >    not_done = .true.
71 >    nmflag  =  .false.
72 >    step = 0
73  
74 <     step = step + 1
74 >    dynamics: do
75 >       if (.not.not_done) exit dynamics
76  
77 +       step = step + 1
78  
79 <     if (checktemp.and.(mod(step,check_temp_steps).eq.0)) then
79 >
80 >       if (checktemp.and.(mod(step,check_temp_steps).eq.0)) then
81            call calc_temp(nose_status%temp)
82 <        if (dabs(nose_status%temp-target_temp).gt.therm_variance) then
83 <           call scale_velocities(target_temp)
84 <        end if
85 <     end if
82 >          if (dabs(nose_status%temp-target_temp).gt.therm_variance) then
83 >             call scale_velocities(target_temp)
84 >          end if
85 >       end if
86  
87 <     call hooverb()
87 >       call nose_hoovera()
88  
89 <     call check(update_nlist)
89 >       call check(update_nlist)
90  
68 –  
69 –     if (do_pot) then
70 –        call calc_forces(update_nlist,nmflag,pe = nose_status%pot_e)
71 –     else
72 –        call calc_forces(update_nlist,nmflag)
73 –     endif
91  
92 <      
93 <     call hoovera()
94 <    
95 <     call evolve_time(step,nose_status,do_pot,not_done)
96 <
92 >       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  
81 – #ifdef MPI
82 –     call mpi_barrier(mpi_comm_world,mpi_err)
83 – #endif
84 –     end do dynamics
98  
99 <   end subroutine nose_hoover_dynamics
99 >       call nose_hooverb()
100  
101 +       call evolve_time(step,nose_status,do_pot,not_done)
102  
89 – subroutine hoovera()
90 –  use velocity, ONLY : remove_cm_momentum,remove_angular_momentum
103  
104 + #ifdef MPI
105 +       call mpi_barrier(mpi_comm_world,mpi_err)
106 + #endif
107 +    end do dynamics
108  
109 <  ! ** CONSTANT-TEMPERATURE MOLECULAR DYNAMICS USING CONSTRAINT.     **
94 <  ! **                                                               **
95 <  ! ** REFERENCES:                                                   **
96 <  ! **                                                               **
97 <  ! ** HOOVER, LADD, AND MORAN, PHYS REV LETT 48, 1818, 1982.        **
98 <  ! ** EVANS, J CHEM PHYS 78, 3297, 1983.                            **
99 <  ! ** BROWN AND CLARKE, MOL PHYS, 51, 1243, 1984.                   **
100 <  ! **                                                               **
101 <  ! ** ROUTINES SUPPLIED:                                            **
102 <  ! **                                                               **
103 <  ! ** SUBROUTINE HOOVERA ( )                                        **
104 <  ! **    FIRST PART OF MOVE WITH VELOCITY CONSTRAINTS APPLIED.      **
105 <  ! ** SUBROUTINE HOOVERB ( DT )                                     **
106 <  ! **    SECOND PART OF MOVE.                                       **
109 >  end subroutine nose_hoover_dynamics
110  
108 –
111  
112 <  integer                 ::  i, dim
113 <  real( kind = DP)        :: ke, v2, instatemp, chi
114 <  real( kind = DP)        :: kebar
113 <
114 <  real( kind = DP ), dimension(ndim) :: v_dim_i
112 > ! 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  
116 –  real(kind=DP) :: ke_local
116  
118 –  !     units for time are femtosec  (10^-15 sec)
119 –  !     units for distance are angstroms (10^-10 m)
120 –  !     units for velocity are angstroms femtosec^-1
121 –  !     units for mass are a.m.u. (1.661 * 10^-27 kg)
122 –  !     units for force are kcal mol^-1 angstrom^-1
123 –  !
124 –  !     converter will put the final terms into angstroms.
125 –  !       or angstrom/femtosecond.
117  
127 – !  converter = 1.0d0/2.390664d3
118  
119 +    integer                 ::  i, dim
120 +    real( kind = DP)        :: ke, v2, instatemp, chi
121 +    real( kind = DP)        :: kebar
122  
123 <  ! first calculate the temperature using the unconstrained velocities:
123 >    real( kind = DP ), dimension(ndim) :: v_dim_i
124  
125 <  ke    = 0.0E0_DP
133 <  ke_local = 0.0E0_DP
125 >    real(kind=DP) :: ke_local
126  
127 <  DO i = 1, nlocal
128 <     do dim = 1, ndim
129 <        v_dim_i(dim) = v(dim,i) + (dt2*f(dim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
130 <        ke_local = ke_local + mass(i)* &
131 <             (v_dim_i(dim)*v_dim_i(dim))
132 <     enddo
133 <  enddo
134 < #ifdef MPI  
143 <  call mpi_allreduce(ke_local,ke,1,mpi_double_precision, &
144 <       mpi_sum,mpi_comm_world,mpi_err)
145 < #else
146 <  ke = ke_local
147 < #endif
127 >    !     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  
136 +    !  converter = 1.0d0/2.390664d3
137  
150 –  ke = ke * AMUA2FS2_TO_KCALMOL
151 –  kebar = kcal_to_kj*1000.0e0_DP*ke / real(3*(natoms - 1),DP) / NAVOGADRO
152 –  instatemp = kebar/k_boltz
153 –  
154 –  !write(*,*) 'before, insta = ', instatemp
138  
139 <  chi = sqrt(target_temp / instatemp)
139 >    ! call nose_hoover chains
140 >    call apply_NHC_par()
141  
142 <  ! ** CALCULATE THE CONSTRAINED VELOCITIES AT TIME T+DT/2 **
142 >    ! advance velocities from t to t + dt/2
143 >    do i = 1, nlocal
144 >       v(1:ndim,i) = v(1:ndim,i) + &
145 >            (dt2*f(1:ndim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
146 >    end do
147  
148 <  do i = 1, nlocal
149 <     do dim = 1, ndim
150 <        v(dim,i) = v(dim,i) * (2.0E0_DP * chi - 1.0E0_DP) + &
151 <             chi*dt*(f(dim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
164 <     enddo
165 <  enddo
148 >    ! 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  
167 –  if (sim_type  == "cluster") then
168 –     call remove_cm_momentum()
169 –     call remove_angular_momentum()
170 –  endif
153  
154 < end subroutine hoovera
173 <
174 <      
175 <      
176 < subroutine hooverb( )
154 >  end subroutine nose_hoovera
155  
156  
157 <  ! *******************************************************************
158 <  ! ** SECOND PART OF THE CONSTANT TEMPERATURE ALGORITHM             **
159 <  ! **                                                               **
160 <  ! ** THIS ADVANCES THE POSITIONS FROM T TO T + DT.                 **
183 <  ! *******************************************************************
184 <
185 <
186 <  integer i, dim
187 <  !     units for time are femtosec  (10^-15 sec)
188 <  !     units for distance are angstroms (10^-10 m)
189 <  !     units for velocity are angstroms femtosec^-1
190 <  !     units for mass are a.m.u. (1.661 * 10^-27 kg)
191 <  !     units for force are kcal mol^-1 angstrom^-1
157 > ! 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  
162 <  do i = 1, nlocal
163 <     do dim = 1, ndim
164 <        q(dim,i) = q(dim,i) + dt * v(dim,i)
165 <     enddo
166 <  enddo
162 >    !     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  
168 < end subroutine hooverb
168 >    do i = 1, nlocal
169 >       v(1:ndim,i) = v(1:ndim,i) + &
170 >            (dt2*f(1:ndim,i)/mass(i))*KCALMOL_TO_AMUA2FS2
171 >    end do
172  
173 +    ! apply nose-hoover thermostat
174 +    call apply_NHC_par()
175  
176 +    if (sim_type  == "cluster") then
177 +       call remove_cm_momentum()
178 +       call remove_angular_momentum()
179 +    endif
180 +  end subroutine nose_hooverb
181 +
182 +
183 + ! 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 +
187 +  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 +
197 +    integer ::  i
198 +
199 +    call calc_temp(system_temp,system_ke)
200 +    
201 +    G_NkT = natoms*ndim*K_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_temp*2 - 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 = natoms*ndim*K_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)*aa*aa &
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) = (scale*scale*ke -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)*aa*aa &
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   end module dynamics_nose_hoover

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines