Skip to content

Commit

Permalink
Revise stress tensor output (resolves #583)
Browse files Browse the repository at this point in the history
  • Loading branch information
@mkrack
mkrack committed Oct 15, 2020
1 parent 6cd7276 commit f9fe711
Show file tree
Hide file tree
Showing 7 changed files with 179 additions and 166 deletions.
32 changes: 15 additions & 17 deletions src/force_env_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,8 +73,7 @@ MODULE force_env_methods
use_eip_force, use_embed, use_fist_force, use_mixed_force, use_prog_name, use_pwdft_force, &
use_qmmm, use_qmmmx, use_qs_force
USE force_env_utils, ONLY: rescale_forces,&
write_forces,&
write_stress_tensor
write_forces
USE force_fields_util, ONLY: get_generic_info
USE fp_methods, ONLY: fp_eval
USE fparser, ONLY: EvalErrType,&
Expand Down Expand Up @@ -163,7 +162,8 @@ MODULE force_env_methods
USE restraint, ONLY: restraint_control
USE scine_utils, ONLY: write_scine
USE string_utilities, ONLY: compress
USE virial_methods, ONLY: write_stress_components
USE virial_methods, ONLY: write_stress_tensor,&
write_stress_tensor_components
USE virial_types, ONLY: cp_virial,&
sym_virial,&
virial_create,&
Expand Down Expand Up @@ -216,10 +216,8 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
nfixed_atoms_total, nkind, &
output_unit, print_forces, print_grrm, &
print_scine
LOGICAL :: calculate_forces, &
calculate_stress_tensor, &
energy_consistency, eval_ef, &
linres_run, my_skip
LOGICAL :: calculate_forces, calculate_stress_tensor, energy_consistency, eval_ef, &
linres_run, my_skip, print_components
REAL(KIND=dp) :: checksum, e_pot, sum_energy, &
sum_pv_virial, sum_stress_tensor
REAL(KIND=dp), DIMENSION(3) :: grand_total_force, total_force
Expand Down Expand Up @@ -392,12 +390,14 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
extension=".stress_tensor")
IF (output_unit > 0) THEN
CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%NDIGITS", &
i_val=ndigits)
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, ndigits, virial%pv_numer)
IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use == use_qs_force)) THEN
CALL write_stress_components(virial, output_unit)
CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%COMPONENTS", &
l_val=print_components)
IF (print_components) THEN
IF ((.NOT. virial%pv_numer) .AND. (force_env%in_use == use_qs_force)) THEN
CALL write_stress_tensor_components(virial, output_unit, cell)
END IF
END IF
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
END IF
CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
"PRINT%STRESS_TENSOR")
Expand Down Expand Up @@ -553,8 +553,8 @@ SUBROUTINE force_env_calc_num_pressure(force_env, dx)
REAL(kind=dp), PARAMETER :: default_dx = 0.001_dp
INTEGER :: i, ip, iq, j, k, natom, ncore, ndigits, &
nshell, output_unit
INTEGER :: i, ip, iq, j, k, natom, ncore, nshell, &
output_unit
REAL(KIND=dp) :: dx_w
REAL(KIND=dp), DIMENSION(2) :: numer_energy
REAL(KIND=dp), DIMENSION(3) :: s
Expand Down Expand Up @@ -700,9 +700,7 @@ SUBROUTINE force_env_calc_num_pressure(force_env, dx)
IF (output_unit > 0) THEN
IF (globenv%run_type_id == debug_run) THEN
CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%NDIGITS", &
i_val=ndigits)
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, ndigits, virial%pv_numer)
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
END IF
WRITE (output_unit, '(/,A,/)') ' **************************** '// &
'NUMERICAL STRESS END *****************************'
Expand Down
97 changes: 9 additions & 88 deletions src/force_env_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,14 +26,11 @@ MODULE force_env_utils
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
USE mathlib, ONLY: det_3x3,&
jacobi
USE molecule_kind_list_types, ONLY: molecule_kind_list_type
USE molecule_list_types, ONLY: molecule_list_type
USE molecule_types, ONLY: global_constraint_type
USE particle_list_types, ONLY: particle_list_type
USE particle_types, ONLY: update_particle_set
USE physcon, ONLY: pascal
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand All @@ -45,7 +42,6 @@ MODULE force_env_utils
PUBLIC :: force_env_shake, &
force_env_rattle, &
rescale_forces, &
write_stress_tensor, &
write_forces

CONTAINS
Expand Down Expand Up @@ -400,98 +396,23 @@ SUBROUTINE rescale_forces(force_env)
CALL timestop(handle)
END SUBROUTINE rescale_forces

! **************************************************************************************************
!> \brief Variable precision output of the stress tensor
!>
!> \param pv_virial ...
!> \param output_unit ...
!> \param cell ...
!> \param ndigits ...
!> \param numerical ...
!> \author MK (26.08.2010)
! **************************************************************************************************
SUBROUTINE write_stress_tensor(pv_virial, output_unit, cell, ndigits, numerical)

REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: pv_virial
INTEGER, INTENT(IN) :: output_unit
TYPE(cell_type), POINTER :: cell
INTEGER, INTENT(IN) :: ndigits
LOGICAL, INTENT(IN) :: numerical

CHARACTER(LEN=15) :: fmtstr3
CHARACTER(LEN=16) :: fmtstr4
CHARACTER(LEN=22) :: fmtstr2
CHARACTER(LEN=27) :: fmtstr5
CHARACTER(LEN=31) :: fmtstr1
INTEGER :: n
REAL(KIND=dp), DIMENSION(3) :: eigval
REAL(KIND=dp), DIMENSION(3, 3) :: eigvec, stress_tensor

IF (output_unit > 0) THEN
CPASSERT(ASSOCIATED(cell))
stress_tensor(:, :) = pv_virial(:, :)/cell%deth*pascal*1.0E-9_dp
n = MIN(MAX(1, ndigits), 20)
fmtstr1 = "(/,T2,A,/,/,T13,A1,2( X,A1))"
WRITE (UNIT=fmtstr1(22:23), FMT="(I2)") n + 7
fmtstr2 = "(T3,A,T5,3(1X,F . ))"
WRITE (UNIT=fmtstr2(16:17), FMT="(I2)") n + 7
WRITE (UNIT=fmtstr2(19:20), FMT="(I2)") n
fmtstr3 = "(/,T3,A,F . )"
WRITE (UNIT=fmtstr3(10:11), FMT="(I2)") n + 8
WRITE (UNIT=fmtstr3(13:14), FMT="(I2)") n
IF (numerical) THEN
WRITE (UNIT=output_unit, FMT=fmtstr1) &
"NUMERICAL STRESS TENSOR [GPa]", "X", "Y", "Z"
ELSE
WRITE (UNIT=output_unit, FMT=fmtstr1) &
"STRESS TENSOR [GPa]", "X", "Y", "Z"
END IF
WRITE (UNIT=output_unit, FMT=fmtstr2) "X", stress_tensor(1, 1:3)
WRITE (UNIT=output_unit, FMT=fmtstr2) "Y", stress_tensor(2, 1:3)
WRITE (UNIT=output_unit, FMT=fmtstr2) "Z", stress_tensor(3, 1:3)
fmtstr4 = "(/,T3,A,ES . )"
WRITE (UNIT=fmtstr4(11:12), FMT="(I2)") n + 8
WRITE (UNIT=fmtstr4(14:15), FMT="(I2)") n
WRITE (UNIT=output_unit, FMT=fmtstr4) &
"1/3 Trace(stress tensor): ", (stress_tensor(1, 1) + &
stress_tensor(2, 2) + &
stress_tensor(3, 3))/3.0_dp, &
"Det(stress tensor) : ", det_3x3(stress_tensor(:, 1), &
stress_tensor(:, 2), &
stress_tensor(:, 3))
eigval(:) = 0.0_dp
eigvec(:, :) = 0.0_dp
CALL jacobi(stress_tensor, eigval, eigvec)
fmtstr5 = "(/,/,T2,A,/,/,T5,3F . ,/)"
WRITE (UNIT=fmtstr5(20:21), FMT="(I2)") n + 8
WRITE (UNIT=fmtstr5(23:24), FMT="(I2)") n
WRITE (UNIT=output_unit, FMT=fmtstr5) &
"EIGENVECTORS AND EIGENVALUES OF THE STRESS TENSOR", &
eigval(1:3)
WRITE (UNIT=output_unit, FMT=fmtstr2) " ", eigvec(1, 1:3)
WRITE (UNIT=output_unit, FMT=fmtstr2) " ", eigvec(2, 1:3)
WRITE (UNIT=output_unit, FMT=fmtstr2) " ", eigvec(3, 1:3)
END IF

END SUBROUTINE write_stress_tensor

! **************************************************************************************************
!> \brief Write forces
!>
!> \param particles ...
!> \param output_unit ...
!> \param iw ...
!> \param label ...
!> \param ndigits ...
!> \param total_force ...
!> \param grand_total_force ...
!> \param zero_force_core_shell_atom ...
!> \author MK (06.09.2010)
! **************************************************************************************************
SUBROUTINE write_forces(particles, output_unit, label, ndigits, total_force, &
SUBROUTINE write_forces(particles, iw, label, ndigits, total_force, &
grand_total_force, zero_force_core_shell_atom)

TYPE(particle_list_type), POINTER :: particles
INTEGER, INTENT(IN) :: output_unit
INTEGER, INTENT(IN) :: iw
CHARACTER(LEN=*), INTENT(IN) :: label
INTEGER, INTENT(IN) :: ndigits
REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: total_force
Expand All @@ -506,7 +427,7 @@ SUBROUTINE write_forces(particles, output_unit, label, ndigits, total_force, &
LOGICAL :: zero_force
REAL(KIND=dp), DIMENSION(3) :: f

IF (output_unit > 0) THEN
IF (iw > 0) THEN
CPASSERT(ASSOCIATED(particles))
n = MIN(MAX(1, ndigits), 20)
fmtstr1 = "(/,T2,A,/,/,T2,A,T11,A,T18,A,T35,A1,2( X,A1))"
Expand All @@ -522,7 +443,7 @@ SUBROUTINE write_forces(particles, output_unit, label, ndigits, total_force, &
ELSE
zero_force = .FALSE.
END IF
WRITE (UNIT=output_unit, FMT=fmtstr1) &
WRITE (UNIT=iw, FMT=fmtstr1) &
label//" FORCES in [a.u.]", "# Atom", "Kind", "Element", "X", "Y", "Z"
total_force(1:3) = 0.0_dp
DO iparticle = 1, particles%n_els
Expand All @@ -537,18 +458,18 @@ SUBROUTINE write_forces(particles, output_unit, label, ndigits, total_force, &
ELSE
f(1:3) = particles%els(iparticle)%f(1:3)
END IF
WRITE (UNIT=output_unit, FMT=fmtstr2) &
WRITE (UNIT=iw, FMT=fmtstr2) &
i, ikind, particles%els(iparticle)%atomic_kind%element_symbol, f(1:3)
total_force(1:3) = total_force(1:3) + f(1:3)
END DO
WRITE (UNIT=output_unit, FMT=fmtstr3) &
WRITE (UNIT=iw, FMT=fmtstr3) &
"SUM OF "//label//" FORCES", total_force(1:3), SQRT(SUM(total_force(:)**2))
END IF

IF (PRESENT(grand_total_force)) THEN
grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3)
WRITE (UNIT=output_unit, FMT="(A)") ""
WRITE (UNIT=output_unit, FMT=fmtstr3) &
WRITE (UNIT=iw, FMT="(A)") ""
WRITE (UNIT=iw, FMT=fmtstr3) &
"GRAND TOTAL FORCE", grand_total_force(1:3), SQRT(SUM(grand_total_force(:)**2))
END IF

Expand Down
13 changes: 6 additions & 7 deletions src/input_cp2k_force_eval.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -332,7 +332,7 @@ SUBROUTINE create_f_env_print_section(section)
CALL section_release(print_key)

CALL cp_print_key_section_create(print_key, __LOCATION__, "DISTRIBUTION2D", &
description="Controls the printing of the distribution of matrix blocks,...", &
description="Controls the printing of the distribution of matrix blocks, ...", &
print_level=high_print_level, filename="__STD_OUT__")
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)
Expand All @@ -347,12 +347,11 @@ SUBROUTINE create_f_env_print_section(section)
description="Controls the printing of the stress tensor", &
print_level=high_print_level, filename="__STD_OUT__")
CALL keyword_create(keyword, __LOCATION__, &
name="NDIGITS", &
description="Specifies the number of digits used "// &
"for the printing of the stress tensor", &
usage="NDIGITS 6", &
default_i_val=8, &
repeats=.FALSE.)
name="COMPONENTS", &
description="Print all GPW/GAPW components contributing to the stress tensor", &
usage="COMPONENTS", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
Expand Down
2 changes: 1 addition & 1 deletion src/motion/md_energies.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -652,7 +652,7 @@ SUBROUTINE md_write_info_low(simpar, md_ener, qmmm, virial, reftraj, cell, &
(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)') &
WRITE (iw, '(T2,A,T61,F20.6)') &
'MD_INI| Barostat temperature [K]', md_ener%temp_baro
END IF
IF (virial%pv_availability) THEN
Expand Down
4 changes: 2 additions & 2 deletions src/motion/md_vel_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2360,8 +2360,8 @@ SUBROUTINE initialize_cascade(simpar, particle_set, molecule_kinds, md_section)
temperature = cp_unit_from_cp2k(temp, "K")
WRITE (UNIT=iw, FMT="(T2,A)") &
"CASCADE|"
WRITE (UNIT=iw, FMT="(T2,A,T61,F18.2,A2)") &
"CASCADE| Temperature after cascade initialization", temperature, " K"
WRITE (UNIT=iw, FMT="(T2,A,T61,F20.6)") &
"CASCADE| Temperature after cascade initialization [K]", temperature
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
Expand Down
20 changes: 10 additions & 10 deletions src/motion_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,27 +135,27 @@ SUBROUTINE rot_ana(particles, mat, dof, print_section, keep_rotations, mass_weig
IF (PRESENT(inertia)) inertia = Ip_eigval
iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO", extension=".vibLog")
IF (iw > 0) THEN
WRITE (iw, '(/,T2,A)') &
WRITE (UNIT=iw, FMT='(/,T2,A)') &
'ROT| Rotational analysis information'
WRITE (iw, '(T2,A)') &
'ROT| Principal axes and moments of inertia in [a.u.]'
WRITE (iw, '(T2,A,T16,3(1X,I18))') &
WRITE (UNIT=iw, FMT='(T2,A)') &
'ROT| Principal axes and moments of inertia [a.u.]'
WRITE (UNIT=iw, FMT='(T2,A,T14,3(1X,I19))') &
'ROT|', 1, 2, 3
WRITE (iw, '(T2,A,T24,3(1X,ES18.11))') &
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,ES19.11))') &
'ROT| Eigenvalues', Ip_eigval(1:3)
WRITE (iw, '(T2,A,T24,3(1X,F18.12))') &
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
'ROT| x', Ip(1, 1:3)
WRITE (iw, '(T2,A,T24,3(1X,F18.12))') &
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
'ROT| y', Ip(2, 1:3)
WRITE (iw, '(T2,A,T24,3(1X,F18.12))') &
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
'ROT| z', Ip(3, 1:3)
END IF
CALL cp_print_key_finished_output(iw, logger, print_section, "ROTATIONAL_INFO")
iw = cp_print_key_unit_nr(logger, print_section, "ROTATIONAL_INFO/COORDINATES", extension=".vibLog")
IF (iw > 0) THEN
WRITE (iw, '(/,T2,A)') 'ROT| Standard molecule orientation in Angstrom'
WRITE (UNIT=iw, FMT='(/,T2,A)') 'ROT| Standard molecule orientation in Angstrom'
DO iparticle = 1, natoms
WRITE (iw, '(T2,"ROT|",T20,A,T27,3(3X,F15.9))') &
WRITE (UNIT=iw, FMT='(T2,"ROT|",T20,A,T27,3(3X,F15.9))') &
TRIM(particles(iparticle)%atomic_kind%name), &
MATMUL(particles(iparticle)%r, Ip)*angstrom
END DO
Expand Down

0 comments on commit f9fe711

Please sign in to comment.