Add output when SIRIUS is not converging

src/input_cp2k_pwdft.F

@@ -72,6 +72,7 @@ MODULE input_cp2k_pwdft
    SUBROUTINE create_pwdft_section(section)
       TYPE(section_type), POINTER                        :: section
 
+      TYPE(keyword_type), POINTER                        :: keyword
       TYPE(section_type), POINTER                        :: subsection
 
 !     ------------------------------------------------------------------------
@@ -86,6 +87,18 @@ SUBROUTINE create_pwdft_section(section)
                           "non collinear magnetism, Hubbard correction, all exchange functionals "// &
                           "supported by libxc and Van der Waals corrections (libvdwxc).")
 
+      NULLIFY (keyword)
+      CALL keyword_create(keyword, __LOCATION__, &
+                          name='ignore_convergence_failure', &
+                          description="when set to true, calculation will continue irrespectively "// &
+                          "of the convergence status of SIRIUS", &
+                          type_of_var=logical_t, &
+                          repeats=.FALSE., &
+                          default_l_val=.FALSE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       NULLIFY (subsection)
       CALL create_sirius_section(subsection, 'control')
       CALL section_add_subsection(section, subsection)

src/pwdft_environment_types.F

@@ -73,6 +73,7 @@ MODULE pwdft_environment_types
       TYPE(pwdft_energy_type), POINTER                 :: energy => NULL()
       REAL(KIND=dp), DIMENSION(:, :), POINTER          :: forces => NULL()
       REAL(KIND=dp), DIMENSION(3, 3)                   :: stress = 0.0_dp
+      LOGICAL                                          :: ignore_convergence_failure
 !     16 different functionals should be enough
       CHARACTER(len=80), DIMENSION(16)                 :: xc_func = ""
 #if defined(__SIRIUS)

src/sirius_interface.F

@@ -131,7 +131,7 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
       CHARACTER(len=default_string_length)               :: label
       INTEGER                                            :: i, iatom, ibeta, ifun, ikind, iwf, j, l, &
                                                             n, natom, nbeta, nkind, nmesh, &
-                                                            num_mag_dims, sirius_mpi_comm, vdw_func, nu, lu
+                                                            num_mag_dims, sirius_mpi_comm, vdw_func, nu, lu, output_unit
       INTEGER, DIMENSION(:), POINTER                     :: mpi_grid_dims
       INTEGER(KIND=C_INT), DIMENSION(3)                  :: k_grid, k_shift
       INTEGER, DIMENSION(:), POINTER                     :: kk
@@ -166,6 +166,7 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
 
       CPASSERT(ASSOCIATED(pwdft_env))
 
+      output_unit = cp_logger_get_default_io_unit()
       ! create context of simulation
       CALL pwdft_env_get(pwdft_env, para_env=para_env)
       sirius_mpi_comm = para_env%get_handle()
@@ -175,6 +176,8 @@ SUBROUTINE cp_sirius_create_env(pwdft_env)
 
       CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_section, xc_input=xc_section)
 
+      CALL section_vals_val_get(pwdft_section, "ignore_convergence_failure", &
+                                l_val=pwdft_env%ignore_convergence_failure)
       ! cp2k should *have* a function that return all xc_functionals. Doing
       ! manually is prone to errors
 
@@ -675,21 +678,38 @@ SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_
          POINTER                                         :: pwdft_env
       LOGICAL, INTENT(IN)                                :: calculate_forces, calculate_stress_tensor
 
-      INTEGER                                            :: n1, n2
-      LOGICAL                                            :: do_print
+      INTEGER                                            :: iw, n1, n2
+      LOGICAL                                            :: do_print, gs_converged
       REAL(KIND=C_DOUBLE)                                :: etotal
       REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :)  :: cforces
       REAL(KIND=C_DOUBLE), DIMENSION(3, 3)               :: cstress
       REAL(KIND=dp), DIMENSION(3, 3)                     :: stress
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: forces
+      TYPE(cp_logger_type), POINTER                      :: logger
       TYPE(pwdft_energy_type), POINTER                   :: energy
       TYPE(section_vals_type), POINTER                   :: print_section, pwdft_input
       TYPE(sirius_ground_state_handler)                  :: gs_handler
 
       CPASSERT(ASSOCIATED(pwdft_env))
 
+      NULLIFY (logger)
+      logger => cp_get_default_logger()
+      iw = cp_logger_get_default_io_unit(logger)
+
       CALL pwdft_env_get(pwdft_env=pwdft_env, gs_handler=gs_handler)
-      CALL sirius_find_ground_state(gs_handler)
+      CALL sirius_find_ground_state(gs_handler, converged=gs_converged)
+
+      IF (gs_converged) THEN
+         IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS: ground state is converged"
+      ELSE
+         IF (pwdft_env%ignore_convergence_failure) THEN
+            IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS Warning: ground state is not converged"
+         ELSE
+            CPABORT("CP2K/SIRIUS (ground state): SIRIUS did not converge.")
+         END IF
+      END IF
+      IF (iw > 0) CALL m_flush(iw)
+
       CALL pwdft_env_get(pwdft_env=pwdft_env, energy=energy)
       etotal = 0.0_C_DOUBLE
 
@@ -736,7 +756,6 @@ SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_
       IF (do_print) THEN
          CALL cp_sirius_print_results(pwdft_env, print_section)
       END IF
-
    END SUBROUTINE cp_sirius_energy_force
 
 !***************************************************************************************************