4-center Hartree-Fock and ADMM for exchange self-energy in GW+kpoints

src/input_cp2k_mp2.F

@@ -1224,7 +1224,17 @@ SUBROUTINE create_low_scaling(section)
          "means smaller expansion coefficients that leads to a more stable calculation at the price "// &
          "of a slightly worse RI approximation. In case the parameter 0.0 is chosen, ordinary RI is used.", &
          usage="REGULARIZATION_RI 1.0E-4", &
-         default_r_val=1.0E-2_dp)
+         default_r_val=0.0_dp)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create( &
+         keyword, __LOCATION__, &
+         name="EPS_EIGVAL_S", &
+         description="Parameter to reduce the expansion coefficients in RI for periodic GW. Removes all "// &
+         "eigenvectors and eigenvalues of S_PQ(k) that are smaller than EPS_EIGVAL_S. ", &
+         usage="EPS_EIGVAL_S 1.0E-3", &
+         default_r_val=1.0E-5_dp)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
@@ -1251,6 +1261,16 @@ SUBROUTINE create_low_scaling(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create( &
+         keyword, __LOCATION__, &
+         name="REL_CUTOFF_TRUNC_COLOUMB_RI_X", &
+         description="Only active in case TRUNC_COLOUMB_RI_X = True. Normally, relative cutoff = 0.5 is "// &
+         "good choice; still needs to be evaluated for RI schemes. ", &
+         usage="REL_CUTOFF_TRUNC_COLOUMB_RI_X 0.3", &
+         default_r_val=0.5_dp)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create( &
          keyword, __LOCATION__, &
          name="KEEP_QUADRATURE", &

src/mp2.F

@@ -697,7 +697,8 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
       DEALLOCATE (mos_mp2)
 
       ! if necessary reallocate hfx buffer
-      IF (free_hfx_buffer .AND. (.NOT. calc_forces)) THEN
+      IF (free_hfx_buffer .AND. (.NOT. calc_forces) .AND. &
+          (mp2_env%ri_g0w0%do_ri_Sigma_x .OR. .NOT. mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma)) THEN
          CALL timeset(routineN//"_alloc_hfx", handle2)
          DO irep = 1, n_rep_hf
             DO i_thread = 0, n_threads - 1

src/mp2_integrals.F

@@ -313,7 +313,7 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
          CPABORT("SVD not implemented for forces.!")
       END IF
 
-      do_kpoints_from_Gamma = SUM(qs_env%mp2_env%ri_rpa_im_time%kp_grid) > 0
+      do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
          CALL get_qs_env(qs_env=qs_env, &
                          kpoints=kpoints)

src/mp2_ri_2c.F

@@ -260,7 +260,8 @@ SUBROUTINE get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_su
                                     do_kpoints, kpoints, put_mat_KS_env=.FALSE.)
 
             CALL inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_matrix_L_kpoints, &
-                                                            fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints)
+                                                            fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints, &
+                                                            qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S)
 
          ELSE
             IF (calc_forces .AND. (.NOT. do_im_time)) THEN
@@ -1404,15 +1405,18 @@ END SUBROUTINE cholesky_decomp
 !> \param fm_matrix_Sinv_Vtrunc_Sinv ...
 !> \param dimen_RI ...
 !> \param kpoints ...
+!> \param eps_eigval_S ...
 ! **************************************************************************************************
    SUBROUTINE inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_matrix_L_kpoints, &
-                                                         fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints)
+                                                         fm_matrix_Sinv_Vtrunc_Sinv, dimen_RI, kpoints, &
+                                                         eps_eigval_S)
 
       TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER       :: fm_matrix_L_RI_metric, &
                                                             fm_matrix_L_kpoints, &
                                                             fm_matrix_Sinv_Vtrunc_Sinv
       INTEGER, INTENT(IN)                                :: dimen_RI
       TYPE(kpoint_type), POINTER                         :: kpoints
+      REAL(KIND=dp), INTENT(IN)                          :: eps_eigval_S
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'inversion_of_S_and_mult_with_chol_dec_of_V'
       COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
@@ -1449,7 +1453,7 @@ SUBROUTINE inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_
          CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_Vtrunc_tmp, cone, fm_matrix_Sinv_Vtrunc_Sinv(ikp, 1)%matrix)
          CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_Vtrunc_tmp, ione, fm_matrix_Sinv_Vtrunc_Sinv(ikp, 2)%matrix)
 
-         CALL cp_cfm_robust_cholesky(cfm_matrix_S_tmp, cfm_matrix_S_inv_tmp, threshold=1.0E-12_dp, exponent=-0.5_dp)
+         CALL cp_cfm_robust_cholesky(cfm_matrix_S_tmp, cfm_matrix_S_inv_tmp, threshold=eps_eigval_S, exponent=-0.5_dp)
 
          ! get S^(-1) = U^(-H)U^(-1)
          CALL cp_cfm_gemm("C", "N", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_S_tmp, cfm_matrix_S_tmp, &
@@ -1501,7 +1505,8 @@ SUBROUTINE setup_abs_cutoffs_chi_and_trunc_coulomb_potential(qs_env, trunc_coulo
 
       INTEGER                                            :: handle
       INTEGER, DIMENSION(3)                              :: periodic
-      REAL(KIND=dp)                                      :: shortest_dist_cell_planes
+      REAL(KIND=dp)                                      :: rel_cutoff_trunc_coulomb_ri_x, &
+                                                            shortest_dist_cell_planes
       REAL(KIND=dp), DIMENSION(2)                        :: abs_cutoffs_chi_W
       TYPE(cell_type), POINTER                           :: cell
 
@@ -1534,9 +1539,12 @@ SUBROUTINE setup_abs_cutoffs_chi_and_trunc_coulomb_potential(qs_env, trunc_coulo
 
       qs_env%mp2_env%ri_rpa_im_time%abs_cutoffs_chi_W(1:2) = abs_cutoffs_chi_W(1:2)
 
+      rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x
+
       IF (PRESENT(trunc_coulomb)) THEN
          trunc_coulomb%potential_type = do_potential_truncated
-         trunc_coulomb%cutoff_radius = (abs_cutoffs_chi_W(1) + abs_cutoffs_chi_W(2))*0.5_dp
+         trunc_coulomb%cutoff_radius = (abs_cutoffs_chi_W(1) + abs_cutoffs_chi_W(2))* &
+                                       rel_cutoff_trunc_coulomb_ri_x
          trunc_coulomb%filename = "t_c_g.dat"
          ! dummy
          trunc_coulomb%omega = 0.0_dp

src/mp2_setup.F

@@ -238,6 +238,7 @@ SUBROUTINE read_mp2_section(input, mp2_env)
                                 l_val=mp2_env%ri_rpa_im_time%do_im_time_kpoints)
       CALL section_vals_val_get(low_scaling_section, "KPOINTS", &
                                 i_vals=mp2_env%ri_rpa_im_time%kp_grid)
+      mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma = SUM(mp2_env%ri_rpa_im_time%kp_grid) > 0
       CALL section_vals_val_get(low_scaling_section, "KPOINT_WEIGHTS_W", &
                                 i_val=mp2_env%ri_rpa_im_time%kpoint_weights_W_method)
       CALL section_vals_val_get(low_scaling_section, "EXPONENT_TAILORED_WEIGHTS", &
@@ -246,10 +247,14 @@ SUBROUTINE read_mp2_section(input, mp2_env)
                                 r_vals=mp2_env%ri_rpa_im_time%rel_cutoffs_chi_W)
       CALL section_vals_val_get(low_scaling_section, "REGULARIZATION_RI", &
                                 r_val=mp2_env%ri_rpa_im_time%regularization_RI)
+      CALL section_vals_val_get(low_scaling_section, "EPS_EIGVAL_S", &
+                                r_val=mp2_env%ri_rpa_im_time%eps_eigval_S)
       CALL section_vals_val_get(low_scaling_section, "MAKE_CHI_POS_DEFINITE", &
                                 l_val=mp2_env%ri_rpa_im_time%make_chi_pos_definite)
       CALL section_vals_val_get(low_scaling_section, "TRUNC_COLOUMB_RI_X", &
                                 l_val=mp2_env%ri_rpa_im_time%trunc_coulomb_ri_x)
+      CALL section_vals_val_get(low_scaling_section, "REL_CUTOFF_TRUNC_COLOUMB_RI_X", &
+                                r_val=mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x)
       CALL section_vals_val_get(low_scaling_section, "K_MESH_G_FACTOR", &
                                 i_val=mp2_env%ri_rpa_im_time%k_mesh_g_factor)
       CALL section_vals_val_get(low_scaling_section, "GREENS_FUNCTION", &

src/mp2_types.F

@@ -153,9 +153,11 @@ MODULE mp2_types
 
    TYPE ri_rpa_im_time_type
       INTEGER                  :: cut_memory
-      LOGICAL                  :: memory_info, make_chi_pos_definite, trunc_coulomb_ri_x, keep_quad
+      LOGICAL                  :: memory_info, make_chi_pos_definite, trunc_coulomb_ri_x, keep_quad, &
+                                  do_kpoints_from_Gamma
       REAL(KIND=dp)            :: eps_filter, &
-                                  eps_filter_factor, eps_compress, exp_tailored_weights, regularization_RI
+                                  eps_filter_factor, eps_compress, exp_tailored_weights, regularization_RI, &
+                                  eps_eigval_S, rel_cutoff_trunc_coulomb_ri_x
       REAL(KIND=dp), DIMENSION(:), POINTER :: rel_cutoffs_chi_W, tau_tj, tau_wj, tj, wj
       REAL(KIND=dp), DIMENSION(:, :), POINTER :: weights_cos_tf_t_to_w, weights_cos_tf_w_to_t
       REAL(KIND=dp), DIMENSION(2) :: abs_cutoffs_chi_W
@@ -211,7 +213,8 @@ MODULE mp2_types
       INTEGER :: n_kp_in_kp_line, n_special_kp, nkp_self_energy, &
                  nkp_self_energy_special_kp, nkp_self_energy_monkh_pack
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp_special_kp
-      TYPE(dbcsr_p_type), DIMENSION(:), ALLOCATABLE :: mat_exchange_for_kp_from_gamma
+      TYPE(dbcsr_p_type), DIMENSION(:), ALLOCATABLE :: &
+         matrix_sigma_x_minus_vxc, matrix_ks
    END TYPE
 
    TYPE ri_basis_opt

src/qs_band_structure.F

@@ -63,7 +63,7 @@ MODULE qs_band_structure
 
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_band_structure'
 
-   PUBLIC :: calculate_band_structure, calculate_kp_orbitals
+   PUBLIC :: calculate_band_structure, calculate_kp_orbitals, calculate_kpoints_for_bs
 
 ! **************************************************************************************************
 
@@ -310,7 +310,6 @@ SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgenera
          OPTIONAL                                        :: kpgeneral
       INTEGER, INTENT(IN), OPTIONAL                      :: group_size_ext
 
-      INTEGER                                            :: i, ix, iy, iz, npoints
       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
       TYPE(cp_para_env_type), POINTER                    :: para_env
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks, matrix_s
@@ -320,6 +319,46 @@ SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgenera
       TYPE(qs_scf_env_type), POINTER                     :: scf_env
       TYPE(scf_control_type), POINTER                    :: scf_control
 
+      CALL calculate_kpoints_for_bs(kpoint, scheme, group_size_ext, mp_grid, kpgeneral)
+
+      CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
+      kpoint%para_env => para_env
+      CALL cp_para_env_retain(para_env)
+      kpoint%blacs_env_all => blacs_env
+      CALL cp_blacs_env_retain(blacs_env)
+      CALL kpoint_env_initialize(kpoint)
+
+      CALL kpoint_initialize_mos(kpoint, qs_env%mos, nadd)
+      CALL kpoint_initialize_mo_set(kpoint)
+
+      CALL get_qs_env(qs_env, sab_kp=sab_nl, dft_control=dft_control)
+      CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control)
+
+      CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
+                      scf_env=scf_env, scf_control=scf_control)
+      CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE.)
+
+   END SUBROUTINE calculate_kp_orbitals
+
+! **************************************************************************************************
+!> \brief ...
+!> \param kpoint ...
+!> \param scheme ...
+!> \param group_size_ext ...
+!> \param mp_grid ...
+!> \param kpgeneral ...
+! **************************************************************************************************
+   SUBROUTINE calculate_kpoints_for_bs(kpoint, scheme, group_size_ext, mp_grid, kpgeneral)
+
+      TYPE(kpoint_type), POINTER                         :: kpoint
+      CHARACTER(LEN=*), INTENT(IN)                       :: scheme
+      INTEGER, INTENT(IN), OPTIONAL                      :: group_size_ext
+      INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL        :: mp_grid
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
+         OPTIONAL                                        :: kpgeneral
+
+      INTEGER                                            :: i, ix, iy, iz, npoints
+
       CPASSERT(.NOT. ASSOCIATED(kpoint))
 
       CALL kpoint_create(kpoint)
@@ -389,25 +428,6 @@ SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgenera
          CPABORT("Unknown kpoint scheme requested")
       END SELECT
 
-      CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
-      kpoint%para_env => para_env
-      CALL cp_para_env_retain(para_env)
-      kpoint%blacs_env_all => blacs_env
-      CALL cp_blacs_env_retain(blacs_env)
-      CALL kpoint_env_initialize(kpoint)
-
-      CALL kpoint_initialize_mos(kpoint, qs_env%mos, nadd)
-      CALL kpoint_initialize_mo_set(kpoint)
-
-      CALL get_qs_env(qs_env, sab_kp=sab_nl, dft_control=dft_control)
-      CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control)
-
-      CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, &
-                      scf_env=scf_env, scf_control=scf_control)
-      CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE.)
-
-   END SUBROUTINE calculate_kp_orbitals
-
-! **************************************************************************************************
+   END SUBROUTINE calculate_kpoints_for_bs
 
 END MODULE qs_band_structure