(git:34ef472)
extended_system_dynamics.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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  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  CALL timeset(routinen, handle)
73  map_info => nhc%map_info
74 
75  ! Compute the kinetic energy of the barostat
76  CALL ke_region_baro(map_info, npt, group)
77 
78  ! Calculate forces on the Nose-Hoover Thermostat and apply chains
79  CALL do_nhc(nhc, map_info)
80 
81  ! Now scale the particle velocities
82  CALL vel_rescale_baro(map_info, npt)
83 
84  CALL timestop(handle)
85  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  SUBROUTINE lnhc_particles(nhc, molecule_kind_set, molecule_set, &
106  particle_set, local_molecules, group, shell_adiabatic, &
107  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  CALL timeset(routinen, handle)
128  my_shell_adiabatic = .false.
129  IF (PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic
130  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  local_molecules, molecule_set, group, vel)
135 
136  ! Calculate forces on the Nose-Hoover Thermostat and apply chains
137  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  vel, shell_vel, core_vel)
143 
144  CALL timestop(handle)
145  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  SUBROUTINE lnhc_shells(nhc, atomic_kind_set, particle_set, local_particles, &
164  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  CALL timeset(routinen, handle)
182  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  group, core_particle_set, shell_particle_set, core_vel, shell_vel)
187 
188  ! Calculate forces on the Nose-Hoover Thermostat and apply chains
189  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  shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
194 
195  CALL timestop(handle)
196  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  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  DO n = 1, nhc%loc_num_nhc
213  imap = nhc%map_info%map_index(n)
214  IF (nhc%nvt(1, n)%nkt == 0.0_dp) cycle
215  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  CALL multiple_step_yoshida(nhc)
220 
221  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  SUBROUTINE shell_scale_comv(atomic_kind_set, local_particles, particle_set, &
236  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  nparticle_kind = SIZE(atomic_kind_set)
253 
254  DO iparticle_kind = 1, nparticle_kind
255  atomic_kind => atomic_kind_set(iparticle_kind)
256  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
257  shell_active=is_shell, shell=shell)
258  IF (is_shell) THEN
259  fac_masss = shell%mass_shell/mass
260  fac_massc = shell%mass_core/mass
261  nparticle_local = local_particles%n_el(iparticle_kind)
262  DO iparticle_local = 1, nparticle_local
263  iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
264  shell_index = particle_set(iparticle)%shell_index
265  vs(1:3) = shell_vel(1:3, shell_index)
266  vc(1:3) = core_vel(1:3, shell_index)
267  shell_vel(1, shell_index) = com_vel(1, iparticle) + fac_massc*(vs(1) - vc(1))
268  shell_vel(2, shell_index) = com_vel(2, iparticle) + fac_massc*(vs(2) - vc(2))
269  shell_vel(3, shell_index) = com_vel(3, iparticle) + fac_massc*(vs(3) - vc(3))
270  core_vel(1, shell_index) = com_vel(1, iparticle) + fac_masss*(vc(1) - vs(1))
271  core_vel(2, shell_index) = com_vel(2, iparticle) + fac_masss*(vc(2) - vs(2))
272  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  END SUBROUTINE shell_scale_comv
277 
278 ! **************************************************************************************************
279 !> \brief ...
280 !> \param nhc ...
281 !> \date 14-NOV-2000
282 !> \par History
283 !> none
284 ! **************************************************************************************************
285  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  nx1 = SIZE(nhc%nvt, 1)
294  nx2 = SIZE(nhc%nvt, 2)
295  map_info => nhc%map_info
296  ! perform multiple time stepping using Yoshida
297  ncloop: DO inc = 1, nhc%nc
298  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  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  DO n = 1, nhc%loc_num_nhc
306  IF (nhc%nvt(1, n)%nkt == 0.0_dp) cycle
307  DO inhc = nhc%nhc_len - 1, 1, -1
308  scale = exp(-0.125_dp*nhc%nvt(inhc + 1, n)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact)
309  nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
310  nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v + &
311  nhc%nvt(inhc, n)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact ! shift
312  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  0.5_dp*nhc%nvt(:, :)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact
320 
321  ! now accumulate the scale factor for particle velocities
322  DO n = 1, nhc%loc_num_nhc
323  imap = nhc%map_info%map_index(n)
324  IF (nhc%nvt(1, n)%nkt == 0.0_dp) cycle
325  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  DO n = 1, nhc%loc_num_nhc
331  imap = nhc%map_info%map_index(n)
332  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  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  DO inhc = 1, nhc%nhc_len - 1
339  DO n = 1, nhc%loc_num_nhc
340  IF (nhc%nvt(1, n)%nkt == 0.0_dp) cycle
341  scale = exp(-0.125_dp*nhc%nvt(inhc + 1, n)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact)
342  nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
343  nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v + &
344  nhc%nvt(inhc, n)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact ! shift
345  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  DO n = 1, nhc%loc_num_nhc
350  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  *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  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  END SUBROUTINE multiple_step_yoshida
361 
362 END MODULE extended_system_dynamics
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
extended_system_dynamics
Definition: extended_system_dynamics.F:15
extended_system_dynamics::lnhc_barostat
subroutine, public lnhc_barostat(nhc, npt, group)
...
Definition: extended_system_dynamics.F:61
extended_system_dynamics::shell_scale_comv
subroutine, public shell_scale_comv(atomic_kind_set, local_particles, particle_set, com_vel, shell_vel, core_vel)
...
Definition: extended_system_dynamics.F:237
extended_system_dynamics::lnhc_shells
subroutine, public lnhc_shells(nhc, atomic_kind_set, particle_set, local_particles, group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
...
Definition: extended_system_dynamics.F:165
extended_system_dynamics::lnhc_particles
subroutine, public lnhc_particles(nhc, molecule_kind_set, molecule_set, particle_set, local_molecules, group, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
...
Definition: extended_system_dynamics.F:108
extended_system_types
Lumps all possible extended system variables into one type for easy access and passing.
Definition: extended_system_types.F:17
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition: molecule_kind_types.F:16
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
shell_potential_types
Definition: shell_potential_types.F:11
thermostat_utils
Utilities for thermostats.
Definition: thermostat_utils.F:12
thermostat_utils::vel_rescale_shells
subroutine, public vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
...
Definition: thermostat_utils.F:1295
thermostat_utils::vel_rescale_particles
subroutine, public vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, local_molecules, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
...
Definition: thermostat_utils.F:1111
thermostat_utils::ke_region_shells
subroutine, public ke_region_shells(map_info, particle_set, atomic_kind_set, local_particles, group, core_particle_set, shell_particle_set, core_vel, shell_vel)
...
Definition: thermostat_utils.F:1216
thermostat_utils::ke_region_baro
subroutine, public ke_region_baro(map_info, npt, group)
...
Definition: thermostat_utils.F:922
thermostat_utils::vel_rescale_baro
subroutine, public vel_rescale_baro(map_info, npt)
...
Definition: thermostat_utils.F:951
thermostat_utils::ke_region_particles
subroutine, public ke_region_particles(map_info, particle_set, molecule_kind_set, local_molecules, molecule_set, group, vel)
...
Definition: thermostat_utils.F:980