Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : ! **************************************************************************************************
16 : MODULE md_energies
17 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
18 : USE atomic_kind_types, ONLY: atomic_kind_type,&
19 : get_atomic_kind_set
20 : USE averages_types, ONLY: average_quantities_type,&
21 : compute_averages
22 : USE barostat_types, ONLY: barostat_type
23 : USE barostat_utils, ONLY: print_barostat_status
24 : USE cell_types, ONLY: cell_type,&
25 : get_cell
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_type,&
28 : cp_to_string
29 : USE cp_output_handling, ONLY: cp_p_file,&
30 : cp_print_key_finished_output,&
31 : cp_print_key_should_output,&
32 : cp_print_key_unit_nr
33 : USE cp_subsys_types, ONLY: cp_subsys_get,&
34 : cp_subsys_type
35 : USE cp_units, ONLY: cp_unit_from_cp2k
36 : USE force_env_types, ONLY: force_env_get,&
37 : force_env_type,&
38 : use_mixed_force
39 : USE input_constants, ONLY: npe_f_ensemble,&
40 : npe_i_ensemble,&
41 : nph_uniaxial_damped_ensemble,&
42 : nph_uniaxial_ensemble,&
43 : npt_f_ensemble,&
44 : npt_i_ensemble,&
45 : npt_ia_ensemble,&
46 : reftraj_ensemble
47 : USE input_cp2k_md, ONLY: create_md_section
48 : USE input_enumeration_types, ONLY: enumeration_type
49 : USE input_keyword_types, ONLY: keyword_get,&
50 : keyword_type
51 : USE input_section_types, ONLY: section_get_keyword,&
52 : section_release,&
53 : section_type,&
54 : section_vals_get_subs_vals,&
55 : section_vals_type,&
56 : section_vals_val_get
57 : USE kinds, ONLY: default_string_length,&
58 : dp,&
59 : int_8
60 : USE machine, ONLY: m_flush,&
61 : m_memory,&
62 : m_memory_max
63 : USE md_conserved_quantities, ONLY: calc_nfree_qm,&
64 : compute_conserved_quantity
65 : USE md_ener_types, ONLY: md_ener_type,&
66 : zero_md_ener
67 : USE md_environment_types, ONLY: get_md_env,&
68 : md_environment_type,&
69 : set_md_env
70 : USE message_passing, ONLY: mp_para_env_type
71 : USE motion_utils, ONLY: write_simulation_cell,&
72 : write_stress_tensor_to_file,&
73 : write_trajectory
74 : USE particle_list_types, ONLY: particle_list_type
75 : USE particle_methods, ONLY: write_structure_data
76 : USE physcon, ONLY: angstrom,&
77 : femtoseconds,&
78 : kelvin
79 : USE qmmm_types, ONLY: qmmm_env_type
80 : USE qs_linres_polar_utils, ONLY: write_polarisability_tensor
81 : USE reftraj_types, ONLY: REFTRAJ_EVAL_NONE,&
82 : reftraj_type
83 : USE simpar_types, ONLY: simpar_type
84 : USE thermal_region_types, ONLY: thermal_regions_type
85 : USE thermal_region_utils, ONLY: print_thermal_regions_temperature
86 : USE thermostat_types, ONLY: thermostats_type
87 : USE thermostat_utils, ONLY: print_thermostats_status
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, &
99 : md_ener_reftraj, &
100 : md_write_output, &
101 : sample_memory
102 :
103 : CONTAINS
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 3544 : 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 1772 : 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 1772 : NULLIFY (subsys)
128 1772 : NULLIFY (atomic_kinds, atomic_kind_set, particles, shell_particles)
129 :
130 0 : CPASSERT(ASSOCIATED(md_ener))
131 1772 : CPASSERT(ASSOCIATED(force_env))
132 :
133 1772 : CALL force_env_get(force_env, subsys=subsys)
134 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
135 1772 : shell_particles=shell_particles)
136 1772 : atomic_kind_set => atomic_kinds%els
137 1772 : nkind = SIZE(atomic_kind_set)
138 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
139 1772 : shell_adiabatic=shell_adiabatic)
140 :
141 1772 : md_ener%nfree = simpar%nfree
142 1772 : md_ener%nfree_shell = -HUGE(0)
143 :
144 1772 : IF (shell_adiabatic) THEN
145 132 : md_ener%nfree_shell = 3*(shell_particles%n_els)
146 : END IF
147 :
148 1772 : IF (simpar%temperature_per_kind) THEN
149 108 : ALLOCATE (md_ener%temp_kind(nkind))
150 72 : ALLOCATE (md_ener%ekin_kind(nkind))
151 108 : ALLOCATE (md_ener%nfree_kind(nkind))
152 110 : md_ener%nfree_kind = 0
153 :
154 36 : IF (shell_adiabatic) THEN
155 36 : ALLOCATE (md_ener%temp_shell_kind(nkind))
156 36 : ALLOCATE (md_ener%ekin_shell_kind(nkind))
157 36 : ALLOCATE (md_ener%nfree_shell_kind(nkind))
158 54 : 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 1772 : tshell=shell_adiabatic)
164 1772 : md_ener%epot = 0.0_dp
165 :
166 1772 : 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 84758 : 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 42379 : 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 42379 : NULLIFY (atomic_kinds, atomic_kind_set, force_env, &
191 42379 : particles, subsys, simpar)
192 : CALL get_md_env(md_env=md_env, force_env=force_env, &
193 42379 : simpar=simpar)
194 :
195 : CALL force_env_get(force_env, &
196 42379 : potential_energy=md_ener%epot, subsys=subsys)
197 :
198 42379 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
199 42379 : atomic_kind_set => atomic_kinds%els
200 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
201 42379 : shell_adiabatic=shell_adiabatic)
202 :
203 42379 : tkind = simpar%temperature_per_kind
204 42379 : tshell = shell_adiabatic
205 :
206 42379 : CALL cp_subsys_get(subsys, particles=particles)
207 42379 : natom = particles%n_els
208 :
209 : CALL compute_conserved_quantity(md_env, md_ener, tkind=tkind, &
210 42379 : tshell=tshell, natom=natom)
211 :
212 42379 : 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 282 : 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 282 : CALL zero_md_ener(md_ener, tkind=.FALSE., tshell=.FALSE.)
230 282 : CALL get_md_env(md_env=md_env, force_env=force_env, reftraj=reftraj)
231 :
232 282 : IF (reftraj%info%eval /= REFTRAJ_EVAL_NONE) THEN
233 142 : CALL force_env_get(force_env, potential_energy=md_ener%epot)
234 : ELSE
235 140 : md_ener%epot = reftraj%epot
236 140 : md_ener%delta_epot = (reftraj%epot - reftraj%epot0)/REAL(reftraj%natom, kind=dp)*kelvin
237 : END IF
238 :
239 282 : 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 170756 : 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 42689 : 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 42689 : NULLIFY (logger)
289 85378 : logger => cp_get_default_logger()
290 42689 : CALL timeset(routineN, handle)
291 :
292 : ! Zeroing
293 42689 : hugoniot = 0.0_dp
294 42689 : econs = 0.0_dp
295 42689 : shell_adiabatic = .FALSE.
296 42689 : shell_present = .FALSE.
297 42689 : NULLIFY (motion_section, atomic_kinds, atomic_kind_set, cell, subsys, &
298 42689 : force_env, md_ener, qmmm_env, reftraj, core_particles, particles, &
299 42689 : shell_particles, print_key, root_section, simpar, virial, &
300 42689 : 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 42689 : para_env=para_env, averages=averages, thermal_regions=thermal_regions)
306 :
307 42689 : root_section => force_env%root_section
308 42689 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
309 :
310 42689 : CALL force_env_get(force_env, cell=cell, subsys=subsys, qmmm_env=qmmm_env)
311 :
312 42689 : qmmm = calc_nfree_qm(md_env, md_ener) > 0
313 42689 : is_mixed = (force_env%in_use == use_mixed_force)
314 :
315 42689 : CALL cp_subsys_get(subsys, particles=particles, virial=virial)
316 42689 : nat = particles%n_els
317 42689 : dt = simpar%dt*simpar%dt_fact
318 :
319 : ! Computing the scalar pressure
320 42689 : IF (virial%pv_availability) THEN
321 4550 : pv_scalar = 0._dp
322 18200 : DO i = 1, 3
323 18200 : pv_scalar = pv_scalar + virial%pv_total(i, i)
324 : END DO
325 4550 : pv_scalar = pv_scalar/3._dp/cell%deth
326 4550 : pv_scalar = cp_unit_from_cp2k(pv_scalar, "bar")
327 4550 : pv_xx_nc = virial%pv_total(1, 1)/cell%deth
328 4550 : pv_xx = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
329 : END IF
330 :
331 42689 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
332 42689 : atomic_kind_set => atomic_kinds%els
333 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
334 : shell_present=shell_present, &
335 42689 : shell_adiabatic=shell_adiabatic)
336 :
337 42689 : 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 42689 : my_pos = "APPEND"
341 42689 : my_act = "WRITE"
342 42689 : IF (init .AND. (itimes == 0)) THEN
343 1542 : my_pos = "REWIND"
344 1542 : my_act = "WRITE"
345 : END IF
346 :
347 : ! In case of REFTRAJ ensemble the POS is determined differently..
348 : ! according the REFTRAJ counter
349 42689 : IF (simpar%ensemble == reftraj_ensemble) THEN
350 282 : IF ((reftraj%isnap == reftraj%info%first_snapshot)) THEN
351 32 : my_pos = "REWIND"
352 : END IF
353 : END IF
354 :
355 : ! Performing protocol relevant to the first step of an MD run
356 42689 : IF (init) THEN
357 : ! Computing the Hugoniot for NPH calculations
358 1770 : IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
359 : simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
360 6 : 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 6 : (simpar%v0 - cell%deth)
363 : END IF
364 :
365 1770 : IF (simpar%ensemble == reftraj_ensemble) reftraj%init = init
366 : ELSE
367 : ! Performing protocol for anything beyond the first step of MD
368 40919 : 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 60 : (simpar%v0 - cell%deth)
371 : END IF
372 :
373 40919 : IF (simpar%ensemble == reftraj_ensemble) THEN
374 246 : time = reftraj%time
375 246 : econs = md_ener%delta_epot
376 246 : itimes = reftraj%itimes
377 : ELSE
378 40673 : 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 40919 : my_act)
385 : END IF
386 :
387 : ! Sample memory, if requested
388 42689 : CALL section_vals_val_get(motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
389 42689 : max_memory = 0
390 42689 : IF (print_memory) THEN
391 42689 : 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 42689 : hugoniot, nat, init, logger, motion_section, my_pos, my_act, max_memory)
398 :
399 : ! Real Output driven by the PRINT sections
400 42689 : 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 42497 : extension=".ener", file_position=my_pos, file_action=my_act, is_new_file=new_file)
404 42497 : IF (ene > 0) THEN
405 18228 : 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 750 : 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 1500 : "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 18228 : itimes, time*femtoseconds, md_ener%ekin, md_ener%temp_part, md_ener%epot, md_ener%constant, used_time
413 18228 : CALL m_flush(ene)
414 : END IF
415 42497 : CALL cp_print_key_finished_output(ene, logger, motion_section, "MD%PRINT%ENERGY")
416 :
417 : ! Possibly Print MIXED Energy
418 42497 : 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 342 : middle_name="mix")
422 342 : 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 545 : itimes, time*femtoseconds, md_ener%epot, force_env%mixed_env%energies, md_ener%constant
425 171 : CALL m_flush(ener_mix)
426 : END IF
427 342 : 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 42497 : IF (qmmm) THEN
432 : trsl = cp_print_key_unit_nr(logger, motion_section, "PRINT%TRANSLATION_VECTOR", &
433 1430 : extension=".translation", middle_name="qmmm")
434 1430 : IF (trsl > 0) THEN
435 0 : 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 1430 : "PRINT%TRANSLATION_VECTOR")
439 : END IF
440 :
441 : ! Write Structure data
442 42497 : 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 42497 : 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 42497 : "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 42497 : "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 42497 : "FORCE_MIXING_LABELS", my_pos, my_act, middle_name="fmlabels", extended_xmol_title=.TRUE.)
459 :
460 : ! Print Simulation Cell
461 42497 : CALL write_simulation_cell(cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
462 :
463 : ! Print Thermostats status
464 42497 : CALL print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
465 :
466 : ! Print Barostat status
467 42497 : CALL print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time)
468 :
469 : ! Print Stress Tensor
470 42497 : CALL write_stress_tensor_to_file(virial, cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
471 :
472 : ! Print Polarisability Tensor
473 42497 : IF (ASSOCIATED(force_env%qs_env)) THEN
474 3868 : CALL write_polarisability_tensor(force_env, motion_section, itimes, time*femtoseconds, my_pos, my_act)
475 : END IF
476 :
477 : ! Temperature per Kinds
478 42497 : IF (simpar%temperature_per_kind) THEN
479 : tempkind = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%TEMP_KIND", &
480 900 : extension=".temp", file_position=my_pos, file_action=my_act)
481 900 : IF (tempkind > 0) THEN
482 266 : nkind = SIZE(md_ener%temp_kind)
483 266 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nkind)))//"F20.9)"
484 : fmd = TRIM(fmd)
485 798 : WRITE (tempkind, fmd) itimes, time*femtoseconds, md_ener%temp_kind(1:nkind)
486 266 : CALL m_flush(tempkind)
487 : END IF
488 900 : CALL cp_print_key_finished_output(tempkind, logger, motion_section, "MD%PRINT%TEMP_KIND")
489 : ELSE
490 41597 : print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_KIND")
491 41597 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
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 302 : "Please switch on the keyword MD%TEMP_KIND.")
496 : END IF
497 : END IF
498 : !Thermal Region
499 42497 : CALL print_thermal_regions_temperature(thermal_regions, itimes, time*femtoseconds, my_pos, my_act)
500 :
501 : ! Core/Shell Model
502 42497 : IF (shell_present) THEN
503 2386 : CALL force_env_get(force_env, harmonic_shell=harm_shell)
504 2386 : CALL cp_subsys_get(subsys, shell_particles=shell_particles, core_particles=core_particles)
505 :
506 : ! Print Shell Energy
507 : shene = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%SHELL_ENERGY", &
508 : extension=".shener", file_position=my_pos, file_action=my_act, &
509 2386 : file_form="FORMATTED", is_new_file=new_file)
510 2386 : IF (shene > 0) THEN
511 1106 : IF (new_file) THEN
512 45 : WRITE (shene, '("#",3X,A,3X,A,3X,3(5X,A,5X))') "Step Nr.", "Time[fs]", "Kin.[a.u.]", &
513 90 : "Temp.[K]", "Pot.[a.u.]"
514 : END IF
515 :
516 : WRITE (shene, "(I8,F12.3,F20.9,F20.9,F20.9,F20.9 )") &
517 1106 : itimes, time*femtoseconds, md_ener%ekin_shell, md_ener%temp_shell, harm_shell
518 1106 : CALL m_flush(shene)
519 : END IF
520 2386 : CALL cp_print_key_finished_output(shene, logger, motion_section, "MD%PRINT%SHELL_ENERGY")
521 :
522 : ! Print Shell Coordinates
523 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
524 2386 : "SHELL_TRAJECTORY", my_pos, my_act, "shpos", shell_particles, extended_xmol_title=.TRUE.)
525 :
526 2386 : IF (shell_adiabatic) THEN
527 : ! Print Shell Velocities
528 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
529 2386 : "SHELL_VELOCITIES", my_pos, my_act, "shvel", shell_particles, extended_xmol_title=.TRUE.)
530 :
531 : ! Print Shell Forces
532 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
533 2386 : "SHELL_FORCES", my_pos, my_act, "shfrc", shell_particles, extended_xmol_title=.TRUE.)
534 :
535 : ! Print Core Coordinates
536 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
537 2386 : "CORE_TRAJECTORY", my_pos, my_act, "copos", core_particles, extended_xmol_title=.TRUE.)
538 :
539 : ! Print Core Velocities
540 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
541 2386 : "CORE_VELOCITIES", my_pos, my_act, "covel", core_particles, extended_xmol_title=.TRUE.)
542 :
543 : ! Print Core Forces
544 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
545 2386 : "CORE_FORCES", my_pos, my_act, "cofrc", core_particles, extended_xmol_title=.TRUE.)
546 :
547 : ! Temperature per Kinds
548 2386 : IF (simpar%temperature_per_kind) THEN
549 : tempkind = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%TEMP_SHELL_KIND", &
550 376 : extension=".shtemp", file_position=my_pos, file_action=my_act)
551 376 : IF (tempkind > 0) THEN
552 21 : nkind = SIZE(md_ener%temp_shell_kind)
553 21 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nkind)))//"F20.9)"
554 : fmd = TRIM(fmd)
555 63 : WRITE (tempkind, fmd) itimes, time*femtoseconds, md_ener%temp_shell_kind(1:nkind)
556 21 : CALL m_flush(tempkind)
557 : END IF
558 : CALL cp_print_key_finished_output(tempkind, logger, motion_section, &
559 376 : "MD%PRINT%TEMP_SHELL_KIND")
560 : ELSE
561 2010 : print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_SHELL_KIND")
562 2010 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
563 : CALL cp_warn(__LOCATION__, &
564 : "The print_key MD%PRINT%TEMP_SHELL_KIND has been activated but the "// &
565 : "calculation of the temperature per kind has not been requested. "// &
566 80 : "Please switch on the keyword MD%TEMP_KIND.")
567 : END IF
568 : END IF
569 : END IF
570 : END IF
571 : END IF
572 42689 : init = .FALSE.
573 42689 : CALL set_md_env(md_env, init=init)
574 42689 : CALL timestop(handle)
575 42689 : END SUBROUTINE md_write_output
576 :
577 : ! **************************************************************************************************
578 : !> \brief This routine prints all basic information during MD steps
579 : !> \param simpar ...
580 : !> \param md_ener ...
581 : !> \param qmmm ...
582 : !> \param virial ...
583 : !> \param reftraj ...
584 : !> \param cell ...
585 : !> \param abc ...
586 : !> \param cell_angle ...
587 : !> \param itimes ...
588 : !> \param dt ...
589 : !> \param time ...
590 : !> \param used_time ...
591 : !> \param averages ...
592 : !> \param econs ...
593 : !> \param pv_scalar ...
594 : !> \param pv_xx ...
595 : !> \param hugoniot ...
596 : !> \param nat ...
597 : !> \param init ...
598 : !> \param logger ...
599 : !> \param motion_section ...
600 : !> \param my_pos ...
601 : !> \param my_act ...
602 : !> \param max_memory ...
603 : !> \par History
604 : !> - 10.2008 - Teodoro Laino [tlaino] - University of Zurich
605 : !> Refactoring: split into an independent routine.
606 : !> All output on screen must be included in this function!
607 : !> \author CJM
608 : ! **************************************************************************************************
609 42689 : SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
610 : abc, cell_angle, itimes, dt, time, used_time, averages, econs, pv_scalar, &
611 : pv_xx, hugoniot, nat, init, logger, motion_section, my_pos, my_act, &
612 : max_memory)
613 :
614 : TYPE(simpar_type), POINTER :: simpar
615 : TYPE(md_ener_type), POINTER :: md_ener
616 : LOGICAL, INTENT(IN) :: qmmm
617 : TYPE(virial_type), POINTER :: virial
618 : TYPE(reftraj_type), POINTER :: reftraj
619 : TYPE(cell_type), POINTER :: cell
620 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: abc, cell_angle
621 : INTEGER, POINTER :: itimes
622 : REAL(KIND=dp), INTENT(IN) :: dt
623 : REAL(KIND=dp), POINTER :: time, used_time
624 : TYPE(average_quantities_type), POINTER :: averages
625 : REAL(KIND=dp), INTENT(IN) :: econs, pv_scalar, pv_xx, hugoniot
626 : INTEGER, INTENT(IN) :: nat
627 : LOGICAL, INTENT(IN) :: init
628 : TYPE(cp_logger_type), POINTER :: logger
629 : TYPE(section_vals_type), POINTER :: motion_section
630 : CHARACTER(LEN=default_string_length), INTENT(IN) :: my_pos, my_act
631 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
632 :
633 : INTEGER :: iw
634 : TYPE(enumeration_type), POINTER :: enum
635 : TYPE(keyword_type), POINTER :: keyword
636 : TYPE(section_type), POINTER :: section
637 :
638 42689 : NULLIFY (enum, keyword, section)
639 : ! Print to the screen info about MD
640 : iw = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%PROGRAM_RUN_INFO", &
641 42689 : extension=".mdLog", file_position=my_pos, file_action=my_act)
642 :
643 : ! Performing protocol relevant to the first step of an MD run
644 42689 : IF (iw > 0) THEN
645 18547 : CALL create_md_section(section)
646 18547 : keyword => section_get_keyword(section, "ENSEMBLE")
647 18547 : CALL keyword_get(keyword, enum=enum)
648 18547 : IF (init) THEN
649 : ! Write initial values of quantities of interest
650 813 : WRITE (iw, '(/,T2,A)') 'MD_INI| MD initialization'
651 : WRITE (iw, '(T2,A,T61,E20.12)') &
652 813 : 'MD_INI| Potential energy [hartree]', md_ener%epot
653 813 : IF (simpar%ensemble /= reftraj_ensemble) THEN
654 795 : IF (.NOT. qmmm) THEN
655 : ! NO QM/MM info
656 : WRITE (iw, '(T2,A,T61,E20.12)') &
657 739 : 'MD_INI| Kinetic energy [hartree]', md_ener%ekin
658 : WRITE (iw, '(T2,A,T61,F20.6)') &
659 739 : 'MD_INI| Temperature [K]', md_ener%temp_part
660 : ELSE
661 : WRITE (iw, '(T2,A,T61,E20.12)') &
662 56 : 'MD_INI| Total kinetic energy [hartree]', md_ener%ekin, &
663 112 : 'MD_INI| QM kinetic energy [hartree]', md_ener%ekin_qm
664 : WRITE (iw, '(T2,A,T61,F20.6)') &
665 56 : 'MD_INI| Total temperature [K]', md_ener%temp_part, &
666 112 : 'MD_INI| QM temperature [K]', md_ener%temp_qm
667 : END IF
668 : END IF
669 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
670 : (simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
671 : (simpar%ensemble == npt_i_ensemble) .OR. &
672 : (simpar%ensemble == npt_ia_ensemble) .OR. &
673 : (simpar%ensemble == npt_f_ensemble) .OR. &
674 813 : (simpar%ensemble == npe_i_ensemble) .OR. &
675 : (simpar%ensemble == npe_f_ensemble)) THEN
676 : WRITE (iw, '(T2,A,T61,F20.6)') &
677 86 : 'MD_INI| Barostat temperature [K]', md_ener%temp_baro
678 : END IF
679 813 : IF (virial%pv_availability) THEN
680 : WRITE (iw, '(T2,A,T61,ES20.12)') &
681 142 : 'MD_INI| Pressure [bar]', pv_scalar
682 : END IF
683 813 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
684 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
685 : WRITE (iw, '(T2,A,T61,ES20.12)') &
686 3 : 'MD_INI| Hugoniot constraint [K]', hugoniot
687 : END IF
688 813 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
689 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
690 : WRITE (iw, '(T2,A,T61,E20.12)') &
691 3 : 'MD_INI| Total energy [hartree]', simpar%e0
692 : END IF
693 : WRITE (iw, '(T2,A,T61,ES20.12)') &
694 813 : 'MD_INI| Cell volume [bohr^3]', cell%deth
695 : WRITE (iw, '(T2,A,T61,ES20.12)') &
696 813 : 'MD_INI| Cell volume [ang^3]', cell%deth*angstrom**3
697 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
698 813 : 'MD_INI| Cell lengths [bohr]', abc(1:3)
699 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
700 3252 : 'MD_INI| Cell lengths [ang]', abc(1:3)*angstrom
701 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
702 813 : 'MD_INI| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
703 : ELSE
704 : ! Write seuquential values of quantities of interest
705 17734 : WRITE (iw, '(/,T2,A)') 'MD| '//REPEAT('*', 75)
706 : !MK WRITE (iw, '(T2,A,T61,A20)') &
707 : !MK 'MD| Ensemble type', ADJUSTR(TRIM(enum_i2c(enum, simpar%ensemble)))
708 : WRITE (iw, '(T2,A,T71,I10)') &
709 17734 : 'MD| Step number', itimes
710 17734 : IF (simpar%variable_dt) THEN
711 : WRITE (iw, '(T2,A,T61,F20.6)') &
712 240 : 'MD| Time step [fs]', dt*femtoseconds
713 : END IF
714 : WRITE (iw, '(T2,A,T61,F20.6)') &
715 17734 : 'MD| Time [fs]', time*femtoseconds
716 : WRITE (iw, '(T2,A,T61,E20.12)') &
717 17734 : 'MD| Conserved quantity [hartree]', md_ener%constant
718 17734 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
719 17734 : WRITE (iw, '(T2,A,T47,A,T73,A)') 'MD|', 'Instantaneous', 'Averages'
720 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
721 17734 : 'MD| CPU time per MD step [s]', used_time, averages%avecpu
722 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
723 17734 : 'MD| Energy drift per atom [K] ', econs, averages%econs
724 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
725 17734 : 'MD| Potential energy [hartree]', md_ener%epot, averages%avepot
726 17734 : IF (simpar%ensemble /= reftraj_ensemble) THEN
727 17611 : IF (.NOT. qmmm) THEN
728 : ! No QM/MM info
729 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
730 16972 : 'MD| Kinetic energy [hartree]', md_ener%ekin, averages%avekin
731 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
732 16972 : 'MD| Temperature [K]', md_ener%temp_part, averages%avetemp
733 : ELSE
734 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
735 639 : 'MD| Total kinetic energy [hartree]', md_ener%ekin, averages%avekin
736 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
737 639 : 'MD| QM kinetic energy [hartree]', md_ener%ekin_qm, averages%avekin_qm
738 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
739 639 : 'MD| Total temperature [K]', md_ener%temp_part, averages%avetemp
740 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
741 639 : 'MD| QM temperature [K]', md_ener%temp_qm, averages%avetemp_qm
742 : END IF
743 : END IF
744 17734 : IF (virial%pv_availability) THEN
745 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
746 1852 : 'MD| Pressure [bar]', pv_scalar, averages%avepress
747 : END IF
748 17734 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
749 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
750 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
751 30 : 'MD| P(xx) [bar]', pv_xx, averages%avepxx
752 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
753 30 : 'MD| Hugoniot [K]', hugoniot/3.0_dp/nat*kelvin, averages%avehugoniot/3.0_dp/nat*kelvin
754 : END IF
755 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
756 : (simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
757 : (simpar%ensemble == npt_i_ensemble) .OR. &
758 : (simpar%ensemble == npt_ia_ensemble) .OR. &
759 : (simpar%ensemble == npt_f_ensemble) .OR. &
760 17734 : (simpar%ensemble == npe_i_ensemble) .OR. &
761 : (simpar%ensemble == npe_f_ensemble)) THEN
762 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
763 990 : 'MD| Barostat temperature [K]', md_ener%temp_baro, averages%avetemp_baro
764 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
765 990 : 'MD| Cell volume [bohr^3]', cell%deth, averages%avevol
766 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
767 990 : 'MD| Cell volume [ang^3]', cell%deth*angstrom**3, averages%avevol*angstrom**3
768 990 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
769 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
770 990 : 'MD| Cell lengths [bohr]', abc(1:3)
771 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
772 3960 : 'MD| Cell lengths [ang]', abc(1:3)*angstrom
773 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
774 990 : 'MD| Average cell lengths [bohr]', averages%aveca, averages%avecb, averages%avecc
775 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
776 990 : 'MD| Average cell lengths [ang]', averages%aveca*angstrom, averages%avecb*angstrom, &
777 1980 : averages%avecc*angstrom
778 : END IF
779 17734 : IF ((simpar%ensemble == npt_f_ensemble) .OR. &
780 : (simpar%ensemble == npe_f_ensemble)) THEN
781 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
782 458 : 'MD| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
783 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
784 458 : 'MD| Average cell angles [deg]', averages%aveal, averages%avebe, averages%avega
785 : END IF
786 17734 : IF (simpar%ensemble == reftraj_ensemble) THEN
787 123 : IF (reftraj%info%msd) THEN
788 6 : IF (reftraj%msd%msd_kind) THEN
789 : WRITE (iw, '(/,T2,A,T51,3F10.5)') &
790 24 : 'MD| COM displacement (dx,dy,dz) [bohr]', reftraj%msd%drcom(1:3)
791 : END IF
792 : END IF
793 : END IF
794 17734 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('*', 75)
795 17734 : IF (max_memory /= 0) THEN
796 : WRITE (iw, '(T2,A,T73,I8)') &
797 17734 : 'MD| Estimated peak process memory after this step [MiB]', &
798 35468 : (max_memory + (1024*1024) - 1)/(1024*1024)
799 : END IF
800 : END IF
801 : END IF
802 42689 : CALL section_release(section)
803 : CALL cp_print_key_finished_output(iw, logger, motion_section, &
804 42689 : "MD%PRINT%PROGRAM_RUN_INFO")
805 42689 : END SUBROUTINE md_write_info_low
806 :
807 : ! **************************************************************************************************
808 : !> \brief Samples memory usage
809 : !> \param para_env ...
810 : !> \return ...
811 : !> \note based on what is done in start/cp2k_runs.F
812 : ! **************************************************************************************************
813 52863 : FUNCTION sample_memory(para_env) RESULT(max_memory)
814 : TYPE(mp_para_env_type), POINTER :: para_env
815 : INTEGER(KIND=int_8) :: max_memory
816 :
817 52863 : CALL m_memory()
818 52863 : max_memory = m_memory_max
819 52863 : CALL para_env%max(max_memory)
820 :
821 52863 : END FUNCTION sample_memory
822 :
823 : END MODULE md_energies
|
