RI HFX: improve defaults & input

cp2k · Apr 29, 2020 · 098a9ae · 098a9ae
1 parent c1d5dea
commit 098a9ae
Show file tree

Hide file tree

Showing 3 changed files with 31 additions and 15 deletions.
diff --git a/src/hfx_ri.F b/src/hfx_ri.F
@@ -1372,8 +1372,8 @@ SUBROUTINE hfx_ri_update_ks_Pmat(qs_env, ri_data, ks_matrix, rho_ao, &
       CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_update_ks_Pmat', &
          routineP = moduleN//':'//routineN
 
-      INTEGER :: col, handle, handle2, i_mem, iblk, iblk_filter, iblkrow, ispin, j_mem, nblk, &
-         nblk_filter, row, row_end, row_start, unit_nr, unit_nr_dbcsr
+      INTEGER :: col, handle, handle2, i_mem, iblk, iblk_filter, iblkrow, ispin, j_mem, n_mem, &
+         nblk, nblk_filter, row, row_end, row_start, unit_nr, unit_nr_dbcsr
       INTEGER(int_8)                                     :: flops_ks_max, flops_ri_max, nblks, &
                                                             nflop, nze, nze_3c, nze_3c_1, &
                                                             nze_3c_2, nze_ks, nze_rho
@@ -1452,6 +1452,8 @@ SUBROUTINE hfx_ri_update_ks_Pmat(qs_env, ri_data, ks_matrix, rho_ao, &
 
       CALL mp_sync(para_env%group)
       t1 = m_walltime()
+
+      n_mem = ri_data%n_mem
       DO ispin = 1, nspins
          CALL dbcsr_t_batched_contract_init(ks_t)
          CALL dbcsr_t_batched_contract_init(ri_data%t_2c_int(1, 1))
@@ -1467,7 +1469,7 @@ SUBROUTINE hfx_ri_update_ks_Pmat(qs_env, ri_data, ks_matrix, rho_ao, &
          occ_3c_2 = 0.0_dp
          nze_ks = 0
          occ_ks = 0.0_dp
-         DO i_mem = 1, ri_data%n_mem
+         DO i_mem = 1, n_mem
             CALL dbcsr_t_copy_matrix_to_tensor(rho_ao(ispin, 1)%matrix, rho_ao_tmp)
             CALL dbcsr_t_get_info(rho_ao_t, nfull_total=dims_2c)
             bounds_2c(:, 1) = [ri_data%starts_array_mem(i_mem), ri_data%ends_array_mem(i_mem)]
@@ -1482,7 +1484,7 @@ SUBROUTINE hfx_ri_update_ks_Pmat(qs_env, ri_data, ks_matrix, rho_ao, &
 
             CALL dbcsr_t_batched_contract_init(rho_ao_t)
 
-            DO j_mem = 1, ri_data%n_mem
+            DO j_mem = 1, n_mem
                CALL timeset(routineN//"_Px3C", handle2)
 
                CALL dbcsr_t_get_info(ri_data%t_3c_int_ctr_1(1, 1), nfull_total=dims_3c)
@@ -1608,7 +1610,7 @@ SUBROUTINE hfx_ri_update_ks_Pmat(qs_env, ri_data, ks_matrix, rho_ao, &
                                      dbcsr_scalar(1.0_dp), ks_t, &
                                      contract_1=[1, 2], notcontract_1=[3], &
                                      contract_2=[1, 2], notcontract_2=[3], &
-                                     map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps/ri_data%n_mem, &
+                                     map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps/n_mem, &
                                      bounds_1=bounds_ii, &
                                      bounds_3=bounds_j, &
                                      unit_nr=unit_nr_dbcsr, &

diff --git a/src/hfx_types.F b/src/hfx_types.F
@@ -1116,12 +1116,14 @@ SUBROUTINE hfx_ri_init_read_input(ri_data, ri_section, qs_kind_set, &
 
       INTEGER                                            :: handle
       LOGICAL                                            :: explicit
+      REAL(dp)                                           :: eps_storage_scaling
 
       CALL timeset(routineN, handle)
 
       CALL section_vals_val_get(ri_section, "EPS_FILTER", r_val=ri_data%filter_eps)
       CALL section_vals_val_get(ri_section, "EPS_FILTER_2C", r_val=ri_data%filter_eps_2c)
-      CALL section_vals_val_get(ri_section, "EPS_FILTER_3C", r_val=ri_data%filter_eps_3c)
+      CALL section_vals_val_get(ri_section, "EPS_STORAGE_SCALING", r_val=eps_storage_scaling)
+      ri_data%filter_eps_3c = ri_data%filter_eps*eps_storage_scaling
       CALL section_vals_val_get(ri_section, "EPS_FILTER_MO", r_val=ri_data%filter_eps_mo)
 
       ASSOCIATE (ri_metric=>ri_data%ri_metric, hfx_pot=>ri_data%hfx_pot)
@@ -1154,6 +1156,10 @@ SUBROUTINE hfx_ri_init_read_input(ri_data, ri_section, qs_kind_set, &
       CALL section_vals_val_get(ri_section, "MIN_BLOCK_SIZE_MO", i_val=ri_data%min_bsize_MO)
       CALL section_vals_val_get(ri_section, "MEMORY_CUT", i_val=ri_data%n_mem)
 
+      IF (ri_data%flavor == ri_pmat) THEN
+         ri_data%n_mem = FLOOR(SQRT(REAL(ri_data%n_mem) - 0.1)) + 1
+      ENDIF
+
       ri_data%orb_basis_type = orb_basis_type
       ri_data%ri_basis_type = ri_basis_type
       ri_data%nelectron_total = nelectron_total
@@ -2624,8 +2630,14 @@ SUBROUTINE hfx_print_ri_info(ri_data, hfx_section)
             "HFX_RI_INFO| Minimum block size", ri_data%min_bsize
          WRITE (UNIT=iw, FMT="(T3, A, T78, I3)") &
             "HFX_RI_INFO| MO block size", ri_data%min_bsize_MO
-         WRITE (UNIT=iw, FMT="(T3, A, T79, I2)") &
-            "HFX_RI_INFO| Memory reduction factor", ri_data%n_mem
+         SELECT CASE (ri_data%flavor)
+         CASE (ri_mo)
+            WRITE (UNIT=iw, FMT="(T3, A, T79, I2)") &
+               "HFX_RI_INFO| Memory reduction factor", ri_data%n_mem
+         CASE (ri_pmat)
+            WRITE (UNIT=iw, FMT="(T3, A, T79, I2)") &
+               "HFX_RI_INFO| Memory reduction factor", ri_data%n_mem**2
+         END SELECT
       ENDIF
 
    END SUBROUTINE

diff --git a/src/input_cp2k_hfx.F b/src/input_cp2k_hfx.F
@@ -525,10 +525,11 @@ SUBROUTINE create_hf_ri_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
-      CALL keyword_create(keyword, __LOCATION__, name="EPS_FILTER_3C", &
-                          description="Filter threshold for 3c integrals, limiting memory consumption "// &
-                          "of 3-center integral storage. Can be set looser than EPS_FILTER.", &
-                          default_r_val=1.0E-07_dp)
+      CALL keyword_create(keyword, __LOCATION__, &
+                          name="EPS_STORAGE_SCALING", &
+                          description="Scaling factor to scale EPS_FILTER for storage of 3-center integrals. Storage threshold "// &
+                          "will be EPS_FILTER*EPS_STORAGE_SCALING.", &
+                          default_r_val=10.0_dp)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
@@ -647,7 +648,7 @@ SUBROUTINE create_hf_ri_section(section)
 
       CALL keyword_create(keyword, __LOCATION__, name="MIN_BLOCK_SIZE", &
                           description="Minimum tensor block size.", &
-                          default_i_val=5)
+                          default_i_val=3)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
@@ -658,8 +659,9 @@ SUBROUTINE create_hf_ri_section(section)
       CALL keyword_release(keyword)
 
       CALL keyword_create(keyword, __LOCATION__, name="MEMORY_CUT", &
-                          description="Memory reduction factor.", &
-                          default_i_val=5)
+                          description="Memory reduction factor. For RI_FLAVOR RHO, this should be a square number "// &
+                          "(if not, the next greater square number is used).", &
+                          default_i_val=9)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)