(git:fc9bb57)
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 ELSE IF (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)) THEN
1261 cpabort("Velocitiy softening with constraints is not supported.")
1262 END IF
1263
1264 !backup positions
1265 DO i = 1, SIZE(part)
1266 x0(i, :) = part(i)%r
1267 END DO
1268
1269 DO k = 1, simpar%soften_nsteps
1270
1271 !use normalized velocities as displace direction
1272 DO i = 1, SIZE(part)
1273 n(i, :) = part(i)%v
1274 END DO
1275 n = n/sqrt(sum(n**2))
1276
1277 ! displace system temporarly to calculate forces
1278 DO i = 1, SIZE(part)
1279 part(i)%r = part(i)%r + simpar%soften_delta*n(i, :)
1280 END DO
1281 CALL force_env_calc_energy_force(force_env)
1282
1283 ! calculate velocity update direction F_t
1284 DO i = 1, SIZE(part)
1285 f(i, :) = part(i)%f
1286 END DO
1287 f_t = f - n*sum(n*f)
1288
1289 ! restore positions and update velocities
1290 DO i = 1, SIZE(part)
1291 part(i)%r = x0(i, :)
1292 part(i)%v = part(i)%v + simpar%soften_alpha*f_t(i, :)
1293 END DO
1294
1295 CALL normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1296 END DO
1297
1298 IF (iw > 0) THEN
1299 WRITE (iw, "(A,T71, I10)") " Velocities softening Steps: ", simpar%soften_nsteps
1300 WRITE (iw, "(A,T71, E10.3)") " Velocities softening NORM(F_t): ", sqrt(sum(f_t**2))
1301 END IF
1302 END SUBROUTINE soften_velocities
1303
1304! **************************************************************************************************
1305!> \brief Scale velocities according to temperature and remove rigid body motion.
1306!> \param simpar ...
1307!> \param part ...
1308!> \param force_env ...
1309!> \param md_env ...
1310!> \param is_fixed ...
1311!> \author CJM,MK,Toon Verstraelen <Toon.Verstraelen@gmail.com>, Ole Schuett
1312! **************************************************************************************************
1313 SUBROUTINE normalize_velocities(simpar, part, force_env, md_env, is_fixed)
1314 TYPE(simpar_type), POINTER :: simpar
1315 TYPE(particle_type), DIMENSION(:), POINTER :: part
1316 TYPE(force_env_type), POINTER :: force_env
1317 TYPE(md_environment_type), POINTER :: md_env
1318 INTEGER, DIMENSION(:), INTENT(INOUT) :: is_fixed
1319
1320 REAL(KIND=dp) :: ekin
1321 REAL(KIND=dp), DIMENSION(3) :: rcom, vang, vcom
1322 TYPE(cell_type), POINTER :: cell
1323
1324 NULLIFY (cell)
1325
1326 ! Subtract the vcom
1327 CALL compute_vcom(part, is_fixed, vcom)
1328 CALL subtract_vcom(part, is_fixed, vcom)
1329 ! If requested and the system is not periodic, subtract the angular velocity
1330 CALL force_env_get(force_env, cell=cell)
1331 IF (sum(cell%perd(1:3)) == 0 .AND. simpar%angvel_zero) THEN
1332 CALL compute_rcom(part, is_fixed, rcom)
1333 CALL compute_vang(part, is_fixed, rcom, vang)
1334 CALL subtract_vang(part, is_fixed, rcom, vang)
1335 END IF
1336 ! Rescale the velocities
1337 IF (simpar%do_thermal_region) THEN
1338 CALL rescale_vel_region(part, md_env, simpar)
1339 ELSE
1340 ekin = compute_ekin(part)
1341 CALL rescale_vel(part, simpar, ekin)
1342 END IF
1343 END SUBROUTINE normalize_velocities
1344
1345! **************************************************************************************************
1346!> \brief Computes Ekin, VCOM and Temp for particles
1347!> \param subsys ...
1348!> \param md_ener ...
1349!> \param vsubtract ...
1350!> \par History
1351!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1352! **************************************************************************************************
1353 SUBROUTINE reset_vcom(subsys, md_ener, vsubtract)
1354 TYPE(cp_subsys_type), POINTER :: subsys
1355 TYPE(md_ener_type), POINTER :: md_ener
1356 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vsubtract
1357
1358 CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_vcom'
1359
1360 INTEGER :: atom, handle, iatom, ikind, natom, &
1361 shell_index
1362 INTEGER, DIMENSION(:), POINTER :: atom_list
1363 LOGICAL :: is_shell
1364 REAL(KIND=dp) :: ekin_old, imass_c, imass_s, mass, v2
1365 REAL(KIND=dp), DIMENSION(3) :: tmp, v
1366 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1367 TYPE(atomic_kind_type), POINTER :: atomic_kind
1368 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1369 shell_particles
1370 TYPE(shell_kind_type), POINTER :: shell
1371
1372 NULLIFY (particles, atomic_kind, atomic_kinds, atom_list, shell)
1373 CALL timeset(routinen, handle)
1374
1375 CALL cp_subsys_get(subsys, &
1376 atomic_kinds=atomic_kinds, &
1377 particles=particles, &
1378 shell_particles=shell_particles, &
1379 core_particles=core_particles)
1380
1381 ekin_old = md_ener%ekin
1382 ! Possibly subtract a quantity from all velocities
1383 DO ikind = 1, atomic_kinds%n_els
1384 atomic_kind => atomic_kinds%els(ikind)
1385 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, &
1386 natom=natom, mass=mass, shell_active=is_shell, shell=shell)
1387 IF (is_shell) THEN
1388 tmp = 0.5_dp*vsubtract*mass
1389 imass_s = 1.0_dp/shell%mass_shell
1390 imass_c = 1.0_dp/shell%mass_core
1391 DO iatom = 1, natom
1392 atom = atom_list(iatom)
1393 shell_index = particles%els(atom)%shell_index
1394 shell_particles%els(shell_index)%v = shell_particles%els(shell_index)%v - tmp*imass_s
1395 core_particles%els(shell_index)%v = core_particles%els(shell_index)%v - tmp*imass_c
1396 particles%els(atom)%v = particles%els(atom)%v - vsubtract
1397 END DO
1398 ELSE
1399 DO iatom = 1, natom
1400 atom = atom_list(iatom)
1401 particles%els(atom)%v = particles%els(atom)%v - vsubtract
1402 END DO
1403 END IF
1404 END DO
1405 ! Compute Kinetic Energy and COM Velocity
1406 md_ener%vcom = 0.0_dp
1407 md_ener%total_mass = 0.0_dp
1408 md_ener%ekin = 0.0_dp
1409 DO ikind = 1, atomic_kinds%n_els
1410 atomic_kind => atomic_kinds%els(ikind)
1411 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom)
1412 v2 = 0.0_dp
1413 v = 0.0_dp
1414 DO iatom = 1, natom
1415 atom = atom_list(iatom)
1416 v2 = v2 + sum(particles%els(atom)%v**2)
1417 v(1) = v(1) + particles%els(atom)%v(1)
1418 v(2) = v(2) + particles%els(atom)%v(2)
1419 v(3) = v(3) + particles%els(atom)%v(3)
1420 END DO
1421 md_ener%ekin = md_ener%ekin + 0.5_dp*mass*v2
1422 md_ener%vcom(1) = md_ener%vcom(1) + mass*v(1)
1423 md_ener%vcom(2) = md_ener%vcom(2) + mass*v(2)
1424 md_ener%vcom(3) = md_ener%vcom(3) + mass*v(3)
1425 md_ener%total_mass = md_ener%total_mass + real(natom, kind=dp)*mass
1426 END DO
1427 md_ener%vcom = md_ener%vcom/md_ener%total_mass
1428 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1429 IF (md_ener%nfree /= 0) THEN
1430 md_ener%temp_part = 2.0_dp*md_ener%ekin/real(md_ener%nfree, kind=dp)*kelvin
1431 END IF
1432 CALL timestop(handle)
1433
1434 END SUBROUTINE reset_vcom
1435
1436! **************************************************************************************************
1437!> \brief Scale velocities to get the correct temperature
1438!> \param subsys ...
1439!> \param md_ener ...
1440!> \param temp_expected ...
1441!> \param temp_tol ...
1442!> \param iw ...
1443!> \par History
1444!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1445! **************************************************************************************************
1446 SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw)
1447 TYPE(cp_subsys_type), POINTER :: subsys
1448 TYPE(md_ener_type), POINTER :: md_ener
1449 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1450 INTEGER, INTENT(IN) :: iw
1451
1452 REAL(KIND=dp) :: ekin_old, scale, temp_old
1453
1454 IF (abs(temp_expected - md_ener%temp_part/kelvin) > temp_tol) THEN
1455 scale = 0.0_dp
1456 IF (md_ener%temp_part > 0.0_dp) scale = sqrt((temp_expected/md_ener%temp_part)*kelvin)
1457 ekin_old = md_ener%ekin
1458 temp_old = md_ener%temp_part
1459 md_ener%ekin = 0.0_dp
1460 md_ener%temp_part = 0.0_dp
1461 md_ener%vcom = 0.0_dp
1462 md_ener%total_mass = 0.0_dp
1463
1464 CALL scale_velocity_low(subsys, scale, ireg=0, ekin=md_ener%ekin, vcom=md_ener%vcom)
1465 IF (md_ener%nfree /= 0) THEN
1466 md_ener%temp_part = 2.0_dp*md_ener%ekin/real(md_ener%nfree, kind=dp)*kelvin
1467 END IF
1468 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1469 IF (iw > 0) THEN
1470 WRITE (unit=iw, fmt='(/,T2,A)') &
1471 'MD_VEL| Temperature scaled to requested temperature'
1472 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1473 'MD_VEL| Old temperature [K]', temp_old, &
1474 'MD_VEL| New temperature [K]', md_ener%temp_part
1475 END IF
1476 END IF
1477
1478 END SUBROUTINE scale_velocity
1479
1480! **************************************************************************************************
1481!> \brief Scale velocities of set of regions
1482!> \param md_env ...
1483!> \param subsys ...
1484!> \param md_ener ...
1485!> \param simpar ...
1486!> \param iw ...
1487!> \par author MI
1488! **************************************************************************************************
1489 SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1490
1491 TYPE(md_environment_type), POINTER :: md_env
1492 TYPE(cp_subsys_type), POINTER :: subsys
1493 TYPE(md_ener_type), POINTER :: md_ener
1494 TYPE(simpar_type), POINTER :: simpar
1495 INTEGER, INTENT(IN) :: iw
1496
1497 INTEGER :: ireg, nfree, nfree_done, nregions
1498 REAL(KIND=dp) :: ekin, ekin_old, ekin_total_new, fscale, &
1499 vcom(3), vcom_total(3)
1500 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: temp_new, temp_old
1501 TYPE(particle_list_type), POINTER :: particles
1502 TYPE(particle_type), DIMENSION(:), POINTER :: part
1503 TYPE(thermal_region_type), POINTER :: t_region
1504 TYPE(thermal_regions_type), POINTER :: thermal_regions
1505
1506 NULLIFY (particles, part, thermal_regions, t_region)
1507 CALL cp_subsys_get(subsys, particles=particles)
1508 part => particles%els
1509 CALL get_md_env(md_env, thermal_regions=thermal_regions)
1510
1511 nregions = thermal_regions%nregions
1512 nfree_done = 0
1513 ekin_total_new = 0.0_dp
1514 ekin_old = md_ener%ekin
1515 vcom_total = 0.0_dp
1516 ALLOCATE (temp_new(0:nregions), temp_old(0:nregions))
1517 temp_new = 0.0_dp
1518 temp_old = 0.0_dp
1519 !loop regions
1520 DO ireg = 1, nregions
1521 NULLIFY (t_region)
1522 t_region => thermal_regions%thermal_region(ireg)
1523 nfree = 3*t_region%npart
1524 ekin = compute_ekin(part, ireg)
1525 IF (nfree > 0) t_region%temperature = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1526 temp_old(ireg) = t_region%temperature
1527 IF (t_region%temp_tol > 0.0_dp .AND. &
1528 abs(t_region%temp_expected - t_region%temperature/kelvin) > t_region%temp_tol) THEN
1529 fscale = sqrt((t_region%temp_expected/t_region%temperature)*kelvin)
1530 CALL scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1531 t_region%temperature = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1532 temp_new(ireg) = t_region%temperature
1533 END IF
1534 nfree_done = nfree_done + nfree
1535 ekin_total_new = ekin_total_new + ekin
1536 END DO
1537 nfree = simpar%nfree - nfree_done
1538 ekin = compute_ekin(part, ireg=0)
1539 IF (nfree > 0) thermal_regions%temp_reg0 = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1540 temp_old(0) = thermal_regions%temp_reg0
1541 IF (simpar%temp_tol > 0.0_dp .AND. nfree > 0) THEN
1542 IF (abs(simpar%temp_ext - thermal_regions%temp_reg0/kelvin) > simpar%temp_tol) THEN
1543 fscale = sqrt((simpar%temp_ext/thermal_regions%temp_reg0)*kelvin)
1544 CALL scale_velocity_low(subsys, fscale, 0, ekin, vcom)
1545 thermal_regions%temp_reg0 = 2.0_dp*ekin/real(nfree, kind=dp)*kelvin
1546 temp_new(0) = thermal_regions%temp_reg0
1547 END IF
1548 END IF
1549 ekin_total_new = ekin_total_new + ekin
1550
1551 md_ener%ekin = ekin_total_new
1552 IF (md_ener%nfree /= 0) THEN
1553 md_ener%temp_part = 2.0_dp*md_ener%ekin/real(md_ener%nfree, kind=dp)*kelvin
1554 END IF
1555 md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin
1556 IF (iw > 0) THEN
1557 DO ireg = 0, nregions
1558 IF (temp_new(ireg) > 0.0_dp) THEN
1559 WRITE (unit=iw, fmt='(/,T2,A,I0,A)') &
1560 'MD_VEL| Temperature region ', ireg, ' scaled to requested temperature'
1561 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1562 'MD_VEL| Old temperature [K]', temp_old(ireg), &
1563 'MD_VEL| New temperature [K]', temp_new(ireg)
1564 END IF
1565 END DO
1566 END IF
1567 DEALLOCATE (temp_new, temp_old)
1568
1569 END SUBROUTINE scale_velocity_region
1570
1571! **************************************************************************************************
1572!> \brief Scale velocities for a specific region
1573!> \param subsys ...
1574!> \param fscale ...
1575!> \param ireg ...
1576!> \param ekin ...
1577!> \param vcom ...
1578!> \par author MI
1579! **************************************************************************************************
1580 SUBROUTINE scale_velocity_low(subsys, fscale, ireg, ekin, vcom)
1581
1582 TYPE(cp_subsys_type), POINTER :: subsys
1583 REAL(KIND=dp), INTENT(IN) :: fscale
1584 INTEGER, INTENT(IN) :: ireg
1585 REAL(KIND=dp), INTENT(OUT) :: ekin, vcom(3)
1586
1587 INTEGER :: atom, iatom, ikind, my_ireg, natom, &
1588 shell_index
1589 INTEGER, DIMENSION(:), POINTER :: atom_list
1590 LOGICAL :: is_shell
1591 REAL(KIND=dp) :: imass, mass, tmass, v2
1592 REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1593 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1594 TYPE(atomic_kind_type), POINTER :: atomic_kind
1595 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1596 shell_particles
1597 TYPE(shell_kind_type), POINTER :: shell
1598
1599 NULLIFY (atomic_kinds, particles, shell_particles, core_particles, shell, atom_list)
1600
1601 my_ireg = ireg
1602 ekin = 0.0_dp
1603 tmass = 0.0_dp
1604 vcom = 0.0_dp
1605
1606 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
1607 shell_particles=shell_particles, core_particles=core_particles)
1608
1609 DO ikind = 1, atomic_kinds%n_els
1610 atomic_kind => atomic_kinds%els(ikind)
1611 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, &
1612 natom=natom, shell_active=is_shell, shell=shell)
1613 IF (is_shell) THEN
1614 imass = 1.0_dp/mass
1615 v2 = 0.0_dp
1616 v = 0.0_dp
1617 DO iatom = 1, natom
1618 atom = atom_list(iatom)
1619 !check region
1620 IF (particles%els(atom)%t_region_index /= my_ireg) cycle
1621
1622 particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1623 shell_index = particles%els(atom)%shell_index
1624 vs = shell_particles%els(shell_index)%v
1625 vc = core_particles%els(shell_index)%v
1626 tmp(1) = imass*(vs(1) - vc(1))
1627 tmp(2) = imass*(vs(2) - vc(2))
1628 tmp(3) = imass*(vs(3) - vc(3))
1629
1630 shell_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) + tmp(1)*shell%mass_core
1631 shell_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) + tmp(2)*shell%mass_core
1632 shell_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) + tmp(3)*shell%mass_core
1633
1634 core_particles%els(shell_index)%v(1) = particles%els(atom)%v(1) - tmp(1)*shell%mass_shell
1635 core_particles%els(shell_index)%v(2) = particles%els(atom)%v(2) - tmp(2)*shell%mass_shell
1636 core_particles%els(shell_index)%v(3) = particles%els(atom)%v(3) - tmp(3)*shell%mass_shell
1637
1638 ! kinetic energy and velocity of COM
1639 v2 = v2 + sum(particles%els(atom)%v**2)
1640 v(1) = v(1) + particles%els(atom)%v(1)
1641 v(2) = v(2) + particles%els(atom)%v(2)
1642 v(3) = v(3) + particles%els(atom)%v(3)
1643 tmass = tmass + mass
1644 END DO
1645 ELSE
1646 v2 = 0.0_dp
1647 v = 0.0_dp
1648 DO iatom = 1, natom
1649 atom = atom_list(iatom)
1650 !check region
1651 IF (particles%els(atom)%t_region_index /= my_ireg) cycle
1652
1653 particles%els(atom)%v(:) = fscale*particles%els(atom)%v
1654 ! kinetic energy and velocity of COM
1655 v2 = v2 + sum(particles%els(atom)%v**2)
1656 v(1) = v(1) + particles%els(atom)%v(1)
1657 v(2) = v(2) + particles%els(atom)%v(2)
1658 v(3) = v(3) + particles%els(atom)%v(3)
1659 tmass = tmass + mass
1660 END DO
1661 END IF
1662 ekin = ekin + 0.5_dp*mass*v2
1663 vcom(1) = vcom(1) + mass*v(1)
1664 vcom(2) = vcom(2) + mass*v(2)
1665 vcom(3) = vcom(3) + mass*v(3)
1666
1667 END DO
1668 vcom = vcom/tmass
1669
1670 END SUBROUTINE scale_velocity_low
1671
1672! **************************************************************************************************
1673!> \brief Scale internal motion of CORE-SHELL model to the correct temperature
1674!> \param subsys ...
1675!> \param md_ener ...
1676!> \param temp_expected ...
1677!> \param temp_tol ...
1678!> \param iw ...
1679!> \par History
1680!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
1681! **************************************************************************************************
1682 SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw)
1683 TYPE(cp_subsys_type), POINTER :: subsys
1684 TYPE(md_ener_type), POINTER :: md_ener
1685 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1686 INTEGER, INTENT(IN) :: iw
1687
1688 INTEGER :: atom, iatom, ikind, natom, shell_index
1689 INTEGER, DIMENSION(:), POINTER :: atom_list
1690 LOGICAL :: is_shell
1691 REAL(KIND=dp) :: ekin_shell_old, fac_mass, mass, scale, &
1692 temp_shell_old, v2
1693 REAL(KIND=dp), DIMENSION(3) :: tmp, v, vc, vs
1694 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
1695 TYPE(atomic_kind_type), POINTER :: atomic_kind
1696 TYPE(particle_list_type), POINTER :: core_particles, particles, &
1697 shell_particles
1698 TYPE(shell_kind_type), POINTER :: shell
1699
1700 NULLIFY (atom_list, atomic_kinds, atomic_kind, core_particles, particles, shell_particles, shell)
1701 IF (abs(temp_expected - md_ener%temp_shell/kelvin) > temp_tol) THEN
1702 scale = 0.0_dp
1703 IF (md_ener%temp_shell > epsilon(0.0_dp)) scale = sqrt((temp_expected/md_ener%temp_shell)*kelvin)
1704 ekin_shell_old = md_ener%ekin_shell
1705 temp_shell_old = md_ener%temp_shell
1706 md_ener%ekin_shell = 0.0_dp
1707 md_ener%temp_shell = 0.0_dp
1708
1709 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, shell_particles=shell_particles, &
1710 core_particles=core_particles)
1711
1712 DO ikind = 1, atomic_kinds%n_els
1713 atomic_kind => atomic_kinds%els(ikind)
1714 CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list, mass=mass, natom=natom, &
1715 shell_active=is_shell, shell=shell)
1716 IF (is_shell) THEN
1717 fac_mass = 1.0_dp/mass
1718 v2 = 0.0_dp
1719 DO iatom = 1, natom
1720 atom = atom_list(iatom)
1721 shell_index = particles%els(atom)%shell_index
1722 vs = shell_particles%els(shell_index)%v
1723 vc = core_particles%els(shell_index)%v
1724 v = particles%els(atom)%v
1725 tmp(1) = fac_mass*(vc(1) - vs(1))
1726 tmp(2) = fac_mass*(vc(2) - vs(2))
1727 tmp(3) = fac_mass*(vc(3) - vs(3))
1728
1729 shell_particles%els(shell_index)%v(1) = v(1) - shell%mass_core*scale*tmp(1)
1730 shell_particles%els(shell_index)%v(2) = v(2) - shell%mass_core*scale*tmp(2)
1731 shell_particles%els(shell_index)%v(3) = v(3) - shell%mass_core*scale*tmp(3)
1732
1733 core_particles%els(shell_index)%v(1) = v(1) + shell%mass_shell*scale*tmp(1)
1734 core_particles%els(shell_index)%v(2) = v(2) + shell%mass_shell*scale*tmp(2)
1735 core_particles%els(shell_index)%v(3) = v(3) + shell%mass_shell*scale*tmp(3)
1736
1737 vs = shell_particles%els(shell_index)%v
1738 vc = core_particles%els(shell_index)%v
1739 tmp(1) = vc(1) - vs(1)
1740 tmp(2) = vc(2) - vs(2)
1741 tmp(3) = vc(3) - vs(3)
1742 v2 = v2 + sum(tmp**2)
1743 END DO
1744 md_ener%ekin_shell = md_ener%ekin_shell + 0.5_dp*shell%mass_core*shell%mass_shell*fac_mass*v2
1745 END IF
1746 END DO
1747 IF (md_ener%nfree_shell > 0) THEN
1748 md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/real(md_ener%nfree_shell, kind=dp)*kelvin
1749 END IF
1750 md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell
1751 IF (iw > 0) THEN
1752 WRITE (unit=iw, fmt='(/,T2,A)') &
1753 'MD_VEL| Temperature of shell internal motion scaled to requested temperature'
1754 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1755 'MD_VEL| Old temperature [K]', temp_shell_old, &
1756 'MD_VEL| New temperature [K]', md_ener%temp_shell
1757 END IF
1758 END IF
1759
1760 END SUBROUTINE scale_velocity_internal
1761
1762! **************************************************************************************************
1763!> \brief Scale barostat velocities to get the desired temperature
1764!> \param md_env ...
1765!> \param md_ener ...
1766!> \param temp_expected ...
1767!> \param temp_tol ...
1768!> \param iw ...
1769!> \par History
1770!> MI 02.2008
1771! **************************************************************************************************
1772 SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw)
1773 TYPE(md_environment_type), POINTER :: md_env
1774 TYPE(md_ener_type), POINTER :: md_ener
1775 REAL(KIND=dp), INTENT(IN) :: temp_expected, temp_tol
1776 INTEGER, INTENT(IN) :: iw
1777
1778 INTEGER :: i, j, nfree
1779 REAL(KIND=dp) :: ekin_old, scale, temp_old
1780 TYPE(npt_info_type), POINTER :: npt(:, :)
1781 TYPE(simpar_type), POINTER :: simpar
1782
1783 NULLIFY (npt, simpar)
1784 CALL get_md_env(md_env, simpar=simpar, npt=npt)
1785 IF (abs(temp_expected - md_ener%temp_baro/kelvin) > temp_tol) THEN
1786 scale = 0.0_dp
1787 IF (md_ener%temp_baro > 0.0_dp) scale = sqrt((temp_expected/md_ener%temp_baro)*kelvin)
1788 ekin_old = md_ener%baro_kin
1789 temp_old = md_ener%temp_baro
1790 md_ener%baro_kin = 0.0_dp
1791 md_ener%temp_baro = 0.0_dp
1792 IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npe_i_ensemble &
1793 .OR. simpar%ensemble == npt_ia_ensemble) THEN
1794 npt(1, 1)%v = npt(1, 1)%v*scale
1795 md_ener%baro_kin = 0.5_dp*npt(1, 1)%v**2*npt(1, 1)%mass
1796 ELSE IF (simpar%ensemble == npt_f_ensemble .OR. simpar%ensemble == npe_f_ensemble) THEN
1797 md_ener%baro_kin = 0.0_dp
1798 DO i = 1, 3
1799 DO j = 1, 3
1800 npt(i, j)%v = npt(i, j)%v*scale
1801 md_ener%baro_kin = md_ener%baro_kin + 0.5_dp*npt(i, j)%v**2*npt(i, j)%mass
1802 END DO
1803 END DO
1804 END IF
1805 nfree = SIZE(npt, 1)*SIZE(npt, 2)
1806 md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/real(nfree, dp)*kelvin
1807 IF (iw > 0) THEN
1808 WRITE (unit=iw, fmt='(/,T2,A)') &
1809 'MD_VEL| Temperature of barostat motion scaled to requested temperature'
1810 WRITE (unit=iw, fmt='(T2,A,T61,F20.6)') &
1811 'MD_VEL| Old temperature [K]', temp_old, &
1812 'MD_VEL| New temperature [K]', md_ener%temp_baro
1813 END IF
1814 END IF
1815
1816 END SUBROUTINE scale_velocity_baro
1817
1818! **************************************************************************************************
1819!> \brief Perform all temperature manipulations during a QS MD run.
1820!> \param simpar ...
1821!> \param md_env ...
1822!> \param md_ener ...
1823!> \param force_env ...
1824!> \param logger ...
1825!> \par History
1826!> Creation (15.09.2003,MK)
1827!> adapted to force_env (05.10.2003,fawzi)
1828!> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1829! **************************************************************************************************
1830 SUBROUTINE temperature_control(simpar, md_env, md_ener, force_env, logger)
1831
1832 TYPE(simpar_type), POINTER :: simpar
1833 TYPE(md_environment_type), POINTER :: md_env
1834 TYPE(md_ener_type), POINTER :: md_ener
1835 TYPE(force_env_type), POINTER :: force_env
1836 TYPE(cp_logger_type), POINTER :: logger
1837
1838 CHARACTER(LEN=*), PARAMETER :: routinen = 'temperature_control'
1839
1840 INTEGER :: handle, iw
1841 TYPE(cp_subsys_type), POINTER :: subsys
1842 TYPE(mp_para_env_type), POINTER :: para_env
1843
1844 CALL timeset(routinen, handle)
1845 NULLIFY (subsys, para_env)
1846 cpassert(ASSOCIATED(simpar))
1847 cpassert(ASSOCIATED(md_ener))
1848 cpassert(ASSOCIATED(force_env))
1849 CALL force_env_get(force_env, subsys=subsys, para_env=para_env)
1850 iw = cp_print_key_unit_nr(logger, force_env%root_section, "MOTION%MD%PRINT%PROGRAM_RUN_INFO", &
1851 extension=".mdLog")
1852
1853 ! Control the particle motion
1854 IF (simpar%do_thermal_region) THEN
1855 CALL scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
1856 ELSE
1857 IF (simpar%temp_tol > 0.0_dp) THEN
1858 CALL scale_velocity(subsys, md_ener, simpar%temp_ext, simpar%temp_tol, iw)
1859 END IF
1860 END IF
1861 ! Control the internal core-shell motion
1862 IF (simpar%temp_sh_tol > 0.0_dp) THEN
1863 CALL scale_velocity_internal(subsys, md_ener, simpar%temp_sh_ext, simpar%temp_sh_tol, iw)
1864 END IF
1865 ! Control cell motion
1866 SELECT CASE (simpar%ensemble)
1867 CASE (nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, &
1868 npt_f_ensemble, npt_i_ensemble, npe_f_ensemble, npe_i_ensemble, npt_ia_ensemble)
1869 IF (simpar%temp_baro_tol > 0.0_dp) THEN
1870 CALL scale_velocity_baro(md_env, md_ener, simpar%temp_baro_ext, simpar%temp_baro_tol, iw)
1871 END IF
1872 END SELECT
1873
1874 CALL cp_print_key_finished_output(iw, logger, force_env%root_section, &
1875 "MOTION%MD%PRINT%PROGRAM_RUN_INFO")
1876 CALL timestop(handle)
1877 END SUBROUTINE temperature_control
1878
1879! **************************************************************************************************
1880!> \brief Set to 0 the velocity of the COM along MD runs, if required.
1881!> \param md_ener ...
1882!> \param force_env ...
1883!> \param md_section ...
1884!> \param logger ...
1885!> \par History
1886!> Creation (29.04.2007,MI)
1887!> Cleaned (09.2007) Teodoro Laino [tlaino] - University of Zurich
1888! **************************************************************************************************
1889 SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
1890
1891 TYPE(md_ener_type), POINTER :: md_ener
1892 TYPE(force_env_type), POINTER :: force_env
1893 TYPE(section_vals_type), POINTER :: md_section
1894 TYPE(cp_logger_type), POINTER :: logger
1895
1896 CHARACTER(LEN=*), PARAMETER :: routinen = 'comvel_control'
1897
1898 INTEGER :: handle, iw
1899 LOGICAL :: explicit
1900 REAL(kind=dp) :: comvel_tol, temp_old, vel_com
1901 REAL(kind=dp), DIMENSION(3) :: vcom_old
1902 TYPE(cp_subsys_type), POINTER :: subsys
1903
1904 CALL timeset(routinen, handle)
1905 NULLIFY (subsys)
1906 cpassert(ASSOCIATED(force_env))
1907 CALL force_env_get(force_env, subsys=subsys)
1908
1909 ! Print COMVEL and COM Position
1910 iw = cp_print_key_unit_nr(logger, md_section, "PRINT%CENTER_OF_MASS", extension=".mdLog")
1911 IF (iw > 0) THEN
1912 WRITE (unit=iw, fmt="(/,T2,A)") &
1913 "MD_VEL| Centre of mass motion (COM)"
1914 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
1915 "MD_VEL| VCOM [a.u.]", md_ener%vcom(1:3)
1916 END IF
1917 CALL cp_print_key_finished_output(iw, logger, md_section, "PRINT%CENTER_OF_MASS")
1918
1919 ! If requested rescale COMVEL
1920 CALL section_vals_val_get(md_section, "COMVEL_TOL", explicit=explicit)
1921 IF (explicit) THEN
1922 CALL section_vals_val_get(md_section, "COMVEL_TOL", r_val=comvel_tol)
1923 iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1924 extension=".mdLog")
1925 vel_com = sqrt(md_ener%vcom(1)**2 + md_ener%vcom(2)**2 + md_ener%vcom(3)**2)
1926
1927 ! Subtract the velocity of the COM, if requested
1928 IF (vel_com > comvel_tol) THEN
1929 temp_old = md_ener%temp_part/kelvin
1930 vcom_old = md_ener%vcom
1931 CALL reset_vcom(subsys, md_ener, vsubtract=vcom_old)
1932 CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
1933 IF (iw > 0) THEN
1934 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
1935 "MD_VEL| Old VCOM [a.u.]", vcom_old(1:3)
1936 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
1937 "MD_VEL| New VCOM [a.u.]", md_ener%vcom(1:3)
1938 END IF
1939 END IF
1940 CALL cp_print_key_finished_output(iw, logger, md_section, &
1941 "PRINT%PROGRAM_RUN_INFO")
1942 END IF
1943
1944 CALL timestop(handle)
1945 END SUBROUTINE comvel_control
1946
1947! **************************************************************************************************
1948!> \brief Set to 0 the angular velocity along MD runs, if required.
1949!> \param md_ener ...
1950!> \param force_env ...
1951!> \param md_section ...
1952!> \param logger ...
1953!> \par History
1954!> Creation (10.2009) Teodoro Laino [tlaino]
1955! **************************************************************************************************
1956 SUBROUTINE angvel_control(md_ener, force_env, md_section, logger)
1957
1958 TYPE(md_ener_type), POINTER :: md_ener
1959 TYPE(force_env_type), POINTER :: force_env
1960 TYPE(section_vals_type), POINTER :: md_section
1961 TYPE(cp_logger_type), POINTER :: logger
1962
1963 CHARACTER(LEN=*), PARAMETER :: routinen = 'angvel_control'
1964
1965 INTEGER :: handle, ifixd, imolecule_kind, iw, natoms
1966 INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
1967 LOGICAL :: explicit
1968 REAL(kind=dp) :: angvel_tol, rcom(3), temp_old, vang(3), &
1969 vang_new(3)
1970 TYPE(cell_type), POINTER :: cell
1971 TYPE(cp_subsys_type), POINTER :: subsys
1972 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
1973 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
1974 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1975 TYPE(molecule_kind_type), POINTER :: molecule_kind
1976 TYPE(particle_list_type), POINTER :: particles
1977
1978 CALL timeset(routinen, handle)
1979 ! If requested rescale ANGVEL
1980 CALL section_vals_val_get(md_section, "ANGVEL_TOL", explicit=explicit)
1981 IF (explicit) THEN
1982 NULLIFY (subsys, cell)
1983 cpassert(ASSOCIATED(force_env))
1984 CALL force_env_get(force_env, subsys=subsys, cell=cell)
1985
1986 IF (sum(cell%perd(1:3)) == 0) THEN
1987 CALL section_vals_val_get(md_section, "ANGVEL_TOL", r_val=angvel_tol)
1988 iw = cp_print_key_unit_nr(logger, md_section, "PRINT%PROGRAM_RUN_INFO", &
1989 extension=".mdLog")
1990
1991 CALL cp_subsys_get(subsys, molecule_kinds=molecule_kinds, &
1992 particles=particles)
1993
1994 natoms = SIZE(particles%els)
1995 ! Build a list of all fixed atoms (if any)
1996 ALLOCATE (is_fixed(natoms))
1997
1998 is_fixed = use_perd_none
1999 molecule_kind_set => molecule_kinds%els
2000 DO imolecule_kind = 1, molecule_kinds%n_els
2001 molecule_kind => molecule_kind_set(imolecule_kind)
2002 CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2003 IF (ASSOCIATED(fixd_list)) THEN
2004 DO ifixd = 1, SIZE(fixd_list)
2005 IF (.NOT. fixd_list(ifixd)%restraint%active) THEN
2006 is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2007 END IF
2008 END DO
2009 END IF
2010 END DO
2011
2012 ! If requested and the system is not periodic, subtract the angular velocity
2013 CALL compute_rcom(particles%els, is_fixed, rcom)
2014 CALL compute_vang(particles%els, is_fixed, rcom, vang)
2015 ! SQRT(DOT_PRODUCT(vang,vang))>angvel_tol
2016 IF (dot_product(vang, vang) > (angvel_tol*angvel_tol)) THEN
2017 CALL subtract_vang(particles%els, is_fixed, rcom, vang)
2018
2019 ! Rescale velocities after removal
2020 temp_old = md_ener%temp_part/kelvin
2021 CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
2022 CALL compute_vang(particles%els, is_fixed, rcom, vang_new)
2023 IF (iw > 0) THEN
2024 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
2025 'MD_VEL| Old VANG [a.u.]', vang(1:3)
2026 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,F16.10))") &
2027 'MD_VEL| New VANG [a.u.]', vang_new(1:3)
2028 END IF
2029 END IF
2030
2031 DEALLOCATE (is_fixed)
2032
2033 CALL cp_print_key_finished_output(iw, logger, md_section, &
2034 "PRINT%PROGRAM_RUN_INFO")
2035 END IF
2036 END IF
2037
2038 CALL timestop(handle)
2039 END SUBROUTINE angvel_control
2040
2041! **************************************************************************************************
2042!> \brief Initialize Velocities for MD runs
2043!> \param force_env ...
2044!> \param simpar ...
2045!> \param globenv ...
2046!> \param md_env ...
2047!> \param md_section ...
2048!> \param constraint_section ...
2049!> \param write_binary_restart_file ...
2050!> \par History
2051!> Teodoro Laino - University of Zurich - 09.2007 [tlaino]
2052! **************************************************************************************************
2053 SUBROUTINE setup_velocities(force_env, simpar, globenv, md_env, md_section, &
2054 constraint_section, write_binary_restart_file)
2055
2056 TYPE(force_env_type), POINTER :: force_env
2057 TYPE(simpar_type), POINTER :: simpar
2058 TYPE(global_environment_type), POINTER :: globenv
2059 TYPE(md_environment_type), POINTER :: md_env
2060 TYPE(section_vals_type), POINTER :: md_section, constraint_section
2061 LOGICAL, INTENT(IN) :: write_binary_restart_file
2062
2063 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_velocities'
2064
2065 INTEGER :: handle, nconstraint, nconstraint_fixd
2066 LOGICAL :: apply_cns0, shell_adiabatic, &
2067 shell_present
2068 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
2069 TYPE(cell_type), POINTER :: cell
2070 TYPE(cp_subsys_type), POINTER :: subsys
2071 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2072 TYPE(mp_para_env_type), POINTER :: para_env
2073 TYPE(particle_list_type), POINTER :: core_particles, particles, &
2074 shell_particles
2075 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
2076 shell_particle_set
2077 TYPE(section_vals_type), POINTER :: force_env_section, print_section, &
2078 subsys_section
2079
2080 CALL timeset(routinen, handle)
2081
2082 NULLIFY (atomic_kinds, cell, para_env, subsys, molecule_kinds, core_particles, particles)
2083 NULLIFY (shell_particles, core_particle_set, particle_set, shell_particle_set)
2084 NULLIFY (force_env_section, print_section, subsys_section)
2085
2086 print_section => section_vals_get_subs_vals(md_section, "PRINT")
2087 apply_cns0 = .false.
2088 IF (simpar%constraint) THEN
2089 CALL section_vals_val_get(constraint_section, "CONSTRAINT_INIT", l_val=apply_cns0)
2090 END IF
2091 ! Always initialize velocities and possibly restart them
2092 CALL force_env_get(force_env, subsys=subsys, cell=cell, para_env=para_env, &
2093 force_env_section=force_env_section)
2094 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2095
2096 CALL cp_subsys_get(subsys, &
2097 atomic_kinds=atomic_kinds, &
2098 core_particles=core_particles, &
2099 molecule_kinds=molecule_kinds, &
2100 particles=particles, &
2101 shell_particles=shell_particles)
2102
2103 CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, &
2104 shell_present=shell_present, &
2105 shell_adiabatic=shell_adiabatic)
2106
2107 NULLIFY (core_particle_set)
2108 NULLIFY (particle_set)
2109 NULLIFY (shell_particle_set)
2110 particle_set => particles%els
2111
2112 IF (shell_present .AND. shell_adiabatic) THEN
2113 ! Constraints are not yet implemented for core-shell models generally
2114 CALL get_molecule_kind_set(molecule_kind_set=molecule_kinds%els, &
2115 nconstraint=nconstraint, &
2116 nconstraint_fixd=nconstraint_fixd)
2117 IF (nconstraint - nconstraint_fixd /= 0) THEN
2118 cpabort("Only the fixed atom constraint is implemented for core-shell models")
2119 END IF
2120!MK CPPostcondition(.NOT.simpar%constraint,cp_failure_level,routineP,failure)
2121 cpassert(ASSOCIATED(shell_particles))
2122 cpassert(ASSOCIATED(core_particles))
2123 shell_particle_set => shell_particles%els
2124 core_particle_set => core_particles%els
2125 END IF
2126
2127 CALL initialize_velocities(simpar, &
2128 particle_set, &
2129 molecule_kinds=molecule_kinds, &
2130 force_env=force_env, &
2131 globenv=globenv, &
2132 md_env=md_env, &
2133 label="Velocities initialization", &
2134 print_section=print_section, &
2135 subsys_section=subsys_section, &
2136 shell_present=(shell_present .AND. shell_adiabatic), &
2137 shell_part=shell_particle_set, &
2138 core_part=core_particle_set, &
2139 force_rescaling=.false., &
2140 para_env=para_env, &
2141 write_binary_restart_file=write_binary_restart_file)
2142
2143 ! Apply constraints if required and rescale velocities..
2144 IF (simpar%ensemble /= reftraj_ensemble) THEN
2145 IF (apply_cns0) THEN
2146 CALL force_env_calc_energy_force(force_env, calc_force=.true.)
2147 CALL force_env_shake(force_env, &
2148 shake_tol=simpar%shake_tol, &
2149 log_unit=simpar%info_constraint, &
2150 lagrange_mult=simpar%lagrange_multipliers, &
2151 dump_lm=simpar%dump_lm, &
2152 compold=.true.)
2153 CALL force_env_rattle(force_env, shake_tol=simpar%shake_tol, &
2154 log_unit=simpar%info_constraint, lagrange_mult=simpar%lagrange_multipliers, &
2155 dump_lm=simpar%dump_lm, reset=.true.)
2156 IF (simpar%do_respa) THEN
2157 CALL force_env_calc_energy_force(force_env%sub_force_env(1)%force_env, &
2158 calc_force=.true.)
2159 CALL force_env_shake(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, compold=.true.)
2162 CALL force_env_rattle(force_env%sub_force_env(1)%force_env, &
2163 shake_tol=simpar%shake_tol, log_unit=simpar%info_constraint, &
2164 lagrange_mult=simpar%lagrange_multipliers, dump_lm=simpar%dump_lm, reset=.true.)
2165 END IF
2166 ! Reinitialize velocities rescaling properly after rattle
2167 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
2168 CALL update_subsys(subsys_section, force_env, .false., write_binary_restart_file)
2169 CALL initialize_velocities(simpar, &
2170 particle_set, &
2171 molecule_kinds=molecule_kinds, &
2172 force_env=force_env, &
2173 globenv=globenv, &
2174 md_env=md_env, &
2175 label="Re-Initializing velocities after applying constraints", &
2176 print_section=print_section, &
2177 subsys_section=subsys_section, &
2178 shell_present=(shell_present .AND. shell_adiabatic), &
2179 shell_part=shell_particle_set, &
2180 core_part=core_particle_set, &
2181 force_rescaling=.true., &
2182 para_env=para_env, &
2183 write_binary_restart_file=write_binary_restart_file)
2184 END IF
2185 END IF
2186
2187 ! Perform setup for a cascade run
2188 CALL initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2189
2190 CALL timestop(handle)
2191
2192 END SUBROUTINE setup_velocities
2193
2194! **************************************************************************************************
2195!> \brief Perform the initialization for a cascade run
2196!> \param simpar ...
2197!> \param particle_set ...
2198!> \param molecule_kinds ...
2199!> \param md_section ...
2200!> \date 05.02.2012
2201!> \author Matthias Krack (MK)
2202!> \version 1.0
2203! **************************************************************************************************
2204 SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
2205
2206 TYPE(simpar_type), POINTER :: simpar
2207 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2208 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
2209 TYPE(section_vals_type), POINTER :: md_section
2210
2211 CHARACTER(LEN=*), PARAMETER :: routinen = 'initialize_cascade'
2212
2213 CHARACTER(len=2*default_string_length) :: line
2214 INTEGER :: handle, iatom, ifixd, imolecule_kind, &
2215 iparticle, iw, natom, nparticle
2216 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_index, is_fixed
2217 LOGICAL :: init_cascade, is_ok, no_read_error
2218 REAL(kind=dp) :: ecom, ekin, energy, norm, temp, &
2219 temperature
2220 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: matom, weight
2221 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vatom
2222 REAL(kind=dp), DIMENSION(3) :: vcom
2223 TYPE(atomic_kind_type), POINTER :: atomic_kind
2224 TYPE(cp_logger_type), POINTER :: logger
2225 TYPE(cp_sll_val_type), POINTER :: atom_list
2226 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
2227 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
2228 TYPE(molecule_kind_type), POINTER :: molecule_kind
2229 TYPE(section_vals_type), POINTER :: atom_list_section, cascade_section, &
2230 print_section
2231 TYPE(val_type), POINTER :: val
2232
2233 CALL timeset(routinen, handle)
2234
2235 NULLIFY (atom_list)
2236 NULLIFY (atom_list_section)
2237 NULLIFY (atomic_kind)
2238 NULLIFY (cascade_section)
2239 NULLIFY (fixd_list)
2240 NULLIFY (molecule_kind)
2241 NULLIFY (molecule_kind_set)
2242 NULLIFY (logger)
2243 NULLIFY (val)
2244
2245 logger => cp_get_default_logger()
2246 print_section => section_vals_get_subs_vals(md_section, "PRINT")
2247 iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
2248
2249 cascade_section => section_vals_get_subs_vals(md_section, "CASCADE")
2250 CALL section_vals_val_get(cascade_section, "_SECTION_PARAMETERS_", l_val=init_cascade)
2251
2252 nparticle = SIZE(particle_set)
2253
2254 IF (init_cascade) THEN
2255
2256 CALL section_vals_val_get(cascade_section, "ENERGY", r_val=energy)
2257 IF (energy < 0.0_dp) THEN
2258 cpabort("Error occurred reading &CASCADE section: Negative energy found")
2259 END IF
2260
2261 IF (iw > 0) THEN
2262 ekin = cp_unit_from_cp2k(energy, "keV")
2263 WRITE (unit=iw, fmt="(/,T2,A,T61,F20.6)") &
2264 "CASCADE| Energy [keV]", ekin
2265 WRITE (unit=iw, fmt="(T2,A)") &
2266 "CASCADE|"
2267 END IF
2268
2269 ! Read the atomic velocities given in the input file
2270 atom_list_section => section_vals_get_subs_vals(cascade_section, "ATOM_LIST")
2271 CALL section_vals_val_get(atom_list_section, "_DEFAULT_KEYWORD_", n_rep_val=natom)
2272 CALL section_vals_list_get(atom_list_section, "_DEFAULT_KEYWORD_", list=atom_list)
2273 IF (natom <= 0) THEN
2274 cpabort("Error occurred reading &CASCADE section: No atom list found")
2275 END IF
2276
2277 IF (iw > 0) THEN
2278 WRITE (unit=iw, fmt="(T2,A,T11,A,3(11X,A),9X,A)") &
2279 "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "weight"
2280 END IF
2281
2282 ALLOCATE (atom_index(natom))
2283 ALLOCATE (matom(natom))
2284 ALLOCATE (vatom(3, natom))
2285 ALLOCATE (weight(natom))
2286
2287 DO iatom = 1, natom
2288 is_ok = cp_sll_val_next(atom_list, val)
2289 CALL val_get(val, c_val=line)
2290 ! Read atomic index, velocity vector, and weight
2291 no_read_error = .false.
2292 READ (unit=line, fmt=*, err=999) atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2293 no_read_error = .true.
2294999 IF (.NOT. no_read_error) THEN
2295 cpabort("Error occurred reading &CASCADE section. Last line read <"//trim(line)//">")
2296 END IF
2297 IF ((atom_index(iatom) <= 0) .OR. ((atom_index(iatom) > nparticle))) THEN
2298 cpabort("Error occurred reading &CASCADE section: Invalid atom index found")
2299 END IF
2300 IF (weight(iatom) < 0.0_dp) THEN
2301 cpabort("Error occurred reading &CASCADE section: Negative weight found")
2302 END IF
2303 IF (iw > 0) THEN
2304 WRITE (unit=iw, fmt="(T2,A,I10,4(1X,F14.6))") &
2305 "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), weight(iatom)
2306 END IF
2307 END DO
2308
2309 ! Normalise velocities and weights
2310 norm = 0.0_dp
2311 DO iatom = 1, natom
2312 iparticle = atom_index(iatom)
2313 IF (particle_set(iparticle)%shell_index /= 0) THEN
2314 cpwarn("Warning: The primary knock-on atom is a core-shell atom")
2315 END IF
2316 atomic_kind => particle_set(iparticle)%atomic_kind
2317 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=matom(iatom))
2318 norm = norm + matom(iatom)*weight(iatom)
2319 END DO
2320 weight(:) = matom(:)*weight(:)*energy/norm
2321 DO iatom = 1, natom
2322 norm = norm2(vatom(1:3, iatom))
2323 vatom(1:3, iatom) = vatom(1:3, iatom)/norm
2324 END DO
2325
2326 IF (iw > 0) THEN
2327 WRITE (unit=iw, fmt="(T2,A)") &
2328 "CASCADE|", &
2329 "CASCADE| Normalised velocities and additional kinetic energy [keV]", &
2330 "CASCADE|"
2331 WRITE (unit=iw, fmt="(T2,A,T11,A,3(11X,A),9X,A)") &
2332 "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "E(kin)"
2333 DO iatom = 1, natom
2334 ekin = cp_unit_from_cp2k(weight(iatom), "keV")
2335 WRITE (unit=iw, fmt="(T2,A,I10,4(1X,F14.6))") &
2336 "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), ekin
2337 END DO
2338 END IF
2339
2340 ! Apply velocity modifications
2341 DO iatom = 1, natom
2342 iparticle = atom_index(iatom)
2343 particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + &
2344 sqrt(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3, iatom)
2345 END DO
2346
2347 DEALLOCATE (atom_index)
2348 DEALLOCATE (matom)
2349 DEALLOCATE (vatom)
2350 DEALLOCATE (weight)
2351
2352 IF (iw > 0) THEN
2353 ! Build a list of all fixed atoms (if any)
2354 ALLOCATE (is_fixed(nparticle))
2355 is_fixed = use_perd_none
2356 molecule_kind_set => molecule_kinds%els
2357 DO imolecule_kind = 1, molecule_kinds%n_els
2358 molecule_kind => molecule_kind_set(imolecule_kind)
2359 CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2360 IF (ASSOCIATED(fixd_list)) THEN
2361 DO ifixd = 1, SIZE(fixd_list)
2362 IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2363 END DO
2364 END IF
2365 END DO
2366 ! Compute vcom, ecom and ekin for printout
2367 CALL compute_vcom(particle_set, is_fixed, vcom, ecom)
2368 ekin = compute_ekin(particle_set) - ecom
2369 IF (simpar%nfree == 0) THEN
2370 cpassert(ekin == 0.0_dp)
2371 temp = 0.0_dp
2372 ELSE
2373 temp = 2.0_dp*ekin/real(simpar%nfree, kind=dp)
2374 END IF
2375 temperature = cp_unit_from_cp2k(temp, "K")
2376 WRITE (unit=iw, fmt="(T2,A)") &
2377 "CASCADE|"
2378 WRITE (unit=iw, fmt="(T2,A,T61,F20.6)") &
2379 "CASCADE| Temperature after cascade initialization [K]", temperature
2380 WRITE (unit=iw, fmt="(T2,A,T30,3(1X,ES16.8))") &
2381 "CASCADE| COM velocity", vcom(1:3)
2382!MK ! compute and log rcom and vang if not periodic
2383!MK CALL force_env_get(force_env,cell=cell)
2384!MK IF (SUM(cell%perd(1:3)) == 0) THEN
2385!MK CALL compute_rcom(particle_set,is_fixed,rcom)
2386!MK CALL compute_vang(particle_set,is_fixed,rcom,vang)
2387!MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3)
2388!MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3)
2389!MK END IF
2390 DEALLOCATE (is_fixed)
2391 END IF
2392
2393 END IF
2394
2395 CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
2396
2397 CALL timestop(handle)
2398
2399 END SUBROUTINE initialize_cascade
2400
2401END 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:284
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:1073
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:1251
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:344
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:105
input_constants::npt_i_ensemble
integer, parameter, public npt_i_ensemble
Definition input_constants.F:105
input_constants::nph_uniaxial_damped_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
Definition input_constants.F:105
input_constants::npe_f_ensemble
integer, parameter, public npe_f_ensemble
Definition input_constants.F:105
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition input_constants.F:105
input_constants::md_init_default
integer, parameter, public md_init_default
Definition input_constants.F:127
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition input_constants.F:105
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition input_constants.F:105
input_constants::md_init_vib
integer, parameter, public md_init_vib
Definition input_constants.F:127
input_constants::reftraj_ensemble
integer, parameter, public reftraj_ensemble
Definition input_constants.F:105
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:1164
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:735
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:708
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:1052
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:381
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:2055
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:1957
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:1831
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:1890
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
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:517
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:665
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:767
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:141
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:152
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