(git:b279b6b)
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-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Collection of utilities for setting-up and handle velocities in MD
10 !> runs
11 !> \author CJM
12 !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2008
13 !> reorganization of the original routines/modules
14 !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
15 !> (patch by Marcel Baer)
16 ! **************************************************************************************************
17 MODULE md_vel_utils
18  USE atomic_kind_list_types, ONLY: atomic_kind_list_type
19  USE atomic_kind_types, ONLY: atomic_kind_type,&
20  get_atomic_kind,&
21  get_atomic_kind_set
22  USE bibliography, ONLY: west2006,&
23  cite_reference
24  USE cell_types, ONLY: &
25  cell_type, use_perd_none, use_perd_x, use_perd_xy, use_perd_xyz, use_perd_xz, use_perd_y, &
26  use_perd_yz, use_perd_z
27  USE cp_linked_list_input, ONLY: cp_sll_val_next,&
28  cp_sll_val_type
29  USE cp_log_handling, ONLY: cp_get_default_logger,&
30  cp_logger_type
31  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
32  cp_print_key_unit_nr
33  USE cp_subsys_types, ONLY: cp_subsys_get,&
34  cp_subsys_type
35  USE cp_units, ONLY: cp_unit_from_cp2k
36  USE extended_system_types, ONLY: npt_info_type
37  USE force_env_methods, ONLY: force_env_calc_energy_force
38  USE force_env_types, ONLY: force_env_get,&
39  force_env_type
40  USE force_env_utils, ONLY: force_env_rattle,&
41  force_env_shake
42  USE global_types, ONLY: global_environment_type
43  USE input_constants, ONLY: &
44  md_init_default, md_init_vib, npe_f_ensemble, npe_i_ensemble, &
45  nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
46  npt_ia_ensemble, reftraj_ensemble
47  USE input_cp2k_binary_restarts, ONLY: read_binary_velocities
48  USE input_restart_force_eval, ONLY: update_subsys
49  USE input_section_types, ONLY: section_vals_get,&
50  section_vals_get_subs_vals,&
51  section_vals_list_get,&
52  section_vals_type,&
53  section_vals_val_get
54  USE input_val_types, ONLY: val_get,&
55  val_type
56  USE kinds, ONLY: default_string_length,&
57  dp
58  USE mathconstants, ONLY: pi
59  USE mathlib, ONLY: diamat_all
60  USE md_ener_types, ONLY: md_ener_type
61  USE md_environment_types, ONLY: get_md_env,&
62  md_environment_type
63  USE md_util, ONLY: read_vib_eigs_unformatted
64  USE message_passing, ONLY: mp_para_env_type
65  USE molecule_kind_list_types, ONLY: molecule_kind_list_type
66  USE molecule_kind_types, ONLY: fixd_constraint_type,&
67  get_molecule_kind,&
68  get_molecule_kind_set,&
69  molecule_kind_type
70  USE parallel_rng_types, ONLY: uniform,&
71  rng_stream_type
72  USE particle_list_types, ONLY: particle_list_type
73  USE particle_types, ONLY: particle_type
74  USE physcon, ONLY: kelvin
75  USE shell_opt, ONLY: optimize_shell_core
76  USE shell_potential_types, ONLY: shell_kind_type
77  USE simpar_types, ONLY: simpar_type
78  USE thermal_region_types, ONLY: thermal_region_type,&
79  thermal_regions_type
80 #include "../base/base_uses.f90"
81 
82  IMPLICIT NONE
83 
84  PRIVATE
85  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_vel_utils'
86 
87  PUBLIC :: temperature_control, &
88  comvel_control, &
89  angvel_control, &
90  setup_velocities
91 
92 CONTAINS
93 
94 ! **************************************************************************************************
95 !> \brief compute center of mass position
96 !> *** is only used by initialize_velocities below ***
97 !> \param part ...
98 !> \param is_fixed ...
99 !> \param rcom ...
100 !> \par History
101 !> 2007-11-6: created
102 !> \author Toon Verstraelen <Toon.Verstraelen@gmail.com>
103 ! **************************************************************************************************
104  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.
2284 999 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) &
2301  cpwarn("Warning: The primary knock-on atom is a core-shell atom")
2302  atomic_kind => particle_set(iparticle)%atomic_kind
2303  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=matom(iatom))
2304  norm = norm + matom(iatom)*weight(iatom)
2305  END DO
2306  weight(:) = matom(:)*weight(:)*energy/norm
2307  DO iatom = 1, natom
2308  norm = sqrt(dot_product(vatom(1:3, iatom), vatom(1:3, iatom)))
2309  vatom(1:3, iatom) = vatom(1:3, iatom)/norm
2310  END DO
2311 
2312  IF (iw > 0) THEN
2313  WRITE (unit=iw, fmt="(T2,A)") &
2314  "CASCADE|", &
2315  "CASCADE| Normalised velocities and additional kinetic energy [keV]", &
2316  "CASCADE|"
2317  WRITE (unit=iw, fmt="(T2,A,T11,A,3(11X,A),9X,A)") &
2318  "CASCADE| ", "Atom index", "v(x)", "v(y)", "v(z)", "E(kin)"
2319  DO iatom = 1, natom
2320  ekin = cp_unit_from_cp2k(weight(iatom), "keV")
2321  WRITE (unit=iw, fmt="(T2,A,I10,4(1X,F14.6))") &
2322  "CASCADE| ", atom_index(iatom), vatom(1:3, iatom), ekin
2323  END DO
2324  END IF
2325 
2326  ! Apply velocity modifications
2327  DO iatom = 1, natom
2328  iparticle = atom_index(iatom)
2329  particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + &
2330  sqrt(2.0_dp*weight(iatom)/matom(iatom))*vatom(1:3, iatom)
2331  END DO
2332 
2333  DEALLOCATE (atom_index)
2334  DEALLOCATE (matom)
2335  DEALLOCATE (vatom)
2336  DEALLOCATE (weight)
2337 
2338  IF (iw > 0) THEN
2339  ! Build a list of all fixed atoms (if any)
2340  ALLOCATE (is_fixed(nparticle))
2341  is_fixed = use_perd_none
2342  molecule_kind_set => molecule_kinds%els
2343  DO imolecule_kind = 1, molecule_kinds%n_els
2344  molecule_kind => molecule_kind_set(imolecule_kind)
2345  CALL get_molecule_kind(molecule_kind=molecule_kind, fixd_list=fixd_list)
2346  IF (ASSOCIATED(fixd_list)) THEN
2347  DO ifixd = 1, SIZE(fixd_list)
2348  IF (.NOT. fixd_list(ifixd)%restraint%active) is_fixed(fixd_list(ifixd)%fixd) = fixd_list(ifixd)%itype
2349  END DO
2350  END IF
2351  END DO
2352  ! Compute vcom, ecom and ekin for printout
2353  CALL compute_vcom(particle_set, is_fixed, vcom, ecom)
2354  ekin = compute_ekin(particle_set) - ecom
2355  IF (simpar%nfree == 0) THEN
2356  cpassert(ekin == 0.0_dp)
2357  temp = 0.0_dp
2358  ELSE
2359  temp = 2.0_dp*ekin/real(simpar%nfree, kind=dp)
2360  END IF
2361  temperature = cp_unit_from_cp2k(temp, "K")
2362  WRITE (unit=iw, fmt="(T2,A)") &
2363  "CASCADE|"
2364  WRITE (unit=iw, fmt="(T2,A,T61,F20.6)") &
2365  "CASCADE| Temperature after cascade initialization [K]", temperature
2366  WRITE (unit=iw, fmt="(T2,A,T30,3(1X,ES16.8))") &
2367  "CASCADE| COM velocity", vcom(1:3)
2368 !MK ! compute and log rcom and vang if not periodic
2369 !MK CALL force_env_get(force_env,cell=cell)
2370 !MK IF (SUM(cell%perd(1:3)) == 0) THEN
2371 !MK CALL compute_rcom(particle_set,is_fixed,rcom)
2372 !MK CALL compute_vang(particle_set,is_fixed,rcom,vang)
2373 !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' COM position:',rcom(1:3)
2374 !MK WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )' ) ' Angular velocity:',vang(1:3)
2375 !MK END IF
2376  DEALLOCATE (is_fixed)
2377  END IF
2378 
2379  END IF
2380 
2381  CALL cp_print_key_finished_output(iw, logger, print_section, "PROGRAM_RUN_INFO")
2382 
2383  CALL timestop(handle)
2384 
2385  END SUBROUTINE initialize_cascade
2386 
2387 END 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:28
bibliography::west2006
integer, save, public west2006
Definition: bibliography.F:43
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:803
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:1013
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:193
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)
returns various attributes about the force environment
Definition: force_env_types.F:329
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:69
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:237
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:95
input_constants::npt_i_ensemble
integer, parameter, public npt_i_ensemble
Definition: input_constants.F:95
input_constants::nph_uniaxial_damped_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
Definition: input_constants.F:95
input_constants::npe_f_ensemble
integer, parameter, public npe_f_ensemble
Definition: input_constants.F:95
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition: input_constants.F:95
input_constants::md_init_default
integer, parameter, public md_init_default
Definition: input_constants.F:117
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition: input_constants.F:95
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition: input_constants.F:95
input_constants::md_init_vib
integer, parameter, public md_init_vib
Definition: input_constants.F:117
input_constants::reftraj_ensemble
integer, parameter, public reftraj_ensemble
Definition: input_constants.F:95
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:1149
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:724
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:697
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:1040
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:372
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