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 |
|
|
(dt2*f(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 |
|
|
(dt2*f(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 = natoms*ndim * 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_temp*2 - 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 = natoms*ndim*K_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)*aa*aa & |
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) = (scale*scale*ke -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)*aa*aa & |
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 |