(git:e7e05ae)
md_conserved_quantities.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 computes the conserved quantities for a given md ensemble
10 !> and also kinetic energies, thermo/barostat stuff
11 !> \author gtb, 05.02.2003
12 ! **************************************************************************************************
13 MODULE md_conserved_quantities
14  USE atomic_kind_list_types, ONLY: atomic_kind_list_type
15  USE atomic_kind_types, ONLY: atomic_kind_type,&
16  get_atomic_kind
17  USE barostat_utils, ONLY: get_baro_energies
18  USE cell_types, ONLY: cell_type
19  USE cp_subsys_types, ONLY: cp_subsys_get,&
20  cp_subsys_type
21  USE distribution_1d_types, ONLY: distribution_1d_type
22  USE extended_system_types, ONLY: npt_info_type
23  USE force_env_types, ONLY: force_env_get,&
24  force_env_type
25  USE input_constants, ONLY: &
26  isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
27  nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
28  npt_ia_ensemble, nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, reftraj_ensemble
29  USE input_section_types, ONLY: section_vals_type,&
30  section_vals_val_get
31  USE kinds, ONLY: dp
32  USE mathconstants, ONLY: zero
33  USE md_ener_types, ONLY: md_ener_type,&
34  zero_md_ener
35  USE md_environment_types, ONLY: get_md_env,&
36  md_environment_type,&
37  set_md_env
38  USE message_passing, ONLY: mp_comm_type,&
39  mp_para_env_type
40  USE particle_list_types, ONLY: particle_list_type
41  USE particle_types, ONLY: particle_type
42  USE physcon, ONLY: kelvin
43  USE qmmm_types, ONLY: qmmm_env_type
44  USE qmmm_types_low, ONLY: force_mixing_label_qm_dynamics
45  USE qmmmx_types, ONLY: qmmmx_env_type
46  USE shell_potential_types, ONLY: shell_kind_type
47  USE simpar_types, ONLY: simpar_type
48  USE thermostat_types, ONLY: thermostat_type
49  USE thermostat_utils, ONLY: get_thermostat_energies
50 #include "../base/base_uses.f90"
51 
52  IMPLICIT NONE
53 
54  PRIVATE
55 
56  PUBLIC :: compute_conserved_quantity, calc_nfree_qm
57  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_conserved_quantities'
58 
59 CONTAINS
60 
61 ! **************************************************************************************************
62 !> \brief calculates conserved quantity.
63 !> \param md_env ...
64 !> \param md_ener ...
65 !> \param tkind ...
66 !> \param tshell ...
67 !> \param natom ...
68 !> \par Input Arguments
69 !> md_env is the md_environment
70 !> epot is the total potential energy
71 !> \par Output Arguments
72 !> cons is the conserved quantity
73 !> \par Output Optional Arguments
74 !> cons_rel : relative cons. quantity (to the first md step)
75 !> ekin : kinetic energy of particles
76 !> temp : temperature
77 !> temp_qm : temperature of the QM system in a QM/MM calculation
78 !> \par History
79 !> none
80 !> \author gloria
81 ! **************************************************************************************************
82  SUBROUTINE compute_conserved_quantity(md_env, md_ener, tkind, tshell, &
83  natom)
84  TYPE(md_environment_type), POINTER :: md_env
85  TYPE(md_ener_type), POINTER :: md_ener
86  LOGICAL, INTENT(IN) :: tkind, tshell
87  INTEGER, INTENT(IN) :: natom
88 
89  INTEGER :: ikind, nfree_qm, nkind
90  INTEGER, POINTER :: itimes
91  LOGICAL :: init
92  REAL(kind=dp), POINTER :: constant
93  TYPE(mp_para_env_type), POINTER :: para_env
94  TYPE(simpar_type), POINTER :: simpar
95 
96  NULLIFY (itimes, para_env, simpar)
97 
98  CALL zero_md_ener(md_ener, tkind, tshell)
99 
100  CALL get_md_env(md_env=md_env, &
101  constant=constant, &
102  itimes=itimes, &
103  init=init, &
104  simpar=simpar, &
105  para_env=para_env)
106 
107  CALL get_part_ke(md_env, md_ener, tkind, tshell, para_env)
108 
109  IF (md_ener%nfree /= 0) THEN
110  md_ener%temp_part = 2.0_dp*md_ener%ekin/real(simpar%nfree, kind=dp)
111  md_ener%temp_part = md_ener%temp_part*kelvin
112  END IF
113 
114  nfree_qm = calc_nfree_qm(md_env, md_ener)
115  IF (nfree_qm > 0) THEN
116  md_ener%temp_qm = 2.0_dp*md_ener%ekin_qm/real(nfree_qm, kind=dp)
117  md_ener%temp_qm = md_ener%temp_qm*kelvin
118  END IF
119 
120  IF (md_ener%nfree_shell > 0) THEN
121  md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/real(md_ener%nfree_shell, kind=dp)
122  md_ener%temp_shell = md_ener%temp_shell*kelvin
123  END IF
124 
125  IF (tkind) THEN
126  nkind = SIZE(md_ener%temp_kind)
127  DO ikind = 1, nkind
128  md_ener%temp_kind(ikind) = 2.0_dp* &
129  md_ener%ekin_kind(ikind)/real(md_ener%nfree_kind(ikind), kind=dp)
130  md_ener%temp_kind(ikind) = md_ener%temp_kind(ikind)*kelvin
131  END DO
132  IF (tshell) THEN
133  DO ikind = 1, nkind
134  md_ener%temp_shell_kind(ikind) = 2.0_dp* &
135  md_ener%ekin_shell_kind(ikind)/real(md_ener%nfree_shell_kind(ikind), kind=dp)
136  md_ener%temp_shell_kind(ikind) = md_ener%temp_shell_kind(ikind)*kelvin
137  END DO
138  END IF
139  END IF
140 
141  SELECT CASE (simpar%ensemble)
142  CASE DEFAULT
143  cpabort('Unknown ensemble')
144  CASE (isokin_ensemble)
145  md_ener%constant = md_ener%ekin
146  CASE (reftraj_ensemble) ! no constant of motion available
147  md_ener%constant = md_ener%epot
148  CASE (nve_ensemble)
149  CALL get_econs_nve(md_env, md_ener, para_env)
150  CASE (nvt_ensemble)
151  CALL get_econs_nvt(md_env, md_ener, para_env)
152  CASE (npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
153  CALL get_econs_npt(md_env, md_ener, para_env)
154  md_ener%temp_baro = md_ener%temp_baro*kelvin
155  CASE (nph_uniaxial_ensemble)
156  CALL get_econs_nph_uniaxial(md_env, md_ener)
157  md_ener%temp_baro = md_ener%temp_baro*kelvin
158  CASE (nph_uniaxial_damped_ensemble)
159  CALL get_econs_nph_uniaxial(md_env, md_ener)
160  md_ener%temp_baro = md_ener%temp_baro*kelvin
161  CASE (langevin_ensemble)
162  md_ener%constant = md_ener%ekin + md_ener%epot
163  CASE (npe_f_ensemble, npe_i_ensemble)
164  CALL get_econs_npe(md_env, md_ener, para_env)
165  md_ener%temp_baro = md_ener%temp_baro*kelvin
166  CASE (nvt_adiabatic_ensemble)
167  CALL get_econs_nvt_adiabatic(md_env, md_ener, para_env)
168  END SELECT
169 
170  IF (init) THEN
171  ! If the value was not read from input let's set it at the begin of the MD
172  IF (constant == 0.0_dp) THEN
173  constant = md_ener%constant
174  CALL set_md_env(md_env=md_env, constant=constant)
175  END IF
176  ELSE
177  CALL get_md_env(md_env=md_env, constant=constant)
178  md_ener%delta_cons = (md_ener%constant - constant)/real(natom, kind=dp)*kelvin
179  END IF
180 
181  END SUBROUTINE compute_conserved_quantity
182 
183 ! **************************************************************************************************
184 !> \brief Calculates the number of QM degress of freedom
185 !> \param md_env ...
186 !> \param md_ener ...
187 !> \return ...
188 ! **************************************************************************************************
189  FUNCTION calc_nfree_qm(md_env, md_ener) RESULT(nfree_qm)
190  TYPE(md_environment_type), POINTER :: md_env
191  TYPE(md_ener_type), POINTER :: md_ener
192  INTEGER :: nfree_qm
193 
194  INTEGER :: ip
195  INTEGER, POINTER :: cur_indices(:), cur_labels(:)
196  TYPE(cp_subsys_type), POINTER :: subsys
197  TYPE(force_env_type), POINTER :: force_env
198  TYPE(particle_list_type), POINTER :: particles
199  TYPE(qmmm_env_type), POINTER :: qmmm_env
200  TYPE(qmmmx_env_type), POINTER :: qmmmx_env
201  TYPE(section_vals_type), POINTER :: force_env_section
202 
203  NULLIFY (qmmm_env, qmmmx_env, subsys, particles, force_env, force_env_section)
204  nfree_qm = 0
205 
206  CALL get_md_env(md_env, force_env=force_env)
207  CALL force_env_get(force_env, &
208  subsys=subsys, &
209  qmmm_env=qmmm_env, &
210  qmmmx_env=qmmmx_env, &
211  force_env_section=force_env_section)
212 
213  IF (ASSOCIATED(qmmm_env)) THEN ! conventional QM/MM
214  CALL cp_subsys_get(subsys, particles=particles)
215  ! The degrees of freedom for the quantum part of the system
216  ! are set to 3*Number of QM atoms and to simpar%nfree in case all the MM
217  ! system is treated at QM level (not really QM/MM, just for consistency).
218  ! The degree of freedom will not be correct if 1-3 atoms are treated only
219  ! MM. In this case we should take care of rotations
220  nfree_qm = 3*SIZE(qmmm_env%qm%qm_atom_index)
221  IF (nfree_qm == 3*(particles%n_els)) nfree_qm = md_ener%nfree
222  END IF
223 
224  IF (ASSOCIATED(qmmmx_env)) THEN ! doing force mixing
225  CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
226  CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
227  nfree_qm = 0
228  DO ip = 1, SIZE(cur_indices)
229  IF (cur_labels(ip) >= force_mixing_label_qm_dynamics) THEN ! this is a QM atom
230  nfree_qm = nfree_qm + 3
231  END IF
232  END DO
233  END IF
234 
235  cpassert(.NOT. (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)))
236  END FUNCTION calc_nfree_qm
237 
238 ! **************************************************************************************************
239 !> \brief calculates conserved quantity for nvt ensemble
240 !> \param md_env ...
241 !> \param md_ener ...
242 !> \param para_env ...
243 !> \par History
244 !> none
245 !> \author gloria
246 ! **************************************************************************************************
247  SUBROUTINE get_econs_nve(md_env, md_ener, para_env)
248  TYPE(md_environment_type), POINTER :: md_env
249  TYPE(md_ener_type), INTENT(inout) :: md_ener
250  TYPE(mp_para_env_type), POINTER :: para_env
251 
252  TYPE(force_env_type), POINTER :: force_env
253  TYPE(thermostat_type), POINTER :: thermostat_coeff, thermostat_shell
254 
255  NULLIFY (force_env, thermostat_coeff, thermostat_shell)
256 
257  CALL get_md_env(md_env, force_env=force_env, thermostat_coeff=thermostat_coeff, &
258  thermostat_shell=thermostat_shell)
259  md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell
260 
261  CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
262  md_ener%thermostat_shell_kin, para_env)
263  md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
264 
265  END SUBROUTINE get_econs_nve
266 
267 ! **************************************************************************************************
268 !> \brief calculates conserved quantity for nvt ensemble
269 !> \param md_env ...
270 !> \param md_ener ...
271 !> \param para_env ...
272 !> \par History
273 !> none
274 !> \author gloria
275 ! **************************************************************************************************
276  SUBROUTINE get_econs_nvt_adiabatic(md_env, md_ener, para_env)
277  TYPE(md_environment_type), POINTER :: md_env
278  TYPE(md_ener_type), INTENT(inout) :: md_ener
279  TYPE(mp_para_env_type), POINTER :: para_env
280 
281  TYPE(force_env_type), POINTER :: force_env
282  TYPE(thermostat_type), POINTER :: thermostat_fast, thermostat_slow
283 
284  NULLIFY (force_env, thermostat_fast, thermostat_slow)
285  CALL get_md_env(md_env, force_env=force_env, thermostat_fast=thermostat_fast, &
286  thermostat_slow=thermostat_slow)
287  CALL get_thermostat_energies(thermostat_fast, md_ener%thermostat_fast_pot, &
288  md_ener%thermostat_fast_kin, para_env)
289  md_ener%constant = md_ener%ekin + md_ener%epot + &
290  md_ener%thermostat_fast_kin + md_ener%thermostat_fast_pot
291  CALL get_thermostat_energies(thermostat_slow, md_ener%thermostat_slow_pot, &
292  md_ener%thermostat_slow_kin, para_env)
293  md_ener%constant = md_ener%constant + &
294  md_ener%thermostat_slow_kin + md_ener%thermostat_slow_pot
295 
296  END SUBROUTINE get_econs_nvt_adiabatic
297 
298 ! **************************************************************************************************
299 !> \brief calculates conserved quantity for nvt ensemble
300 !> \param md_env ...
301 !> \param md_ener ...
302 !> \param para_env ...
303 !> \par History
304 !> none
305 !> \author gloria
306 ! **************************************************************************************************
307  SUBROUTINE get_econs_nvt(md_env, md_ener, para_env)
308  TYPE(md_environment_type), POINTER :: md_env
309  TYPE(md_ener_type), INTENT(inout) :: md_ener
310  TYPE(mp_para_env_type), POINTER :: para_env
311 
312  TYPE(force_env_type), POINTER :: force_env
313  TYPE(thermostat_type), POINTER :: thermostat_coeff, thermostat_part, &
314  thermostat_shell
315 
316  NULLIFY (force_env, thermostat_part, thermostat_coeff, thermostat_shell)
317  CALL get_md_env(md_env, force_env=force_env, thermostat_part=thermostat_part, &
318  thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell)
319  CALL get_thermostat_energies(thermostat_part, md_ener%thermostat_part_pot, &
320  md_ener%thermostat_part_kin, para_env)
321  md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell + &
322  md_ener%thermostat_part_kin + md_ener%thermostat_part_pot
323 
324  CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
325  md_ener%thermostat_shell_kin, para_env)
326  md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
327 
328  END SUBROUTINE get_econs_nvt
329 
330 ! **************************************************************************************************
331 !> \brief calculates conserved quantity for npe ensemble
332 !> \param md_env ...
333 !> \param md_ener ...
334 !> \param para_env ...
335 !> \par History
336 !> none
337 !> \author marcella (02-2008)
338 ! **************************************************************************************************
339  SUBROUTINE get_econs_npe(md_env, md_ener, para_env)
340  TYPE(md_environment_type), POINTER :: md_env
341  TYPE(md_ener_type), INTENT(inout) :: md_ener
342  TYPE(mp_para_env_type), POINTER :: para_env
343 
344  INTEGER :: nfree
345  TYPE(cell_type), POINTER :: box
346  TYPE(npt_info_type), POINTER :: npt(:, :)
347  TYPE(simpar_type), POINTER :: simpar
348  TYPE(thermostat_type), POINTER :: thermostat_baro, thermostat_shell
349 
350  NULLIFY (thermostat_baro, thermostat_shell, npt)
351  CALL get_md_env(md_env, thermostat_baro=thermostat_baro, &
352  simpar=simpar, npt=npt, cell=box, &
353  thermostat_shell=thermostat_shell)
354  CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, &
355  md_ener%baro_pot)
356  nfree = SIZE(npt, 1)*SIZE(npt, 2)
357  md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/nfree
358 
359  md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell &
360  + md_ener%baro_kin + md_ener%baro_pot
361 
362  CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
363  md_ener%thermostat_shell_kin, para_env)
364  md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + &
365  md_ener%thermostat_shell_pot
366 
367  END SUBROUTINE get_econs_npe
368 
369 ! **************************************************************************************************
370 !> \brief calculates conserved quantity for npt ensemble
371 !> \param md_env ...
372 !> \param md_ener ...
373 !> \param para_env ...
374 !> \par History
375 !> none
376 !> \author gloria
377 ! **************************************************************************************************
378  SUBROUTINE get_econs_npt(md_env, md_ener, para_env)
379  TYPE(md_environment_type), POINTER :: md_env
380  TYPE(md_ener_type), INTENT(inout) :: md_ener
381  TYPE(mp_para_env_type), POINTER :: para_env
382 
383  INTEGER :: nfree
384  TYPE(cell_type), POINTER :: box
385  TYPE(npt_info_type), POINTER :: npt(:, :)
386  TYPE(simpar_type), POINTER :: simpar
387  TYPE(thermostat_type), POINTER :: thermostat_baro, thermostat_part, &
388  thermostat_shell
389 
390  NULLIFY (thermostat_baro, thermostat_part, thermostat_shell, npt, simpar, box)
391  CALL get_md_env(md_env, thermostat_part=thermostat_part, thermostat_baro=thermostat_baro, &
392  simpar=simpar, npt=npt, cell=box, thermostat_shell=thermostat_shell)
393  CALL get_thermostat_energies(thermostat_part, md_ener%thermostat_part_pot, &
394  md_ener%thermostat_part_kin, para_env)
395  CALL get_thermostat_energies(thermostat_baro, md_ener%thermostat_baro_pot, &
396  md_ener%thermostat_baro_kin, para_env)
397  CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, md_ener%baro_pot)
398  nfree = SIZE(npt, 1)*SIZE(npt, 2)
399  md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/nfree
400 
401  md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell &
402  + md_ener%thermostat_part_kin + md_ener%thermostat_part_pot &
403  + md_ener%thermostat_baro_kin + md_ener%thermostat_baro_pot &
404  + md_ener%baro_kin + md_ener%baro_pot
405 
406  CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
407  md_ener%thermostat_shell_kin, para_env)
408  md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
409 
410  END SUBROUTINE get_econs_npt
411 
412 ! **************************************************************************************************
413 !> \brief calculates conserved quantity for nph_uniaxial
414 !> \param md_env ...
415 !> \param md_ener ...
416 !> \par History
417 !> none
418 !> \author cjm
419 ! **************************************************************************************************
420  SUBROUTINE get_econs_nph_uniaxial(md_env, md_ener)
421  TYPE(md_environment_type), POINTER :: md_env
422  TYPE(md_ener_type), INTENT(inout) :: md_ener
423 
424  TYPE(cell_type), POINTER :: box
425  TYPE(npt_info_type), POINTER :: npt(:, :)
426  TYPE(simpar_type), POINTER :: simpar
427 
428  CALL get_md_env(md_env, simpar=simpar, npt=npt, cell=box)
429 
430  CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, md_ener%baro_pot)
431  md_ener%temp_baro = 2.0_dp*md_ener%baro_kin
432  md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%baro_kin + md_ener%baro_pot
433  END SUBROUTINE get_econs_nph_uniaxial
434 
435 ! **************************************************************************************************
436 !> \brief Calculates kinetic energy of particles
437 !> \param md_env ...
438 !> \param md_ener ...
439 !> \param tkind ...
440 !> \param tshell ...
441 !> \param group ...
442 !> \par History
443 !> none
444 !> \author CJM
445 ! **************************************************************************************************
446  SUBROUTINE get_part_ke(md_env, md_ener, tkind, tshell, group)
447  TYPE(md_environment_type), POINTER :: md_env
448  TYPE(md_ener_type), POINTER :: md_ener
449  LOGICAL, INTENT(IN) :: tkind, tshell
450 
451  CLASS(mp_comm_type), INTENT(IN) :: group
452 
453  INTEGER :: i, iparticle, iparticle_kind, &
454  iparticle_local, nparticle_kind, &
455  nparticle_local, shell_index
456  INTEGER, POINTER :: cur_indices(:), cur_labels(:)
457  LOGICAL :: is_shell
458  REAL(kind=dp) :: ekin_c, ekin_com, ekin_s, mass
459  TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
460  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
461  TYPE(atomic_kind_type), POINTER :: atomic_kind
462  TYPE(cp_subsys_type), POINTER :: subsys
463  TYPE(distribution_1d_type), POINTER :: local_particles
464  TYPE(force_env_type), POINTER :: force_env
465  TYPE(particle_list_type), POINTER :: core_particles, particles, &
466  shell_particles
467  TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
468  shell_particle_set
469  TYPE(qmmm_env_type), POINTER :: qmmm_env
470  TYPE(qmmmx_env_type), POINTER :: qmmmx_env
471  TYPE(section_vals_type), POINTER :: force_env_section
472  TYPE(shell_kind_type), POINTER :: shell
473 
474  NULLIFY (force_env, qmmm_env, qmmmx_env, subsys, force_env_section)
475  CALL get_md_env(md_env, force_env=force_env)
476  CALL force_env_get(force_env, &
477  subsys=subsys, &
478  qmmm_env=qmmm_env, &
479  qmmmx_env=qmmmx_env, &
480  force_env_section=force_env_section)
481 
482  CALL cp_subsys_get(subsys=subsys, &
483  atomic_kinds=atomic_kinds, &
484  local_particles=local_particles, &
485  particles=particles, shell_particles=shell_particles, &
486  core_particles=core_particles)
487 
488  nparticle_kind = atomic_kinds%n_els
489  atomic_kind_set => atomic_kinds%els
490 
491  ekin_s = zero
492  ekin_c = zero
493  ekin_com = zero
494  IF (tkind) THEN
495  md_ener%nfree_kind = 0
496  IF (tshell) THEN
497  md_ener%nfree_shell_kind = 0
498  END IF
499  END IF
500 
501  particle_set => particles%els
502  IF (tshell) THEN
503  shell_particle_set => shell_particles%els
504  core_particle_set => core_particles%els
505  DO iparticle_kind = 1, nparticle_kind
506  atomic_kind => atomic_kind_set(iparticle_kind)
507  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
508  shell_active=is_shell, shell=shell)
509  nparticle_local = local_particles%n_el(iparticle_kind)
510  IF (is_shell) THEN
511  DO iparticle_local = 1, nparticle_local
512  iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
513  shell_index = particle_set(iparticle)%shell_index
514  !ekin
515  ekin_com = 0.5_dp*mass* &
516  (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
517  + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
518  + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
519  !vcom
520  md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
521  md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
522  md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
523  md_ener%total_mass = md_ener%total_mass + mass
524 
525  md_ener%ekin = md_ener%ekin + ekin_com
526  ekin_c = 0.5_dp*shell%mass_core* &
527  (core_particle_set(shell_index)%v(1)*core_particle_set(shell_index)%v(1) &
528  + core_particle_set(shell_index)%v(2)*core_particle_set(shell_index)%v(2) &
529  + core_particle_set(shell_index)%v(3)*core_particle_set(shell_index)%v(3))
530  ekin_s = 0.5_dp*shell%mass_shell* &
531  (shell_particle_set(shell_index)%v(1)*shell_particle_set(shell_index)%v(1) &
532  + shell_particle_set(shell_index)%v(2)*shell_particle_set(shell_index)%v(2) &
533  + shell_particle_set(shell_index)%v(3)*shell_particle_set(shell_index)%v(3))
534  md_ener%ekin_shell = md_ener%ekin_shell + ekin_c + ekin_s - ekin_com
535 
536  IF (tkind) THEN
537  md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
538  md_ener%nfree_kind(iparticle_kind) = md_ener%nfree_kind(iparticle_kind) + 3
539  md_ener%ekin_shell_kind(iparticle_kind) = md_ener%ekin_shell_kind(iparticle_kind) + &
540  ekin_c + ekin_s - ekin_com
541  md_ener%nfree_shell_kind(iparticle_kind) = md_ener%nfree_shell_kind(iparticle_kind) + 3
542  END IF
543 
544  END DO ! iparticle_local
545  ELSE
546  DO iparticle_local = 1, nparticle_local
547  iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
548  ekin_com = 0.5_dp*mass* &
549  (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
550  + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
551  + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
552  !vcom
553  md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
554  md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
555  md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
556  md_ener%total_mass = md_ener%total_mass + mass
557 
558  md_ener%ekin = md_ener%ekin + ekin_com
559  IF (tkind) THEN
560  md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
561  md_ener%nfree_kind(iparticle_kind) = md_ener%nfree_kind(iparticle_kind) + 3
562  END IF
563  END DO ! iparticle_local
564  END IF
565  END DO ! iparticle_kind
566  IF (tkind) THEN
567  CALL group%sum(md_ener%ekin_kind)
568  CALL group%sum(md_ener%nfree_kind)
569  CALL group%sum(md_ener%ekin_shell_kind)
570  CALL group%sum(md_ener%nfree_shell_kind)
571  END IF
572  ! sum all contributions to energy over calculated parts on all processors
573  CALL group%sum(md_ener%ekin_shell)
574  ELSE
575  DO iparticle_kind = 1, nparticle_kind
576  atomic_kind => atomic_kind_set(iparticle_kind)
577  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
578  nparticle_local = local_particles%n_el(iparticle_kind)
579  DO iparticle_local = 1, nparticle_local
580  iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
581  ! ekin
582  ekin_com = 0.5_dp*mass* &
583  (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
584  + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
585  + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
586 
587  !vcom
588  md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
589  md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
590  md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
591  md_ener%total_mass = md_ener%total_mass + mass
592 
593  md_ener%ekin = md_ener%ekin + ekin_com
594  IF (tkind) THEN
595  md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
596  md_ener%nfree_kind(iparticle_kind) = md_ener%nfree_kind(iparticle_kind) + 3
597  END IF
598  END DO
599  END DO ! iparticle_kind
600  IF (tkind) THEN
601  CALL group%sum(md_ener%ekin_kind)
602  CALL group%sum(md_ener%nfree_kind)
603  END IF
604  END IF
605 
606  ! sum all contributions to energy over calculated parts on all processors
607  CALL group%sum(md_ener%ekin)
608  CALL group%sum(md_ener%vcom)
609  CALL group%sum(md_ener%total_mass)
610  md_ener%vcom = md_ener%vcom/md_ener%total_mass
611  !
612  ! Compute the QM/MM kinetic energy
613 
614  IF (ASSOCIATED(qmmm_env)) THEN ! conventional QM/MM
615  DO i = 1, SIZE(qmmm_env%qm%qm_atom_index)
616  iparticle = qmmm_env%qm%qm_atom_index(i)
617  mass = particle_set(iparticle)%atomic_kind%mass
618  md_ener%ekin_qm = md_ener%ekin_qm + 0.5_dp*mass* &
619  (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
620  + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
621  + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
622  END DO
623  END IF
624 
625  IF (ASSOCIATED(qmmmx_env)) THEN ! doing force mixing
626  CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
627  CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
628  DO i = 1, SIZE(cur_indices)
629  IF (cur_labels(i) >= force_mixing_label_qm_dynamics) THEN ! this is a QM atom
630  iparticle = cur_indices(i)
631  mass = particle_set(iparticle)%atomic_kind%mass
632  md_ener%ekin_qm = md_ener%ekin_qm + 0.5_dp*mass* &
633  (particle_set(iparticle)%v(1)*particle_set(iparticle)%v(1) &
634  + particle_set(iparticle)%v(2)*particle_set(iparticle)%v(2) &
635  + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
636  END IF
637  END DO
638  END IF
639 
640  IF (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)) &
641  cpabort("get_part_ke: qmmm bug")
642  END SUBROUTINE get_part_ke
643 
644 ! **************************************************************************************************
645 
646 END MODULE md_conserved_quantities
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
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
barostat_utils
Barostat utils.
Definition: barostat_utils.F:12
barostat_utils::get_baro_energies
subroutine, public get_baro_energies(cell, simpar, npt, baro_kin, baro_pot)
Calculates kinetic energy and potential of barostat.
Definition: barostat_utils.F:57
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
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
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
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_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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::nvt_adiabatic_ensemble
integer, parameter, public nvt_adiabatic_ensemble
Definition: input_constants.F:95
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::isokin_ensemble
integer, parameter, public isokin_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::langevin_ensemble
integer, parameter, public langevin_ensemble
Definition: input_constants.F:95
input_constants::npe_i_ensemble
integer, parameter, public npe_i_ensemble
Definition: input_constants.F:95
input_constants::npt_ia_ensemble
integer, parameter, public npt_ia_ensemble
Definition: input_constants.F:95
input_constants::nve_ensemble
integer, parameter, public nve_ensemble
Definition: input_constants.F:95
input_constants::npt_f_ensemble
integer, parameter, public npt_f_ensemble
Definition: input_constants.F:95
input_constants::reftraj_ensemble
integer, parameter, public reftraj_ensemble
Definition: input_constants.F:95
input_constants::nvt_ensemble
integer, parameter, public nvt_ensemble
Definition: input_constants.F:95
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_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
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::zero
real(kind=dp), parameter, public zero
Definition: mathconstants.F:123
md_conserved_quantities
computes the conserved quantities for a given md ensemble and also kinetic energies,...
Definition: md_conserved_quantities.F:13
md_conserved_quantities::calc_nfree_qm
integer function, public calc_nfree_qm(md_env, md_ener)
Calculates the number of QM degress of freedom.
Definition: md_conserved_quantities.F:190
md_conserved_quantities::compute_conserved_quantity
subroutine, public compute_conserved_quantity(md_env, md_ener, tkind, tshell, natom)
calculates conserved quantity.
Definition: md_conserved_quantities.F:84
md_ener_types
Split md_ener module from md_environment_type.
Definition: md_ener_types.F:12
md_ener_types::zero_md_ener
subroutine, public zero_md_ener(md_ener, tkind, tshell)
initialize to zero energies and temperatures
Definition: md_ener_types.F:113
md_environment_types
Definition: md_environment_types.F:15
md_environment_types::set_md_env
subroutine, public set_md_env(md_env, itimes, constant, cell, simpar, fe_env, force_env, para_env, init, first_time, thermostats, barostat, reftraj, md_ener, averages, thermal_regions, ehrenfest_md)
Set the integrator environment to the correct program.
Definition: md_environment_types.F:291
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
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
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
qmmm_types_low
Definition: qmmm_types_low.F:13
qmmm_types_low::force_mixing_label_qm_dynamics
integer, parameter, public force_mixing_label_qm_dynamics
Definition: qmmm_types_low.F:44
qmmm_types
Basic container type for QM/MM.
Definition: qmmm_types.F:12
qmmmx_types
Basic container type for QM/MM with force mixing.
Definition: qmmmx_types.F:12
shell_potential_types
Definition: shell_potential_types.F:11
simpar_types
Type for storing MD parameters.
Definition: simpar_types.F:14
thermostat_types
Thermostat structure: module containing thermostat available for MD.
Definition: thermostat_types.F:12
thermostat_utils
Utilities for thermostats.
Definition: thermostat_utils.F:12
thermostat_utils::get_thermostat_energies
subroutine, public get_thermostat_energies(thermostat, thermostat_pot, thermostat_kin, para_env, array_pot, array_kin)
Calculates energy associated with a thermostat.
Definition: thermostat_utils.F:1576