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 : !> \brief Collection of utilities for setting-up and handle velocities in MD
10 : !> runs
11 : !> \author CJM
12 : !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008
13 : !> reorganization of the original routines/modules
14 : !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
15 : !> (patch by Marcel Baer)
16 : ! **************************************************************************************************
17 : MODULE md_vel_utils
18 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind,&
21 : get_atomic_kind_set
22 : USE bibliography, ONLY: West2006,&
23 : cite_reference
24 : USE cell_types, ONLY: &
25 : cell_type, use_perd_none, use_perd_x, use_perd_xy, use_perd_xyz, use_perd_xz, use_perd_y, &
26 : use_perd_yz, use_perd_z
27 : USE cp_linked_list_input, ONLY: cp_sll_val_next,&
28 : cp_sll_val_type
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
32 : cp_print_key_unit_nr
33 : USE cp_subsys_types, ONLY: cp_subsys_get,&
34 : cp_subsys_type
35 : USE cp_units, ONLY: cp_unit_from_cp2k
36 : USE extended_system_types, ONLY: npt_info_type
37 : USE force_env_methods, ONLY: force_env_calc_energy_force
38 : USE force_env_types, ONLY: force_env_get,&
39 : force_env_type
40 : USE force_env_utils, ONLY: force_env_rattle,&
41 : force_env_shake
42 : USE global_types, ONLY: global_environment_type
43 : USE input_constants, ONLY: &
44 : md_init_default, md_init_vib, npe_f_ensemble, npe_i_ensemble, &
45 : nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
46 : npt_ia_ensemble, reftraj_ensemble
47 : USE input_cp2k_binary_restarts, ONLY: read_binary_velocities
48 : USE input_restart_force_eval, ONLY: update_subsys
49 : USE input_section_types, ONLY: section_vals_get,&
50 : section_vals_get_subs_vals,&
51 : section_vals_list_get,&
52 : section_vals_type,&
53 : section_vals_val_get
54 : USE input_val_types, ONLY: val_get,&
55 : val_type
56 : USE kinds, ONLY: default_string_length,&
57 : dp
58 : USE mathconstants, ONLY: pi
59 : USE mathlib, ONLY: diamat_all
60 : USE md_ener_types, ONLY: md_ener_type
61 : USE md_environment_types, ONLY: get_md_env,&
62 : md_environment_type
63 : USE md_util, ONLY: read_vib_eigs_unformatted
64 : USE message_passing, ONLY: mp_para_env_type
65 : USE molecule_kind_list_types, ONLY: molecule_kind_list_type
66 : USE molecule_kind_types, ONLY: fixd_constraint_type,&
67 : get_molecule_kind,&
68 : get_molecule_kind_set,&
69 : molecule_kind_type
70 : USE parallel_rng_types, ONLY: UNIFORM,&
71 : rng_stream_type
72 : USE particle_list_types, ONLY: particle_list_type
73 : USE particle_types, ONLY: particle_type
74 : USE physcon, ONLY: kelvin
75 : USE shell_opt, ONLY: optimize_shell_core
76 : USE shell_potential_types, ONLY: shell_kind_type
77 : USE simpar_types, ONLY: simpar_type
78 : USE thermal_region_types, ONLY: thermal_region_type,&
79 : thermal_regions_type
80 : #include "../base/base_uses.f90"
81 :
82 : IMPLICIT NONE
83 :
84 : PRIVATE
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_vel_utils'
86 :
87 : PUBLIC :: temperature_control, &
88 : comvel_control, &
89 : angvel_control, &
90 : setup_velocities
91 :
92 : CONTAINS
93 :
94 : ! **************************************************************************************************
95 : !> \brief compute center of mass position
96 : !> *** is only used by initialize_velocities below ***
97 : !> \param part ...
98 : !> \param is_fixed ...
99 : !> \param rcom ...
100 : !> \par History
101 : !> 2007-11-6: created
102 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
103 : ! **************************************************************************************************
104 105 : SUBROUTINE compute_rcom(part, is_fixed, rcom)
105 : TYPE(particle_type), DIMENSION(:), POINTER :: part
106 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
107 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: rcom
108 :
109 : INTEGER :: i
110 : REAL(KIND=dp) :: denom, mass
111 : TYPE(atomic_kind_type), POINTER :: atomic_kind
112 :
113 105 : rcom(:) = 0.0_dp
114 105 : denom = 0.0_dp
115 1242 : DO i = 1, SIZE(part)
116 1137 : atomic_kind => part(i)%atomic_kind
117 1137 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
118 1242 : SELECT CASE (is_fixed(i))
119 : CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
120 1137 : rcom(1) = rcom(1) + part(i)%r(1)*mass
121 1137 : rcom(2) = rcom(2) + part(i)%r(2)*mass
122 1137 : rcom(3) = rcom(3) + part(i)%r(3)*mass
123 1137 : denom = denom + mass
124 : END SELECT
125 : END DO
126 420 : rcom = rcom/denom
127 :
128 105 : END SUBROUTINE compute_rcom
129 :
130 : ! **************************************************************************************************
131 : !> \brief compute center of mass velocity
132 : !> *** is only used by initialize_velocities below ***
133 : !> \param part ...
134 : !> \param is_fixed ...
135 : !> \param vcom ...
136 : !> \param ecom ...
137 : !> \par History
138 : !> 2007-11-6: created
139 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
140 : ! **************************************************************************************************
141 3116 : SUBROUTINE compute_vcom(part, is_fixed, vcom, ecom)
142 : TYPE(particle_type), DIMENSION(:), POINTER :: part
143 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
144 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: vcom
145 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: ecom
146 :
147 : INTEGER :: i
148 : REAL(KIND=dp) :: denom, mass
149 : TYPE(atomic_kind_type), POINTER :: atomic_kind
150 :
151 3116 : vcom = 0.0_dp
152 3116 : denom = 0.0_dp
153 739597 : DO i = 1, SIZE(part)
154 736481 : atomic_kind => part(i)%atomic_kind
155 736481 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
156 739597 : IF (mass /= 0.0) THEN
157 1465594 : SELECT CASE (is_fixed(i))
158 : CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
159 730439 : vcom(1) = vcom(1) + part(i)%v(1)*mass
160 730439 : vcom(2) = vcom(2) + part(i)%v(2)*mass
161 730439 : vcom(3) = vcom(3) + part(i)%v(3)*mass
162 735155 : denom = denom + mass
163 : END SELECT
164 : END IF
165 : END DO
166 12464 : vcom = vcom/denom
167 3116 : IF (PRESENT(ecom)) THEN
168 4372 : ecom = 0.5_dp*denom*SUM(vcom*vcom)
169 : END IF
170 :
171 3116 : END SUBROUTINE compute_vcom
172 :
173 : ! **************************************************************************************************
174 : !> \brief Copy atom velocities into core and shell velocities
175 : !> *** is only used by initialize_velocities below ***
176 : !> \param part ...
177 : !> \param shell_part ...
178 : !> \param core_part ...
179 : !> \par History
180 : !> 2007-11-6: created
181 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
182 : ! **************************************************************************************************
183 8 : SUBROUTINE clone_core_shell_vel(part, shell_part, core_part)
184 : TYPE(particle_type), DIMENSION(:), POINTER :: part, shell_part, core_part
185 :
186 : INTEGER :: i
187 : LOGICAL :: is_shell
188 : TYPE(atomic_kind_type), POINTER :: atomic_kind
189 :
190 776 : DO i = 1, SIZE(part)
191 768 : atomic_kind => part(i)%atomic_kind
192 768 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_shell)
193 776 : IF (is_shell) THEN
194 6144 : shell_part(part(i)%shell_index)%v(:) = part(i)%v(:)
195 6144 : core_part(part(i)%shell_index)%v(:) = part(i)%v(:)
196 : END IF
197 : END DO
198 :
199 8 : END SUBROUTINE clone_core_shell_vel
200 :
201 : ! **************************************************************************************************
202 : !> \brief Compute the kinetic energy. Does not subtract the center of mass kinetic
203 : !> energy.
204 : !> *** is only used by initialize_velocities below ***
205 : !> \param part ...
206 : !> \param ireg ...
207 : !> \return ...
208 : !> \par History
209 : !> 2007-11-6: created
210 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
211 : ! **************************************************************************************************
212 3382 : FUNCTION compute_ekin(part, ireg) RESULT(ekin)
213 : TYPE(particle_type), DIMENSION(:), POINTER :: part
214 : INTEGER, INTENT(IN), OPTIONAL :: ireg
215 : REAL(KIND=dp) :: ekin
216 :
217 : INTEGER :: i
218 : REAL(KIND=dp) :: mass
219 : TYPE(atomic_kind_type), POINTER :: atomic_kind
220 :
221 3382 : NULLIFY (atomic_kind)
222 3382 : ekin = 0.0_dp
223 3382 : IF (PRESENT(ireg)) THEN
224 13756 : DO i = 1, SIZE(part)
225 13756 : IF (part(i)%t_region_index == ireg) THEN
226 4560 : atomic_kind => part(i)%atomic_kind
227 4560 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
228 18240 : ekin = ekin + 0.5_dp*mass*SUM(part(i)%v(:)*part(i)%v(:))
229 : END IF
230 : END DO
231 : ELSE
232 739403 : DO i = 1, SIZE(part)
233 736289 : atomic_kind => part(i)%atomic_kind
234 736289 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
235 2948270 : ekin = ekin + 0.5_dp*mass*SUM(part(i)%v(:)*part(i)%v(:))
236 : END DO
237 : END IF
238 :
239 3382 : END FUNCTION compute_ekin
240 :
241 : ! **************************************************************************************************
242 : !> \brief Rescale the velocity to mimic the given external kinetic temperature.
243 : !> Optionally also rescale vcom.
244 : !> *** is only used by initialize_velocities below ***
245 : !> \param part ...
246 : !> \param simpar ...
247 : !> \param ekin ...
248 : !> \param vcom ...
249 : !> \param ireg ...
250 : !> \param nfree ...
251 : !> \param temp ...
252 : !> \par History
253 : !> 2007-11-6: created
254 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
255 : ! **************************************************************************************************
256 2051 : SUBROUTINE rescale_vel(part, simpar, ekin, vcom, ireg, nfree, temp)
257 : TYPE(particle_type), DIMENSION(:), POINTER :: part
258 : TYPE(simpar_type), POINTER :: simpar
259 : REAL(KIND=dp), INTENT(INOUT) :: ekin
260 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
261 : OPTIONAL :: vcom
262 : INTEGER, INTENT(IN), OPTIONAL :: ireg, nfree
263 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: temp
264 :
265 : INTEGER :: i, my_ireg, my_nfree
266 : REAL(KIND=dp) :: factor, my_temp
267 :
268 2051 : IF (PRESENT(ireg) .AND. PRESENT(nfree) .AND. PRESENT(temp)) THEN
269 6 : my_ireg = ireg
270 6 : my_nfree = nfree
271 6 : my_temp = temp
272 2045 : ELSEIF (PRESENT(nfree)) THEN
273 0 : my_ireg = 0
274 0 : my_nfree = nfree
275 0 : my_temp = simpar%temp_ext
276 : ELSE
277 2045 : my_ireg = 0
278 2045 : my_nfree = simpar%nfree
279 2045 : my_temp = simpar%temp_ext
280 : END IF
281 2051 : IF (my_nfree /= 0) THEN
282 2039 : factor = my_temp/(2.0_dp*ekin)*REAL(my_nfree, KIND=dp)
283 : ELSE
284 : factor = 0.0_dp
285 : END IF
286 : ! Note:
287 : ! this rescaling is still wrong, it should take the masses into account
288 : ! rescaling is generally not correct, so needs fixing
289 2051 : ekin = ekin*factor
290 2051 : factor = SQRT(factor)
291 2051 : IF (PRESENT(ireg)) THEN
292 582 : DO i = 1, SIZE(part)
293 1158 : IF (part(i)%t_region_index == my_ireg) part(i)%v(:) = factor*part(i)%v(:)
294 : END DO
295 : ELSE
296 439369 : DO i = 1, SIZE(part)
297 1751341 : part(i)%v(:) = factor*part(i)%v(:)
298 : END DO
299 2045 : IF (PRESENT(vcom)) THEN
300 96 : vcom = factor*vcom
301 : END IF
302 : END IF
303 :
304 2051 : END SUBROUTINE rescale_vel
305 :
306 : ! **************************************************************************************************
307 : !> \brief Rescale the velocity of separated regions independently
308 : !> \param part ...
309 : !> \param md_env ...
310 : !> \param simpar ...
311 : !> \par History
312 : !> 2008-11
313 : !> \author MI
314 : ! **************************************************************************************************
315 2 : SUBROUTINE rescale_vel_region(part, md_env, simpar)
316 :
317 : TYPE(particle_type), DIMENSION(:), POINTER :: part
318 : TYPE(md_environment_type), POINTER :: md_env
319 : TYPE(simpar_type), POINTER :: simpar
320 :
321 : INTEGER :: ireg, nfree, nfree0, nfree_done
322 : REAL(KIND=dp) :: ekin, temp
323 : TYPE(thermal_region_type), POINTER :: t_region
324 : TYPE(thermal_regions_type), POINTER :: thermal_regions
325 :
326 2 : NULLIFY (thermal_regions, t_region)
327 :
328 2 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
329 2 : nfree_done = 0
330 6 : DO ireg = 1, thermal_regions%nregions
331 4 : NULLIFY (t_region)
332 4 : t_region => thermal_regions%thermal_region(ireg)
333 4 : nfree = t_region%npart*3
334 4 : ekin = compute_ekin(part, ireg)
335 4 : temp = t_region%temp_expected
336 4 : CALL rescale_vel(part, simpar, ekin, ireg=ireg, nfree=nfree, temp=temp)
337 4 : nfree_done = nfree_done + nfree
338 4 : ekin = compute_ekin(part, ireg)
339 4 : temp = 2.0_dp*ekin/REAL(nfree, dp)*kelvin
340 6 : t_region%temperature = temp
341 : END DO
342 2 : nfree0 = simpar%nfree - nfree_done
343 2 : IF (nfree0 > 0) THEN
344 2 : ekin = compute_ekin(part, 0)
345 2 : CALL rescale_vel(part, simpar, ekin, ireg=0, nfree=nfree0, temp=simpar%temp_ext)
346 2 : ekin = compute_ekin(part, 0)
347 2 : temp = 2.0_dp*ekin/REAL(nfree0, dp)*kelvin
348 2 : thermal_regions%temp_reg0 = temp
349 : END IF
350 2 : END SUBROUTINE rescale_vel_region
351 :
352 : ! **************************************************************************************************
353 : !> \brief subtract center of mass velocity
354 : !> *** is only used by initialize_velocities below ***
355 : !> \param part ...
356 : !> \param is_fixed ...
357 : !> \param vcom ...
358 : !> \par History
359 : !> 2007-11-6: created
360 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
361 : ! **************************************************************************************************
362 2023 : SUBROUTINE subtract_vcom(part, is_fixed, vcom)
363 : TYPE(particle_type), DIMENSION(:), POINTER :: part
364 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
365 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vcom
366 :
367 : INTEGER :: i
368 :
369 437669 : DO i = 1, SIZE(part)
370 2023 : SELECT CASE (is_fixed(i))
371 : CASE (use_perd_x)
372 512 : part(i)%v(2) = part(i)%v(2) - vcom(2)
373 512 : part(i)%v(3) = part(i)%v(3) - vcom(3)
374 : CASE (use_perd_y)
375 512 : part(i)%v(1) = part(i)%v(1) - vcom(1)
376 512 : part(i)%v(3) = part(i)%v(3) - vcom(3)
377 : CASE (use_perd_z)
378 512 : part(i)%v(1) = part(i)%v(1) - vcom(1)
379 512 : part(i)%v(2) = part(i)%v(2) - vcom(2)
380 : CASE (use_perd_xy)
381 512 : part(i)%v(3) = part(i)%v(3) - vcom(3)
382 : CASE (use_perd_xz)
383 0 : part(i)%v(2) = part(i)%v(2) - vcom(2)
384 : CASE (use_perd_yz)
385 0 : part(i)%v(1) = part(i)%v(1) - vcom(1)
386 : CASE (use_perd_none)
387 1727086 : part(i)%v(:) = part(i)%v(:) - vcom(:)
388 : END SELECT
389 : END DO
390 2023 : END SUBROUTINE subtract_vcom
391 :
392 : ! **************************************************************************************************
393 : !> \brief compute the angular velocity
394 : !> *** is only used by initialize_velocities below ***
395 : !> \param part ...
396 : !> \param is_fixed ...
397 : !> \param rcom ...
398 : !> \param vang ...
399 : !> \par History
400 : !> 2007-11-9: created
401 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
402 : ! **************************************************************************************************
403 107 : SUBROUTINE compute_vang(part, is_fixed, rcom, vang)
404 : TYPE(particle_type), DIMENSION(:), POINTER :: part
405 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
406 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rcom
407 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: vang
408 :
409 : INTEGER :: i
410 : REAL(KIND=dp) :: mass, proj
411 : REAL(KIND=dp), DIMENSION(3) :: evals, mang, r
412 : REAL(KIND=dp), DIMENSION(3, 3) :: iner
413 : TYPE(atomic_kind_type), POINTER :: atomic_kind
414 :
415 107 : NULLIFY (atomic_kind)
416 107 : mang(:) = 0.0_dp
417 107 : iner(:, :) = 0.0_dp
418 1272 : DO i = 1, SIZE(part)
419 : ! compute angular momentum and inertia tensor
420 107 : SELECT CASE (is_fixed(i))
421 : CASE (use_perd_x, use_perd_y, use_perd_z, use_perd_xy, use_perd_xz, use_perd_yz, use_perd_none)
422 4660 : r(:) = part(i)%r(:) - rcom(:)
423 1165 : atomic_kind => part(i)%atomic_kind
424 1165 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
425 1165 : mang(1) = mang(1) + mass*(r(2)*part(i)%v(3) - r(3)*part(i)%v(2))
426 1165 : mang(2) = mang(2) + mass*(r(3)*part(i)%v(1) - r(1)*part(i)%v(3))
427 1165 : mang(3) = mang(3) + mass*(r(1)*part(i)%v(2) - r(2)*part(i)%v(1))
428 :
429 1165 : iner(1, 1) = iner(1, 1) + mass*(r(2)*r(2) + r(3)*r(3))
430 1165 : iner(2, 2) = iner(2, 2) + mass*(r(3)*r(3) + r(1)*r(1))
431 1165 : iner(3, 3) = iner(3, 3) + mass*(r(1)*r(1) + r(2)*r(2))
432 :
433 1165 : iner(1, 2) = iner(1, 2) - mass*r(1)*r(2)
434 1165 : iner(2, 3) = iner(2, 3) - mass*r(2)*r(3)
435 2330 : iner(3, 1) = iner(3, 1) - mass*r(3)*r(1)
436 : END SELECT
437 : END DO
438 107 : iner(2, 1) = iner(1, 2)
439 107 : iner(3, 2) = iner(2, 3)
440 107 : iner(1, 3) = iner(3, 1)
441 :
442 : ! Take the safest route, i.e. diagonalize the inertia tensor and solve
443 : ! the angular velocity only with the non-zero eigenvalues. A plain inversion
444 : ! would fail for linear molecules.
445 107 : CALL diamat_all(iner, evals)
446 :
447 107 : vang(:) = 0.0_dp
448 428 : DO i = 1, 3
449 428 : IF (evals(i) > 0.0_dp) THEN
450 1240 : proj = SUM(iner(:, i)*mang)/evals(i)
451 310 : vang(1) = vang(1) + proj*iner(1, i)
452 310 : vang(2) = vang(2) + proj*iner(2, i)
453 310 : vang(3) = vang(3) + proj*iner(3, i)
454 : END IF
455 : END DO
456 :
457 107 : END SUBROUTINE compute_vang
458 :
459 : ! **************************************************************************************************
460 : !> \brief subtract the angular velocity
461 : !> *** is only used by initialize_velocities below ***
462 : !> \param part ...
463 : !> \param is_fixed ...
464 : !> \param rcom ...
465 : !> \param vang ...
466 : !> \par History
467 : !> 2007-11-9: created
468 : !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
469 : ! **************************************************************************************************
470 6 : SUBROUTINE subtract_vang(part, is_fixed, rcom, vang)
471 : TYPE(particle_type), DIMENSION(:), POINTER :: part
472 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
473 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rcom, vang
474 :
475 : INTEGER :: i
476 : REAL(KIND=dp), DIMENSION(3) :: r
477 :
478 52 : DO i = 1, SIZE(part)
479 184 : r(:) = part(i)%r(:) - rcom(:)
480 6 : SELECT CASE (is_fixed(i))
481 : CASE (use_perd_x)
482 0 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
483 0 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
484 : CASE (use_perd_y)
485 0 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
486 0 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
487 : CASE (use_perd_z)
488 0 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
489 0 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
490 : CASE (use_perd_xy)
491 0 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
492 : CASE (use_perd_xz)
493 0 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
494 : CASE (use_perd_yz)
495 0 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
496 : CASE (use_perd_none)
497 46 : part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
498 46 : part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
499 46 : part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
500 : END SELECT
501 : END DO
502 :
503 6 : END SUBROUTINE subtract_vang
504 :
505 : ! **************************************************************************************************
506 : !> \brief Initializes the velocities to the Maxwellian distribution
507 : !> \param simpar ...
508 : !> \param part ...
509 : !> \param force_env ...
510 : !> \param globenv ...
511 : !> \param md_env ...
512 : !> \param molecule_kinds ...
513 : !> \param label ...
514 : !> \param print_section ...
515 : !> \param subsys_section ...
516 : !> \param shell_present ...
517 : !> \param shell_part ...
518 : !> \param core_part ...
519 : !> \param force_rescaling ...
520 : !> \param para_env ...
521 : !> \param write_binary_restart_file ...
522 : !> \par History
523 : !> - is_fixed removed from particle_type
524 : !> - 2007-11-07: Cleanup (TV)
525 : !> - 2007-11-09: Added angvel_zero feature
526 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>
527 : ! **************************************************************************************************
528 1791 : SUBROUTINE initialize_velocities(simpar, &
529 : part, &
530 : force_env, &
531 : globenv, &
532 : md_env, &
533 : molecule_kinds, &
534 : label, &
535 : print_section, &
536 : subsys_section, &
537 : shell_present, &
538 : shell_part, &
539 : core_part, &
540 : force_rescaling, &
541 : para_env, &
542 : write_binary_restart_file)
543 :
544 : TYPE(simpar_type), POINTER :: simpar
545 : TYPE(particle_type), DIMENSION(:), POINTER :: part
546 : TYPE(force_env_type), POINTER :: force_env
547 : TYPE(global_environment_type), POINTER :: globenv
548 : TYPE(md_environment_type), POINTER :: md_env
549 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
550 : CHARACTER(LEN=*), INTENT(IN) :: label
551 : TYPE(section_vals_type), POINTER :: print_section, subsys_section
552 : LOGICAL, INTENT(IN) :: shell_present
553 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_part, core_part
554 : LOGICAL, INTENT(IN) :: force_rescaling
555 : TYPE(mp_para_env_type), POINTER :: para_env
556 : LOGICAL, INTENT(IN) :: write_binary_restart_file
557 :
558 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_velocities'
559 :
560 : INTEGER :: handle, i, ifixd, imolecule_kind, iw, &
561 : natoms
562 : INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
563 : LOGICAL :: success
564 : REAL(KIND=dp) :: ecom, ekin, mass, mass_tot, temp, tmp_r1
565 : REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vcom
566 : TYPE(atomic_kind_type), POINTER :: atomic_kind
567 : TYPE(cell_type), POINTER :: cell
568 : TYPE(cp_logger_type), POINTER :: logger
569 1791 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
570 1791 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
571 : TYPE(molecule_kind_type), POINTER :: molecule_kind
572 : TYPE(section_vals_type), POINTER :: md_section, root_section, vib_section
573 :
574 1791 : CALL timeset(routineN, handle)
575 :
576 : ! Initializing parameters
577 1791 : natoms = SIZE(part)
578 1791 : NULLIFY (atomic_kind, fixd_list, logger, molecule_kind)
579 1791 : NULLIFY (molecule_kind_set)
580 :
581 : ! Logging
582 1791 : logger => cp_get_default_logger()
583 1791 : iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
584 :
585 : ! Build a list of all fixed atoms (if any)
586 5373 : ALLOCATE (is_fixed(natoms))
587 :
588 489589 : is_fixed = use_perd_none
589 1791 : molecule_kind_set => molecule_kinds%els
590 62005 : DO imolecule_kind = 1, molecule_kinds%n_els
591 60214 : molecule_kind => molecule_kind_set(imolecule_kind)
592 60214 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
593 62005 : IF (ASSOCIATED(fixd_list)) THEN
594 5432 : DO ifixd = 1, SIZE(fixd_list)
595 5432 : IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
596 : END DO
597 : END IF
598 : END DO
599 :
600 : ! Compute the total mass when needed
601 1791 : IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
602 : simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
603 : mass_tot = 0.0_dp
604 1006 : DO i = 1, natoms
605 1000 : atomic_kind => part(i)%atomic_kind
606 1000 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
607 1006 : mass_tot = mass_tot + mass
608 : END DO
609 6 : simpar%v_shock = simpar%v_shock*SQRT(mass_tot)
610 : END IF
611 :
612 : CALL read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
613 1791 : shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
614 1791 : IF (.NOT. success) THEN
615 3048 : SELECT CASE (simpar%initialization_method)
616 : CASE (md_init_default)
617 : CALL generate_velocities(simpar, part, force_env, globenv, md_env, shell_present, &
618 1523 : shell_part, core_part, is_fixed, iw)
619 : CASE (md_init_vib)
620 2 : CALL force_env_get(force_env=force_env, root_section=root_section)
621 2 : md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
622 2 : vib_section => section_vals_get_subs_vals(root_section, "VIBRATIONAL_ANALYSIS")
623 : CALL generate_coords_vels_vib(simpar, &
624 : part, &
625 : md_section, &
626 : vib_section, &
627 : force_env, &
628 : globenv, &
629 : shell_present, &
630 : shell_part, &
631 : core_part, &
632 2 : is_fixed)
633 : ! update restart file for the modified coordinates and velocities
634 1527 : CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
635 : END SELECT
636 : END IF
637 :
638 1791 : IF (iw > 0) THEN
639 : WRITE (iw, '(/,T2,A)') &
640 826 : 'MD_VEL| '//TRIM(ADJUSTL(label))
641 : ! Recompute vcom, ecom and ekin for IO
642 826 : CALL compute_vcom(part, is_fixed, vcom, ecom)
643 826 : ekin = compute_ekin(part) - ecom
644 826 : IF (simpar%nfree == 0) THEN
645 6 : CPASSERT(ekin == 0.0_dp)
646 6 : temp = 0.0_dp
647 : ELSE
648 820 : temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
649 : END IF
650 826 : tmp_r1 = cp_unit_from_cp2k(temp, "K")
651 : WRITE (iw, '(T2,A,T61,F20.6)') &
652 826 : 'MD_VEL| Initial temperature [K]', tmp_r1
653 : WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
654 826 : 'MD_VEL| COM velocity', vcom(1:3)
655 :
656 : ! compute and log rcom and vang if not periodic
657 826 : CALL force_env_get(force_env, cell=cell)
658 3304 : IF (SUM(cell%perd(1:3)) == 0) THEN
659 61 : CALL compute_rcom(part, is_fixed, rcom)
660 61 : CALL compute_vang(part, is_fixed, rcom, vang)
661 : WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
662 61 : 'MD_VEL| COM position', rcom(1:3)
663 : WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
664 61 : 'MD_VEL| Angular velocity', vang(1:3)
665 : END IF
666 : END IF
667 :
668 1791 : DEALLOCATE (is_fixed)
669 1791 : CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
670 1791 : CALL timestop(handle)
671 :
672 3582 : END SUBROUTINE initialize_velocities
673 :
674 : ! **************************************************************************************************
675 : !> \brief Read velocities from binary restart file if available
676 : !> \param simpar ...
677 : !> \param part ...
678 : !> \param force_env ...
679 : !> \param md_env ...
680 : !> \param subsys_section ...
681 : !> \param shell_present ...
682 : !> \param shell_part ...
683 : !> \param core_part ...
684 : !> \param force_rescaling ...
685 : !> \param para_env ...
686 : !> \param is_fixed ...
687 : !> \param success ...
688 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>
689 : ! **************************************************************************************************
690 12537 : SUBROUTINE read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
691 1791 : shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
692 : TYPE(simpar_type), POINTER :: simpar
693 : TYPE(particle_type), DIMENSION(:), POINTER :: part
694 : TYPE(force_env_type), POINTER :: force_env
695 : TYPE(md_environment_type), POINTER :: md_env
696 : TYPE(section_vals_type), POINTER :: subsys_section
697 : LOGICAL, INTENT(IN) :: shell_present
698 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_part, core_part
699 : LOGICAL, INTENT(IN) :: force_rescaling
700 : TYPE(mp_para_env_type), POINTER :: para_env
701 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
702 : LOGICAL, INTENT(OUT) :: success
703 :
704 : INTEGER :: i, natoms, nshell, shell_index
705 : LOGICAL :: atomvel_explicit, atomvel_read, corevel_explicit, corevel_read, is_ok, &
706 : rescale_regions, shellvel_explicit, shellvel_read
707 : REAL(KIND=dp) :: ecom, ekin, fac_massc, fac_masss, mass
708 : REAL(KIND=dp), DIMENSION(3) :: vc, vcom, vs
709 1791 : REAL(KIND=dp), DIMENSION(:), POINTER :: vel
710 : TYPE(atomic_kind_type), POINTER :: atomic_kind
711 : TYPE(cp_sll_val_type), POINTER :: atom_list, core_list, shell_list
712 : TYPE(section_vals_type), POINTER :: atomvel_section, corevel_section, &
713 : shellvel_section
714 : TYPE(shell_kind_type), POINTER :: shell
715 : TYPE(thermal_regions_type), POINTER :: thermal_regions
716 : TYPE(val_type), POINTER :: val
717 :
718 : ! Initializing parameters
719 :
720 1791 : success = .FALSE.
721 1791 : natoms = SIZE(part)
722 1791 : atomvel_read = .FALSE.
723 1791 : corevel_read = .FALSE.
724 1791 : shellvel_read = .FALSE.
725 1791 : NULLIFY (vel, atomic_kind, atom_list, core_list, shell_list)
726 1791 : NULLIFY (atomvel_section, shellvel_section, corevel_section)
727 1791 : NULLIFY (shell, thermal_regions, val)
728 :
729 : ! Core-Shell Model
730 1791 : nshell = 0
731 1791 : IF (shell_present) THEN
732 132 : CPASSERT(ASSOCIATED(core_part))
733 132 : CPASSERT(ASSOCIATED(shell_part))
734 132 : nshell = SIZE(shell_part)
735 : END IF
736 :
737 1791 : atomvel_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
738 1791 : shellvel_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
739 1791 : corevel_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
740 :
741 : ! Read or initialize the particle velocities
742 1791 : CALL section_vals_get(atomvel_section, explicit=atomvel_explicit)
743 1791 : CALL section_vals_get(shellvel_section, explicit=shellvel_explicit)
744 1791 : CALL section_vals_get(corevel_section, explicit=corevel_explicit)
745 1791 : CPASSERT(shellvel_explicit .EQV. corevel_explicit)
746 :
747 : CALL read_binary_velocities("", part, force_env%root_section, para_env, &
748 1791 : subsys_section, atomvel_read)
749 : CALL read_binary_velocities("SHELL", shell_part, force_env%root_section, para_env, &
750 1791 : subsys_section, shellvel_read)
751 : CALL read_binary_velocities("CORE", core_part, force_env%root_section, para_env, &
752 1791 : subsys_section, corevel_read)
753 :
754 1791 : IF (.NOT. (atomvel_explicit .OR. atomvel_read)) RETURN
755 266 : success = .TRUE.
756 :
757 266 : IF (.NOT. atomvel_read) THEN
758 : ! Read the atom velocities if explicitly given in the input file
759 220 : CALL section_vals_list_get(atomvel_section, "_DEFAULT_KEYWORD_", list=atom_list)
760 52950 : DO i = 1, natoms
761 52730 : is_ok = cp_sll_val_next(atom_list, val)
762 52730 : CALL val_get(val, r_vals=vel)
763 369330 : part(i)%v = vel
764 : END DO
765 : END IF
766 57412 : DO i = 1, natoms
767 266 : SELECT CASE (is_fixed(i))
768 : CASE (use_perd_x)
769 0 : part(i)%v(1) = 0.0_dp
770 : CASE (use_perd_y)
771 0 : part(i)%v(2) = 0.0_dp
772 : CASE (use_perd_z)
773 0 : part(i)%v(3) = 0.0_dp
774 : CASE (use_perd_xy)
775 0 : part(i)%v(1) = 0.0_dp
776 0 : part(i)%v(2) = 0.0_dp
777 : CASE (use_perd_xz)
778 0 : part(i)%v(1) = 0.0_dp
779 0 : part(i)%v(3) = 0.0_dp
780 : CASE (use_perd_yz)
781 0 : part(i)%v(2) = 0.0_dp
782 0 : part(i)%v(3) = 0.0_dp
783 : CASE (use_perd_xyz)
784 57224 : part(i)%v = 0.0_dp
785 : END SELECT
786 : END DO
787 266 : IF (shell_present) THEN
788 48 : IF (shellvel_explicit) THEN
789 : ! If the atoms positions are given (?) and core and shell velocities are
790 : ! present in the input, read the latter.
791 24 : CALL section_vals_list_get(shellvel_section, "_DEFAULT_KEYWORD_", list=shell_list)
792 24 : CALL section_vals_list_get(corevel_section, "_DEFAULT_KEYWORD_", list=core_list)
793 2328 : DO i = 1, nshell
794 2304 : is_ok = cp_sll_val_next(shell_list, val)
795 2304 : CALL val_get(val, r_vals=vel)
796 16128 : shell_part(i)%v = vel
797 2304 : is_ok = cp_sll_val_next(core_list, val)
798 2304 : CALL val_get(val, r_vals=vel)
799 16152 : core_part(i)%v = vel
800 : END DO
801 : ELSE
802 24 : IF (.NOT. (shellvel_read .AND. corevel_read)) THEN
803 : ! Otherwise, just copy atom velocties into shell and core velocities.
804 8 : CALL clone_core_shell_vel(part, shell_part, core_part)
805 : END IF
806 : END IF
807 : END IF
808 :
809 : ! compute vcom, ecom and ekin
810 266 : CALL compute_vcom(part, is_fixed, vcom, ecom)
811 266 : ekin = compute_ekin(part) - ecom
812 :
813 266 : IF (simpar%do_thermal_region) THEN
814 12 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
815 12 : IF (ASSOCIATED(thermal_regions)) THEN
816 12 : rescale_regions = thermal_regions%force_rescaling
817 : END IF
818 : ELSE
819 : rescale_regions = .FALSE.
820 : END IF
821 266 : IF (simpar%nfree /= 0 .AND. (force_rescaling .OR. rescale_regions)) THEN
822 24 : IF (simpar%do_thermal_region) THEN
823 0 : CALL rescale_vel_region(part, md_env, simpar)
824 : ELSE
825 24 : CALL rescale_vel(part, simpar, ekin, vcom=vcom)
826 : END IF
827 :
828 : ! After rescaling, the core and shell velocities must also adapt.
829 1894 : DO i = 1, natoms
830 1870 : shell_index = part(i)%shell_index
831 2136 : IF (shell_present .AND. shell_index /= 0) THEN
832 0 : atomic_kind => part(i)%atomic_kind
833 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell)
834 0 : fac_masss = shell%mass_shell/mass
835 0 : fac_massc = shell%mass_core/mass
836 0 : vs = shell_part(shell_index)%v
837 0 : vc = core_part(shell_index)%v
838 :
839 0 : shell_part(shell_index)%v(1) = part(i)%v(1) + fac_massc*(vs(1) - vc(1))
840 0 : shell_part(shell_index)%v(2) = part(i)%v(2) + fac_massc*(vs(2) - vc(2))
841 0 : shell_part(shell_index)%v(3) = part(i)%v(3) + fac_massc*(vs(3) - vc(3))
842 0 : core_part(shell_index)%v(1) = part(i)%v(1) + fac_masss*(vc(1) - vs(1))
843 0 : core_part(shell_index)%v(2) = part(i)%v(2) + fac_masss*(vc(2) - vs(2))
844 0 : core_part(shell_index)%v(3) = part(i)%v(3) + fac_masss*(vc(3) - vs(3))
845 : END IF
846 : END DO
847 : END IF
848 1791 : END SUBROUTINE read_input_velocities
849 :
850 : ! **************************************************************************************************
851 : !> \brief Initializing velocities AND positions randomly on all processors, based on vibrational
852 : !> modes of the system, so that the starting coordinates are already very close to
853 : !> canonical ensumble corresponding to temperature of a head bath.
854 : !> \param simpar : MD simulation parameters
855 : !> \param particles : global array of all particles
856 : !> \param md_section : MD input subsection
857 : !> \param vib_section : vibrational analysis input section
858 : !> \param force_env : force environment container
859 : !> \param global_env : global environment container
860 : !> \param shell_present : if core-shell model is used
861 : !> \param shell_particles : global array of all shell particles in shell model
862 : !> \param core_particles : global array of all core particles in shell model
863 : !> \param is_fixed : array of size of total number of atoms, that determines if any
864 : !> cartesian components are fixed
865 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
866 : ! **************************************************************************************************
867 2 : SUBROUTINE generate_coords_vels_vib(simpar, &
868 : particles, &
869 : md_section, &
870 : vib_section, &
871 : force_env, &
872 : global_env, &
873 : shell_present, &
874 : shell_particles, &
875 : core_particles, &
876 2 : is_fixed)
877 : TYPE(simpar_type), POINTER :: simpar
878 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
879 : TYPE(section_vals_type), POINTER :: md_section, vib_section
880 : TYPE(force_env_type), POINTER :: force_env
881 : TYPE(global_environment_type), POINTER :: global_env
882 : LOGICAL, INTENT(IN) :: shell_present
883 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_particles, core_particles
884 : INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
885 :
886 : INTEGER :: dof, fixed_dof, iatom, ii, imode, &
887 : my_dof, natoms, shell_index
888 : REAL(KIND=dp) :: Erand, mass, my_phase, ratio, temperature
889 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, phase, random
890 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dr, eigenvectors
891 : TYPE(atomic_kind_type), POINTER :: atomic_kind
892 : TYPE(mp_para_env_type), POINTER :: para_env
893 2 : TYPE(rng_stream_type), ALLOCATABLE :: random_stream
894 :
895 2 : CALL cite_reference(West2006)
896 2 : natoms = SIZE(particles)
897 2 : temperature = simpar%temp_ext
898 2 : my_dof = 3*natoms
899 6 : ALLOCATE (eigenvalues(my_dof))
900 8 : ALLOCATE (eigenvectors(my_dof, my_dof))
901 4 : ALLOCATE (phase(my_dof))
902 4 : ALLOCATE (random(my_dof))
903 6 : ALLOCATE (dr(3, natoms))
904 2 : CALL force_env_get(force_env=force_env, para_env=para_env)
905 : ! read vibration modes
906 : CALL read_vib_eigs_unformatted(md_section, &
907 : vib_section, &
908 : para_env, &
909 : dof, &
910 : eigenvalues, &
911 2 : eigenvectors)
912 2 : IF (my_dof /= dof) THEN
913 : CALL cp_abort(__LOCATION__, &
914 : "number of degrees of freedom in vibrational analysis data "// &
915 0 : "do not match total number of cartesian degrees of freedom")
916 : END IF
917 : ! read phases
918 2 : CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%PHASE", r_val=my_phase)
919 2 : my_phase = MIN(1.0_dp, my_phase)
920 : ! generate random numbers
921 2 : random_stream = rng_stream_type(name="MD_INIT_VIB", distribution_type=UNIFORM)
922 2 : CALL random_stream%fill(random)
923 2 : IF (my_phase < 0.0_dp) THEN
924 0 : CALL random_stream%fill(phase)
925 : ELSE
926 20 : phase = my_phase
927 : END IF
928 2 : DEALLOCATE (random_stream)
929 :
930 : ! the first three modes are acoustic with zero frequencies,
931 : ! exclude these from considerations
932 2 : my_dof = dof - 3
933 : ! randomly selects energy from distribution about kT, all
934 : ! energies are scaled so that the sum over vibration modes gives
935 : ! exactly my_dof*kT. Note that k = 1.0 in atomic units
936 2 : Erand = 0.0_dp
937 14 : DO imode = 4, dof
938 14 : Erand = Erand - temperature*LOG(1.0_dp - random(imode))
939 : END DO
940 : ! need to take into account of fixed constraints too
941 2 : fixed_dof = 0
942 8 : DO iatom = 1, natoms
943 2 : SELECT CASE (is_fixed(iatom))
944 : CASE (use_perd_x, use_perd_y, use_perd_z)
945 0 : fixed_dof = fixed_dof + 1
946 : CASE (use_perd_xy, use_perd_xz, use_perd_yz)
947 0 : fixed_dof = fixed_dof + 2
948 : CASE (use_perd_xyz)
949 6 : fixed_dof = fixed_dof + 3
950 : END SELECT
951 : END DO
952 2 : my_dof = my_dof - fixed_dof
953 2 : ratio = REAL(my_dof, KIND=dp)*temperature/Erand
954 : ! update velocities AND positions
955 8 : DO iatom = 1, natoms
956 6 : atomic_kind => particles(iatom)%atomic_kind
957 6 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
958 8 : SELECT CASE (is_fixed(iatom))
959 : CASE (use_perd_x)
960 0 : DO ii = 2, 3
961 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
962 0 : eigenvectors, random, phase, dof, ratio)
963 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
964 : eigenvectors, random, phase, dof, &
965 0 : ratio)
966 : END DO
967 : CASE (use_perd_y)
968 0 : DO ii = 1, 3, 2
969 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
970 0 : eigenvectors, random, phase, dof, ratio)
971 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
972 : eigenvectors, random, phase, dof, &
973 0 : ratio)
974 : END DO
975 : CASE (use_perd_z)
976 0 : DO ii = 1, 2
977 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
978 0 : eigenvectors, random, phase, dof, ratio)
979 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
980 : eigenvectors, random, phase, dof, &
981 0 : ratio)
982 : END DO
983 : CASE (use_perd_xy)
984 : dr(3, iatom) = dr_from_vib_data(iatom, 3, mass, temperature, eigenvalues, &
985 0 : eigenvectors, random, phase, dof, ratio)
986 : particles(iatom)%v(3) = dv_from_vib_data(iatom, 3, mass, temperature, &
987 : eigenvectors, random, phase, dof, &
988 0 : ratio)
989 : CASE (use_perd_xz)
990 : dr(2, iatom) = dr_from_vib_data(iatom, 2, mass, temperature, eigenvalues, &
991 0 : eigenvectors, random, phase, dof, ratio)
992 : particles(iatom)%v(2) = dv_from_vib_data(iatom, 2, mass, temperature, &
993 : eigenvectors, random, phase, dof, &
994 0 : ratio)
995 : CASE (use_perd_yz)
996 : dr(1, iatom) = dr_from_vib_data(iatom, 1, mass, temperature, eigenvalues, &
997 0 : eigenvectors, random, phase, dof, ratio)
998 : particles(iatom)%v(1) = dv_from_vib_data(iatom, 1, mass, temperature, &
999 : eigenvectors, random, phase, dof, &
1000 0 : ratio)
1001 : CASE (use_perd_none)
1002 24 : DO ii = 1, 3
1003 : dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
1004 18 : eigenvectors, random, phase, dof, ratio)
1005 : particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
1006 : eigenvectors, random, phase, dof, &
1007 24 : ratio)
1008 : END DO
1009 : END SELECT
1010 : END DO ! iatom
1011 : ! free memory
1012 2 : DEALLOCATE (eigenvalues)
1013 2 : DEALLOCATE (eigenvectors)
1014 2 : DEALLOCATE (phase)
1015 2 : DEALLOCATE (random)
1016 : ! update particle coordinates
1017 8 : DO iatom = 1, natoms
1018 26 : particles(iatom)%r(:) = particles(iatom)%r(:) + dr(:, iatom)
1019 : END DO
1020 : ! update core-shell model coordinates
1021 2 : IF (shell_present) THEN
1022 : ! particles have moved, and for core-shell model this means
1023 : ! the cores and shells must also move by the same amount. The
1024 : ! shell positions will then be optimised if needed
1025 0 : shell_index = particles(iatom)%shell_index
1026 0 : IF (shell_index /= 0) THEN
1027 : core_particles(shell_index)%r(:) = core_particles(shell_index)%r(:) + &
1028 0 : dr(:, iatom)
1029 : shell_particles(shell_index)%r(:) = shell_particles(shell_index)%r(:) + &
1030 0 : dr(:, iatom)
1031 : END IF
1032 : CALL optimize_shell_core(force_env, &
1033 : particles, &
1034 : shell_particles, &
1035 : core_particles, &
1036 0 : global_env)
1037 : END IF
1038 : ! cleanup
1039 2 : DEALLOCATE (dr)
1040 4 : END SUBROUTINE generate_coords_vels_vib
1041 :
1042 : ! **************************************************************************************************
1043 : !> \brief calculates componbent of initial velocity of an atom from vibreational modes
1044 : !> \param iatom : global atomic index
1045 : !> \param icart : cartesian index (1, 2 or 3)
1046 : !> \param mass : atomic mass
1047 : !> \param temperature : target temperature of canonical ensemble
1048 : !> \param eigenvalues : array containing all cartesian vibrational mode eigenvalues (frequencies)
1049 : !> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
1050 : !> (displacements)
1051 : !> \param random : array containing uniform distributed random numbers, must be the size
1052 : !> of 3*natoms. Numbers must be between 0 and 1
1053 : !> \param phase : array containing numbers between 0 and 1 that determines for each
1054 : !> vibration mode the ratio of potential energy vs kinetic energy contribution
1055 : !> to total energy
1056 : !> \param dof : total number of degrees of freedom, = 3*natoms
1057 : !> \param scale : scale to make sure the sum of vibrational modes give the correct energy
1058 : !> \return : outputs icart-th cartesian component of initial position of atom iatom
1059 : !> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
1060 : ! **************************************************************************************************
1061 18 : PURE FUNCTION dr_from_vib_data(iatom, &
1062 : icart, &
1063 : mass, &
1064 : temperature, &
1065 36 : eigenvalues, &
1066 18 : eigenvectors, &
1067 18 : random, &
1068 18 : phase, &
1069 : dof, &
1070 : scale) &
1071 : RESULT(res)
1072 : INTEGER, INTENT(IN) :: iatom, icart
1073 : REAL(KIND=dp), INTENT(IN) :: mass, temperature
1074 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
1075 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenvectors
1076 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: random, phase
1077 : INTEGER, INTENT(IN) :: dof
1078 : REAL(KIND=dp), INTENT(IN) :: scale
1079 : REAL(KIND=dp) :: res
1080 :
1081 : INTEGER :: imode, ind
1082 :
1083 18 : res = 0.0_dp
1084 : ! assuming the eigenvalues are sorted in ascending order, the
1085 : ! first three modes are acoustic with zero frequencies. These are
1086 : ! excluded from considerations, and should have been reflected in
1087 : ! the calculation of scale outside this function
1088 18 : IF (mass > 0.0_dp) THEN
1089 : ! eigenvector rows assumed to be grouped in atomic blocks
1090 18 : ind = (iatom - 1)*3 + icart
1091 126 : DO imode = 4, dof
1092 : res = res + &
1093 : SQRT(-2.0_dp*scale*temperature*LOG(1 - random(imode))/mass)/ &
1094 : eigenvalues(imode)* &
1095 : eigenvectors(ind, imode)* &
1096 126 : COS(2.0_dp*pi*phase(imode))
1097 : END DO
1098 : END IF
1099 18 : END FUNCTION dr_from_vib_data
1100 :
1101 : ! **************************************************************************************************
1102 : !> \brief calculates componbent of initial velocity of an atom from vibreational modes
1103 : !> \param iatom : global atomic index
1104 : !> \param icart : cartesian index (1, 2 or 3)
1105 : !> \param mass : atomic mass
1106 : !> \param temperature : target temperature of canonical ensemble
1107 : !> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
1108 : !> (displacements)
1109 : !> \param random : array containing uniform distributed random numbers, must be the size
1110 : !> of 3*natoms. Numbers must be between 0 and 1
1111 : !> \param phase : array containing numbers between 0 and 1 that determines for each
1112 : !> vibration mode the ratio of potential energy vs kinetic energy contribution
1113 : !> to total energy
1114 : !> \param dof : total number of degrees of freedom, = 3*natoms
1115 : !> \param scale : scale to make sure the sum of vibrational modes give the correct energy
1116 : !> \return : outputs icart-th cartesian component of initial velocity of atom iatom
1117 : !> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
1118 : ! **************************************************************************************************
1119 18 : PURE FUNCTION dv_from_vib_data(iatom, &
1120 : icart, &
1121 : mass, &
1122 : temperature, &
1123 36 : eigenvectors, &
1124 18 : random, &
1125 18 : phase, &
1126 : dof, &
1127 : scale) &
1128 : RESULT(res)
1129 : INTEGER, INTENT(IN) :: iatom, icart
1130 : REAL(KIND=dp), INTENT(IN) :: mass, temperature
1131 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenvectors
1132 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: random, phase
1133 : INTEGER, INTENT(IN) :: dof
1134 : REAL(KIND=dp), INTENT(IN) :: scale
1135 : REAL(KIND=dp) :: res
1136 :
1137 : INTEGER :: imode, ind
1138 :
1139 18 : res = 0.0_dp
1140 : ! assuming the eigenvalues are sorted in ascending order, the
1141 : ! first three modes are acoustic with zero frequencies. These are
1142 : ! excluded from considerations, and should have been reflected in
1143 : ! the calculation of scale outside this function
1144 18 : IF (mass > 0.0_dp) THEN
1145 : ! eigenvector rows assumed to be grouped in atomic blocks
1146 18 : ind = (iatom - 1)*3 + icart
1147 126 : DO imode = 4, dof
1148 : res = res - &
1149 : SQRT(-2.0_dp*scale*temperature*LOG(1 - random(imode))/mass)* &
1150 : eigenvectors(ind, imode)* &
1151 126 : SIN(2.0_dp*pi*phase(imode))
1152 : END DO
1153 : END IF
1154 18 : END FUNCTION dv_from_vib_data
1155 :
1156 : ! **************************************************************************************************
1157 : !> \brief Initializing velocities deterministically on all processors, if not given in input
1158 : !> \param simpar ...
1159 : !> \param part ...
1160 : !> \param force_env ...
1161 : !> \param globenv ...
1162 : !> \param md_env ...
1163 : !> \param shell_present ...
1164 : !> \param shell_part ...
1165 : !> \param core_part ...
1166 : !> \param is_fixed ...
1167 : !> \param iw ...
1168 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1169 : ! **************************************************************************************************
1170 1523 : SUBROUTINE generate_velocities(simpar, part, force_env, globenv, md_env, &
1171 1523 : shell_present, shell_part, core_part, is_fixed, iw)
1172 : TYPE(simpar_type), POINTER :: simpar
1173 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1174 : TYPE(force_env_type), POINTER :: force_env
1175 : TYPE(global_environment_type), POINTER :: globenv
1176 : TYPE(md_environment_type), POINTER :: md_env
1177 : LOGICAL, INTENT(IN) :: shell_present
1178 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_part, core_part
1179 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1180 : INTEGER, INTENT(IN) :: iw
1181 :
1182 : INTEGER :: i, natoms
1183 : REAL(KIND=dp) :: mass
1184 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1185 :
1186 1523 : NULLIFY (atomic_kind)
1187 1523 : natoms = SIZE(part)
1188 :
1189 432169 : DO i = 1, natoms
1190 430646 : atomic_kind => part(i)%atomic_kind
1191 430646 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1192 430646 : part(i)%v(1) = 0.0_dp
1193 430646 : part(i)%v(2) = 0.0_dp
1194 430646 : part(i)%v(3) = 0.0_dp
1195 432169 : IF (mass /= 0.0) THEN
1196 430712 : SELECT CASE (is_fixed(i))
1197 : CASE (use_perd_x)
1198 512 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1199 512 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1200 : CASE (use_perd_y)
1201 512 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1202 512 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1203 : CASE (use_perd_z)
1204 512 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1205 512 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1206 : CASE (use_perd_xy)
1207 512 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1208 : CASE (use_perd_xz)
1209 0 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1210 : CASE (use_perd_yz)
1211 0 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1212 : CASE (use_perd_none)
1213 425034 : part(i)%v(1) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1214 425034 : part(i)%v(2) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1215 855234 : part(i)%v(3) = globenv%gaussian_rng_stream%next()/SQRT(mass)
1216 : END SELECT
1217 : END IF
1218 : END DO
1219 :
1220 1523 : CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1221 1523 : CALL soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
1222 :
1223 : ! Initialize the core and the shell velocity. Atom velocities are just
1224 : ! copied so that the initial relative core-shell velocity is zero.
1225 1523 : IF (shell_present) THEN
1226 84 : CALL optimize_shell_core(force_env, part, shell_part, core_part, globenv)
1227 : END IF
1228 1523 : END SUBROUTINE generate_velocities
1229 :
1230 : ! **************************************************************************************************
1231 : !> \brief Direct velocities along a low-curvature direction in order to
1232 : !> favors MD trajectories to cross rapidly over small energy barriers
1233 : !> into neighboring basins.
1234 : !> \param simpar ...
1235 : !> \param part ...
1236 : !> \param force_env ...
1237 : !> \param md_env ...
1238 : !> \param is_fixed ...
1239 : !> \param iw ...
1240 : !> \author Ole Schuett
1241 : ! **************************************************************************************************
1242 1523 : SUBROUTINE soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
1243 : TYPE(simpar_type), POINTER :: simpar
1244 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1245 : TYPE(force_env_type), POINTER :: force_env
1246 : TYPE(md_environment_type), POINTER :: md_env
1247 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1248 : INTEGER, INTENT(IN) :: iw
1249 :
1250 : INTEGER :: i, k
1251 1498 : REAL(KIND=dp), DIMENSION(SIZE(part), 3) :: F, F_t, N, x0
1252 :
1253 1523 : IF (simpar%soften_nsteps <= 0) RETURN !nothing todo
1254 :
1255 275 : IF (ANY(is_fixed /= use_perd_none)) &
1256 0 : CPABORT("Velocitiy softening with constraints is not supported.")
1257 :
1258 : !backup positions
1259 275 : DO i = 1, SIZE(part)
1260 1025 : x0(i, :) = part(i)%r
1261 : END DO
1262 :
1263 525 : DO k = 1, simpar%soften_nsteps
1264 :
1265 : !use normalized velocities as displace direction
1266 5500 : DO i = 1, SIZE(part)
1267 20500 : N(i, :) = part(i)%v
1268 : END DO
1269 33500 : N = N/SQRT(SUM(N**2))
1270 :
1271 : ! displace system temporarly to calculate forces
1272 5500 : DO i = 1, SIZE(part)
1273 20500 : part(i)%r = part(i)%r + simpar%soften_delta*N(i, :)
1274 : END DO
1275 500 : CALL force_env_calc_energy_force(force_env)
1276 :
1277 : ! calculate velocity update direction F_t
1278 5500 : DO i = 1, SIZE(part)
1279 20500 : F(i, :) = part(i)%f
1280 : END DO
1281 33500 : F_t = F - N*SUM(N*F)
1282 :
1283 : ! restore positions and update velocities
1284 5500 : DO i = 1, SIZE(part)
1285 20000 : part(i)%r = x0(i, :)
1286 20500 : part(i)%v = part(i)%v + simpar%soften_alpha*F_t(i, :)
1287 : END DO
1288 :
1289 525 : CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1290 : END DO
1291 :
1292 25 : IF (iw > 0) THEN
1293 0 : WRITE (iw, "(A,T71, I10)") " Velocities softening Steps: ", simpar%soften_nsteps
1294 0 : WRITE (iw, "(A,T71, E10.3)") " Velocities softening NORM(F_t): ", SQRT(SUM(F_t**2))
1295 : END IF
1296 25 : END SUBROUTINE soften_velocities
1297 :
1298 : ! **************************************************************************************************
1299 : !> \brief Scale velocities according to temperature and remove rigid body motion.
1300 : !> \param simpar ...
1301 : !> \param part ...
1302 : !> \param force_env ...
1303 : !> \param md_env ...
1304 : !> \param is_fixed ...
1305 : !> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1306 : ! **************************************************************************************************
1307 2023 : SUBROUTINE normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1308 : TYPE(simpar_type), POINTER :: simpar
1309 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1310 : TYPE(force_env_type), POINTER :: force_env
1311 : TYPE(md_environment_type), POINTER :: md_env
1312 : INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1313 :
1314 : REAL(KIND=dp) :: ekin
1315 : REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vcom
1316 : TYPE(cell_type), POINTER :: cell
1317 :
1318 2023 : NULLIFY (cell)
1319 :
1320 : ! Subtract the vcom
1321 2023 : CALL compute_vcom(part, is_fixed, vcom)
1322 2023 : CALL subtract_vcom(part, is_fixed, vcom)
1323 : ! If requested and the system is not periodic, subtract the angular velocity
1324 2023 : CALL force_env_get(force_env, cell=cell)
1325 8092 : IF (SUM(cell%perd(1:3)) == 0 .AND. simpar%angvel_zero) THEN
1326 4 : CALL compute_rcom(part, is_fixed, rcom)
1327 4 : CALL compute_vang(part, is_fixed, rcom, vang)
1328 4 : CALL subtract_vang(part, is_fixed, rcom, vang)
1329 : END IF
1330 : ! Rescale the velocities
1331 2023 : IF (simpar%do_thermal_region) THEN
1332 2 : CALL rescale_vel_region(part, md_env, simpar)
1333 : ELSE
1334 2021 : ekin = compute_ekin(part)
1335 2021 : CALL rescale_vel(part, simpar, ekin)
1336 : END IF
1337 2023 : END SUBROUTINE normalize_velocities
1338 :
1339 : ! **************************************************************************************************
1340 : !> \brief Computes Ekin, VCOM and Temp for particles
1341 : !> \param subsys ...
1342 : !> \param md_ener ...
1343 : !> \param vsubtract ...
1344 : !> \par History
1345 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1346 : ! **************************************************************************************************
1347 42 : SUBROUTINE reset_vcom(subsys, md_ener, vsubtract)
1348 : TYPE(cp_subsys_type), POINTER :: subsys
1349 : TYPE(md_ener_type), POINTER :: md_ener
1350 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vsubtract
1351 :
1352 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_vcom'
1353 :
1354 : INTEGER :: atom, handle, iatom, ikind, natom, &
1355 : shell_index
1356 42 : INTEGER, DIMENSION(:), POINTER :: atom_list
1357 : LOGICAL :: is_shell
1358 : REAL(KIND=dp) :: ekin_old, imass_c, imass_s, mass, v2
1359 : REAL(KIND=dp), DIMENSION(3) :: tmp, v
1360 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1361 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1362 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1363 : shell_particles
1364 : TYPE(shell_kind_type), POINTER :: shell
1365 :
1366 42 : NULLIFY (particles, atomic_kind, atomic_kinds, atom_list, shell)
1367 42 : CALL timeset(routineN, handle)
1368 :
1369 : CALL cp_subsys_get(subsys, &
1370 : atomic_kinds=atomic_kinds, &
1371 : particles=particles, &
1372 : shell_particles=shell_particles, &
1373 42 : core_particles=core_particles)
1374 :
1375 42 : ekin_old = md_ener%ekin
1376 : ! Possibly subtract a quantity from all velocities
1377 126 : DO ikind = 1, atomic_kinds%n_els
1378 84 : atomic_kind => atomic_kinds%els(ikind)
1379 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, &
1380 84 : natom=natom, mass=mass, shell_active=is_shell, shell=shell)
1381 126 : IF (is_shell) THEN
1382 336 : tmp = 0.5_dp*vsubtract*mass
1383 84 : imass_s = 1.0_dp/shell%mass_shell
1384 84 : imass_c = 1.0_dp/shell%mass_core
1385 3780 : DO iatom = 1, natom
1386 3696 : atom = atom_list(iatom)
1387 3696 : shell_index = particles%els(atom)%shell_index
1388 14784 : shell_particles%els(shell_index)%v = shell_particles%els(shell_index)%v - tmp*imass_s
1389 14784 : core_particles%els(shell_index)%v = core_particles%els(shell_index)%v - tmp*imass_c
1390 14868 : particles%els(atom)%v = particles%els(atom)%v - vsubtract
1391 : END DO
1392 : ELSE
1393 0 : DO iatom = 1, natom
1394 0 : atom = atom_list(iatom)
1395 0 : particles%els(atom)%v = particles%els(atom)%v - vsubtract
1396 : END DO
1397 : END IF
1398 : END DO
1399 : ! Compute Kinetic Energy and COM Velocity
1400 168 : md_ener%vcom = 0.0_dp
1401 42 : md_ener%total_mass = 0.0_dp
1402 42 : md_ener%ekin = 0.0_dp
1403 126 : DO ikind = 1, atomic_kinds%n_els
1404 84 : atomic_kind => atomic_kinds%els(ikind)
1405 84 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom)
1406 84 : v2 = 0.0_dp
1407 84 : v = 0.0_dp
1408 3780 : DO iatom = 1, natom
1409 3696 : atom = atom_list(iatom)
1410 14784 : v2 = v2 + SUM(particles%els(atom)%v**2)
1411 3696 : v(1) = v(1) + particles%els(atom)%v(1)
1412 3696 : v(2) = v(2) + particles%els(atom)%v(2)
1413 3780 : v(3) = v(3) + particles%els(atom)%v(3)
1414 : END DO
1415 84 : md_ener%ekin = md_ener%ekin + 0.5_dp*mass*v2
1416 84 : md_ener%vcom(1) = md_ener%vcom(1) + mass*v(1)
1417 84 : md_ener%vcom(2) = md_ener%vcom(2) + mass*v(2)
1418 84 : md_ener%vcom(3) = md_ener%vcom(3) + mass*v(3)
1419 210 : md_ener%total_mass = md_ener%total_mass + REAL(natom, KIND=dp)*mass
1420 : END DO
1421 168 : md_ener%vcom = md_ener%vcom/md_ener%total_mass
1422 42 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1423 42 : IF (md_ener%nfree /= 0) THEN
1424 42 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1425 : END IF
1426 42 : CALL timestop(handle)
1427 :
1428 42 : END SUBROUTINE reset_vcom
1429 :
1430 : ! **************************************************************************************************
1431 : !> \brief Scale velocities to get the correct temperature
1432 : !> \param subsys ...
1433 : !> \param md_ener ...
1434 : !> \param temp_expected ...
1435 : !> \param temp_tol ...
1436 : !> \param iw ...
1437 : !> \par History
1438 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1439 : ! **************************************************************************************************
1440 12914 : SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw)
1441 : TYPE(cp_subsys_type), POINTER :: subsys
1442 : TYPE(md_ener_type), POINTER :: md_ener
1443 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1444 : INTEGER, INTENT(IN) :: iw
1445 :
1446 : REAL(KIND=dp) :: ekin_old, scale, temp_old
1447 :
1448 12914 : IF (ABS(temp_expected - md_ener%temp_part/kelvin) > temp_tol) THEN
1449 2640 : scale = 0.0_dp
1450 2640 : IF (md_ener%temp_part > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_part)*kelvin)
1451 2640 : ekin_old = md_ener%ekin
1452 2640 : temp_old = md_ener%temp_part
1453 2640 : md_ener%ekin = 0.0_dp
1454 2640 : md_ener%temp_part = 0.0_dp
1455 10560 : md_ener%vcom = 0.0_dp
1456 2640 : md_ener%total_mass = 0.0_dp
1457 :
1458 2640 : CALL scale_velocity_low(subsys, scale, ireg=0, ekin=md_ener%ekin, vcom=md_ener%vcom)
1459 2640 : IF (md_ener%nfree /= 0) THEN
1460 2640 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1461 : END IF
1462 2640 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1463 2640 : IF (iw > 0) THEN
1464 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1465 1320 : 'MD_VEL| Temperature scaled to requested temperature'
1466 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1467 1320 : 'MD_VEL| Old temperature [K]', temp_old, &
1468 2640 : 'MD_VEL| New temperature [K]', md_ener%temp_part
1469 : END IF
1470 : END IF
1471 :
1472 12914 : END SUBROUTINE scale_velocity
1473 :
1474 : ! **************************************************************************************************
1475 : !> \brief Scale velocities of set of regions
1476 : !> \param md_env ...
1477 : !> \param subsys ...
1478 : !> \param md_ener ...
1479 : !> \param simpar ...
1480 : !> \param iw ...
1481 : !> \par author MI
1482 : ! **************************************************************************************************
1483 96 : SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1484 :
1485 : TYPE(md_environment_type), POINTER :: md_env
1486 : TYPE(cp_subsys_type), POINTER :: subsys
1487 : TYPE(md_ener_type), POINTER :: md_ener
1488 : TYPE(simpar_type), POINTER :: simpar
1489 : INTEGER, INTENT(IN) :: iw
1490 :
1491 : INTEGER :: ireg, nfree, nfree_done, nregions
1492 : REAL(KIND=dp) :: ekin, ekin_old, ekin_total_new, fscale, &
1493 : vcom(3), vcom_total(3)
1494 96 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp_new, temp_old
1495 : TYPE(particle_list_type), POINTER :: particles
1496 : TYPE(particle_type), DIMENSION(:), POINTER :: part
1497 : TYPE(thermal_region_type), POINTER :: t_region
1498 : TYPE(thermal_regions_type), POINTER :: thermal_regions
1499 :
1500 96 : NULLIFY (particles, part, thermal_regions, t_region)
1501 96 : CALL cp_subsys_get(subsys, particles=particles)
1502 96 : part => particles%els
1503 96 : CALL get_md_env(md_env, thermal_regions=thermal_regions)
1504 :
1505 96 : nregions = thermal_regions%nregions
1506 96 : nfree_done = 0
1507 96 : ekin_total_new = 0.0_dp
1508 96 : ekin_old = md_ener%ekin
1509 : vcom_total = 0.0_dp
1510 384 : ALLOCATE (temp_new(0:nregions), temp_old(0:nregions))
1511 352 : temp_new = 0.0_dp
1512 352 : temp_old = 0.0_dp
1513 : !loop regions
1514 256 : DO ireg = 1, nregions
1515 160 : NULLIFY (t_region)
1516 160 : t_region => thermal_regions%thermal_region(ireg)
1517 160 : nfree = 3*t_region%npart
1518 160 : ekin = compute_ekin(part, ireg)
1519 160 : IF (nfree > 0) t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1520 160 : temp_old(ireg) = t_region%temperature
1521 160 : IF (t_region%temp_tol > 0.0_dp .AND. &
1522 : ABS(t_region%temp_expected - t_region%temperature/kelvin) > t_region%temp_tol) THEN
1523 2 : fscale = SQRT((t_region%temp_expected/t_region%temperature)*kelvin)
1524 2 : CALL scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1525 2 : t_region%temperature = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1526 2 : temp_new(ireg) = t_region%temperature
1527 : END IF
1528 160 : nfree_done = nfree_done + nfree
1529 256 : ekin_total_new = ekin_total_new + ekin
1530 : END DO
1531 96 : nfree = simpar%nfree - nfree_done
1532 96 : ekin = compute_ekin(part, ireg=0)
1533 96 : IF (nfree > 0) thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1534 96 : temp_old(0) = thermal_regions%temp_reg0
1535 96 : IF (simpar%temp_tol > 0.0_dp .AND. nfree > 0) THEN
1536 0 : IF (ABS(simpar%temp_ext - thermal_regions%temp_reg0/kelvin) > simpar%temp_tol) THEN
1537 0 : fscale = SQRT((simpar%temp_ext/thermal_regions%temp_reg0)*kelvin)
1538 0 : CALL scale_velocity_low(subsys, fscale, 0, ekin, vcom)
1539 0 : thermal_regions%temp_reg0 = 2.0_dp*ekin/REAL(nfree, KIND=dp)*kelvin
1540 0 : temp_new(0) = thermal_regions%temp_reg0
1541 : END IF
1542 : END IF
1543 96 : ekin_total_new = ekin_total_new + ekin
1544 :
1545 96 : md_ener%ekin = ekin_total_new
1546 96 : IF (md_ener%nfree /= 0) THEN
1547 96 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
1548 : END IF
1549 96 : md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1550 96 : IF (iw > 0) THEN
1551 176 : DO ireg = 0, nregions
1552 176 : IF (temp_new(ireg) > 0.0_dp) THEN
1553 : WRITE (UNIT=iw, FMT='(/,T2,A,I0,A)') &
1554 1 : 'MD_VEL| Temperature region ', ireg, ' scaled to requested temperature'
1555 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1556 1 : 'MD_VEL| Old temperature [K]', temp_old(ireg), &
1557 2 : 'MD_VEL| New temperature [K]', temp_new(ireg)
1558 : END IF
1559 : END DO
1560 : END IF
1561 96 : DEALLOCATE (temp_new, temp_old)
1562 :
1563 96 : END SUBROUTINE scale_velocity_region
1564 :
1565 : ! **************************************************************************************************
1566 : !> \brief Scale velocities for a specific region
1567 : !> \param subsys ...
1568 : !> \param fscale ...
1569 : !> \param ireg ...
1570 : !> \param ekin ...
1571 : !> \param vcom ...
1572 : !> \par author MI
1573 : ! **************************************************************************************************
1574 2642 : SUBROUTINE scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1575 :
1576 : TYPE(cp_subsys_type), POINTER :: subsys
1577 : REAL(KIND=dp), INTENT(IN) :: fscale
1578 : INTEGER, INTENT(IN) :: ireg
1579 : REAL(KIND=dp), INTENT(OUT) :: ekin, vcom(3)
1580 :
1581 : INTEGER :: atom, iatom, ikind, my_ireg, natom, &
1582 : shell_index
1583 2642 : INTEGER, DIMENSION(:), POINTER :: atom_list
1584 : LOGICAL :: is_shell
1585 : REAL(KIND=dp) :: imass, mass, tmass, v2
1586 : REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1587 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1588 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1589 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1590 : shell_particles
1591 : TYPE(shell_kind_type), POINTER :: shell
1592 :
1593 2642 : NULLIFY (atomic_kinds, particles, shell_particles, core_particles, shell, atom_list)
1594 :
1595 2642 : my_ireg = ireg
1596 2642 : ekin = 0.0_dp
1597 2642 : tmass = 0.0_dp
1598 2642 : vcom = 0.0_dp
1599 :
1600 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
1601 2642 : shell_particles=shell_particles, core_particles=core_particles)
1602 :
1603 9082 : DO ikind = 1, atomic_kinds%n_els
1604 6440 : atomic_kind => atomic_kinds%els(ikind)
1605 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, &
1606 6440 : natom=natom, shell_active=is_shell, shell=shell)
1607 6440 : IF (is_shell) THEN
1608 124 : imass = 1.0_dp/mass
1609 124 : v2 = 0.0_dp
1610 124 : v = 0.0_dp
1611 5740 : DO iatom = 1, natom
1612 5616 : atom = atom_list(iatom)
1613 : !check region
1614 5616 : IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
1615 :
1616 21080 : particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1617 5270 : shell_index = particles%els(atom)%shell_index
1618 21080 : vs = shell_particles%els(shell_index)%v
1619 21080 : vc = core_particles%els(shell_index)%v
1620 5270 : tmp(1) = imass*(vs(1) - vc(1))
1621 5270 : tmp(2) = imass*(vs(2) - vc(2))
1622 5270 : tmp(3) = imass*(vs(3) - vc(3))
1623 :
1624 5270 : shell_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) + tmp(1)*shell%mass_core
1625 5270 : shell_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) + tmp(2)*shell%mass_core
1626 5270 : shell_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) + tmp(3)*shell%mass_core
1627 :
1628 5270 : core_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) - tmp(1)*shell%mass_shell
1629 5270 : core_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) - tmp(2)*shell%mass_shell
1630 5270 : core_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) - tmp(3)*shell%mass_shell
1631 :
1632 : ! kinetic energy and velocity of COM
1633 21080 : v2 = v2 + SUM(particles%els(atom)%v**2)
1634 5270 : v(1) = v(1) + particles%els(atom)%v(1)
1635 5270 : v(2) = v(2) + particles%els(atom)%v(2)
1636 5270 : v(3) = v(3) + particles%els(atom)%v(3)
1637 5740 : tmass = tmass + mass
1638 : END DO
1639 : ELSE
1640 6316 : v2 = 0.0_dp
1641 6316 : v = 0.0_dp
1642 22826 : DO iatom = 1, natom
1643 16510 : atom = atom_list(iatom)
1644 : !check region
1645 16510 : IF (particles%els(atom)%t_region_index /= my_ireg) CYCLE
1646 :
1647 66040 : particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1648 : ! kinetic energy and velocity of COM
1649 66040 : v2 = v2 + SUM(particles%els(atom)%v**2)
1650 16510 : v(1) = v(1) + particles%els(atom)%v(1)
1651 16510 : v(2) = v(2) + particles%els(atom)%v(2)
1652 16510 : v(3) = v(3) + particles%els(atom)%v(3)
1653 22826 : tmass = tmass + mass
1654 : END DO
1655 : END IF
1656 6440 : ekin = ekin + 0.5_dp*mass*v2
1657 6440 : vcom(1) = vcom(1) + mass*v(1)
1658 6440 : vcom(2) = vcom(2) + mass*v(2)
1659 15522 : vcom(3) = vcom(3) + mass*v(3)
1660 :
1661 : END DO
1662 10568 : vcom = vcom/tmass
1663 :
1664 2642 : END SUBROUTINE scale_velocity_low
1665 :
1666 : ! **************************************************************************************************
1667 : !> \brief Scale internal motion of CORE-SHELL model to the correct temperature
1668 : !> \param subsys ...
1669 : !> \param md_ener ...
1670 : !> \param temp_expected ...
1671 : !> \param temp_tol ...
1672 : !> \param iw ...
1673 : !> \par History
1674 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1675 : ! **************************************************************************************************
1676 1060 : SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw)
1677 : TYPE(cp_subsys_type), POINTER :: subsys
1678 : TYPE(md_ener_type), POINTER :: md_ener
1679 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1680 : INTEGER, INTENT(IN) :: iw
1681 :
1682 : INTEGER :: atom, iatom, ikind, natom, shell_index
1683 1060 : INTEGER, DIMENSION(:), POINTER :: atom_list
1684 : LOGICAL :: is_shell
1685 : REAL(KIND=dp) :: ekin_shell_old, fac_mass, mass, scale, &
1686 : temp_shell_old, v2
1687 : REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1688 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1689 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1690 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
1691 : shell_particles
1692 : TYPE(shell_kind_type), POINTER :: shell
1693 :
1694 1060 : NULLIFY (atom_list, atomic_kinds, atomic_kind, core_particles, particles, shell_particles, shell)
1695 1060 : IF (ABS(temp_expected - md_ener%temp_shell/kelvin) > temp_tol) THEN
1696 80 : scale = 0.0_dp
1697 80 : IF (md_ener%temp_shell > EPSILON(0.0_dp)) scale = SQRT((temp_expected/md_ener%temp_shell)*kelvin)
1698 80 : ekin_shell_old = md_ener%ekin_shell
1699 80 : temp_shell_old = md_ener%temp_shell
1700 80 : md_ener%ekin_shell = 0.0_dp
1701 80 : md_ener%temp_shell = 0.0_dp
1702 :
1703 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, shell_particles=shell_particles, &
1704 80 : core_particles=core_particles)
1705 :
1706 240 : DO ikind = 1, atomic_kinds%n_els
1707 160 : atomic_kind => atomic_kinds%els(ikind)
1708 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom, &
1709 160 : shell_active=is_shell, shell=shell)
1710 240 : IF (is_shell) THEN
1711 160 : fac_mass = 1.0_dp/mass
1712 160 : v2 = 0.0_dp
1713 776 : DO iatom = 1, natom
1714 616 : atom = atom_list(iatom)
1715 616 : shell_index = particles%els(atom)%shell_index
1716 2464 : vs = shell_particles%els(shell_index)%v
1717 2464 : vc = core_particles%els(shell_index)%v
1718 2464 : v = particles%els(atom)%v
1719 616 : tmp(1) = fac_mass*(vc(1) - vs(1))
1720 616 : tmp(2) = fac_mass*(vc(2) - vs(2))
1721 616 : tmp(3) = fac_mass*(vc(3) - vs(3))
1722 :
1723 616 : shell_particles%els(shell_index)%v(1) = v(1) - shell%mass_core*scale*tmp(1)
1724 616 : shell_particles%els(shell_index)%v(2) = v(2) - shell%mass_core*scale*tmp(2)
1725 616 : shell_particles%els(shell_index)%v(3) = v(3) - shell%mass_core*scale*tmp(3)
1726 :
1727 616 : core_particles%els(shell_index)%v(1) = v(1) + shell%mass_shell*scale*tmp(1)
1728 616 : core_particles%els(shell_index)%v(2) = v(2) + shell%mass_shell*scale*tmp(2)
1729 616 : core_particles%els(shell_index)%v(3) = v(3) + shell%mass_shell*scale*tmp(3)
1730 :
1731 2464 : vs = shell_particles%els(shell_index)%v
1732 2464 : vc = core_particles%els(shell_index)%v
1733 616 : tmp(1) = vc(1) - vs(1)
1734 616 : tmp(2) = vc(2) - vs(2)
1735 616 : tmp(3) = vc(3) - vs(3)
1736 2624 : v2 = v2 + SUM(tmp**2)
1737 : END DO
1738 160 : md_ener%ekin_shell = md_ener%ekin_shell + 0.5_dp*shell%mass_core*shell%mass_shell*fac_mass*v2
1739 : END IF
1740 : END DO
1741 80 : IF (md_ener%nfree_shell > 0) THEN
1742 80 : md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell, KIND=dp)*kelvin
1743 : END IF
1744 80 : md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell
1745 80 : IF (iw > 0) THEN
1746 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1747 40 : 'MD_VEL| Temperature of shell internal motion scaled to requested temperature'
1748 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1749 40 : 'MD_VEL| Old temperature [K]', temp_shell_old, &
1750 80 : 'MD_VEL| New temperature [K]', md_ener%temp_shell
1751 : END IF
1752 : END IF
1753 :
1754 1060 : END SUBROUTINE scale_velocity_internal
1755 :
1756 : ! **************************************************************************************************
1757 : !> \brief Scale barostat velocities to get the desired temperature
1758 : !> \param md_env ...
1759 : !> \param md_ener ...
1760 : !> \param temp_expected ...
1761 : !> \param temp_tol ...
1762 : !> \param iw ...
1763 : !> \par History
1764 : !> MI 02.2008
1765 : ! **************************************************************************************************
1766 40 : SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw)
1767 : TYPE(md_environment_type), POINTER :: md_env
1768 : TYPE(md_ener_type), POINTER :: md_ener
1769 : REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1770 : INTEGER, INTENT(IN) :: iw
1771 :
1772 : INTEGER :: i, j, nfree
1773 : REAL(KIND=dp) :: ekin_old, scale, temp_old
1774 40 : TYPE(npt_info_type), POINTER :: npt(:, :)
1775 : TYPE(simpar_type), POINTER :: simpar
1776 :
1777 40 : NULLIFY (npt, simpar)
1778 40 : CALL get_md_env(md_env, simpar=simpar, npt=npt)
1779 40 : IF (ABS(temp_expected - md_ener%temp_baro/kelvin) > temp_tol) THEN
1780 2 : scale = 0.0_dp
1781 2 : IF (md_ener%temp_baro > 0.0_dp) scale = SQRT((temp_expected/md_ener%temp_baro)*kelvin)
1782 2 : ekin_old = md_ener%baro_kin
1783 2 : temp_old = md_ener%temp_baro
1784 2 : md_ener%baro_kin = 0.0_dp
1785 2 : md_ener%temp_baro = 0.0_dp
1786 : IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npe_i_ensemble &
1787 2 : .OR. simpar%ensemble == npt_ia_ensemble) THEN
1788 0 : npt(1, 1)%v = npt(1, 1)%v*scale
1789 0 : md_ener%baro_kin = 0.5_dp*npt(1, 1)%v**2*npt(1, 1)%mass
1790 : ELSE IF (simpar%ensemble == npt_f_ensemble .OR. simpar%ensemble == npe_f_ensemble) THEN
1791 2 : md_ener%baro_kin = 0.0_dp
1792 8 : DO i = 1, 3
1793 26 : DO j = 1, 3
1794 18 : npt(i, j)%v = npt(i, j)%v*scale
1795 24 : md_ener%baro_kin = md_ener%baro_kin + 0.5_dp*npt(i, j)%v**2*npt(i, j)%mass
1796 : END DO
1797 : END DO
1798 : END IF
1799 2 : nfree = SIZE(npt, 1)*SIZE(npt, 2)
1800 2 : md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/REAL(nfree, dp)*kelvin
1801 2 : IF (iw > 0) THEN
1802 : WRITE (UNIT=iw, FMT='(/,T2,A)') &
1803 1 : 'MD_VEL| Temperature of barostat motion scaled to requested temperature'
1804 : WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
1805 1 : 'MD_VEL| Old temperature [K]', temp_old, &
1806 2 : 'MD_VEL| New temperature [K]', md_ener%temp_baro
1807 : END IF
1808 : END IF
1809 :
1810 40 : END SUBROUTINE scale_velocity_baro
1811 :
1812 : ! **************************************************************************************************
1813 : !> \brief Perform all temperature manipulations during a QS MD run.
1814 : !> \param simpar ...
1815 : !> \param md_env ...
1816 : !> \param md_ener ...
1817 : !> \param force_env ...
1818 : !> \param logger ...
1819 : !> \par History
1820 : !> Creation (15.09.2003,MK)
1821 : !> adapted to force_env (05.10.2003,fawzi)
1822 : !> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1823 : ! **************************************************************************************************
1824 40312 : SUBROUTINE temperature_control(simpar, md_env, md_ener, force_env, logger)
1825 :
1826 : TYPE(simpar_type), POINTER :: simpar
1827 : TYPE(md_environment_type), POINTER :: md_env
1828 : TYPE(md_ener_type), POINTER :: md_ener
1829 : TYPE(force_env_type), POINTER :: force_env
1830 : TYPE(cp_logger_type), POINTER :: logger
1831 :
1832 : CHARACTER(LEN=*), PARAMETER :: routineN = 'temperature_control'
1833 :
1834 : INTEGER :: handle, iw
1835 : TYPE(cp_subsys_type), POINTER :: subsys
1836 : TYPE(mp_para_env_type), POINTER :: para_env
1837 :
1838 40312 : CALL timeset(routineN, handle)
1839 40312 : NULLIFY (subsys, para_env)
1840 40312 : CPASSERT(ASSOCIATED(simpar))
1841 40312 : CPASSERT(ASSOCIATED(md_ener))
1842 40312 : CPASSERT(ASSOCIATED(force_env))
1843 40312 : CALL force_env_get(force_env, subsys=subsys, para_env=para_env)
1844 : iw = cp_print_key_unit_nr(logger, force_env%root_section, "MOTION%MD%PRINT%PROGRAM_RUN_INFO", &
1845 40312 : extension=".mdLog")
1846 :
1847 : ! Control the particle motion
1848 40312 : IF (simpar%do_thermal_region) THEN
1849 96 : CALL scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1850 : ELSE
1851 40216 : IF (simpar%temp_tol > 0.0_dp) THEN
1852 12870 : CALL scale_velocity(subsys, md_ener, simpar%temp_ext, simpar%temp_tol, iw)
1853 : END IF
1854 : END IF
1855 : ! Control the internal core-shell motion
1856 40312 : IF (simpar%temp_sh_tol > 0.0_dp) THEN
1857 1060 : CALL scale_velocity_internal(subsys, md_ener, simpar%temp_sh_ext, simpar%temp_sh_tol, iw)
1858 : END IF
1859 : ! Control cell motion
1860 42832 : SELECT CASE (simpar%ensemble)
1861 : CASE (nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
1862 : npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble, npt_ia_ensemble)
1863 40312 : IF (simpar%temp_baro_tol > 0.0_dp) THEN
1864 40 : CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw)
1865 : END IF
1866 : END SELECT
1867 :
1868 : CALL cp_print_key_finished_output(iw, logger, force_env%root_section, &
1869 40312 : "MOTION%MD%PRINT%PROGRAM_RUN_INFO")
1870 40312 : CALL timestop(handle)
1871 40312 : END SUBROUTINE temperature_control
1872 :
1873 : ! **************************************************************************************************
1874 : !> \brief Set to 0 the velocity of the COM along MD runs, if required.
1875 : !> \param md_ener ...
1876 : !> \param force_env ...
1877 : !> \param md_section ...
1878 : !> \param logger ...
1879 : !> \par History
1880 : !> Creation (29.04.2007,MI)
1881 : !> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1882 : ! **************************************************************************************************
1883 80624 : SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
1884 :
1885 : TYPE(md_ener_type), POINTER :: md_ener
1886 : TYPE(force_env_type), POINTER :: force_env
1887 : TYPE(section_vals_type), POINTER :: md_section
1888 : TYPE(cp_logger_type), POINTER :: logger
1889 :
1890 : CHARACTER(LEN=*), PARAMETER :: routineN = 'comvel_control'
1891 :
1892 : INTEGER :: handle, iw
1893 : LOGICAL :: explicit
1894 : REAL(KIND=dp) :: comvel_tol, temp_old, vel_com
1895 : REAL(KIND=dp), DIMENSION(3) :: vcom_old
1896 : TYPE(cp_subsys_type), POINTER :: subsys
1897 :
1898 40312 : CALL timeset(routineN, handle)
1899 40312 : NULLIFY (subsys)
1900 40312 : CPASSERT(ASSOCIATED(force_env))
1901 40312 : CALL force_env_get(force_env, subsys=subsys)
1902 :
1903 : ! Print COMVEL and COM Position
1904 40312 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%CENTER_OF_MASS", extension=".mdLog")
1905 40312 : IF (iw > 0) THEN
1906 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1907 7858 : "MD_VEL| Centre of mass motion (COM)"
1908 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1909 7858 : "MD_VEL| VCOM [a.u.]", md_ener%vcom(1:3)
1910 : END IF
1911 40312 : CALL cp_print_key_finished_output(iw, logger, md_section, "PRINT%CENTER_OF_MASS")
1912 :
1913 : ! If requested rescale COMVEL
1914 40312 : CALL section_vals_val_get(md_section, "COMVEL_TOL", explicit=explicit)
1915 40312 : IF (explicit) THEN
1916 826 : CALL section_vals_val_get(md_section, "COMVEL_TOL", r_val=comvel_tol)
1917 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1918 826 : extension=".mdLog")
1919 826 : vel_com = SQRT(md_ener%vcom(1)**2 + md_ener%vcom(2)**2 + md_ener%vcom(3)**2)
1920 :
1921 : ! Subtract the velocity of the COM, if requested
1922 826 : IF (vel_com > comvel_tol) THEN
1923 42 : temp_old = md_ener%temp_part/kelvin
1924 168 : vcom_old = md_ener%vcom
1925 42 : CALL reset_vcom(subsys, md_ener, vsubtract=vcom_old)
1926 42 : CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
1927 42 : IF (iw > 0) THEN
1928 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1929 21 : "MD_VEL| Old VCOM [a.u.]", vcom_old(1:3)
1930 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
1931 21 : "MD_VEL| New VCOM [a.u.]", md_ener%vcom(1:3)
1932 : END IF
1933 : END IF
1934 : CALL cp_print_key_finished_output(iw, logger, md_section, &
1935 826 : "PRINT%PROGRAM_RUN_INFO")
1936 : END IF
1937 :
1938 40312 : CALL timestop(handle)
1939 40312 : END SUBROUTINE comvel_control
1940 :
1941 : ! **************************************************************************************************
1942 : !> \brief Set to 0 the angular velocity along MD runs, if required.
1943 : !> \param md_ener ...
1944 : !> \param force_env ...
1945 : !> \param md_section ...
1946 : !> \param logger ...
1947 : !> \par History
1948 : !> Creation (10.2009) Teodoro Laino [tlaino]
1949 : ! **************************************************************************************************
1950 40312 : SUBROUTINE angvel_control(md_ener, force_env, md_section, logger)
1951 :
1952 : TYPE(md_ener_type), POINTER :: md_ener
1953 : TYPE(force_env_type), POINTER :: force_env
1954 : TYPE(section_vals_type), POINTER :: md_section
1955 : TYPE(cp_logger_type), POINTER :: logger
1956 :
1957 : CHARACTER(LEN=*), PARAMETER :: routineN = 'angvel_control'
1958 :
1959 : INTEGER :: handle, ifixd, imolecule_kind, iw, natoms
1960 40312 : INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
1961 : LOGICAL :: explicit
1962 : REAL(KIND=dp) :: angvel_tol, rcom(3), temp_old, vang(3), &
1963 : vang_new(3)
1964 : TYPE(cell_type), POINTER :: cell
1965 : TYPE(cp_subsys_type), POINTER :: subsys
1966 40312 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
1967 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1968 40312 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1969 : TYPE(molecule_kind_type), POINTER :: molecule_kind
1970 : TYPE(particle_list_type), POINTER :: particles
1971 :
1972 40312 : CALL timeset(routineN, handle)
1973 : ! If requested rescale ANGVEL
1974 40312 : CALL section_vals_val_get(md_section, "ANGVEL_TOL", explicit=explicit)
1975 40312 : IF (explicit) THEN
1976 40 : NULLIFY (subsys, cell)
1977 40 : CPASSERT(ASSOCIATED(force_env))
1978 40 : CALL force_env_get(force_env, subsys=subsys, cell=cell)
1979 :
1980 160 : IF (SUM(cell%perd(1:3)) == 0) THEN
1981 40 : CALL section_vals_val_get(md_section, "ANGVEL_TOL", r_val=angvel_tol)
1982 : iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1983 40 : extension=".mdLog")
1984 :
1985 : CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, &
1986 40 : particles=particles)
1987 :
1988 40 : natoms = SIZE(particles%els)
1989 : ! Build a list of all fixed atoms (if any)
1990 120 : ALLOCATE (is_fixed(natoms))
1991 :
1992 600 : is_fixed = use_perd_none
1993 40 : molecule_kind_set => molecule_kinds%els
1994 600 : DO imolecule_kind = 1, molecule_kinds%n_els
1995 560 : molecule_kind => molecule_kind_set(imolecule_kind)
1996 560 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
1997 600 : IF (ASSOCIATED(fixd_list)) THEN
1998 0 : DO ifixd = 1, SIZE(fixd_list)
1999 0 : IF (.NOT. fixd_list(ifixd)%restraint%active) &
2000 0 : is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2001 : END DO
2002 : END IF
2003 : END DO
2004 :
2005 : ! If requested and the system is not periodic, subtract the angular velocity
2006 40 : CALL compute_rcom(particles%els, is_fixed, rcom)
2007 40 : CALL compute_vang(particles%els, is_fixed, rcom, vang)
2008 : ! SQRT(DOT_PRODUCT(vang,vang))>angvel_tol
2009 160 : IF (DOT_PRODUCT(vang, vang) > (angvel_tol*angvel_tol)) THEN
2010 2 : CALL subtract_vang(particles%els, is_fixed, rcom, vang)
2011 :
2012 : ! Rescale velocities after removal
2013 2 : temp_old = md_ener%temp_part/kelvin
2014 2 : CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
2015 2 : CALL compute_vang(particles%els, is_fixed, rcom, vang_new)
2016 2 : IF (iw > 0) THEN
2017 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
2018 1 : 'MD_VEL| Old VANG [a.u.]', vang(1:3)
2019 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
2020 1 : 'MD_VEL| New VANG [a.u.]', vang_new(1:3)
2021 : END IF
2022 : END IF
2023 :
2024 40 : DEALLOCATE (is_fixed)
2025 :
2026 : CALL cp_print_key_finished_output(iw, logger, md_section, &
2027 80 : "PRINT%PROGRAM_RUN_INFO")
2028 : END IF
2029 : END IF
2030 :
2031 40312 : CALL timestop(handle)
2032 80624 : END SUBROUTINE angvel_control
2033 :
2034 : ! **************************************************************************************************
2035 : !> \brief Initialize Velocities for MD runs
2036 : !> \param force_env ...
2037 : !> \param simpar ...
2038 : !> \param globenv ...
2039 : !> \param md_env ...
2040 : !> \param md_section ...
2041 : !> \param constraint_section ...
2042 : !> \param write_binary_restart_file ...
2043 : !> \par History
2044 : !> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
2045 : ! **************************************************************************************************
2046 3534 : SUBROUTINE setup_velocities(force_env, simpar, globenv, md_env, md_section, &
2047 : constraint_section, write_binary_restart_file)
2048 :
2049 : TYPE(force_env_type), POINTER :: force_env
2050 : TYPE(simpar_type), POINTER :: simpar
2051 : TYPE(global_environment_type), POINTER :: globenv
2052 : TYPE(md_environment_type), POINTER :: md_env
2053 : TYPE(section_vals_type), POINTER :: md_section, constraint_section
2054 : LOGICAL, INTENT(IN) :: write_binary_restart_file
2055 :
2056 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_velocities'
2057 :
2058 : INTEGER :: handle, nconstraint, nconstraint_fixd
2059 : LOGICAL :: apply_cns0, shell_adiabatic, &
2060 : shell_present
2061 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
2062 : TYPE(cell_type), POINTER :: cell
2063 : TYPE(cp_subsys_type), POINTER :: subsys
2064 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2065 : TYPE(mp_para_env_type), POINTER :: para_env
2066 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
2067 : shell_particles
2068 1767 : TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
2069 1767 : shell_particle_set
2070 : TYPE(section_vals_type), POINTER :: force_env_section, print_section, &
2071 : subsys_section
2072 :
2073 1767 : CALL timeset(routineN, handle)
2074 :
2075 1767 : NULLIFY (atomic_kinds, cell, para_env, subsys, molecule_kinds, core_particles, particles)
2076 1767 : NULLIFY (shell_particles, core_particle_set, particle_set, shell_particle_set)
2077 1767 : NULLIFY (force_env_section, print_section, subsys_section)
2078 :
2079 1767 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
2080 1767 : apply_cns0 = .FALSE.
2081 1767 : IF (simpar%constraint) THEN
2082 310 : CALL section_vals_val_get(constraint_section, "CONSTRAINT_INIT", l_val=apply_cns0)
2083 : END IF
2084 : ! Always initialize velocities and possibly restart them
2085 : CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env, &
2086 1767 : force_env_section=force_env_section)
2087 1767 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2088 :
2089 : CALL cp_subsys_get(subsys, &
2090 : atomic_kinds=atomic_kinds, &
2091 : core_particles=core_particles, &
2092 : molecule_kinds=molecule_kinds, &
2093 : particles=particles, &
2094 1767 : shell_particles=shell_particles)
2095 :
2096 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, &
2097 : shell_present=shell_present, &
2098 1767 : shell_adiabatic=shell_adiabatic)
2099 :
2100 1767 : NULLIFY (core_particle_set)
2101 : NULLIFY (particle_set)
2102 1767 : NULLIFY (shell_particle_set)
2103 1767 : particle_set => particles%els
2104 :
2105 1767 : IF (shell_present .AND. shell_adiabatic) THEN
2106 : ! Constraints are not yet implemented for core-shell models generally
2107 : CALL get_molecule_kind_set(molecule_kind_set=molecule_kinds%els, &
2108 : nconstraint=nconstraint, &
2109 132 : nconstraint_fixd=nconstraint_fixd)
2110 132 : IF (nconstraint - nconstraint_fixd /= 0) &
2111 0 : CPABORT("Only the fixed atom constraint is implemented for core-shell models")
2112 : !MK CPPostcondition(.NOT.simpar%constraint,cp_failure_level,routineP,failure)
2113 132 : CPASSERT(ASSOCIATED(shell_particles))
2114 132 : CPASSERT(ASSOCIATED(core_particles))
2115 132 : shell_particle_set => shell_particles%els
2116 132 : core_particle_set => core_particles%els
2117 : END IF
2118 :
2119 : CALL initialize_velocities(simpar, &
2120 : particle_set, &
2121 : molecule_kinds=molecule_kinds, &
2122 : force_env=force_env, &
2123 : globenv=globenv, &
2124 : md_env=md_env, &
2125 : label="Velocities initialization", &
2126 : print_section=print_section, &
2127 : subsys_section=subsys_section, &
2128 : shell_present=(shell_present .AND. shell_adiabatic), &
2129 : shell_part=shell_particle_set, &
2130 : core_part=core_particle_set, &
2131 : force_rescaling=.FALSE., &
2132 : para_env=para_env, &
2133 3402 : write_binary_restart_file=write_binary_restart_file)
2134 :
2135 : ! Apply constraints if required and rescale velocities..
2136 1767 : IF (simpar%ensemble /= reftraj_ensemble) THEN
2137 1729 : IF (apply_cns0) THEN
2138 24 : CALL force_env_calc_energy_force(force_env, calc_force=.TRUE.)
2139 : CALL force_env_shake(force_env, &
2140 : shake_tol=simpar%shake_tol, &
2141 : log_unit=simpar%info_constraint, &
2142 : lagrange_mult=simpar%lagrange_multipliers, &
2143 : dump_lm=simpar%dump_lm, &
2144 24 : compold=.TRUE.)
2145 : CALL force_env_rattle(force_env, shake_tol=simpar%shake_tol, &
2146 : log_unit=simpar%info_constraint, lagrange_mult=simpar%lagrange_multipliers, &
2147 24 : dump_lm=simpar%dump_lm, reset=.TRUE.)
2148 24 : IF (simpar%do_respa) THEN
2149 : CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
2150 0 : calc_force=.TRUE.)
2151 : CALL force_env_shake(force_env%sub_force_env(1)%force_env, &
2152 : shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2153 0 : lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, compold=.TRUE.)
2154 : CALL force_env_rattle(force_env%sub_force_env(1)%force_env, &
2155 : shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2156 0 : lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, reset=.TRUE.)
2157 : END IF
2158 : ! Reinitialize velocities rescaling properly after rattle
2159 24 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2160 24 : CALL update_subsys(subsys_section, force_env, .FALSE., write_binary_restart_file)
2161 : CALL initialize_velocities(simpar, &
2162 : particle_set, &
2163 : molecule_kinds=molecule_kinds, &
2164 : force_env=force_env, &
2165 : globenv=globenv, &
2166 : md_env=md_env, &
2167 : label="Re-Initializing velocities after applying constraints", &
2168 : print_section=print_section, &
2169 : subsys_section=subsys_section, &
2170 : shell_present=(shell_present .AND. shell_adiabatic), &
2171 : shell_part=shell_particle_set, &
2172 : core_part=core_particle_set, &
2173 : force_rescaling=.TRUE., &
2174 : para_env=para_env, &
2175 48 : write_binary_restart_file=write_binary_restart_file)
2176 : END IF
2177 : END IF
2178 :
2179 : ! Perform setup for a cascade run
2180 1767 : CALL initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2181 :
2182 1767 : CALL timestop(handle)
2183 :
2184 1767 : END SUBROUTINE setup_velocities
2185 :
2186 : ! **************************************************************************************************
2187 : !> \brief Perform the initialization for a cascade run
2188 : !> \param simpar ...
2189 : !> \param particle_set ...
2190 : !> \param molecule_kinds ...
2191 : !> \param md_section ...
2192 : !> \date 05.02.2012
2193 : !> \author Matthias Krack (MK)
2194 : !> \version 1.0
2195 : ! **************************************************************************************************
2196 1767 : SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2197 :
2198 : TYPE(simpar_type), POINTER :: simpar
2199 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2200 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2201 : TYPE(section_vals_type), POINTER :: md_section
2202 :
2203 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_cascade'
2204 :
2205 : CHARACTER(len=2*default_string_length) :: line
2206 : INTEGER :: handle, iatom, ifixd, imolecule_kind, &
2207 : iparticle, iw, natom, nparticle
2208 1767 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_index, is_fixed
2209 : LOGICAL :: init_cascade, is_ok, no_read_error
2210 : REAL(KIND=dp) :: ecom, ekin, energy, norm, temp, &
2211 : temperature
2212 1767 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: matom, weight
2213 1767 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vatom
2214 : REAL(KIND=dp), DIMENSION(3) :: vcom
2215 : TYPE(atomic_kind_type), POINTER :: atomic_kind
2216 : TYPE(cp_logger_type), POINTER :: logger
2217 : TYPE(cp_sll_val_type), POINTER :: atom_list
2218 1767 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
2219 1767 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
2220 : TYPE(molecule_kind_type), POINTER :: molecule_kind
2221 : TYPE(section_vals_type), POINTER :: atom_list_section, cascade_section, &
2222 : print_section
2223 : TYPE(val_type), POINTER :: val
2224 :
2225 1767 : CALL timeset(routineN, handle)
2226 :
2227 1767 : NULLIFY (atom_list)
2228 1767 : NULLIFY (atom_list_section)
2229 1767 : NULLIFY (atomic_kind)
2230 1767 : NULLIFY (cascade_section)
2231 1767 : NULLIFY (fixd_list)
2232 1767 : NULLIFY (molecule_kind)
2233 1767 : NULLIFY (molecule_kind_set)
2234 1767 : NULLIFY (logger)
2235 1767 : NULLIFY (val)
2236 :
2237 1767 : logger => cp_get_default_logger()
2238 1767 : print_section => section_vals_get_subs_vals(md_section, "PRINT")
2239 1767 : iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
2240 :
2241 1767 : cascade_section => section_vals_get_subs_vals(md_section, "CASCADE")
2242 1767 : CALL section_vals_val_get(cascade_section, "_SECTION_PARAMETERS_", l_val=init_cascade)
2243 :
2244 1767 : nparticle = SIZE(particle_set)
2245 :
2246 1767 : IF (init_cascade) THEN
2247 :
2248 2 : CALL section_vals_val_get(cascade_section, "ENERGY", r_val=energy)
2249 2 : IF (energy < 0.0_dp) &
2250 0 : CPABORT("Error occurred reading &CASCADE section: Negative energy found")
2251 :
2252 2 : IF (iw > 0) THEN
2253 1 : ekin = cp_unit_from_cp2k(energy, "keV")
2254 : WRITE (UNIT=iw, FMT="(/,T2,A,T61,F20.6)") &
2255 1 : "CASCADE| Energy [keV]", ekin
2256 : WRITE (UNIT=iw, FMT="(T2,A)") &
2257 1 : "CASCADE|"
2258 : END IF
2259 :
2260 : ! Read the atomic velocities given in the input file
2261 2 : atom_list_section => section_vals_get_subs_vals(cascade_section, "ATOM_LIST")
2262 2 : CALL section_vals_val_get(atom_list_section, "_DEFAULT_KEYWORD_", n_rep_val=natom)
2263 2 : CALL section_vals_list_get(atom_list_section, "_DEFAULT_KEYWORD_", list=atom_list)
2264 2 : IF (natom <= 0) &
2265 0 : CPABORT("Error occurred reading &CASCADE section: No atom list found")
2266 :
2267 2 : IF (iw > 0) THEN
2268 : WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
2269 1 : "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "weight"
2270 : END IF
2271 :
2272 6 : ALLOCATE (atom_index(natom))
2273 6 : ALLOCATE (matom(natom))
2274 6 : ALLOCATE (vatom(3, natom))
2275 4 : ALLOCATE (weight(natom))
2276 :
2277 8 : DO iatom = 1, natom
2278 6 : is_ok = cp_sll_val_next(atom_list, val)
2279 6 : CALL val_get(val, c_val=line)
2280 : ! Read atomic index, velocity vector, and weight
2281 6 : no_read_error = .FALSE.
2282 6 : READ (UNIT=line, FMT=*, ERR=999) atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2283 : no_read_error = .TRUE.
2284 : 999 IF (.NOT. no_read_error) &
2285 0 : CPABORT("Error occurred reading &CASCADE section. Last line read <"//TRIM(line)//">")
2286 6 : IF ((atom_index(iatom) <= 0) .OR. ((atom_index(iatom) > nparticle))) &
2287 0 : CPABORT("Error occurred reading &CASCADE section: Invalid atom index found")
2288 6 : IF (weight(iatom) < 0.0_dp) &
2289 0 : CPABORT("Error occurred reading &CASCADE section: Negative weight found")
2290 8 : IF (iw > 0) THEN
2291 : WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
2292 3 : "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2293 : END IF
2294 : END DO
2295 :
2296 : ! Normalise velocities and weights
2297 : norm = 0.0_dp
2298 8 : DO iatom = 1, natom
2299 6 : iparticle = atom_index(iatom)
2300 6 : IF (particle_set(iparticle)%shell_index /= 0) THEN
2301 0 : CPWARN("Warning: The primary knock-on atom is a core-shell atom")
2302 : END IF
2303 6 : atomic_kind => particle_set(iparticle)%atomic_kind
2304 6 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=matom(iatom))
2305 8 : norm = norm + matom(iatom)*weight(iatom)
2306 : END DO
2307 8 : weight(:) = matom(:)*weight(:)*energy/norm
2308 8 : DO iatom = 1, natom
2309 24 : norm = SQRT(DOT_PRODUCT(vatom(1:3, iatom), vatom(1:3, iatom)))
2310 26 : vatom(1:3, iatom) = vatom(1:3, iatom)/norm
2311 : END DO
2312 :
2313 2 : IF (iw > 0) THEN
2314 : WRITE (UNIT=iw, FMT="(T2,A)") &
2315 1 : "CASCADE|", &
2316 1 : "CASCADE| Normalised velocities and additional kinetic energy [keV]", &
2317 2 : "CASCADE|"
2318 : WRITE (UNIT=iw, FMT="(T2,A,T11,A,3(11X,A),9X,A)") &
2319 1 : "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "E(kin)"
2320 4 : DO iatom = 1, natom
2321 3 : ekin = cp_unit_from_cp2k(weight(iatom), "keV")
2322 : WRITE (UNIT=iw, FMT="(T2,A,I10,4(1X,F14.6))") &
2323 4 : "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), ekin
2324 : END DO
2325 : END IF
2326 :
2327 : ! Apply velocity modifications
2328 8 : DO iatom = 1, natom
2329 6 : iparticle = atom_index(iatom)
2330 : particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + &
2331 26 : SQRT(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3, iatom)
2332 : END DO
2333 :
2334 2 : DEALLOCATE (atom_index)
2335 2 : DEALLOCATE (matom)
2336 2 : DEALLOCATE (vatom)
2337 2 : DEALLOCATE (weight)
2338 :
2339 6 : IF (iw > 0) THEN
2340 : ! Build a list of all fixed atoms (if any)
2341 3 : ALLOCATE (is_fixed(nparticle))
2342 97 : is_fixed = use_perd_none
2343 1 : molecule_kind_set => molecule_kinds%els
2344 2 : DO imolecule_kind = 1, molecule_kinds%n_els
2345 1 : molecule_kind => molecule_kind_set(imolecule_kind)
2346 1 : CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2347 2 : IF (ASSOCIATED(fixd_list)) THEN
2348 0 : DO ifixd = 1, SIZE(fixd_list)
2349 0 : IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2350 : END DO
2351 : END IF
2352 : END DO
2353 : ! Compute vcom, ecom and ekin for printout
2354 1 : CALL compute_vcom(particle_set, is_fixed, vcom, ecom)
2355 1 : ekin = compute_ekin(particle_set) - ecom
2356 1 : IF (simpar%nfree == 0) THEN
2357 0 : CPASSERT(ekin == 0.0_dp)
2358 0 : temp = 0.0_dp
2359 : ELSE
2360 1 : temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
2361 : END IF
2362 1 : temperature = cp_unit_from_cp2k(temp, "K")
2363 : WRITE (UNIT=iw, FMT="(T2,A)") &
2364 1 : "CASCADE|"
2365 : WRITE (UNIT=iw, FMT="(T2,A,T61,F20.6)") &
2366 1 : "CASCADE| Temperature after cascade initialization [K]", temperature
2367 : WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,ES16.8))") &
2368 1 : "CASCADE| COM velocity", vcom(1:3)
2369 : !MK ! compute and log rcom and vang if not periodic
2370 : !MK CALL force_env_get(force_env,cell=cell)
2371 : !MK IF (SUM(cell%perd(1:3)) == 0) THEN
2372 : !MK CALL compute_rcom(particle_set,is_fixed,rcom)
2373 : !MK CALL compute_vang(particle_set,is_fixed,rcom,vang)
2374 : !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3)
2375 : !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3)
2376 : !MK END IF
2377 2 : DEALLOCATE (is_fixed)
2378 : END IF
2379 :
2380 : END IF
2381 :
2382 1767 : CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
2383 :
2384 1767 : CALL timestop(handle)
2385 :
2386 3534 : END SUBROUTINE initialize_cascade
2387 :
2388 : END MODULE md_vel_utils
|
