Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2023 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 41085 : 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 41085 : NULLIFY (itimes, para_env, simpar)
97 :
98 41085 : 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 41085 : para_env=para_env)
106 :
107 41085 : CALL get_part_ke(md_env, md_ener, tkind, tshell, para_env)
108 :
109 41085 : IF (md_ener%nfree /= 0) THEN
110 40831 : md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(simpar%nfree, KIND=dp)
111 40831 : md_ener%temp_part = md_ener%temp_part*kelvin
112 : END IF
113 :
114 41085 : nfree_qm = calc_nfree_qm(md_env, md_ener)
115 41085 : IF (nfree_qm > 0) THEN
116 1384 : md_ener%temp_qm = 2.0_dp*md_ener%ekin_qm/REAL(nfree_qm, KIND=dp)
117 1384 : md_ener%temp_qm = md_ener%temp_qm*kelvin
118 : END IF
119 :
120 41085 : IF (md_ener%nfree_shell > 0) THEN
121 2420 : md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell, KIND=dp)
122 2420 : md_ener%temp_shell = md_ener%temp_shell*kelvin
123 : END IF
124 :
125 41085 : IF (tkind) THEN
126 494 : nkind = SIZE(md_ener%temp_kind)
127 1490 : DO ikind = 1, nkind
128 : md_ener%temp_kind(ikind) = 2.0_dp* &
129 996 : md_ener%ekin_kind(ikind)/REAL(md_ener%nfree_kind(ikind), KIND=dp)
130 1490 : md_ener%temp_kind(ikind) = md_ener%temp_kind(ikind)*kelvin
131 : END DO
132 494 : IF (tshell) THEN
133 1134 : DO ikind = 1, nkind
134 : md_ener%temp_shell_kind(ikind) = 2.0_dp* &
135 756 : md_ener%ekin_shell_kind(ikind)/REAL(md_ener%nfree_shell_kind(ikind), KIND=dp)
136 1134 : md_ener%temp_shell_kind(ikind) = md_ener%temp_shell_kind(ikind)*kelvin
137 : END DO
138 : END IF
139 : END IF
140 :
141 41085 : SELECT CASE (simpar%ensemble)
142 : CASE DEFAULT
143 0 : CPABORT('Unknown ensemble')
144 : CASE (isokin_ensemble)
145 8 : md_ener%constant = md_ener%ekin
146 : CASE (reftraj_ensemble) ! no constant of motion available
147 0 : md_ener%constant = md_ener%epot
148 : CASE (nve_ensemble)
149 30831 : CALL get_econs_nve(md_env, md_ener, para_env)
150 : CASE (nvt_ensemble)
151 7714 : CALL get_econs_nvt(md_env, md_ener, para_env)
152 : CASE (npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
153 1908 : CALL get_econs_npt(md_env, md_ener, para_env)
154 1908 : md_ener%temp_baro = md_ener%temp_baro*kelvin
155 : CASE (nph_uniaxial_ensemble)
156 44 : CALL get_econs_nph_uniaxial(md_env, md_ener)
157 44 : md_ener%temp_baro = md_ener%temp_baro*kelvin
158 : CASE (nph_uniaxial_damped_ensemble)
159 22 : CALL get_econs_nph_uniaxial(md_env, md_ener)
160 22 : md_ener%temp_baro = md_ener%temp_baro*kelvin
161 : CASE (langevin_ensemble)
162 264 : md_ener%constant = md_ener%ekin + md_ener%epot
163 : CASE (npe_f_ensemble, npe_i_ensemble)
164 294 : CALL get_econs_npe(md_env, md_ener, para_env)
165 294 : md_ener%temp_baro = md_ener%temp_baro*kelvin
166 : CASE (nvt_adiabatic_ensemble)
167 41085 : CALL get_econs_nvt_adiabatic(md_env, md_ener, para_env)
168 : END SELECT
169 :
170 41085 : IF (init) THEN
171 : ! If the value was not read from input let's set it at the begin of the MD
172 1697 : IF (constant == 0.0_dp) THEN
173 1509 : constant = md_ener%constant
174 1509 : CALL set_md_env(md_env=md_env, constant=constant)
175 : END IF
176 : ELSE
177 39388 : CALL get_md_env(md_env=md_env, constant=constant)
178 39388 : md_ener%delta_cons = (md_ener%constant - constant)/REAL(natom, KIND=dp)*kelvin
179 : END IF
180 :
181 41085 : 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 82457 : 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 82457 : 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 82457 : NULLIFY (qmmm_env, qmmmx_env, subsys, particles, force_env, force_env_section)
204 82457 : nfree_qm = 0
205 :
206 82457 : 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 82457 : force_env_section=force_env_section)
212 :
213 82457 : IF (ASSOCIATED(qmmm_env)) THEN ! conventional QM/MM
214 2744 : 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 2744 : nfree_qm = 3*SIZE(qmmm_env%qm%qm_atom_index)
221 2744 : IF (nfree_qm == 3*(particles%n_els)) nfree_qm = md_ener%nfree
222 : END IF
223 :
224 82457 : IF (ASSOCIATED(qmmmx_env)) THEN ! doing force mixing
225 64 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
226 64 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
227 64 : nfree_qm = 0
228 3520 : DO ip = 1, SIZE(cur_indices)
229 3520 : IF (cur_labels(ip) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom
230 1002 : nfree_qm = nfree_qm + 3
231 : END IF
232 : END DO
233 : END IF
234 :
235 82457 : CPASSERT(.NOT. (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)))
236 82457 : 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 30831 : 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 30831 : NULLIFY (force_env, thermostat_coeff, thermostat_shell)
256 :
257 : CALL get_md_env(md_env, force_env=force_env, thermostat_coeff=thermostat_coeff, &
258 30831 : thermostat_shell=thermostat_shell)
259 30831 : 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 30831 : md_ener%thermostat_shell_kin, para_env)
263 30831 : md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
264 :
265 30831 : 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 0 : 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 0 : NULLIFY (force_env, thermostat_fast, thermostat_slow)
285 : CALL get_md_env(md_env, force_env=force_env, thermostat_fast=thermostat_fast, &
286 0 : thermostat_slow=thermostat_slow)
287 : CALL get_thermostat_energies(thermostat_fast, md_ener%thermostat_fast_pot, &
288 0 : md_ener%thermostat_fast_kin, para_env)
289 : md_ener%constant = md_ener%ekin + md_ener%epot + &
290 0 : md_ener%thermostat_fast_kin + md_ener%thermostat_fast_pot
291 : CALL get_thermostat_energies(thermostat_slow, md_ener%thermostat_slow_pot, &
292 0 : md_ener%thermostat_slow_kin, para_env)
293 : md_ener%constant = md_ener%constant + &
294 0 : md_ener%thermostat_slow_kin + md_ener%thermostat_slow_pot
295 :
296 0 : 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 7714 : 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 7714 : 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 7714 : thermostat_coeff=thermostat_coeff, thermostat_shell=thermostat_shell)
319 : CALL get_thermostat_energies(thermostat_part, md_ener%thermostat_part_pot, &
320 7714 : md_ener%thermostat_part_kin, para_env)
321 : md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%ekin_shell + &
322 7714 : md_ener%thermostat_part_kin + md_ener%thermostat_part_pot
323 :
324 : CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
325 7714 : md_ener%thermostat_shell_kin, para_env)
326 7714 : md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
327 :
328 7714 : 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 294 : 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 294 : TYPE(npt_info_type), POINTER :: npt(:, :)
347 : TYPE(simpar_type), POINTER :: simpar
348 : TYPE(thermostat_type), POINTER :: thermostat_baro, thermostat_shell
349 :
350 294 : 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 294 : thermostat_shell=thermostat_shell)
354 : CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, &
355 294 : md_ener%baro_pot)
356 294 : nfree = SIZE(npt, 1)*SIZE(npt, 2)
357 294 : 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 294 : + md_ener%baro_kin + md_ener%baro_pot
361 :
362 : CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
363 294 : md_ener%thermostat_shell_kin, para_env)
364 : md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + &
365 294 : md_ener%thermostat_shell_pot
366 :
367 294 : 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 1908 : 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 1908 : 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 1908 : 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 1908 : simpar=simpar, npt=npt, cell=box, thermostat_shell=thermostat_shell)
393 : CALL get_thermostat_energies(thermostat_part, md_ener%thermostat_part_pot, &
394 1908 : md_ener%thermostat_part_kin, para_env)
395 : CALL get_thermostat_energies(thermostat_baro, md_ener%thermostat_baro_pot, &
396 1908 : md_ener%thermostat_baro_kin, para_env)
397 1908 : CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, md_ener%baro_pot)
398 1908 : nfree = SIZE(npt, 1)*SIZE(npt, 2)
399 1908 : 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 1908 : + md_ener%baro_kin + md_ener%baro_pot
405 :
406 : CALL get_thermostat_energies(thermostat_shell, md_ener%thermostat_shell_pot, &
407 1908 : md_ener%thermostat_shell_kin, para_env)
408 1908 : md_ener%constant = md_ener%constant + md_ener%thermostat_shell_kin + md_ener%thermostat_shell_pot
409 :
410 1908 : 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 66 : 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 66 : TYPE(npt_info_type), POINTER :: npt(:, :)
426 : TYPE(simpar_type), POINTER :: simpar
427 :
428 66 : CALL get_md_env(md_env, simpar=simpar, npt=npt, cell=box)
429 :
430 66 : CALL get_baro_energies(box, simpar, npt, md_ener%baro_kin, md_ener%baro_pot)
431 66 : md_ener%temp_baro = 2.0_dp*md_ener%baro_kin
432 66 : md_ener%constant = md_ener%ekin + md_ener%epot + md_ener%baro_kin + md_ener%baro_pot
433 66 : 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 41085 : 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 41085 : 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 41085 : 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 41085 : TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
468 41085 : 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 41085 : NULLIFY (force_env, qmmm_env, qmmmx_env, subsys, force_env_section)
475 41085 : 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 41085 : 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 41085 : core_particles=core_particles)
487 :
488 41085 : nparticle_kind = atomic_kinds%n_els
489 41085 : atomic_kind_set => atomic_kinds%els
490 :
491 41085 : ekin_s = zero
492 41085 : ekin_c = zero
493 41085 : ekin_com = zero
494 41085 : IF (tkind) THEN
495 1490 : md_ener%nfree_kind = 0
496 494 : IF (tshell) THEN
497 1134 : md_ener%nfree_shell_kind = 0
498 : END IF
499 : END IF
500 :
501 41085 : particle_set => particles%els
502 41085 : IF (tshell) THEN
503 2420 : shell_particle_set => shell_particles%els
504 2420 : core_particle_set => core_particles%els
505 7238 : DO iparticle_kind = 1, nparticle_kind
506 4818 : atomic_kind => atomic_kind_set(iparticle_kind)
507 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
508 4818 : shell_active=is_shell, shell=shell)
509 4818 : nparticle_local = local_particles%n_el(iparticle_kind)
510 7238 : IF (is_shell) THEN
511 124659 : DO iparticle_local = 1, nparticle_local
512 119883 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
513 119883 : 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 119883 : + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
519 : !vcom
520 119883 : md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
521 119883 : md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
522 119883 : md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
523 119883 : md_ener%total_mass = md_ener%total_mass + mass
524 :
525 119883 : 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 119883 : + 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 119883 : + shell_particle_set(shell_index)%v(3)*shell_particle_set(shell_index)%v(3))
534 119883 : md_ener%ekin_shell = md_ener%ekin_shell + ekin_c + ekin_s - ekin_com
535 :
536 124659 : IF (tkind) THEN
537 18144 : md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
538 18144 : 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 18144 : ekin_c + ekin_s - ekin_com
541 18144 : 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 714 : DO iparticle_local = 1, nparticle_local
547 672 : 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 672 : + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
552 : !vcom
553 672 : md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
554 672 : md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
555 672 : md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
556 672 : md_ener%total_mass = md_ener%total_mass + mass
557 :
558 672 : md_ener%ekin = md_ener%ekin + ekin_com
559 714 : IF (tkind) THEN
560 0 : md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
561 0 : 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 2420 : IF (tkind) THEN
567 1890 : CALL group%sum(md_ener%ekin_kind)
568 1890 : CALL group%sum(md_ener%nfree_kind)
569 1890 : CALL group%sum(md_ener%ekin_shell_kind)
570 1890 : CALL group%sum(md_ener%nfree_shell_kind)
571 : END IF
572 : ! sum all contributions to energy over calculated parts on all processors
573 2420 : CALL group%sum(md_ener%ekin_shell)
574 : ELSE
575 143360 : DO iparticle_kind = 1, nparticle_kind
576 104695 : atomic_kind => atomic_kind_set(iparticle_kind)
577 104695 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
578 104695 : nparticle_local = local_particles%n_el(iparticle_kind)
579 3108878 : DO iparticle_local = 1, nparticle_local
580 2965518 : 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 2965518 : + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
586 :
587 : !vcom
588 2965518 : md_ener%vcom(1) = md_ener%vcom(1) + particle_set(iparticle)%v(1)*mass
589 2965518 : md_ener%vcom(2) = md_ener%vcom(2) + particle_set(iparticle)%v(2)*mass
590 2965518 : md_ener%vcom(3) = md_ener%vcom(3) + particle_set(iparticle)%v(3)*mass
591 2965518 : md_ener%total_mass = md_ener%total_mass + mass
592 :
593 2965518 : md_ener%ekin = md_ener%ekin + ekin_com
594 3070213 : IF (tkind) THEN
595 2270 : md_ener%ekin_kind(iparticle_kind) = md_ener%ekin_kind(iparticle_kind) + ekin_com
596 2270 : 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 38665 : IF (tkind) THEN
601 596 : CALL group%sum(md_ener%ekin_kind)
602 596 : 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 41085 : CALL group%sum(md_ener%ekin)
608 41085 : CALL group%sum(md_ener%vcom)
609 41085 : CALL group%sum(md_ener%total_mass)
610 164340 : md_ener%vcom = md_ener%vcom/md_ener%total_mass
611 : !
612 : ! Compute the QM/MM kinetic energy
613 :
614 41085 : IF (ASSOCIATED(qmmm_env)) THEN ! conventional QM/MM
615 7928 : DO i = 1, SIZE(qmmm_env%qm%qm_atom_index)
616 6556 : iparticle = qmmm_env%qm%qm_atom_index(i)
617 6556 : 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 7928 : + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
622 : END DO
623 : END IF
624 :
625 41085 : IF (ASSOCIATED(qmmmx_env)) THEN ! doing force mixing
626 12 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%INDICES", i_vals=cur_indices)
627 12 : CALL section_vals_val_get(force_env_section, "QMMM%FORCE_MIXING%RESTART_INFO%LABELS", i_vals=cur_labels)
628 1530 : DO i = 1, SIZE(cur_indices)
629 1530 : IF (cur_labels(i) >= force_mixing_label_QM_dynamics) THEN ! this is a QM atom
630 318 : iparticle = cur_indices(i)
631 318 : 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 318 : + particle_set(iparticle)%v(3)*particle_set(iparticle)%v(3))
636 : END IF
637 : END DO
638 : END IF
639 :
640 41085 : IF (ASSOCIATED(qmmm_env) .AND. ASSOCIATED(qmmmx_env)) &
641 0 : CPABORT("get_part_ke: qmmm bug")
642 41085 : END SUBROUTINE get_part_ke
643 :
644 : ! **************************************************************************************************
645 :
646 : END MODULE md_conserved_quantities
|
