mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-27 21:55:16 -04:00
Allow to specify the physical unit for energy and stress
This affects only the output quantities tagged with ENERGY| and STRESS|.
This commit is contained in:
Matthias Krack
parent
e8c0fe6317
commit
097ac98af0
7 changed files with 172 additions and 106 deletions
|
|
@ -65,6 +65,7 @@ MODULE energy_corrections
|
|||
USE cp_result_methods, ONLY: cp_results_erase,&
|
||||
put_results
|
||||
USE cp_result_types, ONLY: cp_result_type
|
||||
USE cp_units, ONLY: cp_unit_from_cp2k
|
||||
USE distribution_1d_types, ONLY: distribution_1d_type
|
||||
USE distribution_2d_types, ONLY: distribution_2d_type
|
||||
USE dm_ls_scf, ONLY: calculate_w_matrix_ls,&
|
||||
|
|
@ -625,6 +626,7 @@ CONTAINS
|
|||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_dc_build_ks_matrix_force'
|
||||
|
||||
CHARACTER(LEN=default_string_length) :: unit_string
|
||||
INTEGER :: handle, i, iounit, ispin, natom, nspins
|
||||
LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
|
||||
use_virial
|
||||
|
|
@ -1111,8 +1113,9 @@ CONTAINS
|
|||
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
|
||||
'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
|
||||
|
||||
CALL write_stress_tensor_components(virdeb, iounit, cell)
|
||||
CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, .FALSE.)
|
||||
unit_string = "GPa" ! old default
|
||||
CALL write_stress_tensor_components(virdeb, iounit, cell, unit_string)
|
||||
CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, unit_string, .FALSE.)
|
||||
|
||||
END BLOCK
|
||||
END IF
|
||||
|
|
@ -1706,6 +1709,7 @@ CONTAINS
|
|||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_build_ks_matrix_force'
|
||||
|
||||
CHARACTER(LEN=default_string_length) :: unit_string
|
||||
INTEGER :: handle, i, iounit, ispin, natom, nspins
|
||||
LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
|
||||
use_virial
|
||||
|
|
@ -1770,7 +1774,10 @@ CONTAINS
|
|||
nspins = dft_control%nspins
|
||||
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
|
||||
|
||||
fconv = 1.0E-9_dp*pascal/cell%deth
|
||||
! Conversion factor a.u. -> GPa
|
||||
unit_string = "GPa"
|
||||
fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, TRIM(unit_string))
|
||||
|
||||
IF (debug_stress .AND. use_virial) THEN
|
||||
sttot = virial%pv_virial
|
||||
END IF
|
||||
|
|
@ -2352,8 +2359,8 @@ CONTAINS
|
|||
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
|
||||
'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
|
||||
|
||||
CALL write_stress_tensor_components(virdeb, iounit, cell)
|
||||
CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, .FALSE.)
|
||||
CALL write_stress_tensor_components(virdeb, iounit, cell, unit_string)
|
||||
CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, unit_string, .FALSE.)
|
||||
|
||||
END BLOCK
|
||||
END IF
|
||||
|
|
|
|||
|
|
@ -55,6 +55,7 @@ MODULE force_env_methods
|
|||
cp_subsys_p_type,&
|
||||
cp_subsys_set,&
|
||||
cp_subsys_type
|
||||
USE cp_units, ONLY: cp_unit_from_cp2k
|
||||
USE eip_environment_types, ONLY: eip_environment_type
|
||||
USE eip_silicon, ONLY: eip_bazant,&
|
||||
eip_lenosky
|
||||
|
|
@ -207,7 +208,7 @@ CONTAINS
|
|||
LOGICAL :: calculate_forces, calculate_stress_tensor, energy_consistency, eval_ef, &
|
||||
linres_run, my_skip, print_components
|
||||
REAL(KIND=dp) :: checksum, e_entropy, e_gap, e_pot, &
|
||||
sum_energy
|
||||
fconv, sum_energy
|
||||
REAL(KIND=dp), DIMENSION(3) :: grand_total_force, total_force
|
||||
TYPE(atprop_type), POINTER :: atprop_env
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
|
|
@ -334,19 +335,25 @@ CONTAINS
|
|||
|
||||
CALL force_env_get(force_env, potential_energy=e_pot)
|
||||
|
||||
! Print always Energy in the same format for all methods
|
||||
! Print energy always in the same format for all methods
|
||||
output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO", &
|
||||
extension=".Log")
|
||||
IF (output_unit > 0) THEN
|
||||
WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) energy [a.u.]: ",T55,F26.15)') &
|
||||
ADJUSTR(TRIM(use_prog_name(force_env%in_use))), e_pot
|
||||
IF (e_gap .GT. -.1_dp) THEN
|
||||
WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) gap [a.u.]: ",T55,F26.15)') &
|
||||
ADJUSTR(TRIM(use_prog_name(force_env%in_use))), e_gap
|
||||
CALL section_vals_val_get(force_env%force_env_section, "PRINT%PROGRAM_RUN_INFO%ENERGY_UNIT", &
|
||||
c_val=unit_string)
|
||||
fconv = cp_unit_from_cp2k(1.0_dp, TRIM(ADJUSTL(unit_string)))
|
||||
WRITE (UNIT=output_unit, FMT="(/,T2,A,T55,F26.15)") &
|
||||
"ENERGY| Total FORCE_EVAL ( "//TRIM(ADJUSTL(use_prog_name(force_env%in_use)))// &
|
||||
" ) energy ["//TRIM(ADJUSTL(unit_string))//"]", e_pot*fconv
|
||||
IF (e_gap > -0.1_dp) THEN
|
||||
WRITE (UNIT=output_unit, FMT="(/,T2,A,T55,F26.15)") &
|
||||
"ENERGY| Total FORCE_EVAL ( "//TRIM(ADJUSTL(use_prog_name(force_env%in_use)))// &
|
||||
" ) gap ["//TRIM(ADJUSTL(unit_string))//"]", e_gap*fconv
|
||||
END IF
|
||||
IF (e_entropy .GT. -0.1_dp) THEN
|
||||
WRITE (output_unit, '(/,T2,"ENERGY| Total FORCE_EVAL ( ",A," ) free energy [a.u.]: ",T55,F26.15)') &
|
||||
ADJUSTR(TRIM(use_prog_name(force_env%in_use))), e_pot - e_entropy
|
||||
IF (e_entropy > -0.1_dp) THEN
|
||||
WRITE (UNIT=output_unit, FMT="(/,T2,A,T55,F26.15)") &
|
||||
"ENERGY| Total FORCE_EVAL ( "//TRIM(ADJUSTL(use_prog_name(force_env%in_use)))// &
|
||||
" ) free energy ["//TRIM(ADJUSTL(unit_string))//"]", (e_pot - e_entropy)*fconv
|
||||
END IF
|
||||
END IF
|
||||
CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
|
||||
|
|
@ -368,7 +375,7 @@ CONTAINS
|
|||
! Variable precision output of the forces
|
||||
CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%NDIGITS", &
|
||||
i_val=ndigits)
|
||||
CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%UNIT", &
|
||||
CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%FORCE_UNIT", &
|
||||
c_val=unit_string)
|
||||
IF (ASSOCIATED(core_particles) .OR. ASSOCIATED(shell_particles)) THEN
|
||||
CALL write_forces(particles, print_forces, "ATOMIC", ndigits, unit_string, total_force, &
|
||||
|
|
@ -400,12 +407,14 @@ CONTAINS
|
|||
IF (output_unit > 0) THEN
|
||||
CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%COMPONENTS", &
|
||||
l_val=print_components)
|
||||
CALL section_vals_val_get(force_env%force_env_section, "PRINT%STRESS_TENSOR%STRESS_UNIT", &
|
||||
c_val=unit_string)
|
||||
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)
|
||||
CALL write_stress_tensor_components(virial, output_unit, cell, unit_string)
|
||||
END IF
|
||||
END IF
|
||||
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
|
||||
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
|
||||
END IF
|
||||
CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
|
||||
"PRINT%STRESS_TENSOR")
|
||||
|
|
@ -496,6 +505,7 @@ CONTAINS
|
|||
|
||||
REAL(kind=dp), PARAMETER :: default_dx = 0.001_dp
|
||||
|
||||
CHARACTER(LEN=default_string_length) :: unit_string
|
||||
INTEGER :: i, ip, iq, j, k, natom, ncore, nshell, &
|
||||
output_unit, symmetry_id
|
||||
REAL(KIND=dp) :: dx_w
|
||||
|
|
@ -536,8 +546,8 @@ CONTAINS
|
|||
output_unit = cp_print_key_unit_nr(logger, force_env%force_env_section, "PRINT%STRESS_TENSOR", &
|
||||
extension=".stress_tensor")
|
||||
IF (output_unit > 0) THEN
|
||||
WRITE (output_unit, '(/A,A/)') ' **************************** ', &
|
||||
'NUMERICAL STRESS ********************************'
|
||||
WRITE (output_unit, "(/A,A/)") " **************************** ", &
|
||||
"NUMERICAL STRESS ********************************"
|
||||
END IF
|
||||
|
||||
! Save all original particle positions
|
||||
|
|
@ -649,10 +659,12 @@ CONTAINS
|
|||
|
||||
IF (output_unit > 0) THEN
|
||||
IF (globenv%run_type_id == debug_run) THEN
|
||||
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, virial%pv_numer)
|
||||
CALL section_vals_val_get(force_env%force_env_section, "PRINT%FORCES%FORCE_UNIT", &
|
||||
c_val=unit_string)
|
||||
CALL write_stress_tensor(virial%pv_virial, output_unit, cell, unit_string, virial%pv_numer)
|
||||
END IF
|
||||
WRITE (output_unit, '(/,A,/)') ' **************************** '// &
|
||||
'NUMERICAL STRESS END *****************************'
|
||||
WRITE (output_unit, "(/,A,/)") " **************************** "// &
|
||||
"NUMERICAL STRESS END *****************************"
|
||||
END IF
|
||||
|
||||
CALL cp_print_key_finished_output(output_unit, logger, force_env%force_env_section, &
|
||||
|
|
|
|||
|
|
@ -420,7 +420,7 @@ CONTAINS
|
|||
INTEGER, INTENT(IN) :: iw
|
||||
CHARACTER(LEN=*), INTENT(IN) :: label
|
||||
INTEGER, INTENT(IN) :: ndigits
|
||||
CHARACTER(LEN=default_string_length) :: unit_string
|
||||
CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: total_force
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
|
||||
OPTIONAL :: grand_total_force
|
||||
|
|
|
|||
|
|
@ -296,17 +296,29 @@ CONTAINS
|
|||
|
||||
CPASSERT(.NOT. ASSOCIATED(section))
|
||||
|
||||
CALL section_create(section, __LOCATION__, name="PRINT", &
|
||||
CALL section_create(section, __LOCATION__, &
|
||||
name="PRINT", &
|
||||
description="Properties that you want to output and that are common to all methods", &
|
||||
n_keywords=0, n_subsections=5, repeats=.FALSE.)
|
||||
n_keywords=0, n_subsections=10, repeats=.FALSE.)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "PROGRAM_RUN_INFO", &
|
||||
description="Controls the printing of basic information generated by force_eval", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="PROGRAM_RUN_INFO", &
|
||||
description="Controls the printing of basic information generated by FORCE_EVAL", &
|
||||
print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
|
||||
CALL keyword_create(keyword, __LOCATION__, &
|
||||
name="ENERGY_UNIT", &
|
||||
description="Specifies the physical unit used for the printing of the total energy. "// &
|
||||
"Note that the meaningfulness of the unit is not checked.", &
|
||||
usage="ENERGY_UNIT eV", &
|
||||
default_c_val="hartree", &
|
||||
repeats=.FALSE.)
|
||||
CALL section_add_keyword(print_key, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "FORCES", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="FORCES", &
|
||||
description="Controls the printing of the forces after each force evaluation", &
|
||||
print_level=high_print_level, filename="__STD_OUT__")
|
||||
CALL keyword_create(keyword, __LOCATION__, &
|
||||
|
|
@ -319,10 +331,11 @@ CONTAINS
|
|||
CALL section_add_keyword(print_key, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
CALL keyword_create(keyword, __LOCATION__, &
|
||||
name="UNIT", &
|
||||
name="FORCE_UNIT", &
|
||||
variants=(/"UNIT"/), & ! add old keyword name for backward compatibility
|
||||
description="Specifies the physical unit used for the printing of the forces. "// &
|
||||
"Note that the meaningfulness of the unit is not checked.", &
|
||||
usage="UNIT eV/angstrom", &
|
||||
usage="FORCE_UNIT eV/angstrom", &
|
||||
default_c_val="hartree/bohr", &
|
||||
repeats=.FALSE.)
|
||||
CALL section_add_keyword(print_key, keyword)
|
||||
|
|
@ -330,38 +343,43 @@ CONTAINS
|
|||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create( &
|
||||
print_key, __LOCATION__, "GRID_INFORMATION", &
|
||||
description="Controls the printing of information regarding the PW and RS grid structures.", &
|
||||
print_level=medium_print_level, filename="__STD_OUT__")
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="GRID_INFORMATION", &
|
||||
description="Controls the printing of information regarding the PW and RS grid structures.", &
|
||||
print_level=medium_print_level, filename="__STD_OUT__")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "TOTAL_NUMBERS", &
|
||||
description="Controls the printing of the total number of atoms, kinds,...", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="TOTAL_NUMBERS", &
|
||||
description="Controls the printing of the total number of atoms, kinds, ...", &
|
||||
print_level=low_print_level, filename="__STD_OUT__")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "DISTRIBUTION", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="DISTRIBUTION", &
|
||||
description="Controls the printing of the distribution of molecules, atoms, ...", &
|
||||
print_level=high_print_level, filename="__STD_OUT__")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "DISTRIBUTION2D", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="DISTRIBUTION2D", &
|
||||
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)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "DISTRIBUTION1D", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="DISTRIBUTION1D", &
|
||||
description="Each node prints out its distribution info ...", &
|
||||
print_level=high_print_level, filename="__STD_OUT__")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "STRESS_TENSOR", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="STRESS_TENSOR", &
|
||||
description="Controls the printing of the stress tensor", &
|
||||
print_level=high_print_level, filename="__STD_OUT__")
|
||||
CALL keyword_create(keyword, __LOCATION__, &
|
||||
|
|
@ -372,16 +390,27 @@ CONTAINS
|
|||
lone_keyword_l_val=.TRUE.)
|
||||
CALL section_add_keyword(print_key, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
CALL keyword_create(keyword, __LOCATION__, &
|
||||
name="STRESS_UNIT", &
|
||||
description="Specifies the physical unit used for the printing of the stress tensor. "// &
|
||||
"Note that the meaningfulness of the unit is not checked.", &
|
||||
usage="STRESS_UNIT kbar", &
|
||||
default_c_val="GPa", &
|
||||
repeats=.FALSE.)
|
||||
CALL section_add_keyword(print_key, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "GRRM", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="GRRM", &
|
||||
description="Controls the printing of the GRRM interface file", &
|
||||
print_level=debug_print_level + 1, filename="CP2K_GRRM")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "SCINE", &
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, &
|
||||
name="SCINE", &
|
||||
description="Controls the printing of the SCINE interface file", &
|
||||
print_level=debug_print_level + 1, filename="")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
|
|
|
|||
|
|
@ -250,7 +250,7 @@ CONTAINS
|
|||
CALL section_release(print_key)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "PROGRAM_RUN_INFO", &
|
||||
description="controls the printing of tests output", &
|
||||
description="Controls the printing of the test output", &
|
||||
print_level=silent_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
|
@ -791,7 +791,7 @@ CONTAINS
|
|||
CALL section_release(sub_section)
|
||||
|
||||
CALL cp_print_key_section_create(print_key, __LOCATION__, "PROGRAM_RUN_INFO", &
|
||||
description="controls the printing of FARMING specific output", &
|
||||
description="Controls the printing of FARMING specific output", &
|
||||
print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
|
||||
CALL section_add_subsection(section, print_key)
|
||||
CALL section_release(print_key)
|
||||
|
|
|
|||
|
|
@ -17,15 +17,16 @@ MODULE virial_methods
|
|||
USE cell_types, ONLY: cell_type
|
||||
USE cp_subsys_types, ONLY: cp_subsys_get,&
|
||||
cp_subsys_type
|
||||
USE cp_units, ONLY: cp_unit_from_cp2k
|
||||
USE distribution_1d_types, ONLY: distribution_1d_type
|
||||
USE kinds, ONLY: dp
|
||||
USE kinds, ONLY: default_string_length,&
|
||||
dp
|
||||
USE mathlib, ONLY: det_3x3,&
|
||||
jacobi
|
||||
USE message_passing, ONLY: mp_comm_type,&
|
||||
mp_para_env_type
|
||||
USE particle_list_types, ONLY: particle_list_type
|
||||
USE particle_types, ONLY: particle_type
|
||||
USE physcon, ONLY: pascal
|
||||
USE virial_types, ONLY: virial_type
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
|
|
@ -48,6 +49,7 @@ CONTAINS
|
|||
!> \author Teodoro Laino [tlaino] - 03.2008 - Zurich University
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE virial_update(virial, subsys, para_env)
|
||||
|
||||
TYPE(virial_type), INTENT(INOUT) :: virial
|
||||
TYPE(cp_subsys_type), POINTER :: subsys
|
||||
TYPE(mp_para_env_type), POINTER :: para_env
|
||||
|
|
@ -78,6 +80,7 @@ CONTAINS
|
|||
! **************************************************************************************************
|
||||
SUBROUTINE virial_evaluate(atomic_kind_set, particle_set, local_particles, &
|
||||
virial, igroup)
|
||||
|
||||
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
|
||||
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
|
||||
TYPE(distribution_1d_type), POINTER :: local_particles
|
||||
|
|
@ -85,7 +88,7 @@ CONTAINS
|
|||
|
||||
CLASS(mp_comm_type), INTENT(IN) :: igroup
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'virial_evaluate'
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = "virial_evaluate"
|
||||
|
||||
INTEGER :: handle, i, iparticle, iparticle_kind, &
|
||||
iparticle_local, j, nparticle_kind, &
|
||||
|
|
@ -155,76 +158,79 @@ CONTAINS
|
|||
!> \param virial ...
|
||||
!> \param iw ...
|
||||
!> \param cell ...
|
||||
!> \param unit_string ...
|
||||
!> \par History
|
||||
!> - Revised virial components (14.10.2020, MK)
|
||||
!> \author JGH
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE write_stress_tensor_components(virial, iw, cell)
|
||||
SUBROUTINE write_stress_tensor_components(virial, iw, cell, unit_string)
|
||||
|
||||
TYPE(virial_type), INTENT(IN) :: virial
|
||||
INTEGER, INTENT(IN) :: iw
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: fmt = '(T2,A,T41,2(1X,ES19.11))'
|
||||
CHARACTER(LEN=*), PARAMETER :: fmt = "(T2,A,T41,2(1X,ES19.11))"
|
||||
|
||||
REAL(KIND=dp) :: fconv
|
||||
REAL(KIND=dp), DIMENSION(3, 3) :: pv
|
||||
|
||||
IF (iw > 0) THEN
|
||||
CPASSERT(ASSOCIATED(cell))
|
||||
! Conversion factor to GPa
|
||||
fconv = 1.0E-9_dp*pascal/cell%deth
|
||||
WRITE (UNIT=iw, FMT='(/,T2,A)') &
|
||||
'STRESS| Stress tensor components (GPW/GAPW) [GPa]'
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T52,A,T70,A)') &
|
||||
'STRESS|', '1/3 Trace', 'Determinant'
|
||||
! Conversion factor
|
||||
fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, TRIM(ADJUSTL(unit_string)))
|
||||
WRITE (UNIT=iw, FMT="(/,T2,A)") &
|
||||
"STRESS| Stress tensor components (GPW/GAPW) ["//TRIM(ADJUSTL(unit_string))//"]"
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T52,A,T70,A)") &
|
||||
"STRESS|", "1/3 Trace", "Determinant"
|
||||
pv = fconv*virial%pv_overlap
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Overlap', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Overlap", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_ekinetic
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Kinetic energy', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Kinetic energy", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_ppl
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Local pseudopotential/core', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Local pseudopotential/core", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_ppnl
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Nonlocal pseudopotential', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Nonlocal pseudopotential", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_ecore_overlap
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Core charge overlap', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Core charge overlap", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_ehartree
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Hartree', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Hartree", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_exc
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Exchange-correlation', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Exchange-correlation", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_exx
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Exact exchange (EXX)', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Exact exchange (EXX)", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_vdw
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| vdW correction', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| vdW correction", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_mp2
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Moller-Plesset (MP2)', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Moller-Plesset (MP2)", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_nlcc
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Nonlinear core correction (NLCC)', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Nonlinear core correction (NLCC)", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_gapw
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Local atomic parts (GAPW)', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Local atomic parts (GAPW)", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_lrigpw
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Resolution-of-the-identity (LRI)', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Resolution-of-the-identity (LRI)", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*(virial%pv_overlap + virial%pv_ekinetic + virial%pv_ppl + virial%pv_ppnl + &
|
||||
virial%pv_ecore_overlap + virial%pv_ehartree + virial%pv_exc + &
|
||||
virial%pv_exx + virial%pv_vdw + virial%pv_mp2 + virial%pv_nlcc + &
|
||||
virial%pv_gapw + virial%pv_lrigpw)
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Sum of components', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Sum of components", one_third_sum_diag(pv), det_3x3(pv)
|
||||
pv = fconv*virial%pv_virial
|
||||
WRITE (UNIT=iw, FMT=fmt) &
|
||||
'STRESS| Total', one_third_sum_diag(pv), det_3x3(pv)
|
||||
"STRESS| Total", one_third_sum_diag(pv), det_3x3(pv)
|
||||
END IF
|
||||
|
||||
END SUBROUTINE write_stress_tensor_components
|
||||
|
|
@ -247,68 +253,74 @@ CONTAINS
|
|||
!> \param pv_virial ...
|
||||
!> \param iw ...
|
||||
!> \param cell ...
|
||||
!> \param unit_string ...
|
||||
!> \param numerical ...
|
||||
!> \author MK (26.08.2010)
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE write_stress_tensor(pv_virial, iw, cell, numerical)
|
||||
SUBROUTINE write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
|
||||
|
||||
REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: pv_virial
|
||||
INTEGER, INTENT(IN) :: iw
|
||||
TYPE(cell_type), POINTER :: cell
|
||||
CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
|
||||
LOGICAL, INTENT(IN) :: numerical
|
||||
|
||||
REAL(KIND=dp) :: fconv
|
||||
REAL(KIND=dp), DIMENSION(3) :: eigval
|
||||
REAL(KIND=dp), DIMENSION(3, 3) :: eigvec, stress_tensor
|
||||
|
||||
IF (iw > 0) THEN
|
||||
CPASSERT(ASSOCIATED(cell))
|
||||
stress_tensor(:, :) = pv_virial(:, :)/cell%deth*pascal*1.0E-9_dp
|
||||
! Conversion factor
|
||||
fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, TRIM(ADJUSTL(unit_string)))
|
||||
stress_tensor(:, :) = fconv*pv_virial(:, :)
|
||||
IF (numerical) THEN
|
||||
WRITE (UNIT=iw, FMT='(/,T2,A)') &
|
||||
'STRESS| Numerical stress tensor [GPa]'
|
||||
WRITE (UNIT=iw, FMT="(/,T2,A)") &
|
||||
"STRESS| Numerical stress tensor ["//TRIM(ADJUSTL(unit_string))//"]"
|
||||
ELSE
|
||||
WRITE (UNIT=iw, FMT='(/,T2,A)') &
|
||||
'STRESS| Analytical stress tensor [GPa]'
|
||||
WRITE (UNIT=iw, FMT="(/,T2,A)") &
|
||||
"STRESS| Analytical stress tensor ["//TRIM(ADJUSTL(unit_string))//"]"
|
||||
END IF
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T14,3(19X,A1))') &
|
||||
'STRESS|', 'x', 'y', 'z'
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,ES19.11))') &
|
||||
'STRESS| x', stress_tensor(1, 1:3)
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,ES19.11))') &
|
||||
'STRESS| y', stress_tensor(2, 1:3)
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,ES19.11))') &
|
||||
'STRESS| z', stress_tensor(3, 1:3)
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T61,ES20.11)') &
|
||||
'STRESS| 1/3 Trace', (stress_tensor(1, 1) + &
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T14,3(19X,A1))") &
|
||||
"STRESS|", "x", "y", "z"
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
|
||||
"STRESS| x", stress_tensor(1, 1:3)
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
|
||||
"STRESS| y", stress_tensor(2, 1:3)
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
|
||||
"STRESS| z", stress_tensor(3, 1:3)
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.11)") &
|
||||
"STRESS| 1/3 Trace", (stress_tensor(1, 1) + &
|
||||
stress_tensor(2, 2) + &
|
||||
stress_tensor(3, 3))/3.0_dp
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T61,ES20.11)') &
|
||||
'STRESS| Determinant', det_3x3(stress_tensor(1:3, 1), &
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.11)") &
|
||||
"STRESS| Determinant", det_3x3(stress_tensor(1:3, 1), &
|
||||
stress_tensor(1:3, 2), &
|
||||
stress_tensor(1:3, 3))
|
||||
eigval(:) = 0.0_dp
|
||||
eigvec(:, :) = 0.0_dp
|
||||
CALL jacobi(stress_tensor, eigval, eigvec)
|
||||
IF (numerical) THEN
|
||||
WRITE (UNIT=iw, FMT='(/,T2,A)') &
|
||||
'STRESS| Eigenvectors and eigenvalues of the numerical stress tensor [GPa]'
|
||||
WRITE (UNIT=iw, FMT="(/,T2,A)") &
|
||||
"STRESS| Eigenvectors and eigenvalues of the numerical stress tensor ["// &
|
||||
TRIM(ADJUSTL(unit_string))//"]"
|
||||
ELSE
|
||||
WRITE (UNIT=iw, FMT='(/,T2,A)') &
|
||||
'STRESS| Eigenvectors and eigenvalues of the analytical stress tensor [GPa]'
|
||||
WRITE (UNIT=iw, FMT="(/,T2,A)") &
|
||||
"STRESS| Eigenvectors and eigenvalues of the analytical stress tensor ["// &
|
||||
TRIM(ADJUSTL(unit_string))//"]"
|
||||
END IF
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T14,3(1X,I19))') &
|
||||
'STRESS|', 1, 2, 3
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,ES19.11))') &
|
||||
'STRESS| Eigenvalues', eigval(1:3)
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
|
||||
'STRESS| x', eigvec(1, 1:3)
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
|
||||
'STRESS| y', eigvec(2, 1:3)
|
||||
WRITE (UNIT=iw, FMT='(T2,A,T21,3(1X,F19.12))') &
|
||||
'STRESS| z', eigvec(3, 1:3)
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T14,3(1X,I19))") &
|
||||
"STRESS|", 1, 2, 3
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
|
||||
"STRESS| Eigenvalues", eigval(1:3)
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
|
||||
"STRESS| x", eigvec(1, 1:3)
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
|
||||
"STRESS| y", eigvec(2, 1:3)
|
||||
WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
|
||||
"STRESS| z", eigvec(3, 1:3)
|
||||
END IF
|
||||
|
||||
END SUBROUTINE write_stress_tensor
|
||||
|
||||
END MODULE virial_methods
|
||||
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@
|
|||
|
||||
&FORCE_EVAL
|
||||
METHOD Quickstep
|
||||
STRESS_TENSOR ANALYTICAL
|
||||
STRESS_TENSOR analytical
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME BASIS_SET
|
||||
POTENTIAL_FILE_NAME "GTH_SOC_POTENTIALS" # quotes needed to allow for this comment
|
||||
|
|
@ -62,9 +62,15 @@
|
|||
&END DFT
|
||||
&PRINT
|
||||
&FORCES
|
||||
FORCE_UNIT eV/Angstrom
|
||||
NDIGITS 6
|
||||
UNIT "eV/Angstrom"
|
||||
&END FORCES
|
||||
&PROGRAM_RUN_INFO
|
||||
ENERGY_UNIT eV
|
||||
&END PROGRAM_RUN_INFO
|
||||
&STRESS_TENSOR
|
||||
STRESS_UNIT kbar
|
||||
&END STRESS_TENSOR
|
||||
&END PRINT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue