Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> CJM 20-Feb-2001: Now npt_ifo is allocated to zero when not used
11 : !> CJM 11-apr-2001: adding routines to thermostat ao_type
12 : !> CJM 02-Aug-2003: renamed
13 : !> \author CJM
14 : ! **************************************************************************************************
15 : MODULE extended_system_dynamics
16 :
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind
19 : USE distribution_1d_types, ONLY: distribution_1d_type
20 : USE extended_system_types, ONLY: lnhc_parameters_type,&
21 : map_info_type,&
22 : npt_info_type
23 : USE kinds, ONLY: dp
24 : USE message_passing, ONLY: mp_comm_type
25 : USE molecule_kind_types, ONLY: molecule_kind_type
26 : USE molecule_types, ONLY: molecule_type
27 : USE particle_types, ONLY: particle_type
28 : USE shell_potential_types, ONLY: shell_kind_type
29 : USE thermostat_utils, ONLY: ke_region_baro,&
30 : ke_region_particles,&
31 : ke_region_shells,&
32 : vel_rescale_baro,&
33 : vel_rescale_particles,&
34 : vel_rescale_shells
35 : #include "../../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
41 : PUBLIC :: shell_scale_comv, &
42 : lnhc_particles, &
43 : lnhc_barostat, &
44 : lnhc_shells
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'extended_system_dynamics'
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief ...
52 : !> \param nhc ...
53 : !> \param npt ...
54 : !> \param group ...
55 : !> \date 13-DEC-2000
56 : !> \par History
57 : !> none
58 : !> \author CJM
59 : ! **************************************************************************************************
60 3276 : SUBROUTINE lnhc_barostat(nhc, npt, group)
61 :
62 : TYPE(lnhc_parameters_type), POINTER :: nhc
63 : TYPE(npt_info_type), DIMENSION(:, :), &
64 : INTENT(INOUT) :: npt
65 : TYPE(mp_comm_type), INTENT(IN) :: group
66 :
67 : CHARACTER(len=*), PARAMETER :: routineN = 'lnhc_barostat'
68 :
69 : INTEGER :: handle
70 : TYPE(map_info_type), POINTER :: map_info
71 :
72 3276 : CALL timeset(routineN, handle)
73 3276 : map_info => nhc%map_info
74 :
75 : ! Compute the kinetic energy of the barostat
76 3276 : CALL ke_region_baro(map_info, npt, group)
77 :
78 : ! Calculate forces on the Nose-Hoover Thermostat and apply chains
79 3276 : CALL do_nhc(nhc, map_info)
80 :
81 : ! Now scale the particle velocities
82 3276 : CALL vel_rescale_baro(map_info, npt)
83 :
84 3276 : CALL timestop(handle)
85 3276 : END SUBROUTINE lnhc_barostat
86 :
87 : ! **************************************************************************************************
88 : !> \brief ...
89 : !> \param nhc ...
90 : !> \param molecule_kind_set ...
91 : !> \param molecule_set ...
92 : !> \param particle_set ...
93 : !> \param local_molecules ...
94 : !> \param group ...
95 : !> \param shell_adiabatic ...
96 : !> \param shell_particle_set ...
97 : !> \param core_particle_set ...
98 : !> \param vel ...
99 : !> \param shell_vel ...
100 : !> \param core_vel ...
101 : !> \date 14-NOV-2000
102 : !> \par History
103 : !> none
104 : ! **************************************************************************************************
105 26536 : SUBROUTINE lnhc_particles(nhc, molecule_kind_set, molecule_set, &
106 : particle_set, local_molecules, group, shell_adiabatic, &
107 13268 : shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
108 :
109 : TYPE(lnhc_parameters_type), POINTER :: nhc
110 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
111 : TYPE(molecule_type), POINTER :: molecule_set(:)
112 : TYPE(particle_type), POINTER :: particle_set(:)
113 : TYPE(distribution_1d_type), POINTER :: local_molecules
114 : TYPE(mp_comm_type), INTENT(IN) :: group
115 : LOGICAL, INTENT(IN), OPTIONAL :: shell_adiabatic
116 : TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
117 : core_particle_set(:)
118 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
119 : core_vel(:, :)
120 :
121 : CHARACTER(len=*), PARAMETER :: routineN = 'lnhc_particles'
122 :
123 : INTEGER :: handle
124 : LOGICAL :: my_shell_adiabatic
125 : TYPE(map_info_type), POINTER :: map_info
126 :
127 13268 : CALL timeset(routineN, handle)
128 13268 : my_shell_adiabatic = .FALSE.
129 13268 : IF (PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic
130 13268 : map_info => nhc%map_info
131 :
132 : ! Compute the kinetic energy for the region to thermostat
133 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
134 19902 : local_molecules, molecule_set, group, vel)
135 :
136 : ! Calculate forces on the Nose-Hoover Thermostat and apply chains
137 13268 : CALL do_nhc(nhc, map_info)
138 :
139 : ! Now scale the particle velocities
140 : CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, &
141 : local_molecules, my_shell_adiabatic, shell_particle_set, core_particle_set, &
142 43778 : vel, shell_vel, core_vel)
143 :
144 13268 : CALL timestop(handle)
145 13268 : END SUBROUTINE lnhc_particles
146 :
147 : ! **************************************************************************************************
148 : !> \brief ...
149 : !> \param nhc ...
150 : !> \param atomic_kind_set ...
151 : !> \param particle_set ...
152 : !> \param local_particles ...
153 : !> \param group ...
154 : !> \param shell_particle_set ...
155 : !> \param core_particle_set ...
156 : !> \param vel ...
157 : !> \param shell_vel ...
158 : !> \param core_vel ...
159 : !> \date 14-NOV-2000
160 : !> \par History
161 : !> none
162 : ! **************************************************************************************************
163 2720 : SUBROUTINE lnhc_shells(nhc, atomic_kind_set, particle_set, local_particles, &
164 1360 : group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
165 :
166 : TYPE(lnhc_parameters_type), POINTER :: nhc
167 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
168 : TYPE(particle_type), POINTER :: particle_set(:)
169 : TYPE(distribution_1d_type), POINTER :: local_particles
170 : TYPE(mp_comm_type), INTENT(IN) :: group
171 : TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
172 : core_particle_set(:)
173 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), &
174 : core_vel(:, :)
175 :
176 : CHARACTER(len=*), PARAMETER :: routineN = 'lnhc_shells'
177 :
178 : INTEGER :: handle
179 : TYPE(map_info_type), POINTER :: map_info
180 :
181 1360 : CALL timeset(routineN, handle)
182 1360 : map_info => nhc%map_info
183 :
184 : ! Compute the kinetic energy of the region to thermostat
185 : CALL ke_region_shells(map_info, particle_set, atomic_kind_set, local_particles, &
186 2720 : group, core_particle_set, shell_particle_set, core_vel, shell_vel)
187 :
188 : ! Calculate forces on the Nose-Hoover Thermostat and apply chains
189 1360 : CALL do_nhc(nhc, map_info)
190 :
191 : ! Now scale the particle velocities
192 : CALL vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, &
193 3400 : shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
194 :
195 1360 : CALL timestop(handle)
196 1360 : END SUBROUTINE lnhc_shells
197 :
198 : ! **************************************************************************************************
199 : !> \brief ...
200 : !> \param nhc ...
201 : !> \param map_info ...
202 : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
203 : ! **************************************************************************************************
204 17904 : SUBROUTINE do_nhc(nhc, map_info)
205 : TYPE(lnhc_parameters_type), POINTER :: nhc
206 : TYPE(map_info_type), POINTER :: map_info
207 :
208 : INTEGER :: imap, n
209 :
210 : ! Force on the first bead in every thermostat chain
211 :
212 619332 : DO n = 1, nhc%loc_num_nhc
213 601428 : imap = nhc%map_info%map_index(n)
214 601428 : IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
215 619332 : nhc%nvt(1, n)%f = (map_info%s_kin(imap) - nhc%nvt(1, n)%nkt)/nhc%nvt(1, n)%mass
216 : END DO
217 :
218 : ! Perform multiple time stepping using Yoshida
219 17904 : CALL multiple_step_yoshida(nhc)
220 :
221 17904 : END SUBROUTINE do_nhc
222 :
223 : ! **************************************************************************************************
224 : !> \brief ...
225 : !> \param atomic_kind_set ...
226 : !> \param local_particles ...
227 : !> \param particle_set ...
228 : !> \param com_vel ...
229 : !> \param shell_vel ...
230 : !> \param core_vel ...
231 : !> \date 14-NOV-2000
232 : !> \par History
233 : !> none
234 : ! **************************************************************************************************
235 80 : SUBROUTINE shell_scale_comv(atomic_kind_set, local_particles, particle_set, &
236 80 : com_vel, shell_vel, core_vel)
237 :
238 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
239 : TYPE(distribution_1d_type), POINTER :: local_particles
240 : TYPE(particle_type), POINTER :: particle_set(:)
241 : REAL(KIND=dp), INTENT(IN) :: com_vel(:, :)
242 : REAL(KIND=dp), INTENT(INOUT) :: shell_vel(:, :), core_vel(:, :)
243 :
244 : INTEGER :: iparticle, iparticle_kind, &
245 : iparticle_local, nparticle_kind, &
246 : nparticle_local, shell_index
247 : LOGICAL :: is_shell
248 : REAL(KIND=dp) :: fac_massc, fac_masss, mass, vc(3), vs(3)
249 : TYPE(atomic_kind_type), POINTER :: atomic_kind
250 : TYPE(shell_kind_type), POINTER :: shell
251 :
252 80 : nparticle_kind = SIZE(atomic_kind_set)
253 :
254 240 : DO iparticle_kind = 1, nparticle_kind
255 160 : atomic_kind => atomic_kind_set(iparticle_kind)
256 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
257 160 : shell_active=is_shell, shell=shell)
258 240 : IF (is_shell) THEN
259 160 : fac_masss = shell%mass_shell/mass
260 160 : fac_massc = shell%mass_core/mass
261 160 : nparticle_local = local_particles%n_el(iparticle_kind)
262 4000 : DO iparticle_local = 1, nparticle_local
263 3840 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
264 3840 : shell_index = particle_set(iparticle)%shell_index
265 15360 : vs(1:3) = shell_vel(1:3, shell_index)
266 15360 : vc(1:3) = core_vel(1:3, shell_index)
267 3840 : shell_vel(1, shell_index) = com_vel(1, iparticle) + fac_massc*(vs(1) - vc(1))
268 3840 : shell_vel(2, shell_index) = com_vel(2, iparticle) + fac_massc*(vs(2) - vc(2))
269 3840 : shell_vel(3, shell_index) = com_vel(3, iparticle) + fac_massc*(vs(3) - vc(3))
270 3840 : core_vel(1, shell_index) = com_vel(1, iparticle) + fac_masss*(vc(1) - vs(1))
271 3840 : core_vel(2, shell_index) = com_vel(2, iparticle) + fac_masss*(vc(2) - vs(2))
272 4000 : core_vel(3, shell_index) = com_vel(3, iparticle) + fac_masss*(vc(3) - vs(3))
273 : END DO
274 : END IF ! is_shell
275 : END DO ! iparticle_kind
276 80 : END SUBROUTINE shell_scale_comv
277 :
278 : ! **************************************************************************************************
279 : !> \brief ...
280 : !> \param nhc ...
281 : !> \date 14-NOV-2000
282 : !> \par History
283 : !> none
284 : ! **************************************************************************************************
285 17904 : SUBROUTINE multiple_step_yoshida(nhc)
286 :
287 : TYPE(lnhc_parameters_type), POINTER :: nhc
288 :
289 : INTEGER :: imap, inc, inhc, iyosh, n, nx1, nx2
290 : REAL(KIND=dp) :: scale
291 : TYPE(map_info_type), POINTER :: map_info
292 :
293 17904 : nx1 = SIZE(nhc%nvt, 1)
294 17904 : nx2 = SIZE(nhc%nvt, 2)
295 17904 : map_info => nhc%map_info
296 : ! perform multiple time stepping using Yoshida
297 53712 : NCLOOP: DO inc = 1, nhc%nc
298 161136 : YOSH: DO iyosh = 1, nhc%nyosh
299 :
300 : ! update velocity on the last thermostat in the chain ! O1
301 : nhc%nvt(nhc%nhc_len, :)%v = nhc%nvt(nhc%nhc_len, :)%v + &
302 3715992 : nhc%nvt(nhc%nhc_len, :)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact
303 :
304 : ! update velocity of other thermostats on chain (from nhc_len-1 to 1) ! O2
305 3715992 : DO n = 1, nhc%loc_num_nhc
306 3608568 : IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
307 13485504 : DO inhc = nhc%nhc_len - 1, 1, -1
308 9813672 : scale = EXP(-0.125_dp*nhc%nvt(inhc + 1, n)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact)
309 9813672 : nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
310 : nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v + &
311 9813672 : nhc%nvt(inhc, n)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact ! shift
312 13422240 : nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
313 : END DO
314 : END DO
315 :
316 : ! the core of the operator ----- START------
317 : ! update nhc positions
318 : nhc%nvt(:, :)%eta = nhc%nvt(:, :)%eta + &
319 17226552 : 0.5_dp*nhc%nvt(:, :)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact
320 :
321 : ! now accumulate the scale factor for particle velocities
322 3715992 : DO n = 1, nhc%loc_num_nhc
323 3608568 : imap = nhc%map_info%map_index(n)
324 3608568 : IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
325 3715992 : map_info%v_scale(imap) = map_info%v_scale(imap)*EXP(-0.5_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact*nhc%nvt(1, n)%v)
326 : END DO
327 : ! the core of the operator ------ END ------
328 :
329 : ! update the force on first thermostat again (since particle velocities changed)
330 3715992 : DO n = 1, nhc%loc_num_nhc
331 3608568 : imap = nhc%map_info%map_index(n)
332 3608568 : IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
333 : nhc%nvt(1, n)%f = (map_info%s_kin(imap)*map_info%v_scale(imap)* &
334 3715992 : map_info%v_scale(imap) - nhc%nvt(1, n)%nkt)/nhc%nvt(1, n)%mass
335 : END DO
336 :
337 : ! update velocity of other thermostats on chain (from 1 to nhc_len-1) ! O2
338 373128 : DO inhc = 1, nhc%nhc_len - 1
339 10167696 : DO n = 1, nhc%loc_num_nhc
340 9901992 : IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
341 9813672 : scale = EXP(-0.125_dp*nhc%nvt(inhc + 1, n)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact)
342 9813672 : nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
343 : nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v + &
344 9813672 : nhc%nvt(inhc, n)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact ! shift
345 10167696 : nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
346 : END DO
347 :
348 : ! updating the forces on all the thermostats
349 10275120 : DO n = 1, nhc%loc_num_nhc
350 9901992 : IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
351 : nhc%nvt(inhc + 1, n)%f = (nhc%nvt(inhc, n)%mass*nhc%nvt(inhc, n)%v &
352 10167696 : *nhc%nvt(inhc, n)%v - nhc%nvt(inhc + 1, n)%nkt)/nhc%nvt(inhc + 1, n)%mass
353 : END DO
354 : END DO
355 : ! update velocity on last thermostat ! O1
356 : nhc%nvt(nhc%nhc_len, :)%v = nhc%nvt(nhc%nhc_len, :)%v + &
357 3751800 : nhc%nvt(nhc%nhc_len, :)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact
358 : END DO YOSH
359 : END DO NCLOOP
360 17904 : END SUBROUTINE multiple_step_yoshida
361 :
362 : END MODULE extended_system_dynamics
|
