(git:1f9fd2c)
Loading...
Searching...
No Matches
md_vel_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \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! **************************************************************************************************
17MODULE 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_transform_input_cartesian, cell_type, use_perd_none, use_perd_x, use_perd_xy, &
26 use_perd_xyz, use_perd_xz, use_perd_y, 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
92CONTAINS
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 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 rcom(:) = 0.0_dp
114 denom = 0.0_dp
115 DO i = 1, SIZE(part)
116 atomic_kind => part(i)%atomic_kind
117 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
118 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 rcom(1) = rcom(1) + part(i)%r(1)*mass
121 rcom(2) = rcom(2) + part(i)%r(2)*mass
122 rcom(3) = rcom(3) + part(i)%r(3)*mass
123 denom = denom + mass
124 END SELECT
125 END DO
126 rcom = rcom/denom
127
128 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 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 vcom = 0.0_dp
152 denom = 0.0_dp
153 DO i = 1, SIZE(part)
154 atomic_kind => part(i)%atomic_kind
155 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
156 IF (mass /= 0.0) THEN
157 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 vcom(1) = vcom(1) + part(i)%v(1)*mass
160 vcom(2) = vcom(2) + part(i)%v(2)*mass
161 vcom(3) = vcom(3) + part(i)%v(3)*mass
162 denom = denom + mass
163 END SELECT
164 END IF
165 END DO
166 vcom = vcom/denom
167 IF (PRESENT(ecom)) THEN
168 ecom = 0.5_dp*denom*sum(vcom*vcom)
169 END IF
170
171 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 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 DO i = 1, SIZE(part)
191 atomic_kind => part(i)%atomic_kind
192 CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_shell)
193 IF (is_shell) THEN
194 shell_part(part(i)%shell_index)%v(:) = part(i)%v(:)
195 core_part(part(i)%shell_index)%v(:) = part(i)%v(:)
196 END IF
197 END DO
198
199 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 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 NULLIFY (atomic_kind)
222 ekin = 0.0_dp
223 IF (PRESENT(ireg)) THEN
224 DO i = 1, SIZE(part)
225 IF (part(i)%t_region_index == ireg) THEN
226 atomic_kind => part(i)%atomic_kind
227 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
228 ekin = ekin + 0.5_dp*mass*sum(part(i)%v(:)*part(i)%v(:))
229 END IF
230 END DO
231 ELSE
232 DO i = 1, SIZE(part)
233 atomic_kind => part(i)%atomic_kind
234 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
235 ekin = ekin + 0.5_dp*mass*sum(part(i)%v(:)*part(i)%v(:))
236 END DO
237 END IF
238
239 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 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 IF (PRESENT(ireg) .AND. PRESENT(nfree) .AND. PRESENT(temp)) THEN
269 my_ireg = ireg
270 my_nfree = nfree
271 my_temp = temp
272 ELSEIF (PRESENT(nfree)) THEN
273 my_ireg = 0
274 my_nfree = nfree
275 my_temp = simpar%temp_ext
276 ELSE
277 my_ireg = 0
278 my_nfree = simpar%nfree
279 my_temp = simpar%temp_ext
280 END IF
281 IF (my_nfree /= 0) THEN
282 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 ekin = ekin*factor
290 factor = sqrt(factor)
291 IF (PRESENT(ireg)) THEN
292 DO i = 1, SIZE(part)
293 IF (part(i)%t_region_index == my_ireg) part(i)%v(:) = factor*part(i)%v(:)
294 END DO
295 ELSE
296 DO i = 1, SIZE(part)
297 part(i)%v(:) = factor*part(i)%v(:)
298 END DO
299 IF (PRESENT(vcom)) THEN
300 vcom = factor*vcom
301 END IF
302 END IF
303
304 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 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 NULLIFY (thermal_regions, t_region)
327
328 CALL get_md_env(md_env, thermal_regions=thermal_regions)
329 nfree_done = 0
330 DO ireg = 1, thermal_regions%nregions
331 NULLIFY (t_region)
332 t_region => thermal_regions%thermal_region(ireg)
333 nfree = t_region%npart*3
334 ekin = compute_ekin(part, ireg)
335 temp = t_region%temp_expected
336 CALL rescale_vel(part, simpar, ekin, ireg=ireg, nfree=nfree, temp=temp)
337 nfree_done = nfree_done + nfree
338 ekin = compute_ekin(part, ireg)
339 temp = 2.0_dp*ekin/real(nfree, dp)*kelvin
340 t_region%temperature = temp
341 END DO
342 nfree0 = simpar%nfree - nfree_done
343 IF (nfree0 > 0) THEN
344 ekin = compute_ekin(part, 0)
345 CALL rescale_vel(part, simpar, ekin, ireg=0, nfree=nfree0, temp=simpar%temp_ext)
346 ekin = compute_ekin(part, 0)
347 temp = 2.0_dp*ekin/real(nfree0, dp)*kelvin
348 thermal_regions%temp_reg0 = temp
349 END IF
350 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 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 DO i = 1, SIZE(part)
370 SELECT CASE (is_fixed(i))
371 CASE (use_perd_x)
372 part(i)%v(2) = part(i)%v(2) - vcom(2)
373 part(i)%v(3) = part(i)%v(3) - vcom(3)
374 CASE (use_perd_y)
375 part(i)%v(1) = part(i)%v(1) - vcom(1)
376 part(i)%v(3) = part(i)%v(3) - vcom(3)
377 CASE (use_perd_z)
378 part(i)%v(1) = part(i)%v(1) - vcom(1)
379 part(i)%v(2) = part(i)%v(2) - vcom(2)
380 CASE (use_perd_xy)
381 part(i)%v(3) = part(i)%v(3) - vcom(3)
382 CASE (use_perd_xz)
383 part(i)%v(2) = part(i)%v(2) - vcom(2)
384 CASE (use_perd_yz)
385 part(i)%v(1) = part(i)%v(1) - vcom(1)
386 CASE (use_perd_none)
387 part(i)%v(:) = part(i)%v(:) - vcom(:)
388 END SELECT
389 END DO
390 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 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 NULLIFY (atomic_kind)
416 mang(:) = 0.0_dp
417 iner(:, :) = 0.0_dp
418 DO i = 1, SIZE(part)
419 ! compute angular momentum and inertia tensor
420 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 r(:) = part(i)%r(:) - rcom(:)
423 atomic_kind => part(i)%atomic_kind
424 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
425 mang(1) = mang(1) + mass*(r(2)*part(i)%v(3) - r(3)*part(i)%v(2))
426 mang(2) = mang(2) + mass*(r(3)*part(i)%v(1) - r(1)*part(i)%v(3))
427 mang(3) = mang(3) + mass*(r(1)*part(i)%v(2) - r(2)*part(i)%v(1))
428
429 iner(1, 1) = iner(1, 1) + mass*(r(2)*r(2) + r(3)*r(3))
430 iner(2, 2) = iner(2, 2) + mass*(r(3)*r(3) + r(1)*r(1))
431 iner(3, 3) = iner(3, 3) + mass*(r(1)*r(1) + r(2)*r(2))
432
433 iner(1, 2) = iner(1, 2) - mass*r(1)*r(2)
434 iner(2, 3) = iner(2, 3) - mass*r(2)*r(3)
435 iner(3, 1) = iner(3, 1) - mass*r(3)*r(1)
436 END SELECT
437 END DO
438 iner(2, 1) = iner(1, 2)
439 iner(3, 2) = iner(2, 3)
440 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 CALL diamat_all(iner, evals)
446
447 vang(:) = 0.0_dp
448 DO i = 1, 3
449 IF (evals(i) > 0.0_dp) THEN
450 proj = sum(iner(:, i)*mang)/evals(i)
451 vang(1) = vang(1) + proj*iner(1, i)
452 vang(2) = vang(2) + proj*iner(2, i)
453 vang(3) = vang(3) + proj*iner(3, i)
454 END IF
455 END DO
456
457 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 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 DO i = 1, SIZE(part)
479 r(:) = part(i)%r(:) - rcom(:)
480 SELECT CASE (is_fixed(i))
481 CASE (use_perd_x)
482 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
483 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
484 CASE (use_perd_y)
485 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
486 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
487 CASE (use_perd_z)
488 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
489 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
490 CASE (use_perd_xy)
491 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
492 CASE (use_perd_xz)
493 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
494 CASE (use_perd_yz)
495 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
496 CASE (use_perd_none)
497 part(i)%v(1) = part(i)%v(1) - (vang(2)*r(3) - vang(3)*r(2))
498 part(i)%v(2) = part(i)%v(2) - (vang(3)*r(1) - vang(1)*r(3))
499 part(i)%v(3) = part(i)%v(3) - (vang(1)*r(2) - vang(2)*r(1))
500 END SELECT
501 END DO
502
503 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 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 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
570 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 CALL timeset(routinen, handle)
575
576 ! Initializing parameters
577 natoms = SIZE(part)
578 NULLIFY (atomic_kind, fixd_list, logger, molecule_kind)
579 NULLIFY (molecule_kind_set)
580
581 ! Logging
582 logger => cp_get_default_logger()
583 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 ALLOCATE (is_fixed(natoms))
587
588 is_fixed = use_perd_none
589 molecule_kind_set => molecule_kinds%els
590 DO imolecule_kind = 1, molecule_kinds%n_els
591 molecule_kind => molecule_kind_set(imolecule_kind)
592 CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
593 IF (ASSOCIATED(fixd_list)) THEN
594 DO ifixd = 1, SIZE(fixd_list)
595 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 IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
602 simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
603 mass_tot = 0.0_dp
604 DO i = 1, natoms
605 atomic_kind => part(i)%atomic_kind
606 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
607 mass_tot = mass_tot + mass
608 END DO
609 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 shell_present, shell_part, core_part, force_rescaling, para_env, is_fixed, success)
614 IF (.NOT. success) THEN
615 SELECT CASE (simpar%initialization_method)
616 CASE (md_init_default)
617 CALL generate_velocities(simpar, part, force_env, globenv, md_env, shell_present, &
618 shell_part, core_part, is_fixed, iw)
619 CASE (md_init_vib)
620 CALL force_env_get(force_env=force_env, root_section=root_section)
621 md_section => section_vals_get_subs_vals(root_section, "MOTION%MD")
622 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 is_fixed)
633 ! update restart file for the modified coordinates and velocities
634 CALL update_subsys(subsys_section, force_env, .false., write_binary_restart_file)
635 END SELECT
636 END IF
637
638 IF (iw > 0) THEN
639 WRITE (iw, '(/,T2,A)') &
640 'MD_VEL| '//trim(adjustl(label))
641 ! Recompute vcom, ecom and ekin for IO
642 CALL compute_vcom(part, is_fixed, vcom, ecom)
643 ekin = compute_ekin(part) - ecom
644 IF (simpar%nfree == 0) THEN
645 cpassert(ekin == 0.0_dp)
646 temp = 0.0_dp
647 ELSE
648 temp = 2.0_dp*ekin/real(simpar%nfree, kind=dp)
649 END IF
650 tmp_r1 = cp_unit_from_cp2k(temp, "K")
651 WRITE (iw, '(T2,A,T61,F20.6)') &
652 'MD_VEL| Initial temperature [K]', tmp_r1
653 WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
654 'MD_VEL| COM velocity', vcom(1:3)
655
656 ! compute and log rcom and vang if not periodic
657 CALL force_env_get(force_env, cell=cell)
658 IF (sum(cell%perd(1:3)) == 0) THEN
659 CALL compute_rcom(part, is_fixed, rcom)
660 CALL compute_vang(part, is_fixed, rcom, vang)
661 WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
662 'MD_VEL| COM position', rcom(1:3)
663 WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
664 'MD_VEL| Angular velocity', vang(1:3)
665 END IF
666 END IF
667
668 DEALLOCATE (is_fixed)
669 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
670 CALL timestop(handle)
671
672 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 SUBROUTINE read_input_velocities(simpar, part, force_env, md_env, subsys_section, &
691 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 REAL(KIND=dp), DIMENSION(:), POINTER :: vel
710 TYPE(atomic_kind_type), POINTER :: atomic_kind
711 TYPE(cell_type), POINTER :: cell
712 TYPE(cp_sll_val_type), POINTER :: atom_list, core_list, shell_list
713 TYPE(section_vals_type), POINTER :: atomvel_section, corevel_section, &
714 shellvel_section
715 TYPE(shell_kind_type), POINTER :: shell
716 TYPE(thermal_regions_type), POINTER :: thermal_regions
717 TYPE(val_type), POINTER :: val
718
719! Initializing parameters
720
721 success = .false.
722 natoms = SIZE(part)
723 atomvel_read = .false.
724 corevel_read = .false.
725 shellvel_read = .false.
726 NULLIFY (vel, atomic_kind, atom_list, core_list, shell_list)
727 NULLIFY (atomvel_section, shellvel_section, corevel_section)
728 NULLIFY (cell, shell, thermal_regions, val)
729 CALL force_env_get(force_env, cell=cell)
730
731 ! Core-Shell Model
732 nshell = 0
733 IF (shell_present) THEN
734 cpassert(ASSOCIATED(core_part))
735 cpassert(ASSOCIATED(shell_part))
736 nshell = SIZE(shell_part)
737 END IF
738
739 atomvel_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
740 shellvel_section => section_vals_get_subs_vals(subsys_section, "SHELL_VELOCITY")
741 corevel_section => section_vals_get_subs_vals(subsys_section, "CORE_VELOCITY")
742
743 ! Read or initialize the particle velocities
744 CALL section_vals_get(atomvel_section, explicit=atomvel_explicit)
745 CALL section_vals_get(shellvel_section, explicit=shellvel_explicit)
746 CALL section_vals_get(corevel_section, explicit=corevel_explicit)
747 cpassert(shellvel_explicit .EQV. corevel_explicit)
748
749 CALL read_binary_velocities("", part, force_env%root_section, para_env, &
750 subsys_section, atomvel_read, cell)
751 CALL read_binary_velocities("SHELL", shell_part, force_env%root_section, para_env, &
752 subsys_section, shellvel_read, cell)
753 CALL read_binary_velocities("CORE", core_part, force_env%root_section, para_env, &
754 subsys_section, corevel_read, cell)
755
756 IF (.NOT. (atomvel_explicit .OR. atomvel_read)) RETURN
757 success = .true.
758
759 IF (.NOT. atomvel_read) THEN
760 ! Read the atom velocities if explicitly given in the input file
761 CALL section_vals_list_get(atomvel_section, "_DEFAULT_KEYWORD_", list=atom_list)
762 DO i = 1, natoms
763 is_ok = cp_sll_val_next(atom_list, val)
764 CALL val_get(val, r_vals=vel)
765 part(i)%v = vel
766 CALL cell_transform_input_cartesian(cell, part(i)%v)
767 END DO
768 END IF
769 DO i = 1, natoms
770 SELECT CASE (is_fixed(i))
771 CASE (use_perd_x)
772 part(i)%v(1) = 0.0_dp
773 CASE (use_perd_y)
774 part(i)%v(2) = 0.0_dp
775 CASE (use_perd_z)
776 part(i)%v(3) = 0.0_dp
777 CASE (use_perd_xy)
778 part(i)%v(1) = 0.0_dp
779 part(i)%v(2) = 0.0_dp
780 CASE (use_perd_xz)
781 part(i)%v(1) = 0.0_dp
782 part(i)%v(3) = 0.0_dp
783 CASE (use_perd_yz)
784 part(i)%v(2) = 0.0_dp
785 part(i)%v(3) = 0.0_dp
786 CASE (use_perd_xyz)
787 part(i)%v = 0.0_dp
788 END SELECT
789 END DO
790 IF (shell_present) THEN
791 IF (shellvel_explicit) THEN
792 ! If the atoms positions are given (?) and core and shell velocities are
793 ! present in the input, read the latter.
794 CALL section_vals_list_get(shellvel_section, "_DEFAULT_KEYWORD_", list=shell_list)
795 CALL section_vals_list_get(corevel_section, "_DEFAULT_KEYWORD_", list=core_list)
796 DO i = 1, nshell
797 is_ok = cp_sll_val_next(shell_list, val)
798 CALL val_get(val, r_vals=vel)
799 shell_part(i)%v = vel
800 CALL cell_transform_input_cartesian(cell, shell_part(i)%v)
801 is_ok = cp_sll_val_next(core_list, val)
802 CALL val_get(val, r_vals=vel)
803 core_part(i)%v = vel
804 CALL cell_transform_input_cartesian(cell, core_part(i)%v)
805 END DO
806 ELSE
807 IF (.NOT. (shellvel_read .AND. corevel_read)) THEN
808 ! Otherwise, just copy atom velocties into shell and core velocities.
809 CALL clone_core_shell_vel(part, shell_part, core_part)
810 END IF
811 END IF
812 END IF
813
814 ! compute vcom, ecom and ekin
815 CALL compute_vcom(part, is_fixed, vcom, ecom)
816 ekin = compute_ekin(part) - ecom
817
818 IF (simpar%do_thermal_region) THEN
819 CALL get_md_env(md_env, thermal_regions=thermal_regions)
820 IF (ASSOCIATED(thermal_regions)) THEN
821 rescale_regions = thermal_regions%force_rescaling
822 END IF
823 ELSE
824 rescale_regions = .false.
825 END IF
826 IF (simpar%nfree /= 0 .AND. (force_rescaling .OR. rescale_regions)) THEN
827 IF (simpar%do_thermal_region) THEN
828 CALL rescale_vel_region(part, md_env, simpar)
829 ELSE
830 CALL rescale_vel(part, simpar, ekin, vcom=vcom)
831 END IF
832
833 ! After rescaling, the core and shell velocities must also adapt.
834 DO i = 1, natoms
835 shell_index = part(i)%shell_index
836 IF (shell_present .AND. shell_index /= 0) THEN
837 atomic_kind => part(i)%atomic_kind
838 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, shell=shell)
839 fac_masss = shell%mass_shell/mass
840 fac_massc = shell%mass_core/mass
841 vs = shell_part(shell_index)%v
842 vc = core_part(shell_index)%v
843
844 shell_part(shell_index)%v(1) = part(i)%v(1) + fac_massc*(vs(1) - vc(1))
845 shell_part(shell_index)%v(2) = part(i)%v(2) + fac_massc*(vs(2) - vc(2))
846 shell_part(shell_index)%v(3) = part(i)%v(3) + fac_massc*(vs(3) - vc(3))
847 core_part(shell_index)%v(1) = part(i)%v(1) + fac_masss*(vc(1) - vs(1))
848 core_part(shell_index)%v(2) = part(i)%v(2) + fac_masss*(vc(2) - vs(2))
849 core_part(shell_index)%v(3) = part(i)%v(3) + fac_masss*(vc(3) - vs(3))
850 END IF
851 END DO
852 END IF
853 END SUBROUTINE read_input_velocities
854
855! **************************************************************************************************
856!> \brief Initializing velocities AND positions randomly on all processors, based on vibrational
857!> modes of the system, so that the starting coordinates are already very close to
858!> canonical ensumble corresponding to temperature of a head bath.
859!> \param simpar : MD simulation parameters
860!> \param particles : global array of all particles
861!> \param md_section : MD input subsection
862!> \param vib_section : vibrational analysis input section
863!> \param force_env : force environment container
864!> \param global_env : global environment container
865!> \param shell_present : if core-shell model is used
866!> \param shell_particles : global array of all shell particles in shell model
867!> \param core_particles : global array of all core particles in shell model
868!> \param is_fixed : array of size of total number of atoms, that determines if any
869!> cartesian components are fixed
870!> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
871! **************************************************************************************************
872 SUBROUTINE generate_coords_vels_vib(simpar, &
873 particles, &
874 md_section, &
875 vib_section, &
876 force_env, &
877 global_env, &
878 shell_present, &
879 shell_particles, &
880 core_particles, &
881 is_fixed)
882 TYPE(simpar_type), POINTER :: simpar
883 TYPE(particle_type), DIMENSION(:), POINTER :: particles
884 TYPE(section_vals_type), POINTER :: md_section, vib_section
885 TYPE(force_env_type), POINTER :: force_env
886 TYPE(global_environment_type), POINTER :: global_env
887 LOGICAL, INTENT(IN) :: shell_present
888 TYPE(particle_type), DIMENSION(:), POINTER :: shell_particles, core_particles
889 INTEGER, DIMENSION(:), INTENT(IN) :: is_fixed
890
891 INTEGER :: dof, fixed_dof, iatom, ii, imode, &
892 my_dof, natoms, shell_index
893 REAL(KIND=dp) :: erand, mass, my_phase, ratio, temperature
894 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, phase, random
895 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dr, eigenvectors
896 TYPE(atomic_kind_type), POINTER :: atomic_kind
897 TYPE(mp_para_env_type), POINTER :: para_env
898 TYPE(rng_stream_type), ALLOCATABLE :: random_stream
899
900 CALL cite_reference(west2006)
901 natoms = SIZE(particles)
902 temperature = simpar%temp_ext
903 my_dof = 3*natoms
904 ALLOCATE (eigenvalues(my_dof))
905 ALLOCATE (eigenvectors(my_dof, my_dof))
906 ALLOCATE (phase(my_dof))
907 ALLOCATE (random(my_dof))
908 ALLOCATE (dr(3, natoms))
909 CALL force_env_get(force_env=force_env, para_env=para_env)
910 ! read vibration modes
911 CALL read_vib_eigs_unformatted(md_section, &
912 vib_section, &
913 para_env, &
914 dof, &
915 eigenvalues, &
916 eigenvectors)
917 IF (my_dof /= dof) THEN
918 CALL cp_abort(__location__, &
919 "number of degrees of freedom in vibrational analysis data "// &
920 "do not match total number of cartesian degrees of freedom")
921 END IF
922 ! read phases
923 CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%PHASE", r_val=my_phase)
924 my_phase = min(1.0_dp, my_phase)
925 ! generate random numbers
926 random_stream = rng_stream_type(name="MD_INIT_VIB", distribution_type=uniform)
927 CALL random_stream%fill(random)
928 IF (my_phase < 0.0_dp) THEN
929 CALL random_stream%fill(phase)
930 ELSE
931 phase = my_phase
932 END IF
933 DEALLOCATE (random_stream)
934
935 ! the first three modes are acoustic with zero frequencies,
936 ! exclude these from considerations
937 my_dof = dof - 3
938 ! randomly selects energy from distribution about kT, all
939 ! energies are scaled so that the sum over vibration modes gives
940 ! exactly my_dof*kT. Note that k = 1.0 in atomic units
941 erand = 0.0_dp
942 DO imode = 4, dof
943 erand = erand - temperature*log(1.0_dp - random(imode))
944 END DO
945 ! need to take into account of fixed constraints too
946 fixed_dof = 0
947 DO iatom = 1, natoms
948 SELECT CASE (is_fixed(iatom))
949 CASE (use_perd_x, use_perd_y, use_perd_z)
950 fixed_dof = fixed_dof + 1
951 CASE (use_perd_xy, use_perd_xz, use_perd_yz)
952 fixed_dof = fixed_dof + 2
953 CASE (use_perd_xyz)
954 fixed_dof = fixed_dof + 3
955 END SELECT
956 END DO
957 my_dof = my_dof - fixed_dof
958 ratio = real(my_dof, kind=dp)*temperature/erand
959 ! update velocities AND positions
960 DO iatom = 1, natoms
961 atomic_kind => particles(iatom)%atomic_kind
962 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
963 SELECT CASE (is_fixed(iatom))
964 CASE (use_perd_x)
965 DO ii = 2, 3
966 dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
967 eigenvectors, random, phase, dof, ratio)
968 particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
969 eigenvectors, random, phase, dof, &
970 ratio)
971 END DO
972 CASE (use_perd_y)
973 DO ii = 1, 3, 2
974 dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
975 eigenvectors, random, phase, dof, ratio)
976 particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
977 eigenvectors, random, phase, dof, &
978 ratio)
979 END DO
980 CASE (use_perd_z)
981 DO ii = 1, 2
982 dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
983 eigenvectors, random, phase, dof, ratio)
984 particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
985 eigenvectors, random, phase, dof, &
986 ratio)
987 END DO
988 CASE (use_perd_xy)
989 dr(3, iatom) = dr_from_vib_data(iatom, 3, mass, temperature, eigenvalues, &
990 eigenvectors, random, phase, dof, ratio)
991 particles(iatom)%v(3) = dv_from_vib_data(iatom, 3, mass, temperature, &
992 eigenvectors, random, phase, dof, &
993 ratio)
994 CASE (use_perd_xz)
995 dr(2, iatom) = dr_from_vib_data(iatom, 2, mass, temperature, eigenvalues, &
996 eigenvectors, random, phase, dof, ratio)
997 particles(iatom)%v(2) = dv_from_vib_data(iatom, 2, mass, temperature, &
998 eigenvectors, random, phase, dof, &
999 ratio)
1000 CASE (use_perd_yz)
1001 dr(1, iatom) = dr_from_vib_data(iatom, 1, mass, temperature, eigenvalues, &
1002 eigenvectors, random, phase, dof, ratio)
1003 particles(iatom)%v(1) = dv_from_vib_data(iatom, 1, mass, temperature, &
1004 eigenvectors, random, phase, dof, &
1005 ratio)
1006 CASE (use_perd_none)
1007 DO ii = 1, 3
1008 dr(ii, iatom) = dr_from_vib_data(iatom, ii, mass, temperature, eigenvalues, &
1009 eigenvectors, random, phase, dof, ratio)
1010 particles(iatom)%v(ii) = dv_from_vib_data(iatom, ii, mass, temperature, &
1011 eigenvectors, random, phase, dof, &
1012 ratio)
1013 END DO
1014 END SELECT
1015 END DO ! iatom
1016 ! free memory
1017 DEALLOCATE (eigenvalues)
1018 DEALLOCATE (eigenvectors)
1019 DEALLOCATE (phase)
1020 DEALLOCATE (random)
1021 ! update particle coordinates
1022 DO iatom = 1, natoms
1023 particles(iatom)%r(:) = particles(iatom)%r(:) + dr(:, iatom)
1024 END DO
1025 ! update core-shell model coordinates
1026 IF (shell_present) THEN
1027 ! particles have moved, and for core-shell model this means
1028 ! the cores and shells must also move by the same amount. The
1029 ! shell positions will then be optimised if needed
1030 shell_index = particles(iatom)%shell_index
1031 IF (shell_index /= 0) THEN
1032 core_particles(shell_index)%r(:) = core_particles(shell_index)%r(:) + &
1033 dr(:, iatom)
1034 shell_particles(shell_index)%r(:) = shell_particles(shell_index)%r(:) + &
1035 dr(:, iatom)
1036 END IF
1037 CALL optimize_shell_core(force_env, &
1038 particles, &
1039 shell_particles, &
1040 core_particles, &
1041 global_env)
1042 END IF
1043 ! cleanup
1044 DEALLOCATE (dr)
1045 END SUBROUTINE generate_coords_vels_vib
1046
1047! **************************************************************************************************
1048!> \brief calculates componbent of initial velocity of an atom from vibreational modes
1049!> \param iatom : global atomic index
1050!> \param icart : cartesian index (1, 2 or 3)
1051!> \param mass : atomic mass
1052!> \param temperature : target temperature of canonical ensemble
1053!> \param eigenvalues : array containing all cartesian vibrational mode eigenvalues (frequencies)
1054!> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
1055!> (displacements)
1056!> \param random : array containing uniform distributed random numbers, must be the size
1057!> of 3*natoms. Numbers must be between 0 and 1
1058!> \param phase : array containing numbers between 0 and 1 that determines for each
1059!> vibration mode the ratio of potential energy vs kinetic energy contribution
1060!> to total energy
1061!> \param dof : total number of degrees of freedom, = 3*natoms
1062!> \param scale : scale to make sure the sum of vibrational modes give the correct energy
1063!> \return : outputs icart-th cartesian component of initial position of atom iatom
1064!> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
1065! **************************************************************************************************
1066 PURE FUNCTION dr_from_vib_data(iatom, &
1067 icart, &
1068 mass, &
1069 temperature, &
1070 eigenvalues, &
1071 eigenvectors, &
1072 random, &
1073 phase, &
1074 dof, &
1075 scale) &
1076 result(res)
1077 INTEGER, INTENT(IN) :: iatom, icart
1078 REAL(KIND=dp), INTENT(IN) :: mass, temperature
1079 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
1080 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenvectors
1081 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: random, phase
1082 INTEGER, INTENT(IN) :: dof
1083 REAL(KIND=dp), INTENT(IN) :: scale
1084 REAL(KIND=dp) :: res
1085
1086 INTEGER :: imode, ind
1087
1088 res = 0.0_dp
1089 ! assuming the eigenvalues are sorted in ascending order, the
1090 ! first three modes are acoustic with zero frequencies. These are
1091 ! excluded from considerations, and should have been reflected in
1092 ! the calculation of scale outside this function
1093 IF (mass > 0.0_dp) THEN
1094 ! eigenvector rows assumed to be grouped in atomic blocks
1095 ind = (iatom - 1)*3 + icart
1096 DO imode = 4, dof
1097 res = res + &
1098 sqrt(-2.0_dp*scale*temperature*log(1 - random(imode))/mass)/ &
1099 eigenvalues(imode)* &
1100 eigenvectors(ind, imode)* &
1101 cos(2.0_dp*pi*phase(imode))
1102 END DO
1103 END IF
1104 END FUNCTION dr_from_vib_data
1105
1106! **************************************************************************************************
1107!> \brief calculates componbent of initial velocity of an atom from vibreational modes
1108!> \param iatom : global atomic index
1109!> \param icart : cartesian index (1, 2 or 3)
1110!> \param mass : atomic mass
1111!> \param temperature : target temperature of canonical ensemble
1112!> \param eigenvectors : array containing all corresponding vibrational mode eigenvectors
1113!> (displacements)
1114!> \param random : array containing uniform distributed random numbers, must be the size
1115!> of 3*natoms. Numbers must be between 0 and 1
1116!> \param phase : array containing numbers between 0 and 1 that determines for each
1117!> vibration mode the ratio of potential energy vs kinetic energy contribution
1118!> to total energy
1119!> \param dof : total number of degrees of freedom, = 3*natoms
1120!> \param scale : scale to make sure the sum of vibrational modes give the correct energy
1121!> \return : outputs icart-th cartesian component of initial velocity of atom iatom
1122!> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
1123! **************************************************************************************************
1124 PURE FUNCTION dv_from_vib_data(iatom, &
1125 icart, &
1126 mass, &
1127 temperature, &
1128 eigenvectors, &
1129 random, &
1130 phase, &
1131 dof, &
1132 scale) &
1133 result(res)
1134 INTEGER, INTENT(IN) :: iatom, icart
1135 REAL(KIND=dp), INTENT(IN) :: mass, temperature
1136 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eigenvectors
1137 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: random, phase
1138 INTEGER, INTENT(IN) :: dof
1139 REAL(KIND=dp), INTENT(IN) :: scale
1140 REAL(KIND=dp) :: res
1141
1142 INTEGER :: imode, ind
1143
1144 res = 0.0_dp
1145 ! assuming the eigenvalues are sorted in ascending order, the
1146 ! first three modes are acoustic with zero frequencies. These are
1147 ! excluded from considerations, and should have been reflected in
1148 ! the calculation of scale outside this function
1149 IF (mass > 0.0_dp) THEN
1150 ! eigenvector rows assumed to be grouped in atomic blocks
1151 ind = (iatom - 1)*3 + icart
1152 DO imode = 4, dof
1153 res = res - &
1154 sqrt(-2.0_dp*scale*temperature*log(1 - random(imode))/mass)* &
1155 eigenvectors(ind, imode)* &
1156 sin(2.0_dp*pi*phase(imode))
1157 END DO
1158 END IF
1159 END FUNCTION dv_from_vib_data
1160
1161! **************************************************************************************************
1162!> \brief Initializing velocities deterministically on all processors, if not given in input
1163!> \param simpar ...
1164!> \param part ...
1165!> \param force_env ...
1166!> \param globenv ...
1167!> \param md_env ...
1168!> \param shell_present ...
1169!> \param shell_part ...
1170!> \param core_part ...
1171!> \param is_fixed ...
1172!> \param iw ...
1173!> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1174! **************************************************************************************************
1175 SUBROUTINE generate_velocities(simpar, part, force_env, globenv, md_env, &
1176 shell_present, shell_part, core_part, is_fixed, iw)
1177 TYPE(simpar_type), POINTER :: simpar
1178 TYPE(particle_type), DIMENSION(:), POINTER :: part
1179 TYPE(force_env_type), POINTER :: force_env
1180 TYPE(global_environment_type), POINTER :: globenv
1181 TYPE(md_environment_type), POINTER :: md_env
1182 LOGICAL, INTENT(IN) :: shell_present
1183 TYPE(particle_type), DIMENSION(:), POINTER :: shell_part, core_part
1184 INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1185 INTEGER, INTENT(IN) :: iw
1186
1187 INTEGER :: i, natoms
1188 REAL(KIND=dp) :: mass
1189 TYPE(atomic_kind_type), POINTER :: atomic_kind
1190
1191 NULLIFY (atomic_kind)
1192 natoms = SIZE(part)
1193
1194 DO i = 1, natoms
1195 atomic_kind => part(i)%atomic_kind
1196 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1197 part(i)%v(1) = 0.0_dp
1198 part(i)%v(2) = 0.0_dp
1199 part(i)%v(3) = 0.0_dp
1200 IF (mass /= 0.0) THEN
1201 SELECT CASE (is_fixed(i))
1202 CASE (use_perd_x)
1203 part(i)%v(2) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1204 part(i)%v(3) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1205 CASE (use_perd_y)
1206 part(i)%v(1) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1207 part(i)%v(3) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1208 CASE (use_perd_z)
1209 part(i)%v(1) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1210 part(i)%v(2) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1211 CASE (use_perd_xy)
1212 part(i)%v(3) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1213 CASE (use_perd_xz)
1214 part(i)%v(2) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1215 CASE (use_perd_yz)
1216 part(i)%v(1) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1217 CASE (use_perd_none)
1218 part(i)%v(1) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1219 part(i)%v(2) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1220 part(i)%v(3) = globenv%gaussian_rng_stream%next()/sqrt(mass)
1221 END SELECT
1222 END IF
1223 END DO
1224
1225 CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1226 CALL soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
1227
1228 ! Initialize the core and the shell velocity. Atom velocities are just
1229 ! copied so that the initial relative core-shell velocity is zero.
1230 IF (shell_present) THEN
1231 CALL optimize_shell_core(force_env, part, shell_part, core_part, globenv)
1232 END IF
1233 END SUBROUTINE generate_velocities
1234
1235! **************************************************************************************************
1236!> \brief Direct velocities along a low-curvature direction in order to
1237!> favors MD trajectories to cross rapidly over small energy barriers
1238!> into neighboring basins.
1239!> \param simpar ...
1240!> \param part ...
1241!> \param force_env ...
1242!> \param md_env ...
1243!> \param is_fixed ...
1244!> \param iw ...
1245!> \author Ole Schuett
1246! **************************************************************************************************
1247 SUBROUTINE soften_velocities(simpar, part, force_env, md_env, is_fixed, iw)
1248 TYPE(simpar_type), POINTER :: simpar
1249 TYPE(particle_type), DIMENSION(:), POINTER :: part
1250 TYPE(force_env_type), POINTER :: force_env
1251 TYPE(md_environment_type), POINTER :: md_env
1252 INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1253 INTEGER, INTENT(IN) :: iw
1254
1255 INTEGER :: i, k
1256 REAL(KIND=dp), DIMENSION(SIZE(part), 3) :: f, f_t, n, x0
1257
1258 IF (simpar%soften_nsteps <= 0) RETURN !nothing todo
1259
1260 IF (any(is_fixed /= use_perd_none)) &
1261 cpabort("Velocitiy softening with constraints is not supported.")
1262
1263 !backup positions
1264 DO i = 1, SIZE(part)
1265 x0(i, :) = part(i)%r
1266 END DO
1267
1268 DO k = 1, simpar%soften_nsteps
1269
1270 !use normalized velocities as displace direction
1271 DO i = 1, SIZE(part)
1272 n(i, :) = part(i)%v
1273 END DO
1274 n = n/sqrt(sum(n**2))
1275
1276 ! displace system temporarly to calculate forces
1277 DO i = 1, SIZE(part)
1278 part(i)%r = part(i)%r + simpar%soften_delta*n(i, :)
1279 END DO
1280 CALL force_env_calc_energy_force(force_env)
1281
1282 ! calculate velocity update direction F_t
1283 DO i = 1, SIZE(part)
1284 f(i, :) = part(i)%f
1285 END DO
1286 f_t = f - n*sum(n*f)
1287
1288 ! restore positions and update velocities
1289 DO i = 1, SIZE(part)
1290 part(i)%r = x0(i, :)
1291 part(i)%v = part(i)%v + simpar%soften_alpha*f_t(i, :)
1292 END DO
1293
1294 CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1295 END DO
1296
1297 IF (iw > 0) THEN
1298 WRITE (iw, "(A,T71, I10)") " Velocities softening Steps: ", simpar%soften_nsteps
1299 WRITE (iw, "(A,T71, E10.3)") " Velocities softening NORM(F_t): ", sqrt(sum(f_t**2))
1300 END IF
1301 END SUBROUTINE soften_velocities
1302
1303! **************************************************************************************************
1304!> \brief Scale velocities according to temperature and remove rigid body motion.
1305!> \param simpar ...
1306!> \param part ...
1307!> \param force_env ...
1308!> \param md_env ...
1309!> \param is_fixed ...
1310!> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1311! **************************************************************************************************
1312 SUBROUTINE normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1313 TYPE(simpar_type), POINTER :: simpar
1314 TYPE(particle_type), DIMENSION(:), POINTER :: part
1315 TYPE(force_env_type), POINTER :: force_env
1316 TYPE(md_environment_type), POINTER :: md_env
1317 INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1318
1319 REAL(KIND=dp) :: ekin
1320 REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vcom
1321 TYPE(cell_type), POINTER :: cell
1322
1323 NULLIFY (cell)
1324
1325 ! Subtract the vcom
1326 CALL compute_vcom(part, is_fixed, vcom)
1327 CALL subtract_vcom(part, is_fixed, vcom)
1328 ! If requested and the system is not periodic, subtract the angular velocity
1329 CALL force_env_get(force_env, cell=cell)
1330 IF (sum(cell%perd(1:3)) == 0 .AND. simpar%angvel_zero) THEN
1331 CALL compute_rcom(part, is_fixed, rcom)
1332 CALL compute_vang(part, is_fixed, rcom, vang)
1333 CALL subtract_vang(part, is_fixed, rcom, vang)
1334 END IF
1335 ! Rescale the velocities
1336 IF (simpar%do_thermal_region) THEN
1337 CALL rescale_vel_region(part, md_env, simpar)
1338 ELSE
1339 ekin = compute_ekin(part)
1340 CALL rescale_vel(part, simpar, ekin)
1341 END IF
1342 END SUBROUTINE normalize_velocities
1343
1344! **************************************************************************************************
1345!> \brief Computes Ekin, VCOM and Temp for particles
1346!> \param subsys ...
1347!> \param md_ener ...
1348!> \param vsubtract ...
1349!> \par History
1350!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1351! **************************************************************************************************
1352 SUBROUTINE reset_vcom(subsys, md_ener, vsubtract)
1353 TYPE(cp_subsys_type), POINTER :: subsys
1354 TYPE(md_ener_type), POINTER :: md_ener
1355 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vsubtract
1356
1357 CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_vcom'
1358
1359 INTEGER :: atom, handle, iatom, ikind, natom, &
1360 shell_index
1361 INTEGER, DIMENSION(:), POINTER :: atom_list
1362 LOGICAL :: is_shell
1363 REAL(KIND=dp) :: ekin_old, imass_c, imass_s, mass, v2
1364 REAL(KIND=dp), DIMENSION(3) :: tmp, v
1365 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1366 TYPE(atomic_kind_type), POINTER :: atomic_kind
1367 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1368 shell_particles
1369 TYPE(shell_kind_type), POINTER :: shell
1370
1371 NULLIFY (particles, atomic_kind, atomic_kinds, atom_list, shell)
1372 CALL timeset(routinen, handle)
1373
1374 CALL cp_subsys_get(subsys, &
1375 atomic_kinds=atomic_kinds, &
1376 particles=particles, &
1377 shell_particles=shell_particles, &
1378 core_particles=core_particles)
1379
1380 ekin_old = md_ener%ekin
1381 ! Possibly subtract a quantity from all velocities
1382 DO ikind = 1, atomic_kinds%n_els
1383 atomic_kind => atomic_kinds%els(ikind)
1384 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, &
1385 natom=natom, mass=mass, shell_active=is_shell, shell=shell)
1386 IF (is_shell) THEN
1387 tmp = 0.5_dp*vsubtract*mass
1388 imass_s = 1.0_dp/shell%mass_shell
1389 imass_c = 1.0_dp/shell%mass_core
1390 DO iatom = 1, natom
1391 atom = atom_list(iatom)
1392 shell_index = particles%els(atom)%shell_index
1393 shell_particles%els(shell_index)%v = shell_particles%els(shell_index)%v - tmp*imass_s
1394 core_particles%els(shell_index)%v = core_particles%els(shell_index)%v - tmp*imass_c
1395 particles%els(atom)%v = particles%els(atom)%v - vsubtract
1396 END DO
1397 ELSE
1398 DO iatom = 1, natom
1399 atom = atom_list(iatom)
1400 particles%els(atom)%v = particles%els(atom)%v - vsubtract
1401 END DO
1402 END IF
1403 END DO
1404 ! Compute Kinetic Energy and COM Velocity
1405 md_ener%vcom = 0.0_dp
1406 md_ener%total_mass = 0.0_dp
1407 md_ener%ekin = 0.0_dp
1408 DO ikind = 1, atomic_kinds%n_els
1409 atomic_kind => atomic_kinds%els(ikind)
1410 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom)
1411 v2 = 0.0_dp
1412 v = 0.0_dp
1413 DO iatom = 1, natom
1414 atom = atom_list(iatom)
1415 v2 = v2 + sum(particles%els(atom)%v**2)
1416 v(1) = v(1) + particles%els(atom)%v(1)
1417 v(2) = v(2) + particles%els(atom)%v(2)
1418 v(3) = v(3) + particles%els(atom)%v(3)
1419 END DO
1420 md_ener%ekin = md_ener%ekin + 0.5_dp*mass*v2
1421 md_ener%vcom(1) = md_ener%vcom(1) + mass*v(1)
1422 md_ener%vcom(2) = md_ener%vcom(2) + mass*v(2)
1423 md_ener%vcom(3) = md_ener%vcom(3) + mass*v(3)
1424 md_ener%total_mass = md_ener%total_mass + real(natom, kind=dp)*mass
1425 END DO
1426 md_ener%vcom = md_ener%vcom/md_ener%total_mass
1427 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1428 IF (md_ener%nfree /= 0) THEN
1429 md_ener%temp_part = 2.0_dp*md_ener%ekin/real(md_ener%nfree, kind=dp)*kelvin
1430 END IF
1431 CALL timestop(handle)
1432
1433 END SUBROUTINE reset_vcom
1434
1435! **************************************************************************************************
1436!> \brief Scale velocities to get the correct temperature
1437!> \param subsys ...
1438!> \param md_ener ...
1439!> \param temp_expected ...
1440!> \param temp_tol ...
1441!> \param iw ...
1442!> \par History
1443!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1444! **************************************************************************************************
1445 SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw)
1446 TYPE(cp_subsys_type), POINTER :: subsys
1447 TYPE(md_ener_type), POINTER :: md_ener
1448 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1449 INTEGER, INTENT(IN) :: iw
1450
1451 REAL(KIND=dp) :: ekin_old, scale, temp_old
1452
1453 IF (abs(temp_expected - md_ener%temp_part/kelvin) > temp_tol) THEN
1454 scale = 0.0_dp
1455 IF (md_ener%temp_part > 0.0_dp) scale = sqrt((temp_expected/md_ener%temp_part)*kelvin)
1456 ekin_old = md_ener%ekin
1457 temp_old = md_ener%temp_part
1458 md_ener%ekin = 0.0_dp
1459 md_ener%temp_part = 0.0_dp
1460 md_ener%vcom = 0.0_dp
1461 md_ener%total_mass = 0.0_dp
1462
1463 CALL scale_velocity_low(subsys, scale, ireg=0, ekin=md_ener%ekin, vcom=md_ener%vcom)
1464 IF (md_ener%nfree /= 0) THEN
1465 md_ener%temp_part = 2.0_dp*md_ener%ekin/real(md_ener%nfree, kind=dp)*kelvin
1466 END IF
1467 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1468 IF (iw > 0) THEN
1469 WRITE (unit=iw, fmt='(/,T2,A)') &
1470 'MD_VEL| Temperature scaled to requested temperature'
1471 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1472 'MD_VEL| Old temperature [K]', temp_old, &
1473 'MD_VEL| New temperature [K]', md_ener%temp_part
1474 END IF
1475 END IF
1476
1477 END SUBROUTINE scale_velocity
1478
1479! **************************************************************************************************
1480!> \brief Scale velocities of set of regions
1481!> \param md_env ...
1482!> \param subsys ...
1483!> \param md_ener ...
1484!> \param simpar ...
1485!> \param iw ...
1486!> \par author MI
1487! **************************************************************************************************
1488 SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1489
1490 TYPE(md_environment_type), POINTER :: md_env
1491 TYPE(cp_subsys_type), POINTER :: subsys
1492 TYPE(md_ener_type), POINTER :: md_ener
1493 TYPE(simpar_type), POINTER :: simpar
1494 INTEGER, INTENT(IN) :: iw
1495
1496 INTEGER :: ireg, nfree, nfree_done, nregions
1497 REAL(KIND=dp) :: ekin, ekin_old, ekin_total_new, fscale, &
1498 vcom(3), vcom_total(3)
1499 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp_new, temp_old
1500 TYPE(particle_list_type), POINTER :: particles
1501 TYPE(particle_type), DIMENSION(:), POINTER :: part
1502 TYPE(thermal_region_type), POINTER :: t_region
1503 TYPE(thermal_regions_type), POINTER :: thermal_regions
1504
1505 NULLIFY (particles, part, thermal_regions, t_region)
1506 CALL cp_subsys_get(subsys, particles=particles)
1507 part => particles%els
1508 CALL get_md_env(md_env, thermal_regions=thermal_regions)
1509
1510 nregions = thermal_regions%nregions
1511 nfree_done = 0
1512 ekin_total_new = 0.0_dp
1513 ekin_old = md_ener%ekin
1514 vcom_total = 0.0_dp
1515 ALLOCATE (temp_new(0:nregions), temp_old(0:nregions))
1516 temp_new = 0.0_dp
1517 temp_old = 0.0_dp
1518 !loop regions
1519 DO ireg = 1, nregions
1520 NULLIFY (t_region)
1521 t_region => thermal_regions%thermal_region(ireg)
1522 nfree = 3*t_region%npart
1523 ekin = compute_ekin(part, ireg)
1524 IF (nfree > 0) t_region%temperature = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1525 temp_old(ireg) = t_region%temperature
1526 IF (t_region%temp_tol > 0.0_dp .AND. &
1527 abs(t_region%temp_expected - t_region%temperature/kelvin) > t_region%temp_tol) THEN
1528 fscale = sqrt((t_region%temp_expected/t_region%temperature)*kelvin)
1529 CALL scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1530 t_region%temperature = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1531 temp_new(ireg) = t_region%temperature
1532 END IF
1533 nfree_done = nfree_done + nfree
1534 ekin_total_new = ekin_total_new + ekin
1535 END DO
1536 nfree = simpar%nfree - nfree_done
1537 ekin = compute_ekin(part, ireg=0)
1538 IF (nfree > 0) thermal_regions%temp_reg0 = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1539 temp_old(0) = thermal_regions%temp_reg0
1540 IF (simpar%temp_tol > 0.0_dp .AND. nfree > 0) THEN
1541 IF (abs(simpar%temp_ext - thermal_regions%temp_reg0/kelvin) > simpar%temp_tol) THEN
1542 fscale = sqrt((simpar%temp_ext/thermal_regions%temp_reg0)*kelvin)
1543 CALL scale_velocity_low(subsys, fscale, 0, ekin, vcom)
1544 thermal_regions%temp_reg0 = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1545 temp_new(0) = thermal_regions%temp_reg0
1546 END IF
1547 END IF
1548 ekin_total_new = ekin_total_new + ekin
1549
1550 md_ener%ekin = ekin_total_new
1551 IF (md_ener%nfree /= 0) THEN
1552 md_ener%temp_part = 2.0_dp*md_ener%ekin/real(md_ener%nfree, kind=dp)*kelvin
1553 END IF
1554 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1555 IF (iw > 0) THEN
1556 DO ireg = 0, nregions
1557 IF (temp_new(ireg) > 0.0_dp) THEN
1558 WRITE (unit=iw, fmt='(/,T2,A,I0,A)') &
1559 'MD_VEL| Temperature region ', ireg, ' scaled to requested temperature'
1560 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1561 'MD_VEL| Old temperature [K]', temp_old(ireg), &
1562 'MD_VEL| New temperature [K]', temp_new(ireg)
1563 END IF
1564 END DO
1565 END IF
1566 DEALLOCATE (temp_new, temp_old)
1567
1568 END SUBROUTINE scale_velocity_region
1569
1570! **************************************************************************************************
1571!> \brief Scale velocities for a specific region
1572!> \param subsys ...
1573!> \param fscale ...
1574!> \param ireg ...
1575!> \param ekin ...
1576!> \param vcom ...
1577!> \par author MI
1578! **************************************************************************************************
1579 SUBROUTINE scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1580
1581 TYPE(cp_subsys_type), POINTER :: subsys
1582 REAL(KIND=dp), INTENT(IN) :: fscale
1583 INTEGER, INTENT(IN) :: ireg
1584 REAL(KIND=dp), INTENT(OUT) :: ekin, vcom(3)
1585
1586 INTEGER :: atom, iatom, ikind, my_ireg, natom, &
1587 shell_index
1588 INTEGER, DIMENSION(:), POINTER :: atom_list
1589 LOGICAL :: is_shell
1590 REAL(KIND=dp) :: imass, mass, tmass, v2
1591 REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1592 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1593 TYPE(atomic_kind_type), POINTER :: atomic_kind
1594 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1595 shell_particles
1596 TYPE(shell_kind_type), POINTER :: shell
1597
1598 NULLIFY (atomic_kinds, particles, shell_particles, core_particles, shell, atom_list)
1599
1600 my_ireg = ireg
1601 ekin = 0.0_dp
1602 tmass = 0.0_dp
1603 vcom = 0.0_dp
1604
1605 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
1606 shell_particles=shell_particles, core_particles=core_particles)
1607
1608 DO ikind = 1, atomic_kinds%n_els
1609 atomic_kind => atomic_kinds%els(ikind)
1610 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, &
1611 natom=natom, shell_active=is_shell, shell=shell)
1612 IF (is_shell) THEN
1613 imass = 1.0_dp/mass
1614 v2 = 0.0_dp
1615 v = 0.0_dp
1616 DO iatom = 1, natom
1617 atom = atom_list(iatom)
1618 !check region
1619 IF (particles%els(atom)%t_region_index /= my_ireg) cycle
1620
1621 particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1622 shell_index = particles%els(atom)%shell_index
1623 vs = shell_particles%els(shell_index)%v
1624 vc = core_particles%els(shell_index)%v
1625 tmp(1) = imass*(vs(1) - vc(1))
1626 tmp(2) = imass*(vs(2) - vc(2))
1627 tmp(3) = imass*(vs(3) - vc(3))
1628
1629 shell_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) + tmp(1)*shell%mass_core
1630 shell_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) + tmp(2)*shell%mass_core
1631 shell_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) + tmp(3)*shell%mass_core
1632
1633 core_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) - tmp(1)*shell%mass_shell
1634 core_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) - tmp(2)*shell%mass_shell
1635 core_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) - tmp(3)*shell%mass_shell
1636
1637 ! kinetic energy and velocity of COM
1638 v2 = v2 + sum(particles%els(atom)%v**2)
1639 v(1) = v(1) + particles%els(atom)%v(1)
1640 v(2) = v(2) + particles%els(atom)%v(2)
1641 v(3) = v(3) + particles%els(atom)%v(3)
1642 tmass = tmass + mass
1643 END DO
1644 ELSE
1645 v2 = 0.0_dp
1646 v = 0.0_dp
1647 DO iatom = 1, natom
1648 atom = atom_list(iatom)
1649 !check region
1650 IF (particles%els(atom)%t_region_index /= my_ireg) cycle
1651
1652 particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1653 ! kinetic energy and velocity of COM
1654 v2 = v2 + sum(particles%els(atom)%v**2)
1655 v(1) = v(1) + particles%els(atom)%v(1)
1656 v(2) = v(2) + particles%els(atom)%v(2)
1657 v(3) = v(3) + particles%els(atom)%v(3)
1658 tmass = tmass + mass
1659 END DO
1660 END IF
1661 ekin = ekin + 0.5_dp*mass*v2
1662 vcom(1) = vcom(1) + mass*v(1)
1663 vcom(2) = vcom(2) + mass*v(2)
1664 vcom(3) = vcom(3) + mass*v(3)
1665
1666 END DO
1667 vcom = vcom/tmass
1668
1669 END SUBROUTINE scale_velocity_low
1670
1671! **************************************************************************************************
1672!> \brief Scale internal motion of CORE-SHELL model to the correct temperature
1673!> \param subsys ...
1674!> \param md_ener ...
1675!> \param temp_expected ...
1676!> \param temp_tol ...
1677!> \param iw ...
1678!> \par History
1679!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1680! **************************************************************************************************
1681 SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw)
1682 TYPE(cp_subsys_type), POINTER :: subsys
1683 TYPE(md_ener_type), POINTER :: md_ener
1684 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1685 INTEGER, INTENT(IN) :: iw
1686
1687 INTEGER :: atom, iatom, ikind, natom, shell_index
1688 INTEGER, DIMENSION(:), POINTER :: atom_list
1689 LOGICAL :: is_shell
1690 REAL(KIND=dp) :: ekin_shell_old, fac_mass, mass, scale, &
1691 temp_shell_old, v2
1692 REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1693 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1694 TYPE(atomic_kind_type), POINTER :: atomic_kind
1695 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1696 shell_particles
1697 TYPE(shell_kind_type), POINTER :: shell
1698
1699 NULLIFY (atom_list, atomic_kinds, atomic_kind, core_particles, particles, shell_particles, shell)
1700 IF (abs(temp_expected - md_ener%temp_shell/kelvin) > temp_tol) THEN
1701 scale = 0.0_dp
1702 IF (md_ener%temp_shell > epsilon(0.0_dp)) scale = sqrt((temp_expected/md_ener%temp_shell)*kelvin)
1703 ekin_shell_old = md_ener%ekin_shell
1704 temp_shell_old = md_ener%temp_shell
1705 md_ener%ekin_shell = 0.0_dp
1706 md_ener%temp_shell = 0.0_dp
1707
1708 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, shell_particles=shell_particles, &
1709 core_particles=core_particles)
1710
1711 DO ikind = 1, atomic_kinds%n_els
1712 atomic_kind => atomic_kinds%els(ikind)
1713 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom, &
1714 shell_active=is_shell, shell=shell)
1715 IF (is_shell) THEN
1716 fac_mass = 1.0_dp/mass
1717 v2 = 0.0_dp
1718 DO iatom = 1, natom
1719 atom = atom_list(iatom)
1720 shell_index = particles%els(atom)%shell_index
1721 vs = shell_particles%els(shell_index)%v
1722 vc = core_particles%els(shell_index)%v
1723 v = particles%els(atom)%v
1724 tmp(1) = fac_mass*(vc(1) - vs(1))
1725 tmp(2) = fac_mass*(vc(2) - vs(2))
1726 tmp(3) = fac_mass*(vc(3) - vs(3))
1727
1728 shell_particles%els(shell_index)%v(1) = v(1) - shell%mass_core*scale*tmp(1)
1729 shell_particles%els(shell_index)%v(2) = v(2) - shell%mass_core*scale*tmp(2)
1730 shell_particles%els(shell_index)%v(3) = v(3) - shell%mass_core*scale*tmp(3)
1731
1732 core_particles%els(shell_index)%v(1) = v(1) + shell%mass_shell*scale*tmp(1)
1733 core_particles%els(shell_index)%v(2) = v(2) + shell%mass_shell*scale*tmp(2)
1734 core_particles%els(shell_index)%v(3) = v(3) + shell%mass_shell*scale*tmp(3)
1735
1736 vs = shell_particles%els(shell_index)%v
1737 vc = core_particles%els(shell_index)%v
1738 tmp(1) = vc(1) - vs(1)
1739 tmp(2) = vc(2) - vs(2)
1740 tmp(3) = vc(3) - vs(3)
1741 v2 = v2 + sum(tmp**2)
1742 END DO
1743 md_ener%ekin_shell = md_ener%ekin_shell + 0.5_dp*shell%mass_core*shell%mass_shell*fac_mass*v2
1744 END IF
1745 END DO
1746 IF (md_ener%nfree_shell > 0) THEN
1747 md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/real(md_ener%nfree_shell, kind=dp)*kelvin
1748 END IF
1749 md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell
1750 IF (iw > 0) THEN
1751 WRITE (unit=iw, fmt='(/,T2,A)') &
1752 'MD_VEL| Temperature of shell internal motion scaled to requested temperature'
1753 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1754 'MD_VEL| Old temperature [K]', temp_shell_old, &
1755 'MD_VEL| New temperature [K]', md_ener%temp_shell
1756 END IF
1757 END IF
1758
1759 END SUBROUTINE scale_velocity_internal
1760
1761! **************************************************************************************************
1762!> \brief Scale barostat velocities to get the desired temperature
1763!> \param md_env ...
1764!> \param md_ener ...
1765!> \param temp_expected ...
1766!> \param temp_tol ...
1767!> \param iw ...
1768!> \par History
1769!> MI 02.2008
1770! **************************************************************************************************
1771 SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw)
1772 TYPE(md_environment_type), POINTER :: md_env
1773 TYPE(md_ener_type), POINTER :: md_ener
1774 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1775 INTEGER, INTENT(IN) :: iw
1776
1777 INTEGER :: i, j, nfree
1778 REAL(KIND=dp) :: ekin_old, scale, temp_old
1779 TYPE(npt_info_type), POINTER :: npt(:, :)
1780 TYPE(simpar_type), POINTER :: simpar
1781
1782 NULLIFY (npt, simpar)
1783 CALL get_md_env(md_env, simpar=simpar, npt=npt)
1784 IF (abs(temp_expected - md_ener%temp_baro/kelvin) > temp_tol) THEN
1785 scale = 0.0_dp
1786 IF (md_ener%temp_baro > 0.0_dp) scale = sqrt((temp_expected/md_ener%temp_baro)*kelvin)
1787 ekin_old = md_ener%baro_kin
1788 temp_old = md_ener%temp_baro
1789 md_ener%baro_kin = 0.0_dp
1790 md_ener%temp_baro = 0.0_dp
1791 IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npe_i_ensemble &
1792 .OR. simpar%ensemble == npt_ia_ensemble) THEN
1793 npt(1, 1)%v = npt(1, 1)%v*scale
1794 md_ener%baro_kin = 0.5_dp*npt(1, 1)%v**2*npt(1, 1)%mass
1795 ELSE IF (simpar%ensemble == npt_f_ensemble .OR. simpar%ensemble == npe_f_ensemble) THEN
1796 md_ener%baro_kin = 0.0_dp
1797 DO i = 1, 3
1798 DO j = 1, 3
1799 npt(i, j)%v = npt(i, j)%v*scale
1800 md_ener%baro_kin = md_ener%baro_kin + 0.5_dp*npt(i, j)%v**2*npt(i, j)%mass
1801 END DO
1802 END DO
1803 END IF
1804 nfree = SIZE(npt, 1)*SIZE(npt, 2)
1805 md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/real(nfree, dp)*kelvin
1806 IF (iw > 0) THEN
1807 WRITE (unit=iw, fmt='(/,T2,A)') &
1808 'MD_VEL| Temperature of barostat motion scaled to requested temperature'
1809 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1810 'MD_VEL| Old temperature [K]', temp_old, &
1811 'MD_VEL| New temperature [K]', md_ener%temp_baro
1812 END IF
1813 END IF
1814
1815 END SUBROUTINE scale_velocity_baro
1816
1817! **************************************************************************************************
1818!> \brief Perform all temperature manipulations during a QS MD run.
1819!> \param simpar ...
1820!> \param md_env ...
1821!> \param md_ener ...
1822!> \param force_env ...
1823!> \param logger ...
1824!> \par History
1825!> Creation (15.09.2003,MK)
1826!> adapted to force_env (05.10.2003,fawzi)
1827!> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1828! **************************************************************************************************
1829 SUBROUTINE temperature_control(simpar, md_env, md_ener, force_env, logger)
1830
1831 TYPE(simpar_type), POINTER :: simpar
1832 TYPE(md_environment_type), POINTER :: md_env
1833 TYPE(md_ener_type), POINTER :: md_ener
1834 TYPE(force_env_type), POINTER :: force_env
1835 TYPE(cp_logger_type), POINTER :: logger
1836
1837 CHARACTER(LEN=*), PARAMETER :: routinen = 'temperature_control'
1838
1839 INTEGER :: handle, iw
1840 TYPE(cp_subsys_type), POINTER :: subsys
1841 TYPE(mp_para_env_type), POINTER :: para_env
1842
1843 CALL timeset(routinen, handle)
1844 NULLIFY (subsys, para_env)
1845 cpassert(ASSOCIATED(simpar))
1846 cpassert(ASSOCIATED(md_ener))
1847 cpassert(ASSOCIATED(force_env))
1848 CALL force_env_get(force_env, subsys=subsys, para_env=para_env)
1849 iw = cp_print_key_unit_nr(logger, force_env%root_section, "MOTION%MD%PRINT%PROGRAM_RUN_INFO", &
1850 extension=".mdLog")
1851
1852 ! Control the particle motion
1853 IF (simpar%do_thermal_region) THEN
1854 CALL scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1855 ELSE
1856 IF (simpar%temp_tol > 0.0_dp) THEN
1857 CALL scale_velocity(subsys, md_ener, simpar%temp_ext, simpar%temp_tol, iw)
1858 END IF
1859 END IF
1860 ! Control the internal core-shell motion
1861 IF (simpar%temp_sh_tol > 0.0_dp) THEN
1862 CALL scale_velocity_internal(subsys, md_ener, simpar%temp_sh_ext, simpar%temp_sh_tol, iw)
1863 END IF
1864 ! Control cell motion
1865 SELECT CASE (simpar%ensemble)
1866 CASE (nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
1867 npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble, npt_ia_ensemble)
1868 IF (simpar%temp_baro_tol > 0.0_dp) THEN
1869 CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw)
1870 END IF
1871 END SELECT
1872
1873 CALL cp_print_key_finished_output(iw, logger, force_env%root_section, &
1874 "MOTION%MD%PRINT%PROGRAM_RUN_INFO")
1875 CALL timestop(handle)
1876 END SUBROUTINE temperature_control
1877
1878! **************************************************************************************************
1879!> \brief Set to 0 the velocity of the COM along MD runs, if required.
1880!> \param md_ener ...
1881!> \param force_env ...
1882!> \param md_section ...
1883!> \param logger ...
1884!> \par History
1885!> Creation (29.04.2007,MI)
1886!> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1887! **************************************************************************************************
1888 SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
1889
1890 TYPE(md_ener_type), POINTER :: md_ener
1891 TYPE(force_env_type), POINTER :: force_env
1892 TYPE(section_vals_type), POINTER :: md_section
1893 TYPE(cp_logger_type), POINTER :: logger
1894
1895 CHARACTER(LEN=*), PARAMETER :: routinen = 'comvel_control'
1896
1897 INTEGER :: handle, iw
1898 LOGICAL :: explicit
1899 REAL(kind=dp) :: comvel_tol, temp_old, vel_com
1900 REAL(kind=dp), DIMENSION(3) :: vcom_old
1901 TYPE(cp_subsys_type), POINTER :: subsys
1902
1903 CALL timeset(routinen, handle)
1904 NULLIFY (subsys)
1905 cpassert(ASSOCIATED(force_env))
1906 CALL force_env_get(force_env, subsys=subsys)
1907
1908 ! Print COMVEL and COM Position
1909 iw = cp_print_key_unit_nr(logger, md_section, "PRINT%CENTER_OF_MASS", extension=".mdLog")
1910 IF (iw > 0) THEN
1911 WRITE (unit=iw, fmt="(/,T2,A)") &
1912 "MD_VEL| Centre of mass motion (COM)"
1913 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
1914 "MD_VEL| VCOM [a.u.]", md_ener%vcom(1:3)
1915 END IF
1916 CALL cp_print_key_finished_output(iw, logger, md_section, "PRINT%CENTER_OF_MASS")
1917
1918 ! If requested rescale COMVEL
1919 CALL section_vals_val_get(md_section, "COMVEL_TOL", explicit=explicit)
1920 IF (explicit) THEN
1921 CALL section_vals_val_get(md_section, "COMVEL_TOL", r_val=comvel_tol)
1922 iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1923 extension=".mdLog")
1924 vel_com = sqrt(md_ener%vcom(1)**2 + md_ener%vcom(2)**2 + md_ener%vcom(3)**2)
1925
1926 ! Subtract the velocity of the COM, if requested
1927 IF (vel_com > comvel_tol) THEN
1928 temp_old = md_ener%temp_part/kelvin
1929 vcom_old = md_ener%vcom
1930 CALL reset_vcom(subsys, md_ener, vsubtract=vcom_old)
1931 CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
1932 IF (iw > 0) THEN
1933 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
1934 "MD_VEL| Old VCOM [a.u.]", vcom_old(1:3)
1935 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
1936 "MD_VEL| New VCOM [a.u.]", md_ener%vcom(1:3)
1937 END IF
1938 END IF
1939 CALL cp_print_key_finished_output(iw, logger, md_section, &
1940 "PRINT%PROGRAM_RUN_INFO")
1941 END IF
1942
1943 CALL timestop(handle)
1944 END SUBROUTINE comvel_control
1945
1946! **************************************************************************************************
1947!> \brief Set to 0 the angular velocity along MD runs, if required.
1948!> \param md_ener ...
1949!> \param force_env ...
1950!> \param md_section ...
1951!> \param logger ...
1952!> \par History
1953!> Creation (10.2009) Teodoro Laino [tlaino]
1954! **************************************************************************************************
1955 SUBROUTINE angvel_control(md_ener, force_env, md_section, logger)
1956
1957 TYPE(md_ener_type), POINTER :: md_ener
1958 TYPE(force_env_type), POINTER :: force_env
1959 TYPE(section_vals_type), POINTER :: md_section
1960 TYPE(cp_logger_type), POINTER :: logger
1961
1962 CHARACTER(LEN=*), PARAMETER :: routinen = 'angvel_control'
1963
1964 INTEGER :: handle, ifixd, imolecule_kind, iw, natoms
1965 INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
1966 LOGICAL :: explicit
1967 REAL(kind=dp) :: angvel_tol, rcom(3), temp_old, vang(3), &
1968 vang_new(3)
1969 TYPE(cell_type), POINTER :: cell
1970 TYPE(cp_subsys_type), POINTER :: subsys
1971 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
1972 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1973 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1974 TYPE(molecule_kind_type), POINTER :: molecule_kind
1975 TYPE(particle_list_type), POINTER :: particles
1976
1977 CALL timeset(routinen, handle)
1978 ! If requested rescale ANGVEL
1979 CALL section_vals_val_get(md_section, "ANGVEL_TOL", explicit=explicit)
1980 IF (explicit) THEN
1981 NULLIFY (subsys, cell)
1982 cpassert(ASSOCIATED(force_env))
1983 CALL force_env_get(force_env, subsys=subsys, cell=cell)
1984
1985 IF (sum(cell%perd(1:3)) == 0) THEN
1986 CALL section_vals_val_get(md_section, "ANGVEL_TOL", r_val=angvel_tol)
1987 iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1988 extension=".mdLog")
1989
1990 CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, &
1991 particles=particles)
1992
1993 natoms = SIZE(particles%els)
1994 ! Build a list of all fixed atoms (if any)
1995 ALLOCATE (is_fixed(natoms))
1996
1997 is_fixed = use_perd_none
1998 molecule_kind_set => molecule_kinds%els
1999 DO imolecule_kind = 1, molecule_kinds%n_els
2000 molecule_kind => molecule_kind_set(imolecule_kind)
2001 CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2002 IF (ASSOCIATED(fixd_list)) THEN
2003 DO ifixd = 1, SIZE(fixd_list)
2004 IF (.NOT. fixd_list(ifixd)%restraint%active) &
2005 is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2006 END DO
2007 END IF
2008 END DO
2009
2010 ! If requested and the system is not periodic, subtract the angular velocity
2011 CALL compute_rcom(particles%els, is_fixed, rcom)
2012 CALL compute_vang(particles%els, is_fixed, rcom, vang)
2013 ! SQRT(DOT_PRODUCT(vang,vang))>angvel_tol
2014 IF (dot_product(vang, vang) > (angvel_tol*angvel_tol)) THEN
2015 CALL subtract_vang(particles%els, is_fixed, rcom, vang)
2016
2017 ! Rescale velocities after removal
2018 temp_old = md_ener%temp_part/kelvin
2019 CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
2020 CALL compute_vang(particles%els, is_fixed, rcom, vang_new)
2021 IF (iw > 0) THEN
2022 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
2023 'MD_VEL| Old VANG [a.u.]', vang(1:3)
2024 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
2025 'MD_VEL| New VANG [a.u.]', vang_new(1:3)
2026 END IF
2027 END IF
2028
2029 DEALLOCATE (is_fixed)
2030
2031 CALL cp_print_key_finished_output(iw, logger, md_section, &
2032 "PRINT%PROGRAM_RUN_INFO")
2033 END IF
2034 END IF
2035
2036 CALL timestop(handle)
2037 END SUBROUTINE angvel_control
2038
2039! **************************************************************************************************
2040!> \brief Initialize Velocities for MD runs
2041!> \param force_env ...
2042!> \param simpar ...
2043!> \param globenv ...
2044!> \param md_env ...
2045!> \param md_section ...
2046!> \param constraint_section ...
2047!> \param write_binary_restart_file ...
2048!> \par History
2049!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
2050! **************************************************************************************************
2051 SUBROUTINE setup_velocities(force_env, simpar, globenv, md_env, md_section, &
2052 constraint_section, write_binary_restart_file)
2053
2054 TYPE(force_env_type), POINTER :: force_env
2055 TYPE(simpar_type), POINTER :: simpar
2056 TYPE(global_environment_type), POINTER :: globenv
2057 TYPE(md_environment_type), POINTER :: md_env
2058 TYPE(section_vals_type), POINTER :: md_section, constraint_section
2059 LOGICAL, INTENT(IN) :: write_binary_restart_file
2060
2061 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_velocities'
2062
2063 INTEGER :: handle, nconstraint, nconstraint_fixd
2064 LOGICAL :: apply_cns0, shell_adiabatic, &
2065 shell_present
2066 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
2067 TYPE(cell_type), POINTER :: cell
2068 TYPE(cp_subsys_type), POINTER :: subsys
2069 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2070 TYPE(mp_para_env_type), POINTER :: para_env
2071 TYPE(particle_list_type), POINTER :: core_particles, particles, &
2072 shell_particles
2073 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
2074 shell_particle_set
2075 TYPE(section_vals_type), POINTER :: force_env_section, print_section, &
2076 subsys_section
2077
2078 CALL timeset(routinen, handle)
2079
2080 NULLIFY (atomic_kinds, cell, para_env, subsys, molecule_kinds, core_particles, particles)
2081 NULLIFY (shell_particles, core_particle_set, particle_set, shell_particle_set)
2082 NULLIFY (force_env_section, print_section, subsys_section)
2083
2084 print_section => section_vals_get_subs_vals(md_section, "PRINT")
2085 apply_cns0 = .false.
2086 IF (simpar%constraint) THEN
2087 CALL section_vals_val_get(constraint_section, "CONSTRAINT_INIT", l_val=apply_cns0)
2088 END IF
2089 ! Always initialize velocities and possibly restart them
2090 CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env, &
2091 force_env_section=force_env_section)
2092 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2093
2094 CALL cp_subsys_get(subsys, &
2095 atomic_kinds=atomic_kinds, &
2096 core_particles=core_particles, &
2097 molecule_kinds=molecule_kinds, &
2098 particles=particles, &
2099 shell_particles=shell_particles)
2100
2101 CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, &
2102 shell_present=shell_present, &
2103 shell_adiabatic=shell_adiabatic)
2104
2105 NULLIFY (core_particle_set)
2106 NULLIFY (particle_set)
2107 NULLIFY (shell_particle_set)
2108 particle_set => particles%els
2109
2110 IF (shell_present .AND. shell_adiabatic) THEN
2111 ! Constraints are not yet implemented for core-shell models generally
2112 CALL get_molecule_kind_set(molecule_kind_set=molecule_kinds%els, &
2113 nconstraint=nconstraint, &
2114 nconstraint_fixd=nconstraint_fixd)
2115 IF (nconstraint - nconstraint_fixd /= 0) &
2116 cpabort("Only the fixed atom constraint is implemented for core-shell models")
2117!MK CPPostcondition(.NOT.simpar%constraint,cp_failure_level,routineP,failure)
2118 cpassert(ASSOCIATED(shell_particles))
2119 cpassert(ASSOCIATED(core_particles))
2120 shell_particle_set => shell_particles%els
2121 core_particle_set => core_particles%els
2122 END IF
2123
2124 CALL initialize_velocities(simpar, &
2125 particle_set, &
2126 molecule_kinds=molecule_kinds, &
2127 force_env=force_env, &
2128 globenv=globenv, &
2129 md_env=md_env, &
2130 label="Velocities initialization", &
2131 print_section=print_section, &
2132 subsys_section=subsys_section, &
2133 shell_present=(shell_present .AND. shell_adiabatic), &
2134 shell_part=shell_particle_set, &
2135 core_part=core_particle_set, &
2136 force_rescaling=.false., &
2137 para_env=para_env, &
2138 write_binary_restart_file=write_binary_restart_file)
2139
2140 ! Apply constraints if required and rescale velocities..
2141 IF (simpar%ensemble /= reftraj_ensemble) THEN
2142 IF (apply_cns0) THEN
2143 CALL force_env_calc_energy_force(force_env, calc_force=.true.)
2144 CALL force_env_shake(force_env, &
2145 shake_tol=simpar%shake_tol, &
2146 log_unit=simpar%info_constraint, &
2147 lagrange_mult=simpar%lagrange_multipliers, &
2148 dump_lm=simpar%dump_lm, &
2149 compold=.true.)
2150 CALL force_env_rattle(force_env, shake_tol=simpar%shake_tol, &
2151 log_unit=simpar%info_constraint, lagrange_mult=simpar%lagrange_multipliers, &
2152 dump_lm=simpar%dump_lm, reset=.true.)
2153 IF (simpar%do_respa) THEN
2154 CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
2155 calc_force=.true.)
2156 CALL force_env_shake(force_env%sub_force_env(1)%force_env, &
2157 shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2158 lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, compold=.true.)
2159 CALL force_env_rattle(force_env%sub_force_env(1)%force_env, &
2160 shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2161 lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, reset=.true.)
2162 END IF
2163 ! Reinitialize velocities rescaling properly after rattle
2164 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2165 CALL update_subsys(subsys_section, force_env, .false., write_binary_restart_file)
2166 CALL initialize_velocities(simpar, &
2167 particle_set, &
2168 molecule_kinds=molecule_kinds, &
2169 force_env=force_env, &
2170 globenv=globenv, &
2171 md_env=md_env, &
2172 label="Re-Initializing velocities after applying constraints", &
2173 print_section=print_section, &
2174 subsys_section=subsys_section, &
2175 shell_present=(shell_present .AND. shell_adiabatic), &
2176 shell_part=shell_particle_set, &
2177 core_part=core_particle_set, &
2178 force_rescaling=.true., &
2179 para_env=para_env, &
2180 write_binary_restart_file=write_binary_restart_file)
2181 END IF
2182 END IF
2183
2184 ! Perform setup for a cascade run
2185 CALL initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2186
2187 CALL timestop(handle)
2188
2189 END SUBROUTINE setup_velocities
2190
2191! **************************************************************************************************
2192!> \brief Perform the initialization for a cascade run
2193!> \param simpar ...
2194!> \param particle_set ...
2195!> \param molecule_kinds ...
2196!> \param md_section ...
2197!> \date 05.02.2012
2198!> \author Matthias Krack (MK)
2199!> \version 1.0
2200! **************************************************************************************************
2201 SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2202
2203 TYPE(simpar_type), POINTER :: simpar
2204 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2205 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2206 TYPE(section_vals_type), POINTER :: md_section
2207
2208 CHARACTER(LEN=*), PARAMETER :: routinen = 'initialize_cascade'
2209
2210 CHARACTER(len=2*default_string_length) :: line
2211 INTEGER :: handle, iatom, ifixd, imolecule_kind, &
2212 iparticle, iw, natom, nparticle
2213 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_index, is_fixed
2214 LOGICAL :: init_cascade, is_ok, no_read_error
2215 REAL(kind=dp) :: ecom, ekin, energy, norm, temp, &
2216 temperature
2217 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: matom, weight
2218 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vatom
2219 REAL(kind=dp), DIMENSION(3) :: vcom
2220 TYPE(atomic_kind_type), POINTER :: atomic_kind
2221 TYPE(cp_logger_type), POINTER :: logger
2222 TYPE(cp_sll_val_type), POINTER :: atom_list
2223 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
2224 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
2225 TYPE(molecule_kind_type), POINTER :: molecule_kind
2226 TYPE(section_vals_type), POINTER :: atom_list_section, cascade_section, &
2227 print_section
2228 TYPE(val_type), POINTER :: val
2229
2230 CALL timeset(routinen, handle)
2231
2232 NULLIFY (atom_list)
2233 NULLIFY (atom_list_section)
2234 NULLIFY (atomic_kind)
2235 NULLIFY (cascade_section)
2236 NULLIFY (fixd_list)
2237 NULLIFY (molecule_kind)
2238 NULLIFY (molecule_kind_set)
2239 NULLIFY (logger)
2240 NULLIFY (val)
2241
2242 logger => cp_get_default_logger()
2243 print_section => section_vals_get_subs_vals(md_section, "PRINT")
2244 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
2245
2246 cascade_section => section_vals_get_subs_vals(md_section, "CASCADE")
2247 CALL section_vals_val_get(cascade_section, "_SECTION_PARAMETERS_", l_val=init_cascade)
2248
2249 nparticle = SIZE(particle_set)
2250
2251 IF (init_cascade) THEN
2252
2253 CALL section_vals_val_get(cascade_section, "ENERGY", r_val=energy)
2254 IF (energy < 0.0_dp) &
2255 cpabort("Error occurred reading &CASCADE section: Negative energy found")
2256
2257 IF (iw > 0) THEN
2258 ekin = cp_unit_from_cp2k(energy, "keV")
2259 WRITE (unit=iw, fmt="(/,T2,A,T61,F20.6)") &
2260 "CASCADE| Energy [keV]", ekin
2261 WRITE (unit=iw, fmt="(T2,A)") &
2262 "CASCADE|"
2263 END IF
2264
2265 ! Read the atomic velocities given in the input file
2266 atom_list_section => section_vals_get_subs_vals(cascade_section, "ATOM_LIST")
2267 CALL section_vals_val_get(atom_list_section, "_DEFAULT_KEYWORD_", n_rep_val=natom)
2268 CALL section_vals_list_get(atom_list_section, "_DEFAULT_KEYWORD_", list=atom_list)
2269 IF (natom <= 0) &
2270 cpabort("Error occurred reading &CASCADE section: No atom list found")
2271
2272 IF (iw > 0) THEN
2273 WRITE (unit=iw, fmt="(T2,A,T11,A,3(11X,A),9X,A)") &
2274 "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "weight"
2275 END IF
2276
2277 ALLOCATE (atom_index(natom))
2278 ALLOCATE (matom(natom))
2279 ALLOCATE (vatom(3, natom))
2280 ALLOCATE (weight(natom))
2281
2282 DO iatom = 1, natom
2283 is_ok = cp_sll_val_next(atom_list, val)
2284 CALL val_get(val, c_val=line)
2285 ! Read atomic index, velocity vector, and weight
2286 no_read_error = .false.
2287 READ (unit=line, fmt=*, err=999) atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2288 no_read_error = .true.
2289999 IF (.NOT. no_read_error) &
2290 cpabort("Error occurred reading &CASCADE section. Last line read <"//trim(line)//">")
2291 IF ((atom_index(iatom) <= 0) .OR. ((atom_index(iatom) > nparticle))) &
2292 cpabort("Error occurred reading &CASCADE section: Invalid atom index found")
2293 IF (weight(iatom) < 0.0_dp) &
2294 cpabort("Error occurred reading &CASCADE section: Negative weight found")
2295 IF (iw > 0) THEN
2296 WRITE (unit=iw, fmt="(T2,A,I10,4(1X,F14.6))") &
2297 "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2298 END IF
2299 END DO
2300
2301 ! Normalise velocities and weights
2302 norm = 0.0_dp
2303 DO iatom = 1, natom
2304 iparticle = atom_index(iatom)
2305 IF (particle_set(iparticle)%shell_index /= 0) THEN
2306 cpwarn("Warning: The primary knock-on atom is a core-shell atom")
2307 END IF
2308 atomic_kind => particle_set(iparticle)%atomic_kind
2309 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=matom(iatom))
2310 norm = norm + matom(iatom)*weight(iatom)
2311 END DO
2312 weight(:) = matom(:)*weight(:)*energy/norm
2313 DO iatom = 1, natom
2314 norm = sqrt(dot_product(vatom(1:3, iatom), vatom(1:3, iatom)))
2315 vatom(1:3, iatom) = vatom(1:3, iatom)/norm
2316 END DO
2317
2318 IF (iw > 0) THEN
2319 WRITE (unit=iw, fmt="(T2,A)") &
2320 "CASCADE|", &
2321 "CASCADE| Normalised velocities and additional kinetic energy [keV]", &
2322 "CASCADE|"
2323 WRITE (unit=iw, fmt="(T2,A,T11,A,3(11X,A),9X,A)") &
2324 "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "E(kin)"
2325 DO iatom = 1, natom
2326 ekin = cp_unit_from_cp2k(weight(iatom), "keV")
2327 WRITE (unit=iw, fmt="(T2,A,I10,4(1X,F14.6))") &
2328 "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), ekin
2329 END DO
2330 END IF
2331
2332 ! Apply velocity modifications
2333 DO iatom = 1, natom
2334 iparticle = atom_index(iatom)
2335 particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + &
2336 sqrt(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3, iatom)
2337 END DO
2338
2339 DEALLOCATE (atom_index)
2340 DEALLOCATE (matom)
2341 DEALLOCATE (vatom)
2342 DEALLOCATE (weight)
2343
2344 IF (iw > 0) THEN
2345 ! Build a list of all fixed atoms (if any)
2346 ALLOCATE (is_fixed(nparticle))
2347 is_fixed = use_perd_none
2348 molecule_kind_set => molecule_kinds%els
2349 DO imolecule_kind = 1, molecule_kinds%n_els
2350 molecule_kind => molecule_kind_set(imolecule_kind)
2351 CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2352 IF (ASSOCIATED(fixd_list)) THEN
2353 DO ifixd = 1, SIZE(fixd_list)
2354 IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2355 END DO
2356 END IF
2357 END DO
2358 ! Compute vcom, ecom and ekin for printout
2359 CALL compute_vcom(particle_set, is_fixed, vcom, ecom)
2360 ekin = compute_ekin(particle_set) - ecom
2361 IF (simpar%nfree == 0) THEN
2362 cpassert(ekin == 0.0_dp)
2363 temp = 0.0_dp
2364 ELSE
2365 temp = 2.0_dp*ekin/real(simpar%nfree, kind=dp)
2366 END IF
2367 temperature = cp_unit_from_cp2k(temp, "K")
2368 WRITE (unit=iw, fmt="(T2,A)") &
2369 "CASCADE|"
2370 WRITE (unit=iw, fmt="(T2,A,T61,F20.6)") &
2371 "CASCADE| Temperature after cascade initialization [K]", temperature
2372 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,ES16.8))") &
2373 "CASCADE| COM velocity", vcom(1:3)
2374!MK ! compute and log rcom and vang if not periodic
2375!MK CALL force_env_get(force_env,cell=cell)
2376!MK IF (SUM(cell%perd(1:3)) == 0) THEN
2377!MK CALL compute_rcom(particle_set,is_fixed,rcom)
2378!MK CALL compute_vang(particle_set,is_fixed,rcom,vang)
2379!MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3)
2380!MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3)
2381!MK END IF
2382 DEALLOCATE (is_fixed)
2383 END IF
2384
2385 END IF
2386
2387 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
2388
2389 CALL timestop(handle)
2390
2391 END SUBROUTINE initialize_cascade
2392
2393END MODULE md_vel_utils
atom
Definition atom.F:9
atomic_kind_list_types
represent a simple array based list of the given type
Definition atomic_kind_list_types.F:15
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::west2006
integer, save, public west2006
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
cell_types::use_perd_y
integer, parameter, public use_perd_y
Definition cell_types.F:42
cell_types::use_perd_xz
integer, parameter, public use_perd_xz
Definition cell_types.F:42
cell_types::use_perd_x
integer, parameter, public use_perd_x
Definition cell_types.F:42
cell_types::use_perd_z
integer, parameter, public use_perd_z
Definition cell_types.F:42
cell_types::use_perd_yz
integer, parameter, public use_perd_yz
Definition cell_types.F:42
cell_types::use_perd_none
integer, parameter, public use_perd_none
Definition cell_types.F:42
cell_types::use_perd_xy
integer, parameter, public use_perd_xy
Definition cell_types.F:42
cell_types::cell_transform_input_cartesian
subroutine, public cell_transform_input_cartesian(cell, vector)
Transform a Cartesian real-space vector from the user input cell frame into CP2K's canonical internal...
Definition cell_types.F:261
cp_linked_list_input
Definition cp_linked_list_input.F:9
cp_linked_list_input::cp_sll_val_next
logical function, public cp_sll_val_next(iterator, el_att)
returns true if the actual element is valid (i.e. iterator ont at end) moves the iterator to the next...
Definition cp_linked_list_input.F:1475
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1067
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1239
extended_system_types
Lumps all possible extended system variables into one type for easy access and passing.
Definition extended_system_types.F:17
force_env_methods
Interface for the force calculations.
Definition force_env_methods.F:17
force_env_methods::force_env_calc_energy_force
recursive subroutine, public force_env_calc_energy_force(force_env, calc_force, consistent_energies, skip_external_control, eval_energy_forces, require_consistent_energy_force, linres, calc_stress_tensor)
Interface routine for force and energy calculations.
Definition force_env_methods.F:218
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
force_env_utils
Util force_env module.
Definition force_env_utils.F:12
force_env_utils::force_env_shake
subroutine, public force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, pos, vel, compold, reset)
perform shake (enforcing of constraints)
Definition force_env_utils.F:74
force_env_utils::force_env_rattle
subroutine, public force_env_rattle(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, vel, reset)
perform rattle (enforcing of constraints on velocities) This routine can be easily adapted to perform...
Definition force_env_utils.F:242
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::nph_uniaxial_ensemble
integer, parameter, public nph_uniaxial_ensemble
Definition input_constants.F:102
input_constants::npt_i_ensemble
integer, parameter, public npt_i_ensemble
Definition input_constants.F:102
input_constants::nph_uniaxial_damped_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
Definition input_constants.F:102
input_constants::npe_f_ensemble
integer, parameter, public npe_f_ensemble
Definition input_constants.F:102
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition input_constants.F:102
input_constants::md_init_default
integer, parameter, public md_init_default
Definition input_constants.F:124
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition input_constants.F:102
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition input_constants.F:102
input_constants::md_init_vib
integer, parameter, public md_init_vib
Definition input_constants.F:124
input_constants::reftraj_ensemble
integer, parameter, public reftraj_ensemble
Definition input_constants.F:102
input_cp2k_binary_restarts
Routines to read the binary restart file of CP2K.
Definition input_cp2k_binary_restarts.F:15
input_cp2k_binary_restarts::read_binary_velocities
subroutine, public read_binary_velocities(prefix, particle_set, root_section, para_env, subsys_section, binary_file_read, cell)
Read the input section &VELOCITY, &CORE_VELOCITY, or &SHELL_VELOCITY from an external file written in...
Definition input_cp2k_binary_restarts.F:506
input_restart_force_eval
Definition input_restart_force_eval.F:8
input_restart_force_eval::update_subsys
subroutine, public update_subsys(subsys_section, force_env, skip_vel_section, write_binary_restart_file)
Updates the subsys section of the input file.
Definition input_restart_force_eval.F:315
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_list_get
subroutine, public section_vals_list_get(section_vals, keyword_name, i_rep_section, list)
returns the requested list
Definition input_section_types.F:1156
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
input_val_types
a wrapper for basic fortran types.
Definition input_val_types.F:14
input_val_types::val_get
subroutine, public val_get(val, has_l, has_i, has_r, has_lc, has_c, l_val, l_vals, i_val, i_vals, r_val, r_vals, c_val, c_vals, len_c, type_of_var, enum)
returns the stored values
Definition input_val_types.F:334
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:380
md_ener_types
Split md_ener module from md_environment_type.
Definition md_ener_types.F:12
md_environment_types
Definition md_environment_types.F:15
md_environment_types::get_md_env
subroutine, public get_md_env(md_env, itimes, constant, used_time, cell, simpar, npt, force_env, para_env, reftraj, t, init, first_time, fe_env, thermostats, barostat, thermostat_coeff, thermostat_part, thermostat_shell, thermostat_baro, thermostat_fast, thermostat_slow, md_ener, averages, thermal_regions, ehrenfest_md)
get components of MD environment type
Definition md_environment_types.F:195
md_util
Utilities for Molecular Dynamics.
Definition md_util.F:12
md_util::read_vib_eigs_unformatted
subroutine, public read_vib_eigs_unformatted(md_section, vib_section, para_env, dof, eigenvalues, eigenvectors)
read eigenvalues and eigenvectors of Hessian from vibrational analysis results, for use of initialisi...
Definition md_util.F:172
md_vel_utils
Collection of utilities for setting-up and handle velocities in MD runs.
Definition md_vel_utils.F:17
md_vel_utils::setup_velocities
subroutine, public setup_velocities(force_env, simpar, globenv, md_env, md_section, constraint_section, write_binary_restart_file)
Initialize Velocities for MD runs.
Definition md_vel_utils.F:2053
md_vel_utils::angvel_control
subroutine, public angvel_control(md_ener, force_env, md_section, logger)
Set to 0 the angular velocity along MD runs, if required.
Definition md_vel_utils.F:1956
md_vel_utils::temperature_control
subroutine, public temperature_control(simpar, md_env, md_ener, force_env, logger)
Perform all temperature manipulations during a QS MD run.
Definition md_vel_utils.F:1830
md_vel_utils::comvel_control
subroutine, public comvel_control(md_ener, force_env, md_section, logger)
Set to 0 the velocity of the COM along MD runs, if required.
Definition md_vel_utils.F:1889
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:512
molecule_kind_types::get_molecule_kind_set
subroutine, public get_molecule_kind_set(molecule_kind_set, maxatom, natom, nbond, nbend, nub, ntorsion, nimpr, nopbend, nconstraint, nconstraint_fixd, nmolecule, nrestraints)
Get informations about a molecule kind set.
Definition molecule_kind_types.F:660
parallel_rng_types
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
Definition parallel_rng_types.F:51
parallel_rng_types::uniform
integer, parameter, public uniform
Definition parallel_rng_types.F:73
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::kelvin
real(kind=dp), parameter, public kelvin
Definition physcon.F:165
shell_opt
Definition shell_opt.F:12
shell_opt::optimize_shell_core
subroutine, public optimize_shell_core(force_env, particle_set, shell_particle_set, core_particle_set, globenv, tmp, check)
Optimize shell-core positions along an MD run.
Definition shell_opt.F:64
shell_potential_types
Definition shell_potential_types.F:11
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
thermal_region_types
Thermal regions type: to initialize and control the temperature of different regions.
Definition thermal_region_types.F:15
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_linked_list_input::cp_sll_val_type
represent a single linked list that stores pointers to the elements
Definition cp_linked_list_input.F:341
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
extended_system_types::npt_info_type
Definition extended_system_types.F:42
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
input_val_types::val_type
a type to have a wrapper that stores any basic fortran type
Definition input_val_types.F:63
md_ener_types::md_ener_type
Definition md_ener_types.F:22
md_environment_types::md_environment_type
Definition md_environment_types.F:53
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:140
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:151
parallel_rng_types::rng_stream_type
Definition parallel_rng_types.F:157
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30
thermal_region_types::thermal_region_type
Definition thermal_region_types.F:40
thermal_region_types::thermal_regions_type
Definition thermal_region_types.F:31