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 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,&
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_type
82 : USE simpar_types, ONLY: simpar_type
83 : USE thermal_region_types, ONLY: thermal_regions_type
84 : USE thermal_region_utils, ONLY: print_thermal_regions_temperature
85 : USE thermostat_types, ONLY: thermostats_type
86 : USE thermostat_utils, ONLY: print_thermostats_status
87 : USE virial_types, ONLY: virial_type
88 : #include "../base/base_uses.f90"
89 :
90 : IMPLICIT NONE
91 :
92 : PRIVATE
93 :
94 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_energies'
95 :
96 : PUBLIC :: initialize_md_ener, &
97 : md_energy, &
98 : md_ener_reftraj, &
99 : md_write_output, &
100 : sample_memory
101 :
102 : CONTAINS
103 :
104 : ! **************************************************************************************************
105 : !> \brief ...
106 : !> \param md_ener ...
107 : !> \param force_env ...
108 : !> \param simpar ...
109 : !> \par History
110 : !> - 10-2007 created
111 : !> \author MI
112 : ! **************************************************************************************************
113 3458 : SUBROUTINE initialize_md_ener(md_ener, force_env, simpar)
114 :
115 : TYPE(md_ener_type), POINTER :: md_ener
116 : TYPE(force_env_type), POINTER :: force_env
117 : TYPE(simpar_type), POINTER :: simpar
118 :
119 : INTEGER :: nkind
120 : LOGICAL :: shell_adiabatic
121 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
122 1729 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
123 : TYPE(cp_subsys_type), POINTER :: subsys
124 : TYPE(particle_list_type), POINTER :: particles, shell_particles
125 :
126 1729 : NULLIFY (subsys)
127 1729 : NULLIFY (atomic_kinds, atomic_kind_set, particles, shell_particles)
128 :
129 0 : CPASSERT(ASSOCIATED(md_ener))
130 1729 : CPASSERT(ASSOCIATED(force_env))
131 :
132 1729 : CALL force_env_get(force_env, subsys=subsys)
133 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
134 1729 : shell_particles=shell_particles)
135 1729 : atomic_kind_set => atomic_kinds%els
136 1729 : nkind = SIZE(atomic_kind_set)
137 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
138 1729 : shell_adiabatic=shell_adiabatic)
139 :
140 1729 : md_ener%nfree = simpar%nfree
141 1729 : md_ener%nfree_shell = -HUGE(0)
142 :
143 1729 : IF (shell_adiabatic) THEN
144 132 : md_ener%nfree_shell = 3*(shell_particles%n_els)
145 : END IF
146 :
147 1729 : IF (simpar%temperature_per_kind) THEN
148 84 : ALLOCATE (md_ener%temp_kind(nkind))
149 84 : ALLOCATE (md_ener%ekin_kind(nkind))
150 84 : ALLOCATE (md_ener%nfree_kind(nkind))
151 86 : md_ener%nfree_kind = 0
152 :
153 28 : IF (shell_adiabatic) THEN
154 54 : ALLOCATE (md_ener%temp_shell_kind(nkind))
155 54 : ALLOCATE (md_ener%ekin_shell_kind(nkind))
156 54 : ALLOCATE (md_ener%nfree_shell_kind(nkind))
157 54 : md_ener%nfree_shell_kind = 0
158 : END IF
159 :
160 : END IF
161 : CALL zero_md_ener(md_ener, tkind=simpar%temperature_per_kind, &
162 1729 : tshell=shell_adiabatic)
163 1729 : md_ener%epot = 0.0_dp
164 :
165 1729 : END SUBROUTINE initialize_md_ener
166 :
167 : ! **************************************************************************************************
168 : !> \brief ...
169 : !> \param md_env ...
170 : !> \param md_ener ...
171 : !> \par History
172 : !> - 10-2007 created
173 : !> \author MI
174 : ! **************************************************************************************************
175 82170 : SUBROUTINE md_energy(md_env, md_ener)
176 :
177 : TYPE(md_environment_type), POINTER :: md_env
178 : TYPE(md_ener_type), POINTER :: md_ener
179 :
180 : INTEGER :: natom
181 : LOGICAL :: shell_adiabatic, tkind, tshell
182 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
183 41085 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
184 : TYPE(cp_subsys_type), POINTER :: subsys
185 : TYPE(force_env_type), POINTER :: force_env
186 : TYPE(particle_list_type), POINTER :: particles
187 : TYPE(simpar_type), POINTER :: simpar
188 :
189 41085 : NULLIFY (atomic_kinds, atomic_kind_set, force_env, &
190 41085 : particles, subsys, simpar)
191 : CALL get_md_env(md_env=md_env, force_env=force_env, &
192 41085 : simpar=simpar)
193 :
194 : CALL force_env_get(force_env, &
195 41085 : potential_energy=md_ener%epot, subsys=subsys)
196 :
197 41085 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
198 41085 : atomic_kind_set => atomic_kinds%els
199 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
200 41085 : shell_adiabatic=shell_adiabatic)
201 :
202 41085 : tkind = simpar%temperature_per_kind
203 41085 : tshell = shell_adiabatic
204 :
205 41085 : CALL cp_subsys_get(subsys, particles=particles)
206 41085 : natom = particles%n_els
207 :
208 : CALL compute_conserved_quantity(md_env, md_ener, tkind=tkind, &
209 41085 : tshell=tshell, natom=natom)
210 :
211 41085 : END SUBROUTINE md_energy
212 :
213 : ! **************************************************************************************************
214 : !> \brief ...
215 : !> \param md_env ...
216 : !> \param md_ener ...
217 : !> \par History
218 : !> - 10.2007 created
219 : !> \author MI
220 : ! **************************************************************************************************
221 262 : SUBROUTINE md_ener_reftraj(md_env, md_ener)
222 : TYPE(md_environment_type), POINTER :: md_env
223 : TYPE(md_ener_type), POINTER :: md_ener
224 :
225 : TYPE(force_env_type), POINTER :: force_env
226 : TYPE(reftraj_type), POINTER :: reftraj
227 :
228 262 : CALL zero_md_ener(md_ener, tkind=.FALSE., tshell=.FALSE.)
229 262 : CALL get_md_env(md_env=md_env, force_env=force_env, reftraj=reftraj)
230 :
231 262 : IF (reftraj%info%eval_ef) THEN
232 124 : CALL force_env_get(force_env, potential_energy=md_ener%epot)
233 : ELSE
234 138 : md_ener%epot = reftraj%epot
235 138 : md_ener%delta_epot = (reftraj%epot - reftraj%epot0)/REAL(reftraj%natom, kind=dp)*kelvin
236 : END IF
237 :
238 262 : END SUBROUTINE md_ener_reftraj
239 :
240 : ! **************************************************************************************************
241 : !> \brief This routine computes the conserved quantity, temperature
242 : !> and things like that and prints them out
243 : !> \param md_env ...
244 : !> \par History
245 : !> - New MD data are appended to the old data (15.09.2003,MK)
246 : !> - 02.2008 - Teodoro Laino [tlaino] - University of Zurich
247 : !> Cleaning code and collecting the many commons routines..
248 : !> \author CJM
249 : ! **************************************************************************************************
250 165488 : SUBROUTINE md_write_output(md_env)
251 :
252 : TYPE(md_environment_type), POINTER :: md_env
253 :
254 : CHARACTER(len=*), PARAMETER :: routineN = 'md_write_output'
255 :
256 : CHARACTER(LEN=default_string_length) :: fmd, my_act, my_pos
257 : INTEGER :: ene, ener_mix, handle, i, nat, nkind, &
258 : shene, tempkind, trsl
259 : INTEGER(KIND=int_8) :: max_memory
260 : INTEGER, POINTER :: itimes
261 : LOGICAL :: init, is_mixed, new_file, print_memory, &
262 : qmmm, shell_adiabatic, shell_present
263 : REAL(dp) :: abc(3), cell_angle(3), dt, econs, &
264 : pv_scalar, pv_xx, pv_xx_nc
265 : REAL(KIND=dp) :: harm_shell, hugoniot
266 : REAL(KIND=dp), POINTER :: time, used_time
267 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
268 41372 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
269 : TYPE(average_quantities_type), POINTER :: averages
270 : TYPE(barostat_type), POINTER :: barostat
271 : TYPE(cell_type), POINTER :: cell
272 : TYPE(cp_logger_type), POINTER :: logger
273 : TYPE(cp_subsys_type), POINTER :: subsys
274 : TYPE(force_env_type), POINTER :: force_env
275 : TYPE(md_ener_type), POINTER :: md_ener
276 : TYPE(mp_para_env_type), POINTER :: para_env
277 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
278 : shell_particles
279 : TYPE(qmmm_env_type), POINTER :: qmmm_env
280 : TYPE(reftraj_type), POINTER :: reftraj
281 : TYPE(section_vals_type), POINTER :: motion_section, print_key, root_section
282 : TYPE(simpar_type), POINTER :: simpar
283 : TYPE(thermal_regions_type), POINTER :: thermal_regions
284 : TYPE(thermostats_type), POINTER :: thermostats
285 : TYPE(virial_type), POINTER :: virial
286 :
287 41372 : NULLIFY (logger)
288 82744 : logger => cp_get_default_logger()
289 41372 : CALL timeset(routineN, handle)
290 :
291 : ! Zeroing
292 41372 : hugoniot = 0.0_dp
293 41372 : econs = 0.0_dp
294 41372 : shell_adiabatic = .FALSE.
295 41372 : shell_present = .FALSE.
296 41372 : NULLIFY (motion_section, atomic_kinds, atomic_kind_set, cell, subsys, &
297 41372 : force_env, md_ener, qmmm_env, reftraj, core_particles, particles, &
298 41372 : shell_particles, print_key, root_section, simpar, virial, &
299 41372 : thermostats, thermal_regions)
300 :
301 : CALL get_md_env(md_env=md_env, itimes=itimes, t=time, used_time=used_time, &
302 : simpar=simpar, force_env=force_env, init=init, md_ener=md_ener, &
303 : reftraj=reftraj, thermostats=thermostats, barostat=barostat, &
304 41372 : para_env=para_env, averages=averages, thermal_regions=thermal_regions)
305 :
306 41372 : root_section => force_env%root_section
307 41372 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
308 :
309 41372 : CALL force_env_get(force_env, cell=cell, subsys=subsys, qmmm_env=qmmm_env)
310 :
311 41372 : qmmm = calc_nfree_qm(md_env, md_ener) > 0
312 41372 : is_mixed = (force_env%in_use == use_mixed_force)
313 :
314 41372 : CALL cp_subsys_get(subsys, particles=particles, virial=virial)
315 41372 : nat = particles%n_els
316 41372 : dt = simpar%dt*simpar%dt_fact
317 :
318 : ! Computing the scalar pressure
319 41372 : IF (virial%pv_availability) THEN
320 4320 : pv_scalar = 0._dp
321 17280 : DO i = 1, 3
322 17280 : pv_scalar = pv_scalar + virial%pv_total(i, i)
323 : END DO
324 4320 : pv_scalar = pv_scalar/3._dp/cell%deth
325 4320 : pv_scalar = cp_unit_from_cp2k(pv_scalar, "bar")
326 4320 : pv_xx_nc = virial%pv_total(1, 1)/cell%deth
327 4320 : pv_xx = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
328 : END IF
329 :
330 41372 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
331 41372 : atomic_kind_set => atomic_kinds%els
332 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
333 : shell_present=shell_present, &
334 41372 : shell_adiabatic=shell_adiabatic)
335 :
336 41372 : CALL get_cell(cell, abc=abc, alpha=cell_angle(3), beta=cell_angle(2), gamma=cell_angle(1))
337 :
338 : ! Determine POS and ACT for I/O
339 41372 : my_pos = "APPEND"
340 41372 : my_act = "WRITE"
341 41372 : IF (init .AND. (itimes == 0)) THEN
342 1488 : my_pos = "REWIND"
343 1488 : my_act = "WRITE"
344 : END IF
345 :
346 : ! In case of REFTRAJ ensemble the POS is determined differently..
347 : ! according the REFTRAJ counter
348 41372 : IF (simpar%ensemble == reftraj_ensemble) THEN
349 262 : IF ((reftraj%isnap == reftraj%info%first_snapshot)) THEN
350 28 : my_pos = "REWIND"
351 : END IF
352 : END IF
353 :
354 : ! Performing protocol relevant to the first step of an MD run
355 41372 : IF (init) THEN
356 : ! Computing the Hugoniot for NPH calculations
357 1727 : IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
358 : simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
359 6 : IF (simpar%e0 == 0._dp) simpar%e0 = md_ener%epot + md_ener%ekin
360 : hugoniot = md_ener%epot + md_ener%ekin - simpar%e0 - 0.5_dp*(pv_xx_nc + simpar%p0)* &
361 6 : (simpar%v0 - cell%deth)
362 : END IF
363 :
364 1727 : IF (simpar%ensemble == reftraj_ensemble) reftraj%init = init
365 : ELSE
366 : ! Performing protocol for anything beyond the first step of MD
367 39645 : IF (simpar%ensemble == nph_uniaxial_ensemble .OR. simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
368 : hugoniot = md_ener%epot + md_ener%ekin - simpar%e0 - 0.5_dp*(pv_xx_nc + simpar%p0)* &
369 60 : (simpar%v0 - cell%deth)
370 : END IF
371 :
372 39645 : IF (simpar%ensemble == reftraj_ensemble) THEN
373 230 : time = reftraj%time
374 230 : econs = md_ener%delta_epot
375 230 : itimes = reftraj%itimes
376 : ELSE
377 39415 : econs = md_ener%delta_cons
378 : END IF
379 :
380 : ! Compute average quantities
381 : CALL compute_averages(averages, force_env, md_ener, cell, virial, pv_scalar, &
382 : pv_xx, used_time, hugoniot, abc, cell_angle, nat, itimes, time, my_pos, &
383 39645 : my_act)
384 : END IF
385 :
386 : ! Sample memory, if requested
387 41372 : CALL section_vals_val_get(motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
388 41372 : max_memory = 0
389 41372 : IF (print_memory) THEN
390 41372 : max_memory = sample_memory(para_env)
391 : END IF
392 :
393 : ! Print md information
394 : CALL md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, abc, &
395 : cell_angle, itimes, dt, time, used_time, averages, econs, pv_scalar, pv_xx, &
396 41372 : hugoniot, nat, init, logger, motion_section, my_pos, my_act, max_memory)
397 :
398 : ! Real Output driven by the PRINT sections
399 41372 : IF ((.NOT. init) .OR. (itimes == 0) .OR. simpar%ensemble == reftraj_ensemble) THEN
400 : ! Print Energy
401 : ene = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%ENERGY", &
402 41165 : extension=".ener", file_position=my_pos, file_action=my_act, is_new_file=new_file)
403 41165 : IF (ene > 0) THEN
404 17955 : IF (new_file) THEN
405 : ! Please change also the corresponding format explanation below
406 : ! keep the constant of motion the true constant of motion !
407 726 : WRITE (ene, '("#",5X,A,10X,A,8X,A,10X,A,12X,A,2(8X,A))') "Step Nr.", "Time[fs]", "Kin.[a.u.]", "Temp[K]", &
408 1452 : "Pot.[a.u.]", "Cons Qty[a.u.]", "UsedTime[s]"
409 : END IF
410 : WRITE (ene, "(I10,F20.6,F20.9,F20.9,F20.9,F20.9,F20.9)") &
411 17955 : itimes, time*femtoseconds, md_ener%ekin, md_ener%temp_part, md_ener%epot, md_ener%constant, used_time
412 17955 : CALL m_flush(ene)
413 : END IF
414 41165 : CALL cp_print_key_finished_output(ene, logger, motion_section, "MD%PRINT%ENERGY")
415 :
416 : ! Possibly Print MIXED Energy
417 41165 : IF (is_mixed) THEN
418 : ener_mix = cp_print_key_unit_nr(logger, motion_section, "PRINT%MIXED_ENERGIES", &
419 : extension=".ener", file_position=my_pos, file_action=my_act, &
420 342 : middle_name="mix")
421 342 : IF (ener_mix > 0) THEN
422 : WRITE (ener_mix, "(I8,F12.3,F20.9,"//cp_to_string(SIZE(force_env%mixed_env%energies))//"F20.9,F20.9)") &
423 545 : itimes, time*femtoseconds, md_ener%epot, force_env%mixed_env%energies, md_ener%constant
424 171 : CALL m_flush(ener_mix)
425 : END IF
426 342 : CALL cp_print_key_finished_output(ener_mix, logger, motion_section, "PRINT%MIXED_ENERGIES")
427 : END IF
428 :
429 : ! Print QMMM translation vector if requested
430 41165 : IF (qmmm) THEN
431 : trsl = cp_print_key_unit_nr(logger, motion_section, "PRINT%TRANSLATION_VECTOR", &
432 1424 : extension=".translation", middle_name="qmmm")
433 1424 : IF (trsl > 0) THEN
434 0 : WRITE (trsl, '(I10,3F15.10)') itimes, qmmm_env%qm%transl_v
435 : END IF
436 : CALL cp_print_key_finished_output(trsl, logger, motion_section, &
437 1424 : "PRINT%TRANSLATION_VECTOR")
438 : END IF
439 :
440 : ! Write Structure data
441 41165 : CALL write_structure_data(particles%els, cell, motion_section)
442 :
443 : ! Print Coordinates
444 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
445 41165 : pos=my_pos, act=my_act, extended_xmol_title=.TRUE.)
446 :
447 : ! Print Velocities
448 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
449 41165 : "VELOCITIES", my_pos, my_act, middle_name="vel", extended_xmol_title=.TRUE.)
450 :
451 : ! Print Force
452 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
453 41165 : "FORCES", my_pos, my_act, middle_name="frc", extended_xmol_title=.TRUE.)
454 :
455 : ! Print Force-Mixing labels
456 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
457 41165 : "FORCE_MIXING_LABELS", my_pos, my_act, middle_name="fmlabels", extended_xmol_title=.TRUE.)
458 :
459 : ! Print Simulation Cell
460 41165 : CALL write_simulation_cell(cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
461 :
462 : ! Print Thermostats status
463 41165 : CALL print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
464 :
465 : ! Print Barostat status
466 41165 : CALL print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time)
467 :
468 : ! Print Stress Tensor
469 41165 : CALL write_stress_tensor(virial, cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
470 :
471 : ! Print Polarisability Tensor
472 41165 : IF (ASSOCIATED(force_env%qs_env)) THEN
473 3594 : CALL write_polarisability_tensor(force_env, motion_section, itimes, time*femtoseconds, my_pos, my_act)
474 : END IF
475 :
476 : ! Temperature per Kinds
477 41165 : IF (simpar%temperature_per_kind) THEN
478 : tempkind = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%TEMP_KIND", &
479 492 : extension=".temp", file_position=my_pos, file_action=my_act)
480 492 : IF (tempkind > 0) THEN
481 242 : nkind = SIZE(md_ener%temp_kind)
482 242 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nkind)))//"F20.9)"
483 : fmd = TRIM(fmd)
484 242 : WRITE (tempkind, fmd) itimes, time*femtoseconds, md_ener%temp_kind(1:nkind)
485 242 : CALL m_flush(tempkind)
486 : END IF
487 492 : CALL cp_print_key_finished_output(tempkind, logger, motion_section, "MD%PRINT%TEMP_KIND")
488 : ELSE
489 40673 : print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_KIND")
490 40673 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) &
491 : CALL cp_warn(__LOCATION__, &
492 : "The print_key MD%PRINT%TEMP_KIND has been activated but the "// &
493 : "calculation of the temperature per kind has not been requested. "// &
494 316 : "Please switch on the keyword MD%TEMP_KIND.")
495 : END IF
496 : !Thermal Region
497 41165 : CALL print_thermal_regions_temperature(thermal_regions, itimes, time*femtoseconds, my_pos, my_act)
498 :
499 : ! Core/Shell Model
500 41165 : IF (shell_present) THEN
501 2386 : CALL force_env_get(force_env, harmonic_shell=harm_shell)
502 2386 : CALL cp_subsys_get(subsys, shell_particles=shell_particles, core_particles=core_particles)
503 :
504 : ! Print Shell Energy
505 : shene = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%SHELL_ENERGY", &
506 : extension=".shener", file_position=my_pos, file_action=my_act, &
507 2386 : file_form="FORMATTED", is_new_file=new_file)
508 2386 : IF (shene > 0) THEN
509 1106 : IF (new_file) THEN
510 45 : WRITE (shene, '("#",3X,A,3X,A,3X,3(5X,A,5X))') "Step Nr.", "Time[fs]", "Kin.[a.u.]", &
511 90 : "Temp.[K]", "Pot.[a.u.]"
512 : END IF
513 :
514 : WRITE (shene, "(I8,F12.3,F20.9,F20.9,F20.9,F20.9 )") &
515 1106 : itimes, time*femtoseconds, md_ener%ekin_shell, md_ener%temp_shell, harm_shell
516 1106 : CALL m_flush(shene)
517 : END IF
518 2386 : CALL cp_print_key_finished_output(shene, logger, motion_section, "MD%PRINT%SHELL_ENERGY")
519 :
520 : ! Print Shell Coordinates
521 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
522 2386 : "SHELL_TRAJECTORY", my_pos, my_act, "shpos", shell_particles, extended_xmol_title=.TRUE.)
523 :
524 2386 : IF (shell_adiabatic) THEN
525 : ! Print Shell Velocities
526 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
527 2386 : "SHELL_VELOCITIES", my_pos, my_act, "shvel", shell_particles, extended_xmol_title=.TRUE.)
528 :
529 : ! Print Shell Forces
530 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
531 2386 : "SHELL_FORCES", my_pos, my_act, "shfrc", shell_particles, extended_xmol_title=.TRUE.)
532 :
533 : ! Print Core Coordinates
534 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
535 2386 : "CORE_TRAJECTORY", my_pos, my_act, "copos", core_particles, extended_xmol_title=.TRUE.)
536 :
537 : ! Print Core Velocities
538 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
539 2386 : "CORE_VELOCITIES", my_pos, my_act, "covel", core_particles, extended_xmol_title=.TRUE.)
540 :
541 : ! Print Core Forces
542 : CALL write_trajectory(force_env, root_section, itimes, time*femtoseconds, dt*femtoseconds, md_ener%epot, &
543 2386 : "CORE_FORCES", my_pos, my_act, "cofrc", core_particles, extended_xmol_title=.TRUE.)
544 :
545 : ! Temperature per Kinds
546 2386 : IF (simpar%temperature_per_kind) THEN
547 : tempkind = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%TEMP_SHELL_KIND", &
548 376 : extension=".shtemp", file_position=my_pos, file_action=my_act)
549 376 : IF (tempkind > 0) THEN
550 21 : nkind = SIZE(md_ener%temp_shell_kind)
551 21 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nkind)))//"F20.9)"
552 : fmd = TRIM(fmd)
553 21 : WRITE (tempkind, fmd) itimes, time*femtoseconds, md_ener%temp_shell_kind(1:nkind)
554 21 : CALL m_flush(tempkind)
555 : END IF
556 : CALL cp_print_key_finished_output(tempkind, logger, motion_section, &
557 376 : "MD%PRINT%TEMP_SHELL_KIND")
558 : ELSE
559 2010 : print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_SHELL_KIND")
560 2010 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) &
561 : CALL cp_warn(__LOCATION__, &
562 : "The print_key MD%PRINT%TEMP_SHELL_KIND has been activated but the "// &
563 : "calculation of the temperature per kind has not been requested. "// &
564 80 : "Please switch on the keyword MD%TEMP_KIND.")
565 : END IF
566 : END IF
567 : END IF
568 : END IF
569 41372 : init = .FALSE.
570 41372 : CALL set_md_env(md_env, init=init)
571 41372 : CALL timestop(handle)
572 41372 : END SUBROUTINE md_write_output
573 :
574 : ! **************************************************************************************************
575 : !> \brief This routine prints all basic information during MD steps
576 : !> \param simpar ...
577 : !> \param md_ener ...
578 : !> \param qmmm ...
579 : !> \param virial ...
580 : !> \param reftraj ...
581 : !> \param cell ...
582 : !> \param abc ...
583 : !> \param cell_angle ...
584 : !> \param itimes ...
585 : !> \param dt ...
586 : !> \param time ...
587 : !> \param used_time ...
588 : !> \param averages ...
589 : !> \param econs ...
590 : !> \param pv_scalar ...
591 : !> \param pv_xx ...
592 : !> \param hugoniot ...
593 : !> \param nat ...
594 : !> \param init ...
595 : !> \param logger ...
596 : !> \param motion_section ...
597 : !> \param my_pos ...
598 : !> \param my_act ...
599 : !> \param max_memory ...
600 : !> \par History
601 : !> - 10.2008 - Teodoro Laino [tlaino] - University of Zurich
602 : !> Refactoring: split into an independent routine.
603 : !> All output on screen must be included in this function!
604 : !> \author CJM
605 : ! **************************************************************************************************
606 41372 : SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
607 : abc, cell_angle, itimes, dt, time, used_time, averages, econs, pv_scalar, &
608 : pv_xx, hugoniot, nat, init, logger, motion_section, my_pos, my_act, &
609 : max_memory)
610 :
611 : TYPE(simpar_type), POINTER :: simpar
612 : TYPE(md_ener_type), POINTER :: md_ener
613 : LOGICAL, INTENT(IN) :: qmmm
614 : TYPE(virial_type), POINTER :: virial
615 : TYPE(reftraj_type), POINTER :: reftraj
616 : TYPE(cell_type), POINTER :: cell
617 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: abc, cell_angle
618 : INTEGER, POINTER :: itimes
619 : REAL(KIND=dp), INTENT(IN) :: dt
620 : REAL(KIND=dp), POINTER :: time, used_time
621 : TYPE(average_quantities_type), POINTER :: averages
622 : REAL(KIND=dp), INTENT(IN) :: econs, pv_scalar, pv_xx, hugoniot
623 : INTEGER, INTENT(IN) :: nat
624 : LOGICAL, INTENT(IN) :: init
625 : TYPE(cp_logger_type), POINTER :: logger
626 : TYPE(section_vals_type), POINTER :: motion_section
627 : CHARACTER(LEN=default_string_length), INTENT(IN) :: my_pos, my_act
628 : INTEGER(KIND=int_8), INTENT(IN) :: max_memory
629 :
630 : INTEGER :: iw
631 : TYPE(enumeration_type), POINTER :: enum
632 : TYPE(keyword_type), POINTER :: keyword
633 : TYPE(section_type), POINTER :: section
634 :
635 41372 : NULLIFY (enum, keyword, section)
636 : ! Print to the screen info about MD
637 : iw = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%PROGRAM_RUN_INFO", &
638 41372 : extension=".mdLog", file_position=my_pos, file_action=my_act)
639 :
640 : ! Performing protocol relevant to the first step of an MD run
641 41372 : IF (iw > 0) THEN
642 18283 : CALL create_md_section(section)
643 18283 : keyword => section_get_keyword(section, "ENSEMBLE")
644 18283 : CALL keyword_get(keyword, enum=enum)
645 18283 : IF (init) THEN
646 : ! Write initial values of quantities of interest
647 793 : WRITE (iw, '(/,T2,A)') 'MD_INI| MD initialization'
648 : WRITE (iw, '(T2,A,T61,E20.12)') &
649 793 : 'MD_INI| Potential energy [hartree]', md_ener%epot
650 793 : IF (simpar%ensemble /= reftraj_ensemble) THEN
651 777 : IF (.NOT. qmmm) THEN
652 : ! NO QM/MM info
653 : WRITE (iw, '(T2,A,T61,E20.12)') &
654 722 : 'MD_INI| Kinetic energy [hartree]', md_ener%ekin
655 : WRITE (iw, '(T2,A,T61,F20.6)') &
656 722 : 'MD_INI| Temperature [K]', md_ener%temp_part
657 : ELSE
658 : WRITE (iw, '(T2,A,T61,E20.12)') &
659 55 : 'MD_INI| Total kinetic energy [hartree]', md_ener%ekin, &
660 110 : 'MD_INI| QM kinetic energy [hartree]', md_ener%ekin_qm
661 : WRITE (iw, '(T2,A,T61,F20.6)') &
662 55 : 'MD_INI| Total temperature [K]', md_ener%temp_part, &
663 110 : 'MD_INI| QM temperature [K]', md_ener%temp_qm
664 : END IF
665 : END IF
666 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
667 : (simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
668 : (simpar%ensemble == npt_i_ensemble) .OR. &
669 : (simpar%ensemble == npt_ia_ensemble) .OR. &
670 : (simpar%ensemble == npt_f_ensemble) .OR. &
671 793 : (simpar%ensemble == npe_i_ensemble) .OR. &
672 : (simpar%ensemble == npe_f_ensemble)) THEN
673 : WRITE (iw, '(T2,A,T61,F20.6)') &
674 79 : 'MD_INI| Barostat temperature [K]', md_ener%temp_baro
675 : END IF
676 793 : IF (virial%pv_availability) THEN
677 : WRITE (iw, '(T2,A,T61,ES20.12)') &
678 143 : 'MD_INI| Pressure [bar]', pv_scalar
679 : END IF
680 793 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
681 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
682 : WRITE (iw, '(T2,A,T61,ES20.12)') &
683 3 : 'MD_INI| Hugoniot constraint [K]', hugoniot
684 : END IF
685 793 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
686 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
687 : WRITE (iw, '(T2,A,T61,E20.12)') &
688 3 : 'MD_INI| Total energy [hartree]', simpar%e0
689 : END IF
690 : WRITE (iw, '(T2,A,T61,ES20.12)') &
691 793 : 'MD_INI| Cell volume [bohr^3]', cell%deth
692 : WRITE (iw, '(T2,A,T61,ES20.12)') &
693 793 : 'MD_INI| Cell volume [ang^3]', cell%deth*angstrom**3
694 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
695 793 : 'MD_INI| Cell lengths [bohr]', abc(1:3)
696 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
697 3172 : 'MD_INI| Cell lengths [ang]', abc(1:3)*angstrom
698 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
699 793 : 'MD_INI| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
700 : ELSE
701 : ! Write seuquential values of quantities of interest
702 17490 : WRITE (iw, '(/,T2,A)') 'MD| '//REPEAT('*', 75)
703 : !MK WRITE (iw, '(T2,A,T61,A20)') &
704 : !MK 'MD| Ensemble type', ADJUSTR(TRIM(enum_i2c(enum, simpar%ensemble)))
705 : WRITE (iw, '(T2,A,T71,I10)') &
706 17490 : 'MD| Step number', itimes
707 17490 : IF (simpar%variable_dt) THEN
708 : WRITE (iw, '(T2,A,T61,F20.6)') &
709 240 : 'MD| Time step [fs]', dt*femtoseconds
710 : END IF
711 : WRITE (iw, '(T2,A,T61,F20.6)') &
712 17490 : 'MD| Time [fs]', time*femtoseconds
713 : WRITE (iw, '(T2,A,T61,E20.12)') &
714 17490 : 'MD| Conserved quantity [hartree]', md_ener%constant
715 17490 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
716 17490 : WRITE (iw, '(T2,A,T47,A,T73,A)') 'MD|', 'Instantaneous', 'Averages'
717 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
718 17490 : 'MD| CPU time per MD step [s]', used_time, averages%avecpu
719 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
720 17490 : 'MD| Energy drift per atom [K] ', econs, averages%econs
721 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
722 17490 : 'MD| Potential energy [hartree]', md_ener%epot, averages%avepot
723 17490 : IF (simpar%ensemble /= reftraj_ensemble) THEN
724 17375 : IF (.NOT. qmmm) THEN
725 : ! No QM/MM info
726 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
727 16738 : 'MD| Kinetic energy [hartree]', md_ener%ekin, averages%avekin
728 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
729 16738 : 'MD| Temperature [K]', md_ener%temp_part, averages%avetemp
730 : ELSE
731 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
732 637 : 'MD| Total kinetic energy [hartree]', md_ener%ekin, averages%avekin
733 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
734 637 : 'MD| QM kinetic energy [hartree]', md_ener%ekin_qm, averages%avekin_qm
735 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
736 637 : 'MD| Total temperature [K]', md_ener%temp_part, averages%avetemp
737 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
738 637 : 'MD| QM temperature [K]', md_ener%temp_qm, averages%avetemp_qm
739 : END IF
740 : END IF
741 17490 : IF (virial%pv_availability) THEN
742 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
743 1916 : 'MD| Pressure [bar]', pv_scalar, averages%avepress
744 : END IF
745 17490 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
746 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
747 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
748 30 : 'MD| P(xx) [bar]', pv_xx, averages%avepxx
749 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
750 30 : 'MD| Hugoniot [K]', hugoniot/3.0_dp/nat*kelvin, averages%avehugoniot/3.0_dp/nat*kelvin
751 : END IF
752 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
753 : (simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
754 : (simpar%ensemble == npt_i_ensemble) .OR. &
755 : (simpar%ensemble == npt_ia_ensemble) .OR. &
756 : (simpar%ensemble == npt_f_ensemble) .OR. &
757 17490 : (simpar%ensemble == npe_i_ensemble) .OR. &
758 : (simpar%ensemble == npe_f_ensemble)) THEN
759 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
760 954 : 'MD| Barostat temperature [K]', md_ener%temp_baro, averages%avetemp_baro
761 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
762 954 : 'MD| Cell volume [bohr^3]', cell%deth, averages%avevol
763 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
764 954 : 'MD| Cell volume [ang^3]', cell%deth*angstrom**3, averages%avevol*angstrom**3
765 954 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
766 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
767 954 : 'MD| Cell lengths [bohr]', abc(1:3)
768 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
769 3816 : 'MD| Cell lengths [ang]', abc(1:3)*angstrom
770 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
771 954 : 'MD| Average cell lengths [bohr]', averages%aveca, averages%avecb, averages%avecc
772 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
773 954 : 'MD| Average cell lengths [ang]', averages%aveca*angstrom, averages%avecb*angstrom, &
774 1908 : averages%avecc*angstrom
775 : END IF
776 17490 : IF ((simpar%ensemble == npt_f_ensemble) .OR. &
777 : (simpar%ensemble == npe_f_ensemble)) THEN
778 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
779 448 : 'MD| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
780 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
781 448 : 'MD| Average cell angles [deg]', averages%aveal, averages%avebe, averages%avega
782 : END IF
783 17490 : IF (simpar%ensemble == reftraj_ensemble) THEN
784 115 : IF (reftraj%info%msd) THEN
785 6 : IF (reftraj%msd%msd_kind) THEN
786 : WRITE (iw, '(/,T2,A,T51,3F10.5)') &
787 6 : 'MD| COM displacement (dx,dy,dz) [bohr]', reftraj%msd%drcom(1:3)
788 : END IF
789 : END IF
790 : END IF
791 17490 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('*', 75)
792 17490 : IF (max_memory .NE. 0) THEN
793 : WRITE (iw, '(T2,A,T73,I8)') &
794 17490 : 'MD| Estimated peak process memory after this step [MiB]', &
795 34980 : (max_memory + (1024*1024) - 1)/(1024*1024)
796 : END IF
797 : END IF
798 : END IF
799 41372 : CALL section_release(section)
800 : CALL cp_print_key_finished_output(iw, logger, motion_section, &
801 41372 : "MD%PRINT%PROGRAM_RUN_INFO")
802 41372 : END SUBROUTINE md_write_info_low
803 :
804 : ! **************************************************************************************************
805 : !> \brief Samples memory usage
806 : !> \param para_env ...
807 : !> \return ...
808 : !> \note based on what is done in start/cp2k_runs.F
809 : ! **************************************************************************************************
810 53726 : FUNCTION sample_memory(para_env) RESULT(max_memory)
811 : TYPE(mp_para_env_type), POINTER :: para_env
812 : INTEGER(KIND=int_8) :: max_memory
813 :
814 53726 : CALL m_memory()
815 53726 : max_memory = m_memory_max
816 53726 : CALL para_env%max(max_memory)
817 :
818 53726 : END FUNCTION sample_memory
819 :
820 : END MODULE md_energies
|
