Skip to content

Commit

Permalink
Revise MD initialization output
Browse files Browse the repository at this point in the history
  • Loading branch information
@mkrack
mkrack committed Oct 14, 2020
1 parent cac8ca0 commit 7bae583
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 38 deletions.
10 changes: 5 additions & 5 deletions src/motion/md_energies.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -627,7 +627,7 @@ SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
CALL keyword_get(keyword, enum=enum)
IF (init) THEN
! Write initial values of quantities of interest
WRITE (iw, '(T2,A)') 'MD_INI| MD initialization'
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
Expand Down Expand Up @@ -681,9 +681,9 @@ SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
'MD_INI| Cell angles [deg]', cell_angle(3), cell_angle(2), cell_angle(1)
ELSE
! Write seuquential values of quantities of interest
!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)') '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)') &
'MD| Step number', itimes
IF (simpar%variable_dt) THEN
Expand Down Expand Up @@ -769,7 +769,7 @@ SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
END IF
END IF
END IF
!MK WRITE (iw, '(T2,A)') 'MD| '//REPEAT('*', 75)
WRITE (iw, '(T2,A)') 'MD| '//REPEAT('*', 75)
END IF
END IF
CALL section_release(section)
Expand Down
78 changes: 45 additions & 33 deletions src/motion/md_vel_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -581,10 +581,6 @@ SUBROUTINE initialize_velocities(simpar, &
! Logging
logger => cp_get_default_logger()
iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".log")
IF (iw > 0) THEN
WRITE (iw, '(T2,A)') &
'MD_VEL| '//TRIM(ADJUSTL(label))
END IF

! Build a list of all fixed atoms (if any)
ALLOCATE (is_fixed(natoms))
Expand Down Expand Up @@ -640,6 +636,8 @@ SUBROUTINE initialize_velocities(simpar, &
END IF

IF (iw > 0) THEN
WRITE (iw, '(/,T2,A)') &
'MD_VEL| '//TRIM(ADJUSTL(label))
! Recompute vcom, ecom and ekin for IO
CALL compute_vcom(part, is_fixed, vcom, ecom)
ekin = compute_ekin(part) - ecom
Expand Down Expand Up @@ -1462,12 +1460,15 @@ SUBROUTINE scale_velocity(subsys, md_ener, temp_expected, temp_tol, iw)
md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
END IF
md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin

IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,T2,A,F10.2,A,F10.2,A)") "Temperature scaled to requested temperature:", &
temp_old, " K ->", md_ener%temp_part, " K"
WRITE (UNIT=iw, FMT='(/,T2,A)') &
'MD_VEL| Temperature scaled to requested temperature'
WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
'MD_VEL| Old temperature [K]', temp_old, &
'MD_VEL| New temperature [K]', md_ener%temp_part
END IF
END IF

END SUBROUTINE scale_velocity

! **************************************************************************************************
Expand Down Expand Up @@ -1546,13 +1547,17 @@ SUBROUTINE scale_velocity_region(md_env, subsys, md_ener, simpar, iw)
md_ener%temp_part = 2.0_dp*md_ener%ekin/REAL(md_ener%nfree, KIND=dp)*kelvin
END IF
md_ener%constant = md_ener%constant - ekin_old + md_ener%ekin

DO ireg = 0, nregions
IF (iw > 0 .AND. temp_new(ireg) > 0.0_dp) THEN
WRITE (UNIT=iw, FMT="(/,T2,A,I5, A,F10.2,A,F10.2,A)") "Temperature region ", ireg, &
" rescaled from:", temp_old(ireg), " K to ", temp_new(ireg), " K"
END IF
END DO
IF (iw > 0) THEN
DO ireg = 0, nregions
IF (temp_new(ireg) > 0.0_dp) THEN
WRITE (UNIT=iw, FMT='(/,T2,A,I0,A)') &
'MD_VEL| Temperature region ', ireg, ' scaled to requested temperature'
WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
'MD_VEL| Old temperature [K]', temp_old(ireg), &
'MD_VEL| New temperature [K]', temp_new(ireg)
END IF
END DO
END IF
DEALLOCATE (temp_new, temp_old)

END SUBROUTINE scale_velocity_region
Expand Down Expand Up @@ -1737,13 +1742,15 @@ SUBROUTINE scale_velocity_internal(subsys, md_ener, temp_expected, temp_tol, iw)
md_ener%temp_shell = 2.0_dp*md_ener%ekin_shell/REAL(md_ener%nfree_shell, KIND=dp)*kelvin
END IF
md_ener%constant = md_ener%constant - ekin_shell_old + md_ener%ekin_shell

IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,T2,A,F10.2,A,F10.2,A)") &
"Temperature shell internal motion scaled to requested temperature:", &
temp_shell_old, " K ->", md_ener%temp_shell, " K"
WRITE (UNIT=iw, FMT='(/,T2,A)') &
'MD_VEL| Temperature of shell internal motion scaled to requested temperature'
WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
'MD_VEL| Old temperature [K]', temp_shell_old, &
'MD_VEL| New temperature [K]', md_ener%temp_shell
END IF
ENDIF
END IF

END SUBROUTINE scale_velocity_internal

! **************************************************************************************************
Expand Down Expand Up @@ -1788,15 +1795,17 @@ SUBROUTINE scale_velocity_baro(md_env, md_ener, temp_expected, temp_tol, iw)
END DO
END DO
END IF

nfree = SIZE(npt, 1)*SIZE(npt, 2)
md_ener%temp_baro = 2.0_dp*md_ener%baro_kin/REAL(nfree, dp)*kelvin
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,T2,A,F10.2,A,F10.2,A)") &
"Temperature of barostat motion scaled to requested temperature:", &
temp_old, " K ->", md_ener%temp_baro, " K"
WRITE (UNIT=iw, FMT='(/,T2,A)') &
'MD_VEL| Temperature of barostat motion scaled to requested temperature'
WRITE (UNIT=iw, FMT='(T2,A,T61,F20.6)') &
'MD_VEL| Old temperature [K]', temp_old, &
'MD_VEL| New temperature [K]', md_ener%temp_baro
END IF
END IF

END SUBROUTINE scale_velocity_baro

! **************************************************************************************************
Expand Down Expand Up @@ -1893,8 +1902,10 @@ SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
! Print COMVEL and COM Position
iw = cp_print_key_unit_nr(logger, md_section, "PRINT%CENTER_OF_MASS", extension=".mdLog")
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,T2,A,(T58,A3,F20.10))") &
"Centre of mass motion (COM):", "x =", md_ener%vcom(1), "y =", md_ener%vcom(2), "z =", md_ener%vcom(3)
WRITE (UNIT=iw, FMT="(/,T2,A)") &
"MD_VEL| Centre of mass motion (COM)"
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
"MD_VEL| VCOM [a.u.]", md_ener%vcom(1:3)
END IF
CALL cp_print_key_finished_output(iw, logger, md_section, "PRINT%CENTER_OF_MASS")

Expand All @@ -1914,9 +1925,9 @@ SUBROUTINE comvel_control(md_ener, force_env, md_section, logger)
CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
"MD| Old VCOM [a.u.] =", vcom_old(1:3)
"MD_VEL| 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)
"MD_VEL| 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 Expand Up @@ -2002,9 +2013,10 @@ SUBROUTINE angvel_control(md_ener, force_env, md_section, logger)
CALL scale_velocity(subsys, md_ener, temp_old, 0.0_dp, iw)
CALL compute_vang(particles%els, is_fixed, rcom, vang_new)
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(T2,'MD| ',A,3F16.10,A)") &
"Old VANG = ", vang(1:3), " a.u.", &
"New VANG = ", vang_new(1:3), " a.u"
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
'MD_VEL| Old VANG [a.u.]', vang(1:3)
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,F16.10))") &
'MD_VEL| New VANG [a.u.]', vang_new(1:3)
END IF
END IF

Expand Down Expand Up @@ -2238,7 +2250,7 @@ SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)

IF (iw > 0) THEN
ekin = cp_unit_from_cp2k(energy, "keV")
WRITE (UNIT=iw, FMT="(T2,A,T61,F20.6)") &
WRITE (UNIT=iw, FMT="(/,T2,A,T61,F20.6)") &
"CASCADE| Energy [keV]", ekin
WRITE (UNIT=iw, FMT="(T2,A)") &
"CASCADE|"
Expand Down Expand Up @@ -2350,8 +2362,8 @@ SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
"CASCADE|"
WRITE (UNIT=iw, FMT="(T2,A,T61,F18.2,A2)") &
"CASCADE| Temperature after cascade initialization", temperature, " K"
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,ES16.8),/)") &
"CASCADE| COM velocity: ", vcom(1:3)
WRITE (UNIT=iw, FMT="(T2,A,T30,3(1X,ES16.8))") &
"CASCADE| COM velocity", vcom(1:3)
!MK ! compute and log rcom and vang if not periodic
!MK CALL force_env_get(force_env,cell=cell)
!MK IF (SUM(cell%perd(1:3)) == 0) THEN
Expand Down

0 comments on commit 7bae583

Please sign in to comment.