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 : !> \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) &
76 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "INIT")
77 :
78 16 : IF (al%tau_nh <= 0.0_dp) THEN
79 : CALL al_OU_step(0.5_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
80 0 : particle_set, local_molecules, local_particles, vel)
81 : IF (debug_this_module) &
82 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post OU")
83 : ELSE
84 : ! quarter step of Langevin using Ornstein-Uhlenbeck
85 : CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
86 24 : particle_set, local_molecules, local_particles, vel)
87 : IF (debug_this_module) &
88 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 1st OU")
89 :
90 : ! Compute the kinetic energy for the region to thermostat for the (T dependent chi step)
91 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
92 24 : local_molecules, molecule_set, group, vel=vel)
93 : ! quarter step of chi, and set vel drag factors for a half step
94 16 : CALL al_NH_quarter_step(al, map_info, set_half_step_vel_factors=.TRUE.)
95 :
96 : ! Now scale the particle velocities for a NH half step
97 : CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, &
98 24 : local_molecules, my_shell_adiabatic, vel=vel)
99 : ! Recompute the kinetic energy for the region to thermostat (for the T dependent chi step)
100 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
101 24 : local_molecules, molecule_set, group, vel=vel)
102 : IF (debug_this_module) &
103 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post rescale_vel")
104 :
105 : ! quarter step of chi
106 16 : CALL al_NH_quarter_step(al, map_info, set_half_step_vel_factors=.FALSE.)
107 :
108 : ! quarter step of Langevin using Ornstein-Uhlenbeck
109 : CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
110 24 : particle_set, local_molecules, local_particles, vel)
111 : IF (debug_this_module) &
112 : CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 2nd OU")
113 : END IF
114 :
115 : ! Recompute the final kinetic energy for the region to thermostat
116 : CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
117 24 : local_molecules, molecule_set, group, vel=vel)
118 :
119 16 : CALL timestop(handle)
120 16 : END SUBROUTINE al_particles
121 :
122 : ! **************************************************************************************************
123 : !> \brief ...
124 : !> \param molecule_kind_set ...
125 : !> \param molecule_set ...
126 : !> \param local_molecules ...
127 : !> \param particle_set ...
128 : !> \param vel ...
129 : !> \param label ...
130 : ! **************************************************************************************************
131 0 : SUBROUTINE dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, label)
132 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
133 : TYPE(molecule_type), POINTER :: molecule_set(:)
134 : TYPE(distribution_1d_type), POINTER :: local_molecules
135 : TYPE(particle_type), POINTER :: particle_set(:)
136 : REAL(dp), OPTIONAL :: vel(:, :)
137 : CHARACTER(len=*) :: label
138 :
139 : INTEGER :: first_atom, ikind, imol, imol_local, &
140 : ipart, last_atom, nmol_local
141 : TYPE(molecule_type), POINTER :: molecule
142 :
143 0 : DO ikind = 1, SIZE(molecule_kind_set)
144 0 : nmol_local = local_molecules%n_el(ikind)
145 0 : DO imol_local = 1, nmol_local
146 0 : imol = local_molecules%list(ikind)%array(imol_local)
147 0 : molecule => molecule_set(imol)
148 0 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
149 0 : DO ipart = first_atom, last_atom
150 0 : IF (PRESENT(vel)) THEN
151 0 : WRITE (unit=*, fmt='("VEL ",A20," IPART ",I6," V ",3F20.10)') TRIM(label), ipart, vel(:, ipart)
152 : ELSE
153 0 : WRITE (unit=*, fmt='("PARTICLE_SET%VEL ",A20," IPART ",I6," V ",3F20.10)') TRIM(label), &
154 0 : ipart, particle_set(ipart)%v(:)
155 : END IF
156 : END DO
157 : END DO
158 : END DO
159 0 : END SUBROUTINE dump_vel
160 :
161 : ! **************************************************************************************************
162 : !> \brief ...
163 : !> \param step ...
164 : !> \param al ...
165 : !> \param force_env ...
166 : !> \param map_info ...
167 : !> \param molecule_kind_set ...
168 : !> \param molecule_set ...
169 : !> \param particle_set ...
170 : !> \param local_molecules ...
171 : !> \param local_particles ...
172 : !> \param vel ...
173 : ! **************************************************************************************************
174 32 : SUBROUTINE al_OU_step(step, al, force_env, map_info, molecule_kind_set, molecule_set, &
175 32 : particle_set, local_molecules, local_particles, vel)
176 : REAL(dp), INTENT(in) :: step
177 : TYPE(al_system_type), POINTER :: al
178 : TYPE(force_env_type), POINTER :: force_env
179 : TYPE(map_info_type), POINTER :: map_info
180 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
181 : TYPE(molecule_type), POINTER :: molecule_set(:)
182 : TYPE(particle_type), POINTER :: particle_set(:)
183 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
184 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :)
185 :
186 : INTEGER :: first_atom, i, ii, ikind, imap, imol, imol_local, ipart, iparticle_kind, &
187 : iparticle_local, jj, last_atom, nmol_local, nparticle, nparticle_kind, nparticle_local
188 : LOGICAL :: check, present_vel
189 : REAL(KIND=dp) :: mass
190 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: w(:, :)
191 : TYPE(atomic_kind_type), POINTER :: atomic_kind
192 : TYPE(molecule_type), POINTER :: molecule
193 :
194 32 : present_vel = PRESENT(vel)
195 :
196 : ![NB] not a big deal, but could this be done once at init time?
197 396776 : DO i = 1, al%loc_num_al
198 396744 : imap = map_info%map_index(i)
199 : ! drag on velocities
200 396776 : IF (al%tau_langevin > 0.0_dp) THEN
201 396744 : map_info%v_scale(imap) = EXP(-step*al%dt/al%tau_langevin)
202 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))
203 : ELSE
204 0 : map_info%v_scale(imap) = 1.0_dp
205 0 : map_info%s_kin(imap) = 0.0_dp
206 : END IF
207 : ! magnitude of random force, not including 1/sqrt(mass) part
208 : END DO
209 :
210 32 : nparticle = SIZE(particle_set)
211 32 : nparticle_kind = SIZE(local_particles%n_el)
212 96 : ALLOCATE (w(3, nparticle))
213 1058016 : w(:, :) = 0.0_dp
214 32 : check = (nparticle_kind <= SIZE(local_particles%n_el) .AND. nparticle_kind <= SIZE(local_particles%list))
215 0 : CPASSERT(check)
216 32 : check = ASSOCIATED(local_particles%local_particle_set)
217 32 : CPASSERT(check)
218 5152 : DO iparticle_kind = 1, nparticle_kind
219 5120 : nparticle_local = local_particles%n_el(iparticle_kind)
220 5120 : check = (nparticle_local <= SIZE(local_particles%list(iparticle_kind)%array))
221 5120 : CPASSERT(check)
222 137400 : DO iparticle_local = 1, nparticle_local
223 132248 : ipart = local_particles%list(iparticle_kind)%array(iparticle_local)
224 132248 : w(1, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
225 132248 : w(2, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
226 137368 : w(3, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
227 : END DO
228 : END DO
229 :
230 32 : CALL fix_atom_control(force_env, w)
231 :
232 32 : ii = 0
233 61320 : DO ikind = 1, SIZE(molecule_kind_set)
234 61288 : nmol_local = local_molecules%n_el(ikind)
235 101200 : DO imol_local = 1, nmol_local
236 39880 : imol = local_molecules%list(ikind)%array(imol_local)
237 39880 : molecule => molecule_set(imol)
238 39880 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
239 233416 : DO ipart = first_atom, last_atom
240 132248 : ii = ii + 1
241 132248 : atomic_kind => particle_set(ipart)%atomic_kind
242 132248 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
243 172128 : IF (present_vel) THEN
244 264496 : DO jj = 1, 3
245 : vel(jj, ipart) = vel(jj, ipart)*map_info%p_scale(jj, ii)%point + &
246 264496 : map_info%p_kin(jj, ii)%point/SQRT(mass)*w(jj, ipart)
247 : END DO
248 : ELSE
249 264496 : DO jj = 1, 3
250 : particle_set(ipart)%v(jj) = particle_set(ipart)%v(jj)*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 : END IF
254 : END DO
255 : END DO
256 : END DO
257 :
258 32 : DEALLOCATE (w)
259 :
260 32 : END SUBROUTINE al_OU_step
261 :
262 : ! **************************************************************************************************
263 : !> \brief ...
264 : !> \param al ...
265 : !> \param map_info ...
266 : !> \param set_half_step_vel_factors ...
267 : !> \author Noam Bernstein [noamb] 02.2012
268 : ! **************************************************************************************************
269 32 : SUBROUTINE al_NH_quarter_step(al, map_info, set_half_step_vel_factors)
270 : TYPE(al_system_type), POINTER :: al
271 : TYPE(map_info_type), POINTER :: map_info
272 : LOGICAL, INTENT(in) :: set_half_step_vel_factors
273 :
274 : INTEGER :: i, imap
275 : REAL(KIND=dp) :: decay, delta_K
276 :
277 : ![NB] how to deal with dt_fact?
278 :
279 396776 : DO i = 1, al%loc_num_al
280 396776 : IF (al%nvt(i)%mass > 0.0_dp) THEN
281 396744 : imap = map_info%map_index(i)
282 396744 : delta_K = 0.5_dp*(map_info%s_kin(imap) - al%nvt(i)%nkt)
283 396744 : al%nvt(i)%chi = al%nvt(i)%chi + 0.5_dp*al%dt*delta_K/al%nvt(i)%mass
284 396744 : IF (set_half_step_vel_factors) THEN
285 198372 : decay = EXP(-0.5_dp*al%dt*al%nvt(i)%chi)
286 198372 : map_info%v_scale(imap) = decay
287 : END IF
288 : ELSE
289 0 : al%nvt(i)%chi = 0.0_dp
290 0 : IF (set_half_step_vel_factors) THEN
291 0 : map_info%v_scale(imap) = 1.0_dp
292 : END IF
293 : END IF
294 : END DO
295 :
296 32 : END SUBROUTINE al_NH_quarter_step
297 :
298 : END MODULE al_system_dynamics
|
