RI HFX: optionally calculate condition number of (P|Q) matrices

src/hfx_ri.F

@@ -9,6 +9,7 @@
 
 MODULE hfx_ri
 
+   USE arnoldi_api,                     ONLY: arnoldi_extremal
    USE atomic_kind_types,               ONLY: atomic_kind_type
    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
                                               gto_basis_set_type
@@ -338,6 +339,8 @@ SUBROUTINE hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_RI, t_2c_int_po
          starts_array_mc_block_int, starts_array_mc_int
       INTEGER, DIMENSION(3)                              :: pcoord, pdims
       INTEGER, DIMENSION(:), POINTER                     :: col_bsize, row_bsize
+      LOGICAL                                            :: converged
+      REAL(dp)                                           :: max_ev, min_ev
       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
       TYPE(dbcsr_t_type)                                 :: t_3c_tmp
@@ -473,6 +476,46 @@ SUBROUTINE hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_RI, t_2c_int_po
          CALL init_interaction_radii_orb_basis(orb_basis, dft_control%qs_control%eps_pgf_orb)
       ENDDO
 
+      IF (ri_data%calc_condnum) THEN
+         CALL arnoldi_extremal(t_2c_int_pot(1), max_ev, min_ev, threshold=ri_data%eps_lanczos, &
+                               max_iter=ri_data%max_iter_lanczos, converged=converged)
+
+         IF (.NOT. converged) THEN
+            CPWARN("Condition number estimate of (P|Q) (HFX potential) is not reliable (not converged).")
+         ENDIF
+
+         IF (ri_data%unit_nr > 0) THEN
+            WRITE (ri_data%unit_nr, '(T2,A)') "2-Norm Condition Number of (P|Q) integrals (HFX potential)"
+            IF (min_ev > 0) THEN
+               WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
+                  "CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", LOG10(max_ev/min_ev)
+            ELSE
+               WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
+                  "CN : max/min ev: ", max_ev, " / ", min_ev, "Log(CN): infinity"
+            ENDIF
+         ENDIF
+
+         IF (.NOT. ri_data%same_op) THEN
+            CALL arnoldi_extremal(t_2c_int_RI(1), max_ev, min_ev, threshold=ri_data%eps_lanczos, &
+                                  max_iter=ri_data%max_iter_lanczos, converged=converged)
+
+            IF (.NOT. converged) THEN
+               CPWARN("Condition number estimate of (P|Q) matrix (RI metric) is not reliable (not converged).")
+            ENDIF
+
+            IF (ri_data%unit_nr > 0) THEN
+               WRITE (ri_data%unit_nr, '(T2,A)') "2-Norm Condition Number of (P|Q) integrals (RI metric)"
+               IF (min_ev > 0) THEN
+                  WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
+                     "CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", LOG10(max_ev/min_ev)
+               ELSE
+                  WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
+                     "CN : max/min ev: ", max_ev, " / ", min_ev, "Log(CN): infinity"
+               ENDIF
+            ENDIF
+         ENDIF
+      ENDIF
+
       CALL timestop(handle)
    END SUBROUTINE
 

src/hfx_types.F

@@ -361,7 +361,7 @@ MODULE hfx_types
                        eps_lanczos, eps_pgf_orb
       INTEGER :: t2c_sqrt_order, max_iter_lanczos, flavor, unit_nr_dbcsr, unit_nr, &
                  min_bsize, min_bsize_MO, t2c_method, nelectron_total
-      LOGICAL :: check_2c_inv
+      LOGICAL :: check_2c_inv, calc_condnum
 
       TYPE(libint_potential_type) :: ri_metric
 
@@ -1133,6 +1133,7 @@ SUBROUTINE hfx_ri_init_read_input(ri_data, ri_section, qs_kind_set, &
       CALL section_vals_val_get(ri_section, "2C_MATRIX_FUNCTIONS", i_val=ri_data%t2c_method)
 
       CALL section_vals_val_get(ri_section, "CHECK_2C_MATRIX", l_val=ri_data%check_2c_inv)
+      CALL section_vals_val_get(ri_section, "CALC_COND_NUM", l_val=ri_data%calc_condnum)
       CALL section_vals_val_get(ri_section, "SQRT_ORDER", i_val=ri_data%t2c_sqrt_order)
       CALL section_vals_val_get(ri_section, "EPS_LANCZOS", r_val=ri_data%eps_lanczos)
       CALL section_vals_val_get(ri_section, "MAX_ITER_LANCZOS", i_val=ri_data%max_iter_lanczos)

src/input_cp2k_hfx.F

@@ -591,6 +591,17 @@ SUBROUTINE create_hf_ri_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create( &
+         keyword, __LOCATION__, &
+         name="CALC_COND_NUM", &
+         variants=(/"CALC_CONDITION_NUMBER"/), &
+         description="Calculate the condition number of integral matrices.", &
+         usage="CALC_COND_NUM", &
+         default_l_val=.FALSE., &
+         lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="SQRT_ORDER", &
                           description="Order of the iteration method for the calculation of "// &
                           "the sqrt of 2-center integral matrix.", &
@@ -612,7 +623,7 @@ SUBROUTINE create_hf_ri_section(section)
 
       CALL keyword_create(keyword, __LOCATION__, name="MAX_ITER_LANCZOS", &
                           description="Maximum number of lanczos iterations.", &
-                          usage="MAX_ITER_LANCZOS ", default_i_val=128)
+                          usage="MAX_ITER_LANCZOS ", default_i_val=500)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 

tests/QS/regtest-hfx-ri/H2O-hfx-periodic-ri-truncated.inp

@@ -28,6 +28,7 @@
                     CUTOFF_RADIUS 1.5
                     EPS_FILTER 1.0E-06
                     EPS_PGF_ORB 1.0E-04
+                    CALC_COND_NUM T
                 &END
                 &INTERACTION_POTENTIAL
                     POTENTIAL_TYPE TRUNCATED

tests/QS/regtest-hfx-ri/Ne-hybrid-periodic-shortrange.inp

@@ -34,6 +34,7 @@
             &END XC_FUNCTIONAL
             &HF
                 &RI
+                    CALC_COND_NUM T
                 &END
                 &SCREENING
                     EPS_SCHWARZ 1.0E-5