diff --git a/src/energy_corrections.F b/src/energy_corrections.F index d5cdb5d659..a625ef2feb 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -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 diff --git a/src/force_env_methods.F b/src/force_env_methods.F index d59f7c2c65..d391989c30 100644 --- a/src/force_env_methods.F +++ b/src/force_env_methods.F @@ -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, & diff --git a/src/force_env_utils.F b/src/force_env_utils.F index b501b71511..fa1f1152bd 100644 --- a/src/force_env_utils.F +++ b/src/force_env_utils.F @@ -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 diff --git a/src/input_cp2k_force_eval.F b/src/input_cp2k_force_eval.F index f89c529bfc..86f0a3c271 100644 --- a/src/input_cp2k_force_eval.F +++ b/src/input_cp2k_force_eval.F @@ -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) diff --git a/src/start/input_cp2k.F b/src/start/input_cp2k.F index 82f0cdd7d6..740a4cc61e 100644 --- a/src/start/input_cp2k.F +++ b/src/start/input_cp2k.F @@ -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) diff --git a/src/virial_methods.F b/src/virial_methods.F index 34ca8e48f9..e935008478 100644 --- a/src/virial_methods.F +++ b/src/virial_methods.F @@ -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 - diff --git a/tests/QS/regtest-gpw-1/H2O-geoopt.inp b/tests/QS/regtest-gpw-1/H2O-geoopt.inp index dfe35aaef6..0ec76260a5 100644 --- a/tests/QS/regtest-gpw-1/H2O-geoopt.inp +++ b/tests/QS/regtest-gpw-1/H2O-geoopt.inp @@ -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