(git:374b731)
Loading...
Searching...
No Matches
md_energies.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 prints all energy info per timestep to the screen or to
10!> user defined output files
11!> \author Joost VandeVondele (copy from md_fist_energies)
12!>
13!> \par History
14!> - New MD data are appended to the old data (15.09.2003,MK)
15! **************************************************************************************************
24 USE cell_types, ONLY: cell_type,&
29 USE cp_output_handling, ONLY: cp_p_file,&
57 USE kinds, ONLY: default_string_length,&
58 dp,&
59 int_8
60 USE machine, ONLY: m_flush,&
61 m_memory,&
65 USE md_ener_types, ONLY: md_ener_type,&
76 USE physcon, ONLY: angstrom,&
78 kelvin
79 USE qmmm_types, ONLY: qmmm_env_type
83 USE simpar_types, ONLY: simpar_type
88 USE virial_types, ONLY: virial_type
89#include "../base/base_uses.f90"
90
91 IMPLICIT NONE
92
93 PRIVATE
94
95 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_energies'
96
97 PUBLIC :: initialize_md_ener, &
98 md_energy, &
102
103CONTAINS
104
105! **************************************************************************************************
106!> \brief ...
107!> \param md_ener ...
108!> \param force_env ...
109!> \param simpar ...
110!> \par History
111!> - 10-2007 created
112!> \author MI
113! **************************************************************************************************
114 SUBROUTINE initialize_md_ener(md_ener, force_env, simpar)
115
116 TYPE(md_ener_type), POINTER :: md_ener
117 TYPE(force_env_type), POINTER :: force_env
118 TYPE(simpar_type), POINTER :: simpar
119
120 INTEGER :: nkind
121 LOGICAL :: shell_adiabatic
122 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
123 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
124 TYPE(cp_subsys_type), POINTER :: subsys
125 TYPE(particle_list_type), POINTER :: particles, shell_particles
126
127 NULLIFY (subsys)
128 NULLIFY (atomic_kinds, atomic_kind_set, particles, shell_particles)
129
130 cpassert(ASSOCIATED(md_ener))
131 cpassert(ASSOCIATED(force_env))
132
133 CALL force_env_get(force_env, subsys=subsys)
134 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
135 shell_particles=shell_particles)
136 atomic_kind_set => atomic_kinds%els
137 nkind = SIZE(atomic_kind_set)
138 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
139 shell_adiabatic=shell_adiabatic)
140
141 md_ener%nfree = simpar%nfree
142 md_ener%nfree_shell = -huge(0)
143
144 IF (shell_adiabatic) THEN
145 md_ener%nfree_shell = 3*(shell_particles%n_els)
146 END IF
147
148 IF (simpar%temperature_per_kind) THEN
149 ALLOCATE (md_ener%temp_kind(nkind))
150 ALLOCATE (md_ener%ekin_kind(nkind))
151 ALLOCATE (md_ener%nfree_kind(nkind))
152 md_ener%nfree_kind = 0
153
154 IF (shell_adiabatic) THEN
155 ALLOCATE (md_ener%temp_shell_kind(nkind))
156 ALLOCATE (md_ener%ekin_shell_kind(nkind))
157 ALLOCATE (md_ener%nfree_shell_kind(nkind))
158 md_ener%nfree_shell_kind = 0
159 END IF
160
161 END IF
162 CALL zero_md_ener(md_ener, tkind=simpar%temperature_per_kind, &
163 tshell=shell_adiabatic)
164 md_ener%epot = 0.0_dp
165
166 END SUBROUTINE initialize_md_ener
167
168! **************************************************************************************************
169!> \brief ...
170!> \param md_env ...
171!> \param md_ener ...
172!> \par History
173!> - 10-2007 created
174!> \author MI
175! **************************************************************************************************
176 SUBROUTINE md_energy(md_env, md_ener)
177
178 TYPE(md_environment_type), POINTER :: md_env
179 TYPE(md_ener_type), POINTER :: md_ener
180
181 INTEGER :: natom
182 LOGICAL :: shell_adiabatic, tkind, tshell
183 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
184 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
185 TYPE(cp_subsys_type), POINTER :: subsys
186 TYPE(force_env_type), POINTER :: force_env
187 TYPE(particle_list_type), POINTER :: particles
188 TYPE(simpar_type), POINTER :: simpar
189
190 NULLIFY (atomic_kinds, atomic_kind_set, force_env, &
191 particles, subsys, simpar)
192 CALL get_md_env(md_env=md_env, force_env=force_env, &
193 simpar=simpar)
194
195 CALL force_env_get(force_env, &
196 potential_energy=md_ener%epot, subsys=subsys)
197
198 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
199 atomic_kind_set => atomic_kinds%els
200 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
201 shell_adiabatic=shell_adiabatic)
202
203 tkind = simpar%temperature_per_kind
204 tshell = shell_adiabatic
205
206 CALL cp_subsys_get(subsys, particles=particles)
207 natom = particles%n_els
208
209 CALL compute_conserved_quantity(md_env, md_ener, tkind=tkind, &
210 tshell=tshell, natom=natom)
211
212 END SUBROUTINE md_energy
213
214! **************************************************************************************************
215!> \brief ...
216!> \param md_env ...
217!> \param md_ener ...
218!> \par History
219!> - 10.2007 created
220!> \author MI
221! **************************************************************************************************
222 SUBROUTINE md_ener_reftraj(md_env, md_ener)
223 TYPE(md_environment_type), POINTER :: md_env
224 TYPE(md_ener_type), POINTER :: md_ener
225
226 TYPE(force_env_type), POINTER :: force_env
227 TYPE(reftraj_type), POINTER :: reftraj
228
229 CALL zero_md_ener(md_ener, tkind=.false., tshell=.false.)
230 CALL get_md_env(md_env=md_env, force_env=force_env, reftraj=reftraj)
231
232 IF (reftraj%info%eval /= reftraj_eval_none) THEN
233 CALL force_env_get(force_env, potential_energy=md_ener%epot)
234 ELSE
235 md_ener%epot = reftraj%epot
236 md_ener%delta_epot = (reftraj%epot - reftraj%epot0)/real(reftraj%natom, kind=dp)*kelvin
237 END IF
238
239 END SUBROUTINE md_ener_reftraj
240
241! **************************************************************************************************
242!> \brief This routine computes the conserved quantity, temperature
243!> and things like that and prints them out
244!> \param md_env ...
245!> \par History
246!> - New MD data are appended to the old data (15.09.2003,MK)
247!> - 02.2008 - Teodoro Laino [tlaino] - University of Zurich
248!> Cleaning code and collecting the many commons routines..
249!> \author CJM
250! **************************************************************************************************
251 SUBROUTINE md_write_output(md_env)
252
253 TYPE(md_environment_type), POINTER :: md_env
254
255 CHARACTER(len=*), PARAMETER :: routinen = 'md_write_output'
256
257 CHARACTER(LEN=default_string_length) :: fmd, my_act, my_pos
258 INTEGER :: ene, ener_mix, handle, i, nat, nkind, &
259 shene, tempkind, trsl
260 INTEGER(KIND=int_8) :: max_memory
261 INTEGER, POINTER :: itimes
262 LOGICAL :: init, is_mixed, new_file, print_memory, &
263 qmmm, shell_adiabatic, shell_present
264 REAL(dp) :: abc(3), cell_angle(3), dt, econs, &
265 pv_scalar, pv_xx, pv_xx_nc
266 REAL(kind=dp) :: harm_shell, hugoniot
267 REAL(kind=dp), POINTER :: time, used_time
268 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
269 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
270 TYPE(average_quantities_type), POINTER :: averages
271 TYPE(barostat_type), POINTER :: barostat
272 TYPE(cell_type), POINTER :: cell
273 TYPE(cp_logger_type), POINTER :: logger
274 TYPE(cp_subsys_type), POINTER :: subsys
275 TYPE(force_env_type), POINTER :: force_env
276 TYPE(md_ener_type), POINTER :: md_ener
277 TYPE(mp_para_env_type), POINTER :: para_env
278 TYPE(particle_list_type), POINTER :: core_particles, particles, &
279 shell_particles
280 TYPE(qmmm_env_type), POINTER :: qmmm_env
281 TYPE(reftraj_type), POINTER :: reftraj
282 TYPE(section_vals_type), POINTER :: motion_section, print_key, root_section
283 TYPE(simpar_type), POINTER :: simpar
284 TYPE(thermal_regions_type), POINTER :: thermal_regions
285 TYPE(thermostats_type), POINTER :: thermostats
286 TYPE(virial_type), POINTER :: virial
287
288 NULLIFY (logger)
289 logger => cp_get_default_logger()
290 CALL timeset(routinen, handle)
291
292 ! Zeroing
293 hugoniot = 0.0_dp
294 econs = 0.0_dp
295 shell_adiabatic = .false.
296 shell_present = .false.
297 NULLIFY (motion_section, atomic_kinds, atomic_kind_set, cell, subsys, &
298 force_env, md_ener, qmmm_env, reftraj, core_particles, particles, &
299 shell_particles, print_key, root_section, simpar, virial, &
300 thermostats, thermal_regions)
301
302 CALL get_md_env(md_env=md_env, itimes=itimes, t=time, used_time=used_time, &
303 simpar=simpar, force_env=force_env, init=init, md_ener=md_ener, &
304 reftraj=reftraj, thermostats=thermostats, barostat=barostat, &
305 para_env=para_env, averages=averages, thermal_regions=thermal_regions)
306
307 root_section => force_env%root_section
308 motion_section => section_vals_get_subs_vals(root_section, "MOTION")
309
310 CALL force_env_get(force_env, cell=cell, subsys=subsys, qmmm_env=qmmm_env)
311
312 qmmm = calc_nfree_qm(md_env, md_ener) > 0
313 is_mixed = (force_env%in_use == use_mixed_force)
314
315 CALL cp_subsys_get(subsys, particles=particles, virial=virial)
316 nat = particles%n_els
317 dt = simpar%dt*simpar%dt_fact
318
319 ! Computing the scalar pressure
320 IF (virial%pv_availability) THEN
321 pv_scalar = 0._dp
322 DO i = 1, 3
323 pv_scalar = pv_scalar + virial%pv_total(i, i)
324 END DO
325 pv_scalar = pv_scalar/3._dp/cell%deth
326 pv_scalar = cp_unit_from_cp2k(pv_scalar, "bar")
327 pv_xx_nc = virial%pv_total(1, 1)/cell%deth
328 pv_xx = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
329 END IF
330
331 CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
332 atomic_kind_set => atomic_kinds%els
333 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
334 shell_present=shell_present, &
335 shell_adiabatic=shell_adiabatic)
336
337 CALL get_cell(cell, abc=abc, alpha=cell_angle(3), beta=cell_angle(2), gamma=cell_angle(1))
338
339 ! Determine POS and ACT for I/O
340 my_pos = "APPEND"
341 my_act = "WRITE"
342 IF (init .AND. (itimes == 0)) THEN
343 my_pos = "REWIND"
344 my_act = "WRITE"
345 END IF
346
347 ! In case of REFTRAJ ensemble the POS is determined differently..
348 ! according the REFTRAJ counter
349 IF (simpar%ensemble == reftraj_ensemble) THEN
350 IF ((reftraj%isnap == reftraj%info%first_snapshot)) THEN
351 my_pos = "REWIND"
352 END IF
353 END IF
354
355 ! Performing protocol relevant to the first step of an MD run
356 IF (init) THEN
357 ! Computing the Hugoniot for NPH calculations
358 IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
359 simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
360 IF (simpar%e0 == 0._dp) simpar%e0 = md_ener%epot + md_ener%ekin
361 hugoniot = md_ener%epot + md_ener%ekin - simpar%e0 - 0.5_dp*(pv_xx_nc + simpar%p0)* &
362 (simpar%v0 - cell%deth)
363 END IF
364
365 IF (simpar%ensemble == reftraj_ensemble) reftraj%init = init
366 ELSE
367 ! Performing protocol for anything beyond the first step of MD
368 IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
369 hugoniot = md_ener%epot + md_ener%ekin - simpar%e0 - 0.5_dp*(pv_xx_nc + simpar%p0)* &
370 (simpar%v0 - cell%deth)
371 END IF
372
373 IF (simpar%ensemble == reftraj_ensemble) THEN
374 time = reftraj%time
375 econs = md_ener%delta_epot
376 itimes = reftraj%itimes
377 ELSE
378 econs = md_ener%delta_cons
379 END IF
380
381 ! Compute average quantities
382 CALL compute_averages(averages, force_env, md_ener, cell, virial, pv_scalar, &
383 pv_xx, used_time, hugoniot, abc, cell_angle, nat, itimes, time, my_pos, &
384 my_act)
385 END IF
386
387 ! Sample memory, if requested
388 CALL section_vals_val_get(motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
389 max_memory = 0
390 IF (print_memory) THEN
391 max_memory = sample_memory(para_env)
392 END IF
393
394 ! Print md information
395 CALL md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, abc, &
396 cell_angle, itimes, dt, time, used_time, averages, econs, pv_scalar, pv_xx, &
397 hugoniot, nat, init, logger, motion_section, my_pos, my_act, max_memory)
398
399 ! Real Output driven by the PRINT sections
400 IF ((.NOT. init) .OR. (itimes == 0) .OR. simpar%ensemble == reftraj_ensemble) THEN
401 ! Print Energy
402 ene = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%ENERGY", &
403 extension=".ener", file_position=my_pos, file_action=my_act, is_new_file=new_file)
404 IF (ene > 0) THEN
405 IF (new_file) THEN
406 ! Please change also the corresponding format explanation below
407 ! keep the constant of motion the true constant of motion !
408 WRITE (ene, '("#",5X,A,10X,A,8X,A,10X,A,12X,A,2(8X,A))') "Step Nr.", "Time[fs]", "Kin.[a.u.]", "Temp[K]", &
409 "Pot.[a.u.]", "Cons Qty[a.u.]", "UsedTime[s]"
410 END IF
411 WRITE (ene, "(I10,F20.6,F20.9,F20.9,F20.9,F20.9,F20.9)") &
412 itimes, time*femtoseconds, md_ener%ekin, md_ener%temp_part, md_ener%epot, md_ener%constant, used_time
413 CALL m_flush(ene)
414 END IF
415 CALL cp_print_key_finished_output(ene, logger, motion_section, "MD%PRINT%ENERGY")
416
417 ! Possibly Print MIXED Energy
418 IF (is_mixed) THEN
419 ener_mix = cp_print_key_unit_nr(logger, motion_section, "PRINT%MIXED_ENERGIES", &
420 extension=".ener", file_position=my_pos, file_action=my_act, &
421 middle_name="mix")
422 IF (ener_mix > 0) THEN
423 WRITE (ener_mix, "(I8,F12.3,F20.9,"//cp_to_string(SIZE(force_env%mixed_env%energies))//"F20.9,F20.9)") &
424 itimes, time*femtoseconds, md_ener%epot, force_env%mixed_env%energies, md_ener%constant
425 CALL m_flush(ener_mix)
426 END IF
427 CALL cp_print_key_finished_output(ener_mix, logger, motion_section, "PRINT%MIXED_ENERGIES")
428 END IF
429
430 ! Print QMMM translation vector if requested
431 IF (qmmm) THEN
432 trsl = cp_print_key_unit_nr(logger, motion_section, "PRINT%TRANSLATION_VECTOR", &
433 extension=".translation", middle_name="qmmm")
434 IF (trsl > 0) THEN
435 WRITE (trsl, '(I10,3F15.10)') itimes, qmmm_env%qm%transl_v
436 END IF
437 CALL cp_print_key_finished_output(trsl, logger, motion_section, &
438 "PRINT%TRANSLATION_VECTOR")
439 END IF
440
441 ! Write Structure data
442 CALL write_structure_data(particles%els, cell, motion_section)
443
444 ! Print Coordinates
445 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
446 pos=my_pos, act=my_act, extended_xmol_title=.true.)
447
448 ! Print Velocities
449 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
450 "VELOCITIES", my_pos, my_act, middle_name="vel", extended_xmol_title=.true.)
451
452 ! Print Force
453 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
454 "FORCES", my_pos, my_act, middle_name="frc", extended_xmol_title=.true.)
455
456 ! Print Force-Mixing labels
457 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
458 "FORCE_MIXING_LABELS", my_pos, my_act, middle_name="fmlabels", extended_xmol_title=.true.)
459
460 ! Print Simulation Cell
461 CALL write_simulation_cell(cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
462
463 ! Print Thermostats status
464 CALL print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
465
466 ! Print Barostat status
467 CALL print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time)
468
469 ! Print Stress Tensor
470 CALL write_stress_tensor(virial, cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
471
472 ! Print Polarisability Tensor
473 IF (ASSOCIATED(force_env%qs_env)) THEN
474 CALL write_polarisability_tensor(force_env, motion_section, itimes, time*femtoseconds, my_pos, my_act)
475 END IF
476
477 ! Temperature per Kinds
478 IF (simpar%temperature_per_kind) THEN
479 tempkind = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%TEMP_KIND", &
480 extension=".temp", file_position=my_pos, file_action=my_act)
481 IF (tempkind > 0) THEN
482 nkind = SIZE(md_ener%temp_kind)
483 fmd = "(I10,F20.3,"//trim(adjustl(cp_to_string(nkind)))//"F20.9)"
484 fmd = trim(fmd)
485 WRITE (tempkind, fmd) itimes, time*femtoseconds, md_ener%temp_kind(1:nkind)
486 CALL m_flush(tempkind)
487 END IF
488 CALL cp_print_key_finished_output(tempkind, logger, motion_section, "MD%PRINT%TEMP_KIND")
489 ELSE
490 print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_KIND")
491 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) &
492 CALL cp_warn(__location__, &
493 "The print_key MD%PRINT%TEMP_KIND has been activated but the "// &
494 "calculation of the temperature per kind has not been requested. "// &
495 "Please switch on the keyword MD%TEMP_KIND.")
496 END IF
497 !Thermal Region
498 CALL print_thermal_regions_temperature(thermal_regions, itimes, time*femtoseconds, my_pos, my_act)
499
500 ! Core/Shell Model
501 IF (shell_present) THEN
502 CALL force_env_get(force_env, harmonic_shell=harm_shell)
503 CALL cp_subsys_get(subsys, shell_particles=shell_particles, core_particles=core_particles)
504
505 ! Print Shell Energy
506 shene = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%SHELL_ENERGY", &
507 extension=".shener", file_position=my_pos, file_action=my_act, &
508 file_form="FORMATTED", is_new_file=new_file)
509 IF (shene > 0) THEN
510 IF (new_file) THEN
511 WRITE (shene, '("#",3X,A,3X,A,3X,3(5X,A,5X))') "Step Nr.", "Time[fs]", "Kin.[a.u.]", &
512 "Temp.[K]", "Pot.[a.u.]"
513 END IF
514
515 WRITE (shene, "(I8,F12.3,F20.9,F20.9,F20.9,F20.9 )") &
516 itimes, time*femtoseconds, md_ener%ekin_shell, md_ener%temp_shell, harm_shell
517 CALL m_flush(shene)
518 END IF
519 CALL cp_print_key_finished_output(shene, logger, motion_section, "MD%PRINT%SHELL_ENERGY")
520
521 ! Print Shell Coordinates
522 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
523 "SHELL_TRAJECTORY", my_pos, my_act, "shpos", shell_particles, extended_xmol_title=.true.)
524
525 IF (shell_adiabatic) THEN
526 ! Print Shell Velocities
527 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
528 "SHELL_VELOCITIES", my_pos, my_act, "shvel", shell_particles, extended_xmol_title=.true.)
529
530 ! Print Shell Forces
531 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
532 "SHELL_FORCES", my_pos, my_act, "shfrc", shell_particles, extended_xmol_title=.true.)
533
534 ! Print Core Coordinates
535 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
536 "CORE_TRAJECTORY", my_pos, my_act, "copos", core_particles, extended_xmol_title=.true.)
537
538 ! Print Core Velocities
539 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
540 "CORE_VELOCITIES", my_pos, my_act, "covel", core_particles, extended_xmol_title=.true.)
541
542 ! Print Core Forces
543 CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
544 "CORE_FORCES", my_pos, my_act, "cofrc", core_particles, extended_xmol_title=.true.)
545
546 ! Temperature per Kinds
547 IF (simpar%temperature_per_kind) THEN
548 tempkind = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%TEMP_SHELL_KIND", &
549 extension=".shtemp", file_position=my_pos, file_action=my_act)
550 IF (tempkind > 0) THEN
551 nkind = SIZE(md_ener%temp_shell_kind)
552 fmd = "(I10,F20.3,"//trim(adjustl(cp_to_string(nkind)))//"F20.9)"
553 fmd = trim(fmd)
554 WRITE (tempkind, fmd) itimes, time*femtoseconds, md_ener%temp_shell_kind(1:nkind)
555 CALL m_flush(tempkind)
556 END IF
557 CALL cp_print_key_finished_output(tempkind, logger, motion_section, &
558 "MD%PRINT%TEMP_SHELL_KIND")
559 ELSE
560 print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_SHELL_KIND")
561 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) &
562 CALL cp_warn(__location__, &
563 "The print_key MD%PRINT%TEMP_SHELL_KIND has been activated but the "// &
564 "calculation of the temperature per kind has not been requested. "// &
565 "Please switch on the keyword MD%TEMP_KIND.")
566 END IF
567 END IF
568 END IF
569 END IF
570 init = .false.
571 CALL set_md_env(md_env, init=init)
572 CALL timestop(handle)
573 END SUBROUTINE md_write_output
574
575! **************************************************************************************************
576!> \brief This routine prints all basic information during MD steps
577!> \param simpar ...
578!> \param md_ener ...
579!> \param qmmm ...
580!> \param virial ...
581!> \param reftraj ...
582!> \param cell ...
583!> \param abc ...
584!> \param cell_angle ...
585!> \param itimes ...
586!> \param dt ...
587!> \param time ...
588!> \param used_time ...
589!> \param averages ...
590!> \param econs ...
591!> \param pv_scalar ...
592!> \param pv_xx ...
593!> \param hugoniot ...
594!> \param nat ...
595!> \param init ...
596!> \param logger ...
597!> \param motion_section ...
598!> \param my_pos ...
599!> \param my_act ...
600!> \param max_memory ...
601!> \par History
602!> - 10.2008 - Teodoro Laino [tlaino] - University of Zurich
603!> Refactoring: split into an independent routine.
604!> All output on screen must be included in this function!
605!> \author CJM
606! **************************************************************************************************
607 SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
608 abc, cell_angle, itimes, dt, time, used_time, averages, econs, pv_scalar, &
609 pv_xx, hugoniot, nat, init, logger, motion_section, my_pos, my_act, &
610 max_memory)
611
612 TYPE(simpar_type), POINTER :: simpar
613 TYPE(md_ener_type), POINTER :: md_ener
614 LOGICAL, INTENT(IN) :: qmmm
615 TYPE(virial_type), POINTER :: virial
616 TYPE(reftraj_type), POINTER :: reftraj
617 TYPE(cell_type), POINTER :: cell
618 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: abc, cell_angle
619 INTEGER, POINTER :: itimes
620 REAL(kind=dp), INTENT(IN) :: dt
621 REAL(kind=dp), POINTER :: time, used_time
622 TYPE(average_quantities_type), POINTER :: averages
623 REAL(kind=dp), INTENT(IN) :: econs, pv_scalar, pv_xx, hugoniot
624 INTEGER, INTENT(IN) :: nat
625 LOGICAL, INTENT(IN) :: init
626 TYPE(cp_logger_type), POINTER :: logger
627 TYPE(section_vals_type), POINTER :: motion_section
628 CHARACTER(LEN=default_string_length), INTENT(IN) :: my_pos, my_act
629 INTEGER(KIND=int_8), INTENT(IN) :: max_memory
630
631 INTEGER :: iw
632 TYPE(enumeration_type), POINTER :: enum
633 TYPE(keyword_type), POINTER :: keyword
634 TYPE(section_type), POINTER :: section
635
636 NULLIFY (enum, keyword, section)
637 ! Print to the screen info about MD
638 iw = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%PROGRAM_RUN_INFO", &
639 extension=".mdLog", file_position=my_pos, file_action=my_act)
640
641 ! Performing protocol relevant to the first step of an MD run
642 IF (iw > 0) THEN
643 CALL create_md_section(section)
644 keyword => section_get_keyword(section, "ENSEMBLE")
645 CALL keyword_get(keyword, enum=enum)
646 IF (init) THEN
647 ! Write initial values of quantities of interest
648 WRITE (iw, '(/,T2,A)') 'MD_INI| MD initialization'
649 WRITE (iw, '(T2,A,T61,E20.12)') &
650 'MD_INI| Potential energy [hartree]', md_ener%epot
651 IF (simpar%ensemble /= reftraj_ensemble) THEN
652 IF (.NOT. qmmm) THEN
653 ! NO QM/MM info
654 WRITE (iw, '(T2,A,T61,E20.12)') &
655 'MD_INI| Kinetic energy [hartree]', md_ener%ekin
656 WRITE (iw, '(T2,A,T61,F20.6)') &
657 'MD_INI| Temperature [K]', md_ener%temp_part
658 ELSE
659 WRITE (iw, '(T2,A,T61,E20.12)') &
660 'MD_INI| Total kinetic energy [hartree]', md_ener%ekin, &
661 'MD_INI| QM kinetic energy [hartree]', md_ener%ekin_qm
662 WRITE (iw, '(T2,A,T61,F20.6)') &
663 'MD_INI| Total temperature [K]', md_ener%temp_part, &
664 'MD_INI| QM temperature [K]', md_ener%temp_qm
665 END IF
666 END IF
667 IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
668 (simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
669 (simpar%ensemble == npt_i_ensemble) .OR. &
670 (simpar%ensemble == npt_ia_ensemble) .OR. &
671 (simpar%ensemble == npt_f_ensemble) .OR. &
672 (simpar%ensemble == npe_i_ensemble) .OR. &
673 (simpar%ensemble == npe_f_ensemble)) THEN
674 WRITE (iw, '(T2,A,T61,F20.6)') &
675 'MD_INI| Barostat temperature [K]', md_ener%temp_baro
676 END IF
677 IF (virial%pv_availability) THEN
678 WRITE (iw, '(T2,A,T61,ES20.12)') &
679 'MD_INI| Pressure [bar]', pv_scalar
680 END IF
681 IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
682 (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
683 WRITE (iw, '(T2,A,T61,ES20.12)') &
684 'MD_INI| Hugoniot constraint [K]', hugoniot
685 END IF
686 IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
687 (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
688 WRITE (iw, '(T2,A,T61,E20.12)') &
689 'MD_INI| Total energy [hartree]', simpar%e0
690 END IF
691 WRITE (iw, '(T2,A,T61,ES20.12)') &
692 'MD_INI| Cell volume [bohr^3]', cell%deth
693 WRITE (iw, '(T2,A,T61,ES20.12)') &
694 'MD_INI| Cell volume [ang^3]', cell%deth*angstrom**3
695 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
696 'MD_INI| Cell lengths [bohr]', abc(1:3)
697 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
698 'MD_INI| Cell lengths [ang]', abc(1:3)*angstrom
699 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
700 'MD_INI| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
701 ELSE
702 ! Write seuquential values of quantities of interest
703 WRITE (iw, '(/,T2,A)') 'MD| '//repeat('*', 75)
704!MK WRITE (iw, '(T2,A,T61,A20)') &
705!MK 'MD| Ensemble type', ADJUSTR(TRIM(enum_i2c(enum, simpar%ensemble)))
706 WRITE (iw, '(T2,A,T71,I10)') &
707 'MD| Step number', itimes
708 IF (simpar%variable_dt) THEN
709 WRITE (iw, '(T2,A,T61,F20.6)') &
710 'MD| Time step [fs]', dt*femtoseconds
711 END IF
712 WRITE (iw, '(T2,A,T61,F20.6)') &
713 'MD| Time [fs]', time*femtoseconds
714 WRITE (iw, '(T2,A,T61,E20.12)') &
715 'MD| Conserved quantity [hartree]', md_ener%constant
716 WRITE (iw, '(T2,A)') 'MD| '//repeat('-', 75)
717 WRITE (iw, '(T2,A,T47,A,T73,A)') 'MD|', 'Instantaneous', 'Averages'
718 WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
719 'MD| CPU time per MD step [s]', used_time, averages%avecpu
720 WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
721 'MD| Energy drift per atom [K] ', econs, averages%econs
722 WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
723 'MD| Potential energy [hartree]', md_ener%epot, averages%avepot
724 IF (simpar%ensemble /= reftraj_ensemble) THEN
725 IF (.NOT. qmmm) THEN
726 ! No QM/MM info
727 WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
728 'MD| Kinetic energy [hartree]', md_ener%ekin, averages%avekin
729 WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
730 'MD| Temperature [K]', md_ener%temp_part, averages%avetemp
731 ELSE
732 WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
733 'MD| Total kinetic energy [hartree]', md_ener%ekin, averages%avekin
734 WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
735 'MD| QM kinetic energy [hartree]', md_ener%ekin_qm, averages%avekin_qm
736 WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
737 'MD| Total temperature [K]', md_ener%temp_part, averages%avetemp
738 WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
739 'MD| QM temperature [K]', md_ener%temp_qm, averages%avetemp_qm
740 END IF
741 END IF
742 IF (virial%pv_availability) THEN
743 WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
744 'MD| Pressure [bar]', pv_scalar, averages%avepress
745 END IF
746 IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
747 (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
748 WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
749 'MD| P(xx) [bar]', pv_xx, averages%avepxx
750 WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
751 'MD| Hugoniot [K]', hugoniot/3.0_dp/nat*kelvin, averages%avehugoniot/3.0_dp/nat*kelvin
752 END IF
753 IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
754 (simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
755 (simpar%ensemble == npt_i_ensemble) .OR. &
756 (simpar%ensemble == npt_ia_ensemble) .OR. &
757 (simpar%ensemble == npt_f_ensemble) .OR. &
758 (simpar%ensemble == npe_i_ensemble) .OR. &
759 (simpar%ensemble == npe_f_ensemble)) THEN
760 WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
761 'MD| Barostat temperature [K]', md_ener%temp_baro, averages%avetemp_baro
762 WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
763 'MD| Cell volume [bohr^3]', cell%deth, averages%avevol
764 WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
765 'MD| Cell volume [ang^3]', cell%deth*angstrom**3, averages%avevol*angstrom**3
766 WRITE (iw, '(T2,A)') 'MD| '//repeat('-', 75)
767 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
768 'MD| Cell lengths [bohr]', abc(1:3)
769 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
770 'MD| Cell lengths [ang]', abc(1:3)*angstrom
771 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
772 'MD| Average cell lengths [bohr]', averages%aveca, averages%avecb, averages%avecc
773 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
774 'MD| Average cell lengths [ang]', averages%aveca*angstrom, averages%avecb*angstrom, &
775 averages%avecc*angstrom
776 END IF
777 IF ((simpar%ensemble == npt_f_ensemble) .OR. &
778 (simpar%ensemble == npe_f_ensemble)) THEN
779 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
780 'MD| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
781 WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
782 'MD| Average cell angles [deg]', averages%aveal, averages%avebe, averages%avega
783 END IF
784 IF (simpar%ensemble == reftraj_ensemble) THEN
785 IF (reftraj%info%msd) THEN
786 IF (reftraj%msd%msd_kind) THEN
787 WRITE (iw, '(/,T2,A,T51,3F10.5)') &
788 'MD| COM displacement (dx,dy,dz) [bohr]', reftraj%msd%drcom(1:3)
789 END IF
790 END IF
791 END IF
792 WRITE (iw, '(T2,A)') 'MD| '//repeat('*', 75)
793 IF (max_memory .NE. 0) THEN
794 WRITE (iw, '(T2,A,T73,I8)') &
795 'MD| Estimated peak process memory after this step [MiB]', &
796 (max_memory + (1024*1024) - 1)/(1024*1024)
797 END IF
798 END IF
799 END IF
800 CALL section_release(section)
801 CALL cp_print_key_finished_output(iw, logger, motion_section, &
802 "MD%PRINT%PROGRAM_RUN_INFO")
803 END SUBROUTINE md_write_info_low
804
805! **************************************************************************************************
806!> \brief Samples memory usage
807!> \param para_env ...
808!> \return ...
809!> \note based on what is done in start/cp2k_runs.F
810! **************************************************************************************************
811 FUNCTION sample_memory(para_env) RESULT(max_memory)
812 TYPE(mp_para_env_type), POINTER :: para_env
813 INTEGER(KIND=int_8) :: max_memory
814
815 CALL m_memory()
816 max_memory = m_memory_max
817 CALL para_env%max(max_memory)
818
819 END FUNCTION sample_memory
820
821END MODULE md_energies
represent a simple array based list of the given type
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Handles the type to compute averages during an MD.
subroutine, public compute_averages(averages, force_env, md_ener, cell, virial, pv_scalar, pv_xx, used_time, hugoniot, abc, cell_angle, nat, itimes, time, my_pos, my_act)
computes the averages
Barostat structure: module containing barostat available for MD.
Barostat utils.
subroutine, public print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time)
Prints status of barostat during an MD run.
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
types that represent a subsys, i.e. a part of the system
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
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1179
Interface for the force calculations.
integer, parameter, public use_mixed_force
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
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public nph_uniaxial_ensemble
integer, parameter, public npt_i_ensemble
integer, parameter, public nph_uniaxial_damped_ensemble
integer, parameter, public npe_f_ensemble
integer, parameter, public npe_i_ensemble
integer, parameter, public npt_ia_ensemble
integer, parameter, public npt_f_ensemble
integer, parameter, public reftraj_ensemble
subroutine, public create_md_section(section)
...
represents an enumeration, i.e. a mapping between integers and strings
represents keywords in an input
subroutine, public keyword_get(keyword, names, usage, description, type_of_var, n_var, default_value, lone_keyword_value, repeats, enum, citations)
...
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
recursive type(keyword_type) function, pointer, public section_get_keyword(section, keyword_name)
returns the requested keyword
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:332
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:106
integer(kind=int_8), save, public m_memory_max
Definition machine.F:93
computes the conserved quantities for a given md ensemble and also kinetic energies,...
integer function, public calc_nfree_qm(md_env, md_ener)
Calculates the number of QM degress of freedom.
subroutine, public compute_conserved_quantity(md_env, md_ener, tkind, tshell, natom)
calculates conserved quantity.
Split md_ener module from md_environment_type.
subroutine, public zero_md_ener(md_ener, tkind, tshell)
initialize to zero energies and temperatures
prints all energy info per timestep to the screen or to user defined output files
Definition md_energies.F:16
subroutine, public md_ener_reftraj(md_env, md_ener)
...
subroutine, public initialize_md_ener(md_ener, force_env, simpar)
...
subroutine, public md_energy(md_env, md_ener)
...
integer(kind=int_8) function, public sample_memory(para_env)
Samples memory usage.
subroutine, public md_write_output(md_env)
This routine computes the conserved quantity, temperature and things like that and prints them out.
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.
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
Interface to the message passing library MPI.
Output Utilities for MOTION_SECTION.
subroutine, public write_stress_tensor(virial, cell, motion_section, itimes, time, pos, act)
Prints the Stress Tensor.
subroutine, public write_simulation_cell(cell, motion_section, itimes, time, pos, act)
Prints the Simulation Cell.
subroutine, public write_trajectory(force_env, root_section, it, time, dtime, etot, pk_name, pos, act, middle_name, particles, extended_xmol_title)
Prints the information controlled by the TRAJECTORY section.
represent a simple array based list of the given type
Define methods related to particle_type.
subroutine, public write_structure_data(particle_set, cell, input_section)
Write structure data requested by a separate structure data input section to the output unit....
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public femtoseconds
Definition physcon.F:153
real(kind=dp), parameter, public kelvin
Definition physcon.F:165
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
Basic container type for QM/MM.
Definition qmmm_types.F:12
Polarizability calculation by dfpt Initialization of the polar_env, Perturbation Hamiltonian by appli...
subroutine, public write_polarisability_tensor(force_env, motion_section, itimes, time, pos, act)
Prints the polarisability tensor to a file during MD runs.
initialization of the reftraj structure used to analyse previously generated trajectories
integer, parameter, public reftraj_eval_none
Type for storing MD parameters.
Thermal regions type: to initialize and control the temperature of different regions.
Setup of regions with different temperature.
subroutine, public print_thermal_regions_temperature(thermal_regions, itimes, time, pos, act)
print_thermal_regions_temperature
Thermostat structure: module containing thermostat available for MD.
Utilities for thermostats.
subroutine, public print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
Prints status of all thermostats during an MD run.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represents a system: atoms, molecules, their pos,vel,...
wrapper to abstract the force evaluation of the various methods
represent a keyword in the input
represent a section of the input file
stores all the informations relevant to an mpi environment
Simulation parameter type for molecular dynamics.