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