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