RPA based SIGMA Functionals (cp2k#3695)

Co-authored-by: Raviraj Chandreshbhai Mandalia <[email protected]>

src/CMakeLists.txt

@@ -739,6 +739,7 @@ list(
   rpa_gw_sigma_x.F
   rpa_im_time.F
   rpa_main.F
+  rpa_sigma_functional.F
   rpa_rse.F
   rpa_im_time_force_types.F
   rpa_im_time_force_methods.F

src/input_constants.F

@@ -1319,5 +1319,11 @@ MODULE input_constants
    INTEGER, PARAMETER, PUBLIC               :: bse_iter_en_cond = 0, &
                                                bse_iter_res_cond = 1, &
                                                bse_iter_both_cond = 2
+   ! SIGMA_functional parameters selection
+   INTEGER, PARAMETER, PUBLIC               :: sigma_none = 0, &
+                                               sigma_PBE_S2 = 1, &
+                                               sigma_PBE0_S1 = 2, &
+                                               sigma_PBE0_S2 = 3, &
+                                               sigma_PBE_S1 = 4
 
 END MODULE input_constants

src/input_cp2k_mp2.F

@@ -37,7 +37,8 @@ MODULE input_cp2k_mp2
       ot_precond_s_inverse, ri_default, ri_rpa_g0w0_crossing_bisection, &
       ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, soc_lda, soc_none, soc_pbe, &
       wfc_mm_style_gemm, wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, &
-      z_solver_sd, rpa_exchange_none, rpa_exchange_axk, rpa_exchange_sosex, G0W0, evGW0, evGW
+      z_solver_sd, rpa_exchange_none, rpa_exchange_axk, rpa_exchange_sosex, G0W0, evGW0, evGW, &
+      sigma_none, sigma_PBE0_S1, sigma_PBE0_S2, sigma_PBE_S1, sigma_PBE_S2
    USE input_cp2k_hfx, ONLY: create_hfx_section
    USE input_cp2k_kpoints, ONLY: create_kpoint_set_section
    USE input_keyword_types, ONLY: keyword_create, &
@@ -519,6 +520,23 @@ SUBROUTINE create_ri_rpa(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create( &
+         keyword, __LOCATION__, &
+         name="SIGMA_FUNCTIONAL", &
+         description="Determine parametrization for sigma-functional", &
+         usage="SIGMA_FUNCTIONAL PBE_S2", &
+         enum_c_vals=s2a("NONE", "PBE0_S1", "PBE0_S2", "PBE_S1", "PBE_S2"), &
+         enum_i_vals=(/sigma_none, sigma_PBE0_S1, sigma_PBE0_S2, sigma_PBE_S1, sigma_PBE_S2/), &
+         enum_desc=s2a("No sigma functional calculation", &
+                       "use parameters based on PBE0 with S1 set.", &
+                       "use parameters based on PBE0 with S2 set.", &
+                       "use parameters based on PBE with S1 set.", &
+                       "use parameters based on PBE with S2 set." &
+                       ), &
+         default_i_val=sigma_none)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="QUADRATURE_POINTS", &
                           variants=(/"RPA_NUM_QUAD_POINTS"/), &
                           description="Number of quadrature points for the numerical integration in the RI-RPA method.", &

src/mp2.F

@@ -55,7 +55,8 @@ MODULE mp2
                                               do_eri_mme,&
                                               rpa_exchange_axk,&
                                               rpa_exchange_none,&
-                                              rpa_exchange_sosex
+                                              rpa_exchange_sosex,&
+                                              sigma_none
    USE input_section_types,             ONLY: section_vals_get,&
                                               section_vals_get_subs_vals,&
                                               section_vals_type
@@ -606,7 +607,6 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
          ! go with ri_rpa_gpw
          CALL mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T, &
                            mos_mp2, para_env, unit_nr, calc_forces, calc_ex, do_ri_rpa=.TRUE.)
-
          ! Scale energy contributions
          Emp2 = Emp2*mp2_env%ri_rpa%scale_rpa
          mp2_env%ri_rpa%ener_exchange = mp2_env%ri_rpa%ener_exchange*mp2_env%ri_rpa%scale_rpa
@@ -665,6 +665,10 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
          IF (.NOT. update_xc_energy) Emp2 = 0.0_dp
 
          IF (unit_nr > 0 .AND. update_xc_energy) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA energy  =   ', Emp2
+         IF (unit_nr > 0 .AND. mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
+            WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Sigma corr. to RI-RPA energy  =   ', &
+               mp2_env%ri_rpa%e_sigma_corr
+         END IF
          IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_axk) THEN
             IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA-AXK energy=', mp2_env%ri_rpa%ener_exchange
          ELSE IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_sosex) THEN
@@ -687,6 +691,10 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
       IF (mp2_env%ri_rpa%do_rse) THEN
          Emp2 = Emp2 + mp2_env%ri_rpa%rse_corr
       END IF
+      IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
+         !WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Sigma corr. to RI-RPA energy  =   ',&
+         Emp2 = Emp2 + mp2_env%ri_rpa%e_sigma_corr
+      END IF
       energy%mp2 = Emp2
       energy%total = energy%total + Emp2
 

src/mp2_gpw.F

@@ -265,7 +265,6 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
       CALL get_cell(cell=cell, periodic=periodic)
 
       IF (do_im_time) THEN
-
          IF (mp2_env%ri_metric%potential_type == ri_default) THEN
             IF (SUM(periodic) == 1 .OR. SUM(periodic) == 3) THEN
                mp2_env%ri_metric%potential_type = do_potential_id

src/mp2_setup.F

@@ -100,6 +100,10 @@ SUBROUTINE read_mp2_section(input, mp2_env)
       END IF
 
       CALL section_vals_val_get(mp2_section, "RI_RPA%_SECTION_PARAMETERS_", l_val=do_rpa)
+
+      !CALL section_vals_val_get(mp2_section, "RI_RPA%SIGMA_FUNCTIONAL",l_val=mp2_env%ri_rpa%do_sigma)
+      CALL section_vals_val_get(mp2_section, "RI_RPA%SIGMA_FUNCTIONAL", i_val=mp2_env%ri_rpa%sigma_param)
+
       IF (do_rpa) THEN
          CALL check_method(mp2_env%method)
          mp2_env%method = ri_rpa_method_gpw
@@ -120,6 +124,7 @@ SUBROUTINE read_mp2_section(input, mp2_env)
                                 l_val=mp2_env%ri_rpa%use_hfx_implementation)
 
       CALL section_vals_val_get(mp2_section, "RI_RPA%RSE", l_val=mp2_env%ri_rpa%do_rse)
+
       CALL section_vals_val_get(mp2_section, "RI_RPA%PRINT_DGEMM_INFO", l_val=mp2_env%ri_rpa%print_dgemm_info)
 
       NULLIFY (gw_section)

src/mp2_types.F

@@ -130,6 +130,7 @@ MODULE mp2_types
                                                             do_admm = .FALSE., &
                                                             do_rse = .FALSE., &
                                                             print_dgemm_info = .FALSE.
+
       ! GCC 8 has an issue with this being an ALLOCATABLE
       TYPE(dbcsr_type), DIMENSION(:), POINTER            :: mo_coeff_o => NULL(), &
                                                             mo_coeff_v => NULL()
@@ -140,6 +141,11 @@ MODULE mp2_types
                                                             rse_corr_diag = 0.0_dp, &
                                                             rse_corr = 0.0_dp, &
                                                             scale_rpa = 0.0_dp
+
+      !LOGICAL                                            :: do_sigma = .FALSE.
+      INTEGER                                            :: sigma_param = 0.0_dp
+      REAL(KIND=dp)                                      :: e_sigma_corr = 0.0_dp
+
    END TYPE ri_rpa_type
 
    TYPE ri_rpa_im_time_type

src/rpa_main.F

@@ -51,6 +51,7 @@ MODULE rpa_main
    USE input_constants,                 ONLY: rpa_exchange_axk,&
                                               rpa_exchange_none,&
                                               rpa_exchange_sosex,&
+                                              sigma_none,&
                                               wfc_mm_style_gemm
    USE kinds,                           ONLY: dp,&
                                               int_8
@@ -101,6 +102,10 @@ MODULE rpa_main
                                               keep_initial_quad
    USE rpa_im_time_force_types,         ONLY: im_time_force_release,&
                                               im_time_force_type
+   USE rpa_sigma_functional,            ONLY: finalize_rpa_sigma,&
+                                              rpa_sigma_create,&
+                                              rpa_sigma_matrix_spectral,&
+                                              rpa_sigma_type
    USE rpa_util,                        ONLY: Q_trace_and_add_unit_matrix,&
                                               alloc_im_time,&
                                               calc_mat_Q,&
@@ -323,7 +328,6 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i
       ELSE
          num_integ_points = mp2_env%ri_rpa%rpa_num_quad_points
          input_num_integ_groups = mp2_env%ri_rpa%rpa_num_integ_groups
-
          IF (my_do_gw .AND. do_minimax_quad) THEN
             IF (num_integ_points > 34) THEN
                IF (unit_nr > 0) &
@@ -478,6 +482,7 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i
       ! not necessary for imaginary time
 
       ALLOCATE (BIb_C_2D(nspins))
+
       IF (.NOT. do_im_time) THEN
 
          ! reorder the local data in such a way to help the next stage of matrix creation
@@ -1153,12 +1158,14 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
       TYPE(im_time_force_type)                           :: force_data
       TYPE(rpa_exchange_work_type)                       :: exchange_work
       TYPE(rpa_grad_type)                                :: rpa_grad
+      TYPE(rpa_sigma_type)                               :: rpa_sigma
       TYPE(two_dim_int_array), ALLOCATABLE, DIMENSION(:) :: t_3c_O_mo_ind
 
       CALL timeset(routineN, handle)
 
       nspins = SIZE(homo)
       nmo = homo(1) + virtual(1)
+
       my_open_shell = (nspins == 2)
 
       do_gw_im_time = my_do_gw .AND. do_im_time
@@ -1255,7 +1262,6 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
          END IF
 
       END IF
-
       IF (do_ic_model) THEN
          ! image charge model only implemented for cubic scaling GW
          CPASSERT(do_gw_im_time)
@@ -1272,7 +1278,6 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
                                do_kpoints_cubic_RPA, e_fermi(1), tj, wj, &
                                weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
                                qs_env%mp2_env%ri_g0w0%regularization_minimax)
-
          !For sos_laplace_mp2 and low-scaling RPA, potentially need to store/retrieve the initial weights
          IF (qs_env%mp2_env%ri_rpa_im_time%keep_quad) THEN
             CALL keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, &
@@ -1298,7 +1303,6 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
       ELSE
          alpha = 4.0_dp
       END IF
-
       IF (my_do_gw) THEN
          CALL allocate_matrices_gw(vec_Sigma_c_gw, color_rpa_group, dimen_nm_gw, &
                                    gw_corr_lev_occ, gw_corr_lev_virt, homo, &
@@ -1325,19 +1329,17 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
       IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_create(rpa_grad, fm_mat_Q(1), &
                                                                    fm_mat_S, homo, virtual, mp2_env, Eigenval(:, 1, :), &
                                                                    unit_nr, do_ri_sos_laplace_mp2)
-
+!doubtful
       IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) &
          CALL exchange_work%create(qs_env, para_env_sub, mat_munu, dimen_RI_red, &
                                    fm_mat_S, fm_mat_Q(1), fm_mat_Q_gemm(1), homo, virtual)
-
       Erpa = 0.0_dp
       IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) e_exchange = 0.0_dp
       first_cycle = .TRUE.
       omega_old = 0.0_dp
       CALL dgemm_counter_init(dgemm_counter, unit_nr, mp2_env%ri_rpa%print_dgemm_info)
 
       DO count_ev_sc_GW = 1, iter_evGW
-
          dbcsr_time = 0.0_dp
          dbcsr_nflop = 0
 
@@ -1351,7 +1353,7 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
          END IF
 
          ! calculate Q_PQ(it)
-         IF (do_im_time) THEN
+         IF (do_im_time) THEN ! not using Imaginary time
 
             IF (.NOT. do_kpoints_cubic_RPA) THEN
                DO ispin = 1, nspins
@@ -1424,8 +1426,11 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
 
          END IF ! do im time
 
-         DO jquad = 1, num_integ_points
+         IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
+            CALL rpa_sigma_create(rpa_sigma, mp2_env%ri_rpa%sigma_param, fm_mat_Q(1), unit_nr, para_env)
+         END IF
 
+         DO jquad = 1, num_integ_points
             IF (MODULO(jquad, num_integ_group) /= color_rpa_group) CYCLE
 
             CALL timeset(routineN//"_RPA_matrix_operations", handle3)
@@ -1453,7 +1458,6 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
                END IF
 
             ELSE
-
                IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A, 1X, I3, 1X, A, 1X, I3)") &
                   "INTEG_INFO| Started with Integration point", jquad, "of", num_integ_points
 
@@ -1515,7 +1519,13 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
                e_exchange = e_exchange + e_exchange_corr*wj(jquad)
             END IF
 
+            ! for developing Sigma functional  closed and open shell are taken cared for
+            IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
+               CALL rpa_sigma_matrix_spectral(rpa_sigma, fm_mat_Q(1), wj(jquad), para_env_RPA)
+            END IF
+
             IF (do_ri_sos_laplace_mp2) THEN
+
                CALL SOS_MP2_postprocessing(fm_mat_Q, Erpa, tau_wj(jquad))
 
                IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, &
@@ -1579,9 +1589,15 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
 
          END DO ! jquad
 
+         IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
+            CALL finalize_rpa_sigma(rpa_sigma, unit_nr, mp2_env%ri_rpa%e_sigma_corr, para_env, do_minimax_quad)
+            IF (do_minimax_quad) mp2_env%ri_rpa%e_sigma_corr = mp2_env%ri_rpa%e_sigma_corr/2.0_dp
+            CALL para_env%sum(mp2_env%ri_rpa%e_sigma_corr)
+         END IF
+
          CALL para_env%sum(Erpa)
 
-         IF (.NOT. do_ri_sos_laplace_mp2) THEN
+         IF (.NOT. (do_ri_sos_laplace_mp2)) THEN
             Erpa = Erpa/(pi*2.0_dp)
             IF (do_minimax_quad) Erpa = Erpa/2.0_dp
          END IF