RI HFX: use different filter_eps for 2c integral tensor

for better convergence of calculation of matrix inverse/sqrt
cp2k · Mar 7, 2020 · a11e60e · a11e60e
1 parent f7de199
commit a11e60e
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 5 deletions.
diff --git a/src/hfx_ri.F b/src/hfx_ri.F
@@ -328,7 +328,7 @@ SUBROUTINE hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_RI, t_2c_int_po
 
       CALL dbcsr_distribution_release(dbcsr_dist)
 
-      CALL build_2c_integrals(t_2c_int_pot, ri_data%filter_eps, qs_env, nl_2c_pot, basis_set_RI, basis_set_RI, &
+      CALL build_2c_integrals(t_2c_int_pot, ri_data%filter_eps_int, qs_env, nl_2c_pot, basis_set_RI, basis_set_RI, &
                               ri_data%hfx_pot)
 
       CALL release_neighbor_list_sets(nl_2c_pot)
@@ -340,8 +340,10 @@ SUBROUTINE hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_RI, t_2c_int_po
                                       dist_2d=dist_2d)
 
          CALL dbcsr_create(t_2c_int_RI(1), template=t_2c_int_pot(1), matrix_type=dbcsr_type_symmetric)
-         CALL build_2c_integrals(t_2c_int_RI, ri_data%filter_eps, qs_env, nl_2c_RI, basis_set_RI, basis_set_RI, &
+
+         CALL build_2c_integrals(t_2c_int_RI, ri_data%filter_eps_int, qs_env, nl_2c_RI, basis_set_RI, basis_set_RI, &
                                  ri_data%ri_metric)
+
          CALL release_neighbor_list_sets(nl_2c_RI)
       ENDIF
 

diff --git a/src/hfx_types.F b/src/hfx_types.F
@@ -361,7 +361,7 @@ MODULE hfx_types
 
    TYPE hfx_ri_type
       ! input parameters
-      REAL(KIND=dp) :: filter_eps, eps_schwarz, &
+      REAL(KIND=dp) :: filter_eps, filter_eps_int, eps_schwarz, &
                        eps_lanczos, eps_pgf_orb
       INTEGER :: t2c_sqrt_order, max_iter_lanczos, flavor, unit_nr, max_bsize
       LOGICAL :: check_2c_inv
@@ -986,6 +986,7 @@ SUBROUTINE hfx_ri_init_read_input(ri_data, x_data, hfx_section, ri_section, qs_k
       ri_data%eps_schwarz = x_data%screening_parameter%eps_schwarz
 
       CALL section_vals_val_get(ri_section, "EPS_FILTER", r_val=ri_data%filter_eps)
+      CALL section_vals_val_get(ri_section, "EPS_FILTER_INT", r_val=ri_data%filter_eps_int)
 
       ASSOCIATE (ri_metric=>ri_data%ri_metric, hfx_pot=>ri_data%hfx_pot)
          CALL section_vals_val_get(ri_section, "RI_METRIC", i_val=ri_metric%potential_type, explicit=explicit)
@@ -1105,7 +1106,7 @@ SUBROUTINE hfx_ri_init(ri_data, qs_kind_set, particle_set, para_env, basis_cntrl
          IF (hfx_pot%potential_type == do_potential_short) THEN
             ! need a more accurate threshold for determining 2-center integral operator range
             ! because stability of matrix inversion/sqrt is sensitive to this
-            CALL erfc_cutoff(ri_data%eps_schwarz, hfx_pot%omega, hfx_pot%cutoff_radius)
+            CALL erfc_cutoff(ri_data%filter_eps_int, hfx_pot%omega, hfx_pot%cutoff_radius)
          ENDIF
          ! determine whether RI metric is same operator as used in HFX
          same_op = ri_metric%potential_type == hfx_pot%potential_type

diff --git a/src/input_cp2k_hfx.F b/src/input_cp2k_hfx.F
@@ -512,11 +512,17 @@ SUBROUTINE create_hf_ri_section(section)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="EPS_FILTER", &
-                          description="Filter threshold for DBCSR tensors.", &
+                          description="Filter threshold for DBCSR tensor contraction.", &
                           default_r_val=1.0E-09_dp)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="EPS_FILTER_INT", &
+                          description="Filter threshold for 2c integrals. Should be tighter than EPS_FILTER.", &
+                          default_r_val=1.0E-10_dp)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="OMEGA", &
                           description="The range parameter for the short range operator (in 1/a0). "// &
                           "Default is OMEGA from INTERACTION_POTENTIAL. ", &