Line data Source code
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 : ! **************************************************************************************************
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_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 3592 : 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 1796 : 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 1796 : NULLIFY (subsys)
128 1796 : NULLIFY (atomic_kinds, atomic_kind_set, particles, shell_particles)
129 :
130 0 : CPASSERT(ASSOCIATED(md_ener))
131 1796 : CPASSERT(ASSOCIATED(force_env))
132 :
133 1796 : CALL force_env_get(force_env, subsys=subsys)
134 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds, particles=particles, &
135 1796 : shell_particles=shell_particles)
136 1796 : atomic_kind_set => atomic_kinds%els
137 1796 : nkind = SIZE(atomic_kind_set)
138 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
139 1796 : shell_adiabatic=shell_adiabatic)
140 :
141 1796 : md_ener%nfree = simpar%nfree
142 1796 : md_ener%nfree_shell = -HUGE(0)
143 :
144 1796 : IF (shell_adiabatic) THEN
145 132 : md_ener%nfree_shell = 3*(shell_particles%n_els)
146 : END IF
147 :
148 1796 : IF (simpar%temperature_per_kind) THEN
149 108 : ALLOCATE (md_ener%temp_kind(nkind))
150 108 : 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 54 : ALLOCATE (md_ener%temp_shell_kind(nkind))
156 54 : ALLOCATE (md_ener%ekin_shell_kind(nkind))
157 54 : 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 1796 : tshell=shell_adiabatic)
164 1796 : md_ener%epot = 0.0_dp
165 :
166 1796 : 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 89078 : 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 44539 : 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 44539 : NULLIFY (atomic_kinds, atomic_kind_set, force_env, &
191 44539 : particles, subsys, simpar)
192 : CALL get_md_env(md_env=md_env, force_env=force_env, &
193 44539 : simpar=simpar)
194 :
195 : CALL force_env_get(force_env, &
196 44539 : potential_energy=md_ener%epot, subsys=subsys)
197 :
198 44539 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
199 44539 : atomic_kind_set => atomic_kinds%els
200 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
201 44539 : shell_adiabatic=shell_adiabatic)
202 :
203 44539 : tkind = simpar%temperature_per_kind
204 44539 : tshell = shell_adiabatic
205 :
206 44539 : CALL cp_subsys_get(subsys, particles=particles)
207 44539 : natom = particles%n_els
208 :
209 : CALL compute_conserved_quantity(md_env, md_ener, tkind=tkind, &
210 44539 : tshell=tshell, natom=natom)
211 :
212 44539 : 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 288 : 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 288 : CALL zero_md_ener(md_ener, tkind=.FALSE., tshell=.FALSE.)
230 288 : CALL get_md_env(md_env=md_env, force_env=force_env, reftraj=reftraj)
231 :
232 288 : IF (reftraj%info%eval /= REFTRAJ_EVAL_NONE) THEN
233 148 : 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 288 : 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 179548 : 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 44887 : 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 44887 : NULLIFY (logger)
289 89774 : logger => cp_get_default_logger()
290 44887 : CALL timeset(routineN, handle)
291 :
292 : ! Zeroing
293 44887 : hugoniot = 0.0_dp
294 44887 : econs = 0.0_dp
295 44887 : shell_adiabatic = .FALSE.
296 44887 : shell_present = .FALSE.
297 44887 : NULLIFY (motion_section, atomic_kinds, atomic_kind_set, cell, subsys, &
298 44887 : force_env, md_ener, qmmm_env, reftraj, core_particles, particles, &
299 44887 : shell_particles, print_key, root_section, simpar, virial, &
300 44887 : 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 44887 : para_env=para_env, averages=averages, thermal_regions=thermal_regions)
306 :
307 44887 : root_section => force_env%root_section
308 44887 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
309 :
310 44887 : CALL force_env_get(force_env, cell=cell, subsys=subsys, qmmm_env=qmmm_env)
311 :
312 44887 : qmmm = calc_nfree_qm(md_env, md_ener) > 0
313 44887 : is_mixed = (force_env%in_use == use_mixed_force)
314 :
315 44887 : CALL cp_subsys_get(subsys, particles=particles, virial=virial)
316 44887 : nat = particles%n_els
317 44887 : dt = simpar%dt*simpar%dt_fact
318 :
319 : ! Computing the scalar pressure
320 44887 : IF (virial%pv_availability) THEN
321 4744 : pv_scalar = 0._dp
322 18976 : DO i = 1, 3
323 18976 : pv_scalar = pv_scalar + virial%pv_total(i, i)
324 : END DO
325 4744 : pv_scalar = pv_scalar/3._dp/cell%deth
326 4744 : pv_scalar = cp_unit_from_cp2k(pv_scalar, "bar")
327 4744 : pv_xx_nc = virial%pv_total(1, 1)/cell%deth
328 4744 : pv_xx = cp_unit_from_cp2k(virial%pv_total(1, 1)/cell%deth, "bar")
329 : END IF
330 :
331 44887 : CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
332 44887 : atomic_kind_set => atomic_kinds%els
333 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
334 : shell_present=shell_present, &
335 44887 : shell_adiabatic=shell_adiabatic)
336 :
337 44887 : 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 44887 : my_pos = "APPEND"
341 44887 : my_act = "WRITE"
342 44887 : IF (init .AND. (itimes == 0)) THEN
343 1514 : my_pos = "REWIND"
344 1514 : my_act = "WRITE"
345 : END IF
346 :
347 : ! In case of REFTRAJ ensemble the POS is determined differently..
348 : ! according the REFTRAJ counter
349 44887 : IF (simpar%ensemble == reftraj_ensemble) THEN
350 288 : 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 44887 : IF (init) THEN
357 : ! Computing the Hugoniot for NPH calculations
358 1794 : 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 1794 : IF (simpar%ensemble == reftraj_ensemble) reftraj%init = init
366 : ELSE
367 : ! Performing protocol for anything beyond the first step of MD
368 43093 : 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 43093 : IF (simpar%ensemble == reftraj_ensemble) THEN
374 250 : time = reftraj%time
375 250 : econs = md_ener%delta_epot
376 250 : itimes = reftraj%itimes
377 : ELSE
378 42843 : 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 43093 : my_act)
385 : END IF
386 :
387 : ! Sample memory, if requested
388 44887 : CALL section_vals_val_get(motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
389 44887 : max_memory = 0
390 44887 : IF (print_memory) THEN
391 44887 : 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 44887 : hugoniot, nat, init, logger, motion_section, my_pos, my_act, max_memory)
398 :
399 : ! Real Output driven by the PRINT sections
400 44887 : 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 44645 : extension=".ener", file_position=my_pos, file_action=my_act, is_new_file=new_file)
404 44645 : IF (ene > 0) THEN
405 18019 : 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 741 : 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 1482 : "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 18019 : itimes, time*femtoseconds, md_ener%ekin, md_ener%temp_part, md_ener%epot, md_ener%constant, used_time
413 18019 : CALL m_flush(ene)
414 : END IF
415 44645 : CALL cp_print_key_finished_output(ene, logger, motion_section, "MD%PRINT%ENERGY")
416 :
417 : ! Possibly Print MIXED Energy
418 44645 : 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 44645 : 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 44645 : 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 44645 : 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 44645 : "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 44645 : "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 44645 : "FORCE_MIXING_LABELS", my_pos, my_act, middle_name="fmlabels", extended_xmol_title=.TRUE.)
459 :
460 : ! Print Simulation Cell
461 44645 : CALL write_simulation_cell(cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
462 :
463 : ! Print Thermostats status
464 44645 : CALL print_thermostats_status(thermostats, para_env, my_pos, my_act, itimes, time)
465 :
466 : ! Print Barostat status
467 44645 : CALL print_barostat_status(barostat, simpar, my_pos, my_act, cell, itimes, time)
468 :
469 : ! Print Stress Tensor
470 44645 : CALL write_stress_tensor(virial, cell, motion_section, itimes, time*femtoseconds, my_pos, my_act)
471 :
472 : ! Print Polarisability Tensor
473 44645 : IF (ASSOCIATED(force_env%qs_env)) THEN
474 3632 : CALL write_polarisability_tensor(force_env, motion_section, itimes, time*femtoseconds, my_pos, my_act)
475 : END IF
476 :
477 : ! Temperature per Kinds
478 44645 : 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 266 : 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 43745 : print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_KIND")
491 43745 : 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 302 : "Please switch on the keyword MD%TEMP_KIND.")
496 : END IF
497 : !Thermal Region
498 44645 : CALL print_thermal_regions_temperature(thermal_regions, itimes, time*femtoseconds, my_pos, my_act)
499 :
500 : ! Core/Shell Model
501 44645 : IF (shell_present) THEN
502 2386 : CALL force_env_get(force_env, harmonic_shell=harm_shell)
503 2386 : 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 2386 : file_form="FORMATTED", is_new_file=new_file)
509 2386 : IF (shene > 0) THEN
510 1106 : IF (new_file) THEN
511 45 : WRITE (shene, '("#",3X,A,3X,A,3X,3(5X,A,5X))') "Step Nr.", "Time[fs]", "Kin.[a.u.]", &
512 90 : "Temp.[K]", "Pot.[a.u.]"
513 : END IF
514 :
515 : WRITE (shene, "(I8,F12.3,F20.9,F20.9,F20.9,F20.9 )") &
516 1106 : itimes, time*femtoseconds, md_ener%ekin_shell, md_ener%temp_shell, harm_shell
517 1106 : CALL m_flush(shene)
518 : END IF
519 2386 : 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 2386 : "SHELL_TRAJECTORY", my_pos, my_act, "shpos", shell_particles, extended_xmol_title=.TRUE.)
524 :
525 2386 : 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 2386 : "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 2386 : "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 2386 : "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 2386 : "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 2386 : "CORE_FORCES", my_pos, my_act, "cofrc", core_particles, extended_xmol_title=.TRUE.)
545 :
546 : ! Temperature per Kinds
547 2386 : IF (simpar%temperature_per_kind) THEN
548 : tempkind = cp_print_key_unit_nr(logger, motion_section, "MD%PRINT%TEMP_SHELL_KIND", &
549 376 : extension=".shtemp", file_position=my_pos, file_action=my_act)
550 376 : IF (tempkind > 0) THEN
551 21 : nkind = SIZE(md_ener%temp_shell_kind)
552 21 : fmd = "(I10,F20.3,"//TRIM(ADJUSTL(cp_to_string(nkind)))//"F20.9)"
553 : fmd = TRIM(fmd)
554 21 : WRITE (tempkind, fmd) itimes, time*femtoseconds, md_ener%temp_shell_kind(1:nkind)
555 21 : CALL m_flush(tempkind)
556 : END IF
557 : CALL cp_print_key_finished_output(tempkind, logger, motion_section, &
558 376 : "MD%PRINT%TEMP_SHELL_KIND")
559 : ELSE
560 2010 : print_key => section_vals_get_subs_vals(motion_section, "MD%PRINT%TEMP_SHELL_KIND")
561 2010 : 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 80 : "Please switch on the keyword MD%TEMP_KIND.")
566 : END IF
567 : END IF
568 : END IF
569 : END IF
570 44887 : init = .FALSE.
571 44887 : CALL set_md_env(md_env, init=init)
572 44887 : CALL timestop(handle)
573 44887 : 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 44887 : 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 44887 : 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 44887 : extension=".mdLog", file_position=my_pos, file_action=my_act)
640 :
641 : ! Performing protocol relevant to the first step of an MD run
642 44887 : IF (iw > 0) THEN
643 18347 : CALL create_md_section(section)
644 18347 : keyword => section_get_keyword(section, "ENSEMBLE")
645 18347 : CALL keyword_get(keyword, enum=enum)
646 18347 : IF (init) THEN
647 : ! Write initial values of quantities of interest
648 809 : WRITE (iw, '(/,T2,A)') 'MD_INI| MD initialization'
649 : WRITE (iw, '(T2,A,T61,E20.12)') &
650 809 : 'MD_INI| Potential energy [hartree]', md_ener%epot
651 809 : IF (simpar%ensemble /= reftraj_ensemble) THEN
652 790 : IF (.NOT. qmmm) THEN
653 : ! NO QM/MM info
654 : WRITE (iw, '(T2,A,T61,E20.12)') &
655 734 : 'MD_INI| Kinetic energy [hartree]', md_ener%ekin
656 : WRITE (iw, '(T2,A,T61,F20.6)') &
657 734 : 'MD_INI| Temperature [K]', md_ener%temp_part
658 : ELSE
659 : WRITE (iw, '(T2,A,T61,E20.12)') &
660 56 : 'MD_INI| Total kinetic energy [hartree]', md_ener%ekin, &
661 112 : 'MD_INI| QM kinetic energy [hartree]', md_ener%ekin_qm
662 : WRITE (iw, '(T2,A,T61,F20.6)') &
663 56 : 'MD_INI| Total temperature [K]', md_ener%temp_part, &
664 112 : '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 809 : (simpar%ensemble == npe_i_ensemble) .OR. &
673 : (simpar%ensemble == npe_f_ensemble)) THEN
674 : WRITE (iw, '(T2,A,T61,F20.6)') &
675 85 : 'MD_INI| Barostat temperature [K]', md_ener%temp_baro
676 : END IF
677 809 : IF (virial%pv_availability) THEN
678 : WRITE (iw, '(T2,A,T61,ES20.12)') &
679 149 : 'MD_INI| Pressure [bar]', pv_scalar
680 : END IF
681 809 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
682 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
683 : WRITE (iw, '(T2,A,T61,ES20.12)') &
684 3 : 'MD_INI| Hugoniot constraint [K]', hugoniot
685 : END IF
686 809 : IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
687 : (simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
688 : WRITE (iw, '(T2,A,T61,E20.12)') &
689 3 : 'MD_INI| Total energy [hartree]', simpar%e0
690 : END IF
691 : WRITE (iw, '(T2,A,T61,ES20.12)') &
692 809 : 'MD_INI| Cell volume [bohr^3]', cell%deth
693 : WRITE (iw, '(T2,A,T61,ES20.12)') &
694 809 : 'MD_INI| Cell volume [ang^3]', cell%deth*angstrom**3
695 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
696 809 : 'MD_INI| Cell lengths [bohr]', abc(1:3)
697 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
698 3236 : 'MD_INI| Cell lengths [ang]', abc(1:3)*angstrom
699 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
700 809 : '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 17538 : 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 17538 : 'MD| Step number', itimes
708 17538 : IF (simpar%variable_dt) THEN
709 : WRITE (iw, '(T2,A,T61,F20.6)') &
710 240 : 'MD| Time step [fs]', dt*femtoseconds
711 : END IF
712 : WRITE (iw, '(T2,A,T61,F20.6)') &
713 17538 : 'MD| Time [fs]', time*femtoseconds
714 : WRITE (iw, '(T2,A,T61,E20.12)') &
715 17538 : 'MD| Conserved quantity [hartree]', md_ener%constant
716 17538 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
717 17538 : WRITE (iw, '(T2,A,T47,A,T73,A)') 'MD|', 'Instantaneous', 'Averages'
718 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
719 17538 : 'MD| CPU time per MD step [s]', used_time, averages%avecpu
720 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
721 17538 : 'MD| Energy drift per atom [K] ', econs, averages%econs
722 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
723 17538 : 'MD| Potential energy [hartree]', md_ener%epot, averages%avepot
724 17538 : IF (simpar%ensemble /= reftraj_ensemble) THEN
725 17413 : IF (.NOT. qmmm) THEN
726 : ! No QM/MM info
727 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
728 16774 : 'MD| Kinetic energy [hartree]', md_ener%ekin, averages%avekin
729 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
730 16774 : 'MD| Temperature [K]', md_ener%temp_part, averages%avetemp
731 : ELSE
732 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
733 639 : 'MD| Total kinetic energy [hartree]', md_ener%ekin, averages%avekin
734 : WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
735 639 : 'MD| QM kinetic energy [hartree]', md_ener%ekin_qm, averages%avekin_qm
736 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
737 639 : 'MD| Total temperature [K]', md_ener%temp_part, averages%avetemp
738 : WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
739 639 : 'MD| QM temperature [K]', md_ener%temp_qm, averages%avetemp_qm
740 : END IF
741 : END IF
742 17538 : IF (virial%pv_availability) THEN
743 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
744 1942 : 'MD| Pressure [bar]', pv_scalar, averages%avepress
745 : END IF
746 17538 : 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 30 : 'MD| P(xx) [bar]', pv_xx, averages%avepxx
750 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
751 30 : '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 17538 : (simpar%ensemble == npe_i_ensemble) .OR. &
759 : (simpar%ensemble == npe_f_ensemble)) THEN
760 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
761 980 : 'MD| Barostat temperature [K]', md_ener%temp_baro, averages%avetemp_baro
762 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
763 980 : 'MD| Cell volume [bohr^3]', cell%deth, averages%avevol
764 : WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
765 980 : 'MD| Cell volume [ang^3]', cell%deth*angstrom**3, averages%avevol*angstrom**3
766 980 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
767 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
768 980 : 'MD| Cell lengths [bohr]', abc(1:3)
769 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
770 3920 : 'MD| Cell lengths [ang]', abc(1:3)*angstrom
771 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
772 980 : 'MD| Average cell lengths [bohr]', averages%aveca, averages%avecb, averages%avecc
773 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
774 980 : 'MD| Average cell lengths [ang]', averages%aveca*angstrom, averages%avecb*angstrom, &
775 1960 : averages%avecc*angstrom
776 : END IF
777 17538 : 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 448 : 'MD| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
781 : WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
782 448 : 'MD| Average cell angles [deg]', averages%aveal, averages%avebe, averages%avega
783 : END IF
784 17538 : IF (simpar%ensemble == reftraj_ensemble) THEN
785 125 : IF (reftraj%info%msd) THEN
786 6 : IF (reftraj%msd%msd_kind) THEN
787 : WRITE (iw, '(/,T2,A,T51,3F10.5)') &
788 6 : 'MD| COM displacement (dx,dy,dz) [bohr]', reftraj%msd%drcom(1:3)
789 : END IF
790 : END IF
791 : END IF
792 17538 : WRITE (iw, '(T2,A)') 'MD| '//REPEAT('*', 75)
793 17538 : IF (max_memory .NE. 0) THEN
794 : WRITE (iw, '(T2,A,T73,I8)') &
795 17538 : 'MD| Estimated peak process memory after this step [MiB]', &
796 35076 : (max_memory + (1024*1024) - 1)/(1024*1024)
797 : END IF
798 : END IF
799 : END IF
800 44887 : CALL section_release(section)
801 : CALL cp_print_key_finished_output(iw, logger, motion_section, &
802 44887 : "MD%PRINT%PROGRAM_RUN_INFO")
803 44887 : 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 58866 : 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 58866 : CALL m_memory()
816 58866 : max_memory = m_memory_max
817 58866 : CALL para_env%max(max_memory)
818 :
819 58866 : END FUNCTION sample_memory
820 :
821 : END MODULE md_energies
|
