Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \author Noam Bernstein [noamb] 02.2012
10 : ! **************************************************************************************************
11 : MODULE al_system_dynamics
12 :
13 : USE al_system_types, ONLY: al_system_type
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE constraint_fxd, ONLY: fix_atom_control
17 : USE distribution_1d_types, ONLY: distribution_1d_type
18 : USE extended_system_types, ONLY: map_info_type
19 : USE force_env_types, ONLY: force_env_type
20 : USE kinds, ONLY: dp
21 : USE message_passing, ONLY: mp_comm_type
22 : USE molecule_kind_types, ONLY: molecule_kind_type
23 : USE molecule_types, ONLY: get_molecule,&
24 : molecule_type
25 : USE particle_types, ONLY: particle_type
26 : USE thermostat_utils, ONLY: ke_region_particles,&
27 : vel_rescale_particles
28 : #include "../../base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
34 : PUBLIC :: al_particles
35 :
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'al_system_dynamics'
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief ...
42 : !> \param al ...
43 : !> \param force_env ...
44 : !> \param molecule_kind_set ...
45 : !> \param molecule_set ...
46 : !> \param particle_set ...
47 : !> \param local_molecules ...
48 : !> \param local_particles ...
49 : !> \param group ...
50 : !> \param vel ...
51 : !> \author Noam Bernstein [noamb] 02.2012
52 : ! **************************************************************************************************
53 32 : SUBROUTINE al_particles(al, force_env, molecule_kind_set, molecule_set, &
54 16 : particle_set, local_molecules, local_particles, group, vel)
55 :
56 : TYPE(al_system_type), POINTER :: al
57 : TYPE(force_env_type), POINTER :: force_env
58 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
59 : TYPE(molecule_type), POINTER :: molecule_set(:)
60 : TYPE(particle_type), POINTER :: particle_set(:)
61 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
62 : TYPE(mp_comm_type), INTENT(IN) :: group
63 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :)
64 :
65 : CHARACTER(len=*), PARAMETER :: routineN = 'al_particles'
66 :
67 : INTEGER :: handle
68 : LOGICAL :: my_shell_adiabatic
69 : TYPE(map_info_type), POINTER :: map_info
70 :
71 16 : CALL timeset(routineN, handle)
72 16 : my_shell_adiabatic = .FALSE.
73 16 : map_info => al%map_info
74 :
75 : IF (debug_this_module) THEN
76 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "INIT")
77 : END IF
78 :
79 16 : IF (al%tau_nh <= 0.0_dp) THEN
80 : CALL al_OU_step(0.5_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
81 0 : particle_set, local_molecules, local_particles, vel)
82 : IF (debug_this_module) THEN
83 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post OU")
84 : END IF
85 : ELSE
86 : ! quarter step of Langevin using Ornstein-Uhlenbeck
87 : CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
88 24 : particle_set, local_molecules, local_particles, vel)
89 : IF (debug_this_module) THEN
90 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 1st OU")
91 : END IF
92 :
93 : ! Compute the kinetic energy for the region to thermostat for the (T dependent chi step)
94 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
95 24 : local_molecules, molecule_set, group, vel=vel)
96 : ! quarter step of chi, and set vel drag factors for a half step
97 16 : CALL al_NH_quarter_step(al, map_info, set_half_step_vel_factors=.TRUE.)
98 :
99 : ! Now scale the particle velocities for a NH half step
100 : CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, &
101 24 : local_molecules, my_shell_adiabatic, vel=vel)
102 : ! Recompute the kinetic energy for the region to thermostat (for the T dependent chi step)
103 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
104 24 : local_molecules, molecule_set, group, vel=vel)
105 : IF (debug_this_module) THEN
106 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post rescale_vel")
107 : END IF
108 :
109 : ! quarter step of chi
110 16 : CALL al_NH_quarter_step(al, map_info, set_half_step_vel_factors=.FALSE.)
111 :
112 : ! quarter step of Langevin using Ornstein-Uhlenbeck
113 : CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
114 24 : particle_set, local_molecules, local_particles, vel)
115 : IF (debug_this_module) THEN
116 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 2nd OU")
117 : END IF
118 : END IF
119 :
120 : ! Recompute the final kinetic energy for the region to thermostat
121 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
122 24 : local_molecules, molecule_set, group, vel=vel)
123 :
124 16 : CALL timestop(handle)
125 16 : END SUBROUTINE al_particles
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param molecule_kind_set ...
130 : !> \param molecule_set ...
131 : !> \param local_molecules ...
132 : !> \param particle_set ...
133 : !> \param vel ...
134 : !> \param label ...
135 : ! **************************************************************************************************
136 0 : SUBROUTINE dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, label)
137 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
138 : TYPE(molecule_type), POINTER :: molecule_set(:)
139 : TYPE(distribution_1d_type), POINTER :: local_molecules
140 : TYPE(particle_type), POINTER :: particle_set(:)
141 : REAL(dp), OPTIONAL :: vel(:, :)
142 : CHARACTER(len=*) :: label
143 :
144 : INTEGER :: first_atom, ikind, imol, imol_local, &
145 : ipart, last_atom, nmol_local
146 : TYPE(molecule_type), POINTER :: molecule
147 :
148 0 : DO ikind = 1, SIZE(molecule_kind_set)
149 0 : nmol_local = local_molecules%n_el(ikind)
150 0 : DO imol_local = 1, nmol_local
151 0 : imol = local_molecules%list(ikind)%array(imol_local)
152 0 : molecule => molecule_set(imol)
153 0 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
154 0 : DO ipart = first_atom, last_atom
155 0 : IF (PRESENT(vel)) THEN
156 0 : WRITE (unit=*, fmt='("VEL ",A20," IPART ",I6," V ",3F20.10)') TRIM(label), ipart, vel(:, ipart)
157 : ELSE
158 0 : WRITE (unit=*, fmt='("PARTICLE_SET%VEL ",A20," IPART ",I6," V ",3F20.10)') TRIM(label), &
159 0 : ipart, particle_set(ipart)%v(:)
160 : END IF
161 : END DO
162 : END DO
163 : END DO
164 0 : END SUBROUTINE dump_vel
165 :
166 : ! **************************************************************************************************
167 : !> \brief ...
168 : !> \param step ...
169 : !> \param al ...
170 : !> \param force_env ...
171 : !> \param map_info ...
172 : !> \param molecule_kind_set ...
173 : !> \param molecule_set ...
174 : !> \param particle_set ...
175 : !> \param local_molecules ...
176 : !> \param local_particles ...
177 : !> \param vel ...
178 : ! **************************************************************************************************
179 32 : SUBROUTINE al_OU_step(step, al, force_env, map_info, molecule_kind_set, molecule_set, &
180 32 : particle_set, local_molecules, local_particles, vel)
181 : REAL(dp), INTENT(in) :: step
182 : TYPE(al_system_type), POINTER :: al
183 : TYPE(force_env_type), POINTER :: force_env
184 : TYPE(map_info_type), POINTER :: map_info
185 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
186 : TYPE(molecule_type), POINTER :: molecule_set(:)
187 : TYPE(particle_type), POINTER :: particle_set(:)
188 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
189 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :)
190 :
191 : INTEGER :: first_atom, i, ii, ikind, imap, imol, imol_local, ipart, iparticle_kind, &
192 : iparticle_local, jj, last_atom, nmol_local, nparticle, nparticle_kind, nparticle_local
193 : LOGICAL :: check, present_vel
194 : REAL(KIND=dp) :: mass
195 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: w
196 : TYPE(atomic_kind_type), POINTER :: atomic_kind
197 : TYPE(molecule_type), POINTER :: molecule
198 :
199 32 : present_vel = PRESENT(vel)
200 :
201 : ![NB] not a big deal, but could this be done once at init time?
202 396776 : DO i = 1, al%loc_num_al
203 396744 : imap = map_info%map_index(i)
204 : ! drag on velocities
205 396776 : IF (al%tau_langevin > 0.0_dp) THEN
206 396744 : map_info%v_scale(imap) = EXP(-step*al%dt/al%tau_langevin)
207 396744 : map_info%s_kin(imap) = SQRT((al%nvt(i)%nkt/al%nvt(i)%degrees_of_freedom)*(1.0_dp - map_info%v_scale(imap)**2))
208 : ELSE
209 0 : map_info%v_scale(imap) = 1.0_dp
210 0 : map_info%s_kin(imap) = 0.0_dp
211 : END IF
212 : ! magnitude of random force, not including 1/sqrt(mass) part
213 : END DO
214 :
215 32 : nparticle = SIZE(particle_set)
216 32 : nparticle_kind = SIZE(local_particles%n_el)
217 96 : ALLOCATE (w(3, nparticle))
218 32 : w(:, :) = 0.0_dp
219 32 : check = (nparticle_kind <= SIZE(local_particles%n_el) .AND. nparticle_kind <= SIZE(local_particles%list))
220 0 : CPASSERT(check)
221 32 : check = ASSOCIATED(local_particles%local_particle_set)
222 32 : CPASSERT(check)
223 5152 : DO iparticle_kind = 1, nparticle_kind
224 5120 : nparticle_local = local_particles%n_el(iparticle_kind)
225 5120 : check = (nparticle_local <= SIZE(local_particles%list(iparticle_kind)%array))
226 5120 : CPASSERT(check)
227 137400 : DO iparticle_local = 1, nparticle_local
228 132248 : ipart = local_particles%list(iparticle_kind)%array(iparticle_local)
229 132248 : w(1, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
230 132248 : w(2, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
231 137368 : w(3, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
232 : END DO
233 : END DO
234 :
235 32 : CALL fix_atom_control(force_env, w)
236 :
237 32 : ii = 0
238 61320 : DO ikind = 1, SIZE(molecule_kind_set)
239 61288 : nmol_local = local_molecules%n_el(ikind)
240 101200 : DO imol_local = 1, nmol_local
241 39880 : imol = local_molecules%list(ikind)%array(imol_local)
242 39880 : molecule => molecule_set(imol)
243 39880 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
244 233416 : DO ipart = first_atom, last_atom
245 132248 : ii = ii + 1
246 132248 : atomic_kind => particle_set(ipart)%atomic_kind
247 132248 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
248 172128 : IF (present_vel) THEN
249 264496 : DO jj = 1, 3
250 : vel(jj, ipart) = vel(jj, ipart)*map_info%p_scale(jj, ii)%point + &
251 264496 : map_info%p_kin(jj, ii)%point/SQRT(mass)*w(jj, ipart)
252 : END DO
253 : ELSE
254 264496 : DO jj = 1, 3
255 : particle_set(ipart)%v(jj) = particle_set(ipart)%v(jj)*map_info%p_scale(jj, ii)%point + &
256 264496 : map_info%p_kin(jj, ii)%point/SQRT(mass)*w(jj, ipart)
257 : END DO
258 : END IF
259 : END DO
260 : END DO
261 : END DO
262 :
263 32 : DEALLOCATE (w)
264 :
265 32 : END SUBROUTINE al_OU_step
266 :
267 : ! **************************************************************************************************
268 : !> \brief ...
269 : !> \param al ...
270 : !> \param map_info ...
271 : !> \param set_half_step_vel_factors ...
272 : !> \author Noam Bernstein [noamb] 02.2012
273 : ! **************************************************************************************************
274 32 : SUBROUTINE al_NH_quarter_step(al, map_info, set_half_step_vel_factors)
275 : TYPE(al_system_type), POINTER :: al
276 : TYPE(map_info_type), POINTER :: map_info
277 : LOGICAL, INTENT(in) :: set_half_step_vel_factors
278 :
279 : INTEGER :: i, imap
280 : REAL(KIND=dp) :: decay, delta_K
281 :
282 : ![NB] how to deal with dt_fact?
283 :
284 396776 : DO i = 1, al%loc_num_al
285 396776 : IF (al%nvt(i)%mass > 0.0_dp) THEN
286 396744 : imap = map_info%map_index(i)
287 396744 : delta_K = 0.5_dp*(map_info%s_kin(imap) - al%nvt(i)%nkt)
288 396744 : al%nvt(i)%chi = al%nvt(i)%chi + 0.5_dp*al%dt*delta_K/al%nvt(i)%mass
289 396744 : IF (set_half_step_vel_factors) THEN
290 198372 : decay = EXP(-0.5_dp*al%dt*al%nvt(i)%chi)
291 198372 : map_info%v_scale(imap) = decay
292 : END IF
293 : ELSE
294 0 : al%nvt(i)%chi = 0.0_dp
295 0 : IF (set_half_step_vel_factors) THEN
296 0 : map_info%v_scale(imap) = 1.0_dp
297 : END IF
298 : END IF
299 : END DO
300 :
301 32 : END SUBROUTINE al_NH_quarter_step
302 :
303 : END MODULE al_system_dynamics
|
