(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
13MODULE 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
59CONTAINS
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
646END 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
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
extended_system_types::npt_info_type
Definition extended_system_types.F:42
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:139
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
md_ener_types::md_ener_type
Definition md_ener_types.F:22
md_environment_types::md_environment_type
Definition md_environment_types.F:53
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
qmmm_types::qmmm_env_type
Definition qmmm_types.F:34
qmmmx_types::qmmmx_env_type
Definition qmmmx_types.F:27
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30
thermostat_types::thermostat_type
Definition thermostat_types.F:78