Skip to content

Commit

Permalink
Revise MD output format
Browse files Browse the repository at this point in the history
  • Loading branch information
@mkrack
mkrack committed Oct 14, 2020
1 parent 3740dfe commit 61b146e
Show file tree
Hide file tree
Showing 4 changed files with 207 additions and 187 deletions.
225 changes: 115 additions & 110 deletions src/motion/md_energies.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,7 @@ MODULE md_energies
npt_i_ensemble,&
reftraj_ensemble
USE input_cp2k_md, ONLY: create_md_section
USE input_enumeration_types, ONLY: enum_i2c,&
enumeration_type
USE input_enumeration_types, ONLY: enumeration_type
USE input_keyword_types, ONLY: keyword_get,&
keyword_type
USE input_section_types, ONLY: section_get_keyword,&
Expand Down Expand Up @@ -626,145 +625,151 @@ SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
CALL create_md_section(section)
keyword => section_get_keyword(section, "ENSEMBLE")
CALL keyword_get(keyword, enum=enum)

IF (init) THEN
! Write initial values of quantities of interest
WRITE (iw, *)
WRITE (iw, '( A )') ' MD_ENERGIES| Initialization proceeding'
WRITE (iw, *)
WRITE (iw, '( )')
WRITE (iw, '( A,A )') ' ******************************** ', &
'GO CP2K GO! **********************************'
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL POTENTIAL ENERGY', &
'[hartree]', '= ', md_ener%epot
WRITE (iw, '(T2,A)') 'MD_INI| MD initialization'
WRITE (iw, '(T2,A,T61,E20.12)') &
'MD_INI| Potential energy [hartree]', md_ener%epot
IF (simpar%ensemble /= reftraj_ensemble) THEN
IF (.NOT. qmmm) THEN
! NO QM/MM info
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL KINETIC ENERGY', &
'[hartree]', '= ', md_ener%ekin
WRITE (iw, '( A,A,T40,A,T60,1(1X,F20.3) )') ' INITIAL TEMPERATURE', &
'[K]', '= ', md_ener%temp_part
WRITE (iw, '(T2,A,T61,E20.12)') &
'MD_INI| Kinetic energy [hartree]', md_ener%ekin
WRITE (iw, '(T2,A,T61,F20.6)') &
'MD_INI| Temperature [K]', md_ener%temp_part
ELSE
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' TOTAL INITIAL KINETIC ENERGY', &
'[hartree]', '= ', md_ener%ekin
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' QM INITIAL KINETIC ENERGY', &
'[hartree]', '= ', md_ener%ekin_qm
WRITE (iw, '( A,A,T40,A,T60,1(1X,F20.3) )') ' TOTAL INITIAL TEMPERATURE', &
'[K]', '= ', md_ener%temp_part
WRITE (iw, '( A,A,T40,A,T60,1(1X,F20.3) )') ' QM INITIAL TEMPERATURE', &
'[K]', '= ', md_ener%temp_qm
WRITE (iw, '(T2,A,T61,E20.12)') &
'MD_INI| Total kinetic energy [hartree]', md_ener%ekin, &
'MD_INI| QM kinetic energy [hartree]', md_ener%ekin_qm
WRITE (iw, '(T2,A,T61,F20.6)') &
'MD_INI| Total temperature [K]', md_ener%temp_part, &
'MD_INI| QM temperature [K]', md_ener%temp_qm
END IF
END IF
IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
simpar%ensemble == nph_uniaxial_damped_ensemble .OR. &
simpar%ensemble == npt_i_ensemble .OR. &
simpar%ensemble == npt_f_ensemble .OR. &
simpar%ensemble == npe_i_ensemble .OR. &
simpar%ensemble == npe_f_ensemble) &
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL BAROSTAT TEMP', &
'[K]', '= ', md_ener%temp_baro
IF (virial%pv_availability) &
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL PRESSURE', &
'[bar]', '= ', pv_scalar
IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
simpar%ensemble == nph_uniaxial_damped_ensemble) &
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL HUGONIOT CONSTRAINT', &
'[K]', '= ', hugoniot
IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
simpar%ensemble == nph_uniaxial_damped_ensemble) &
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL E0', &
'[hartree]', '= ', simpar%e0
WRITE (iw, '( A,A,T40,A,T60,1(1X,E20.12) )') ' INITIAL VOLUME', &
'[bohr^3]', '= ', cell%deth
WRITE (iw, '( A,A,T29,A,T33,3(1X,E15.7) )') ' INITIAL CELL LNTHS', &
'[bohr]', '= ', abc(1), abc(2), abc(3)
WRITE (iw, '( A,A,T29,A,T33,3(1X,E15.7) )') ' INITIAL CELL ANGLS', &
'[deg]', '= ', cell_angle(3), cell_angle(2), cell_angle(1)
WRITE (iw, '( A,A )') ' ******************************** ', &
'GO CP2K GO! **********************************'
IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
(simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
(simpar%ensemble == npt_i_ensemble) .OR. &
(simpar%ensemble == npt_f_ensemble) .OR. &
(simpar%ensemble == npe_i_ensemble) .OR. &
(simpar%ensemble == npe_f_ensemble)) THEN
WRITE (iw, '(T2,A,T61,F20.3)') &
'MD_INI| Barostat temperature [K]', md_ener%temp_baro
END IF
IF (virial%pv_availability) THEN
WRITE (iw, '(T2,A,T61,ES20.12)') &
'MD_INI| Pressure [bar]', pv_scalar
END IF
IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
(simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
WRITE (iw, '(T2,A,T61,ES20.12)') &
'MD_INI| Hugoniot constraint [K]', hugoniot
END IF
IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
(simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
WRITE (iw, '(T2,A,T61,E20.12)') &
'MD_INI| Total energy [hartree]', simpar%e0
END IF
WRITE (iw, '(T2,A,T61,ES20.12)') &
'MD_INI| Cell volume [bohr^3]', cell%deth
WRITE (iw, '(T2,A,T61,ES20.12)') &
'MD_INI| Cell volume [ang^3]', cell%deth*angstrom**3
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD_INI| Cell lengths [bohr]', abc(1:3)
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD_INI| Cell lengths [ang]', abc(1:3)*angstrom
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD_INI| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
ELSE
! Write seuquential values of quantities of interest
WRITE (iw, '(/,T2,A)') REPEAT('*', 79)
WRITE (iw, '(T2,A,T61,A20)') &
'ENSEMBLE TYPE = ', ADJUSTR(TRIM(enum_i2c(enum, simpar%ensemble)))
!MK WRITE (iw, '(T2,A)') 'MD| '//REPEAT('*', 75)
!MK WRITE (iw, '(T2,A,T61,A20)') &
!MK 'MD| Ensemble type', ADJUSTR(TRIM(enum_i2c(enum, simpar%ensemble)))
WRITE (iw, '(T2,A,T71,I10)') &
'STEP NUMBER = ', itimes
'MD| Step number', itimes
IF (simpar%variable_dt) THEN
WRITE (iw, '(T2,A,T60,1X,F20.6)') &
'TIME STEP [fs] = ', dt*femtoseconds
WRITE (iw, '(T2,A,T61,F20.6)') &
'MD| Time step [fs]', dt*femtoseconds
END IF
WRITE (iw, '(T2,A,T60,1X,F20.6)') &
'TIME [fs] = ', time*femtoseconds
WRITE (iw, '(T2,A,T60,1X,E20.12)') &
'CONSERVED QUANTITY [hartree] = ', md_ener%constant
WRITE (iw, '(/,T47,A,T73,A)') 'INSTANTANEOUS', 'AVERAGES'
WRITE (iw, '(T2,A,T39,2(1X,F20.2))') &
'CPU TIME [s] = ', used_time, averages%avecpu
WRITE (iw, '(T2,A,T61,F20.6)') &
'MD| Time [fs]', time*femtoseconds
WRITE (iw, '(T2,A,T61,E20.12)') &
'MD| Conserved quantity [hartree]', md_ener%constant
WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
WRITE (iw, '(T2,A,T47,A,T73,A)') 'MD|', 'Instantaneous', 'Averages'
WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
'MD| CPU time per MD step [s]', used_time, averages%avecpu
WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
'ENERGY DRIFT PER ATOM [K] = ', econs, averages%econs
'MD| Energy drift per atom [K] ', econs, averages%econs
WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
'POTENTIAL ENERGY[hartree] = ', md_ener%epot, averages%avepot
'MD| Potential energy [hartree]', md_ener%epot, averages%avepot
IF (simpar%ensemble /= reftraj_ensemble) THEN
IF (.NOT. qmmm) THEN
! No QM/MM info
WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
'KINETIC ENERGY [hartree] = ', md_ener%ekin, averages%avekin
WRITE (iw, '(T2,A,T39,2(1X,F20.3))') &
'TEMPERATURE [K] = ', md_ener%temp_part, averages%avetemp
'MD| Kinetic energy [hartree]', md_ener%ekin, averages%avekin
WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
'MD| Temperature [K]', md_ener%temp_part, averages%avetemp
ELSE
WRITE (iw, '( A,A,T31,A,T39,2(1X,E20.12) )') ' TOTAL KINETIC ENERGY', &
'[hartree]', '= ', md_ener%ekin, averages%avekin
WRITE (iw, '( A,A,T31,A,T39,2(1X,E20.12) )') ' QM KINETIC ENERGY', &
'[hartree]', '= ', md_ener%ekin_qm, averages%avekin_qm
WRITE (iw, '( A,A,T31,A,T39,2(1X,F20.3) )') ' TOTAL TEMPERATURE', &
'[K]', '= ', md_ener%temp_part, averages%avetemp
WRITE (iw, '( A,A,T31,A,T39,2(1X,F20.3) )') ' QM TEMPERATURE', &
'[K]', '= ', md_ener%temp_qm, averages%avetemp_qm
WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
'MD| Total kinetic energy [hartree]', md_ener%ekin, averages%avekin
WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
'MD| QM kinetic energy [hartree]', md_ener%ekin_qm, averages%avekin_qm
WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
'MD| Total temperature [K]', md_ener%temp_part, averages%avetemp
WRITE (iw, '(T2,A,T39,2(1X,F20.6))') &
'MD| QM temperature [K]', md_ener%temp_qm, averages%avetemp_qm
END IF
END IF
IF (virial%pv_availability) THEN
WRITE (iw, '(T2,A,T39,2(1X,E20.12))') &
'PRESSURE [bar] = ', pv_scalar, averages%avepress
WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
'MD| Pressure [bar]', pv_scalar, averages%avepress
END IF
IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
simpar%ensemble == nph_uniaxial_damped_ensemble) THEN
WRITE (iw, '( A,A,T31,A,T39,2(1X,E20.12) )') ' P_XX', &
'[bar]', '= ', pv_xx, averages%avepxx
WRITE (iw, '( A,A,T31,A,T39,2(1X,E20.12) )') ' HUGONIOT', &
'[K]', '= ', hugoniot/3._dp/nat*kelvin, &
averages%avehugoniot/3._dp/nat*kelvin
IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
(simpar%ensemble == nph_uniaxial_damped_ensemble)) THEN
WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
'MD| P(xx) [bar]', pv_xx, averages%avepxx
WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
'MD| Hugoniot [K]', hugoniot/3.0_dp/nat*kelvin, averages%avehugoniot/3.0_dp/nat*kelvin
END IF
IF (simpar%ensemble == nph_uniaxial_ensemble .OR. &
simpar%ensemble == nph_uniaxial_damped_ensemble .OR. &
simpar%ensemble == npt_i_ensemble .OR. &
simpar%ensemble == npt_f_ensemble .OR. &
simpar%ensemble == npe_i_ensemble .OR. &
simpar%ensemble == npe_f_ensemble) THEN
WRITE (iw, '( A,A,T31,A,T39,2(1X,E20.12) )') ' BAROSTAT TEMP', &
'[K]', '= ', md_ener%temp_baro, averages%avetemp_baro
WRITE (iw, '( A,A,T31,A,T39,2(1X,E20.12) )') ' VOLUME', &
'[bohr^3]', '= ', cell%deth, averages%avevol
WRITE (iw, '( A,A,T31,A,T33,3(1X,E15.7) )') ' CELL LNTHS', &
'[bohr]', '= ', abc(1), abc(2), abc(3)
WRITE (iw, '( A,A,T31,A,T33,3(1X,E15.7) )') ' AVE. CELL LNTHS', &
'[bohr]', '= ', averages%aveca, averages%avecb, &
averages%avecc
IF ((simpar%ensemble == nph_uniaxial_ensemble) .OR. &
(simpar%ensemble == nph_uniaxial_damped_ensemble) .OR. &
(simpar%ensemble == npt_i_ensemble) .OR. &
(simpar%ensemble == npt_f_ensemble) .OR. &
(simpar%ensemble == npe_i_ensemble) .OR. &
(simpar%ensemble == npe_f_ensemble)) THEN
WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
'MD| Barostat temperature [K]', md_ener%temp_baro, averages%avetemp_baro
WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
'MD| Cell volume [bohr^3]', cell%deth, averages%avevol
WRITE (iw, '(T2,A,T39,2(1X,ES20.12))') &
'MD| Cell volume [ang^3]', cell%deth*angstrom**3, averages%avevol*angstrom**3
WRITE (iw, '(T2,A)') 'MD| '//REPEAT('-', 75)
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD| Cell lengths [bohr]', abc(1:3)
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD| Cell lengths [ang]', abc(1:3)*angstrom
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD| Average cell lengths [bohr]', averages%aveca, averages%avecb, averages%avecc
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD| Average cell lengths [ang]', averages%aveca*angstrom, averages%avecb*angstrom, &
averages%avecc*angstrom
END IF
IF (simpar%ensemble == npt_f_ensemble .OR. &
simpar%ensemble == npe_f_ensemble) THEN
WRITE (iw, '( A,A,T31,A,T33,3(1X,E15.7) )') ' CELL ANGLS', &
'[deg]', '= ', cell_angle(3), cell_angle(2), cell_angle(1)
WRITE (iw, '( A,A,T31,A,T33,3(1X,E15.7) )') ' AVE. CELL ANGLS', &
'[deg]', '= ', averages%aveal, averages%avebe, averages%avega
IF ((simpar%ensemble == npt_f_ensemble) .OR. &
(simpar%ensemble == npe_f_ensemble)) THEN
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
WRITE (iw, '(T2,A,T33,3(1X,ES15.8))') &
'MD| Average cell angles [deg]', averages%aveal, averages%avebe, averages%avega
END IF
IF (simpar%ensemble == reftraj_ensemble) THEN
IF (reftraj%info%msd) THEN
IF (reftraj%msd%msd_kind) &
WRITE (iw, '(/,A,T50,3f10.5)') ' COM displacement (dx,dy,dz) [angstrom]: ', &
reftraj%msd%drcom(1:3)*angstrom
IF (reftraj%msd%msd_kind) THEN
WRITE (iw, '(/,T2,A,T51,3F10.5)') &
'MD| COM displacement (dx,dy,dz) [bohr]', reftraj%msd%drcom(1:3)
END IF
END IF
END IF
WRITE (iw, '(T2,A,/)') REPEAT('*', 79)
!MK WRITE (iw, '(T2,A)') 'MD| '//REPEAT('*', 75)
END IF
END IF
CALL section_release(section)
Expand Down
28 changes: 15 additions & 13 deletions src/motion/md_vel_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -557,9 +557,8 @@ SUBROUTINE initialize_velocities(simpar, &

CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_velocities'

CHARACTER(LEN=default_string_length) :: my_format
INTEGER :: handle, i, ifixd, imolecule_kind, iw, &
natoms, num_x
natoms
INTEGER, ALLOCATABLE, DIMENSION(:) :: is_fixed
LOGICAL :: success
REAL(KIND=dp) :: ecom, ekin, mass, mass_tot, temp, tmp_r1
Expand All @@ -583,9 +582,8 @@ SUBROUTINE initialize_velocities(simpar, &
logger => cp_get_default_logger()
iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
IF (iw > 0) THEN
num_x = (79 - LEN_TRIM(ADJUSTL(label)) - 2)/2
WRITE (my_format, '(A,I0,A,I0,A)') '(1X,', num_x, '("*"),1X,A,1X,', num_x, '("*"))'
WRITE (iw, TRIM(my_format)) TRIM(ADJUSTL(label))
WRITE (iw, '(T2,A)') &
'MD_VEL| '//TRIM(ADJUSTL(label))
END IF

! Build a list of all fixed atoms (if any)
Expand Down Expand Up @@ -652,18 +650,21 @@ SUBROUTINE initialize_velocities(simpar, &
temp = 2.0_dp*ekin/REAL(simpar%nfree, KIND=dp)
END IF
tmp_r1 = cp_unit_from_cp2k(temp, "K")
WRITE (iw, '( A, T61, F18.2, A2 )') ' Initial Temperature ', tmp_r1, " K"
WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )') ' COM velocity:', vcom(1), vcom(2), vcom(3)
WRITE (iw, '(T2,A,T61,F20.6)') &
'MD_VEL| Initial temperature [K]', tmp_r1
WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
'MD_VEL| COM velocity', vcom(1:3)

! compute and log rcom and vang if not periodic
CALL force_env_get(force_env, cell=cell)
IF (SUM(cell%perd(1:3)) == 0) THEN
CALL compute_rcom(part, is_fixed, rcom)
CALL compute_vang(part, is_fixed, rcom, vang)
WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )') ' COM position:', rcom(1), rcom(2), rcom(3)
WRITE (iw, '( A, T21, F20.12 , F20.12 , F20.12 )') ' Angular velocity:', vang(1), vang(2), vang(3)
WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
'MD_VEL| COM position', rcom(1:3)
WRITE (iw, '(T2,A,T30,3(1X,F16.10))') &
'MD_VEL| Angular velocity', vang(1:3)
END IF
WRITE (iw, '( 1X,79("*"),/)')
END IF

DEALLOCATE (is_fixed)
Expand Down Expand Up @@ -1912,9 +1913,10 @@ SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
CALL reset_vcom(subsys, md_ener, vsubtract=vcom_old)
CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(T2,'MD| ',A,3F16.10,A)") &
"Old VCOM = ", vcom_old(1:3), " a.u.", &
"New VCOM = ", md_ener%vcom(1:3), " a.u"
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
"MD| Old VCOM [a.u.] =", vcom_old(1:3)
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
"MD| New VCOM [a.u.] =", md_ener%vcom(1:3)
END IF
END IF
CALL cp_print_key_finished_output(iw, logger, md_section, &
Expand Down

0 comments on commit 61b146e

Please sign in to comment.