Allow to specify the physical unit for energy and stress

This affects only the output quantities tagged with ENERGY| and STRESS|.

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 @@ SUBROUTINE ec_dc_build_ks_matrix_force(qs_env, ec_env)
 
       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 @@ SUBROUTINE ec_dc_build_ks_matrix_force(qs_env, ec_env)
             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 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
 
       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 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
       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 @@ SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
             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

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 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
       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 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
 
       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 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
          ! 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 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
             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 @@ SUBROUTINE force_env_calc_num_pressure(force_env, dx)
 
       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 @@ SUBROUTINE force_env_calc_num_pressure(force_env, dx)
       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 @@ SUBROUTINE force_env_calc_num_pressure(force_env, dx)
 
       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, &

src/force_env_utils.F

@@ -420,7 +420,7 @@ SUBROUTINE write_forces(particles, iw, label, ndigits, unit_string, total_force,
       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

src/input_cp2k_force_eval.F

@@ -296,17 +296,29 @@ SUBROUTINE create_f_env_print_section(section)
 
       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,49 +331,55 @@ SUBROUTINE create_f_env_print_section(section)
       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)
       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__, "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 @@ SUBROUTINE create_f_env_print_section(section)
                           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)