diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f6cfc73875..d70e2e5cdf 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/input_constants.F b/src/input_constants.F index b9b0f566d4..f1a6bab8d2 100644 --- a/src/input_constants.F +++ b/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 diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F index 0cc7ed4f9b..2398c55ce6 100644 --- a/src/input_cp2k_mp2.F +++ b/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.", & diff --git a/src/mp2.F b/src/mp2.F index e004521c0d..c33c0516a0 100644 --- a/src/mp2.F +++ b/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 diff --git a/src/mp2_gpw.F b/src/mp2_gpw.F index b034ba3a69..e52f24001e 100644 --- a/src/mp2_gpw.F +++ b/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 diff --git a/src/mp2_setup.F b/src/mp2_setup.F index a50d426e8a..9cb85265ac 100644 --- a/src/mp2_setup.F +++ b/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) diff --git a/src/mp2_types.F b/src/mp2_types.F index f074c96e22..66b0afaf54 100644 --- a/src/mp2_types.F +++ b/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 diff --git a/src/rpa_main.F b/src/rpa_main.F index 94676f1ed8..a09a614bba 100644 --- a/src/rpa_main.F +++ b/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,11 +1329,10 @@ 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. @@ -1337,7 +1340,6 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s 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 diff --git a/src/rpa_sigma_functional.F b/src/rpa_sigma_functional.F new file mode 100644 index 0000000000..0bbe4ab78c --- /dev/null +++ b/src/rpa_sigma_functional.F @@ -0,0 +1,746 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Routines to calculate RI-RPA energy and Sigma correction to the RPA energies +!> using the cubic spline based on eigen values of Q(w). +!> \par History +! ************************************************************************************************** +MODULE rpa_sigma_functional + USE cp_fm_diag, ONLY: choose_eigv_solver + USE cp_fm_struct, ONLY: cp_fm_struct_type + USE cp_fm_types, ONLY: cp_fm_create,& + cp_fm_get_info,& + cp_fm_release,& + cp_fm_to_fm,& + cp_fm_type + USE input_constants, ONLY: sigma_PBE0_S1,& + sigma_PBE0_S2,& + sigma_PBE_S1,& + sigma_PBE_S2,& + sigma_none + USE kinds, ONLY: dp + USE machine, ONLY: m_flush + USE mathconstants, ONLY: pi + USE message_passing, ONLY: mp_comm_type,& + mp_para_env_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_sigma_functional' + + PUBLIC :: rpa_sigma_matrix_spectral, rpa_sigma_create, rpa_sigma_type, finalize_rpa_sigma + + TYPE rpa_sigma_type + PRIVATE + REAL(KIND=dp) :: e_sigma_corr = 0.0_dp + REAL(KIND=dp) :: e_rpa_by_eig_val = 0.0_dp + INTEGER :: sigma_param = 0 + TYPE(cp_fm_type) :: mat_Q_diagonal = cp_fm_type() + TYPE(cp_fm_type) :: fm_evec = cp_fm_type() + REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: sigma_eigenvalue + INTEGER :: dimen_RI_red = 0 + END TYPE + +CONTAINS + +! ************************************************************************************************** +!> \brief ... Collect the Q(w) (fm_mat_Q) matrix to create rpa_sigma a derived type variable. +!> and write out the choosen parametrization for the cubic spline. +!> \param rpa_sigma ... +!> \param sigma_param ... +!> \param fm_mat_Q ... +!> \param unit_nr ... +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE rpa_sigma_create(rpa_sigma, sigma_param, fm_mat_Q, unit_nr, para_env) + + TYPE(rpa_sigma_type), INTENT(OUT) :: rpa_sigma + INTEGER :: sigma_param + TYPE(cp_fm_type) :: fm_mat_Q + INTEGER :: unit_nr + + CLASS(mp_comm_type), INTENT(IN) :: para_env + + TYPE(cp_fm_struct_type), POINTER :: matrix_struct + + ! Getting information about the Q matrix and creating initializing two matrices to pass it to the diagonalising driver. + CALL cp_fm_get_info(fm_mat_Q, matrix_struct=matrix_struct, nrow_global=rpa_sigma%dimen_RI_red) + + ALLOCATE (rpa_sigma%sigma_eigenvalue(rpa_sigma%dimen_RI_red)) + + CALL cp_fm_create(rpa_sigma%fm_evec, matrix_struct) + CALL cp_fm_create(rpa_sigma%mat_Q_diagonal, matrix_struct) + + rpa_sigma%sigma_param = sigma_param + + SELECT CASE (rpa_sigma%sigma_param) + + CASE (sigma_none) + ! There is nothing to do + CASE DEFAULT + CPABORT("Unknown parameterization") + + CASE (sigma_PBE0_S1) + IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") & + "SIGMA_INFO| Sigma eigenvalues parameterized with PBE0_S1 reference" + + CASE (sigma_PBE0_S2) + IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") & + "SIGMA_INFO| Sigma eigenvalues parameterized with PBE0_S2 reference" + + CASE (sigma_PBE_S1) + IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") & + "SIGMA_INFO| Sigma eigenvalues parameterized with PBE_S1 reference" + + CASE (sigma_PBE_S2) + IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") & + "SIGMA_INFO| Sigma eigenvalues parameterized with PBE_S2 reference" + END SELECT + IF (unit_nr > 0) CALL m_flush(unit_nr) + CALL para_env%sync() + + END SUBROUTINE rpa_sigma_create + +! ************************************************************************************************** +!> \brief ... Memory cleanup routine +!> \param rpa_sigma ... +! ************************************************************************************************** + SUBROUTINE rpa_sigma_cleanup(rpa_sigma) + + TYPE(rpa_sigma_type), INTENT(INOUT) :: rpa_sigma + + CALL cp_fm_release(rpa_sigma%mat_Q_diagonal) + CALL cp_fm_release(rpa_sigma%fm_evec) + DEALLOCATE (rpa_sigma%sigma_eigenvalue) + + END SUBROUTINE rpa_sigma_cleanup + +! ************************************************************************************************** +!> \brief ... Diagonalize and store the eigenvalues of fm_mat_Q in rpa_sigma%sigma_eigenvalue. +!> \param rpa_sigma ... +!> \param fm_mat_Q ... +!> \param wj ... +!> \param para_env_RPA ... +! ************************************************************************************************** + SUBROUTINE rpa_sigma_matrix_spectral(rpa_sigma, fm_mat_Q, wj, para_env_RPA) + + TYPE(rpa_sigma_type) :: rpa_sigma + TYPE(cp_fm_type) :: fm_mat_Q + REAL(KIND=dp) :: wj + TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA + + ! copy the Q matrix into the dummy matrix to avoid changing it. + CALL cp_fm_to_fm(fm_mat_Q, rpa_sigma%mat_Q_diagonal) + + !diagonalising driver + CALL choose_eigv_solver(rpa_sigma%mat_Q_diagonal, rpa_sigma%fm_evec, rpa_sigma%sigma_eigenvalue) + + ! Computing the integration to calculate the sigma correction. + CALL compute_e_sigma_corr_by_freq_int(rpa_sigma, wj, para_env_RPA) + + END SUBROUTINE rpa_sigma_matrix_spectral +! ************************************************************************************************** + +! ************************************************************************************************** +!> \brief ... To compute the e_sigma_corr and e_rpa_by_eig_val by freq integration over the eigenvalues of Q(w) +!> e_sigma_corr = - H(sigma) & e_rpa_by_eig_val = log(1+sigma)-sigma +!> \param rpa_sigma ... +!> \param wj ... +!> \param para_env_RPA ... +! ************************************************************************************************** + SUBROUTINE compute_e_sigma_corr_by_freq_int(rpa_sigma, wj, para_env_RPA) + TYPE(rpa_sigma_type), INTENT(INOUT) :: rpa_sigma + REAL(KIND=dp), INTENT(IN) :: wj + TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA + + INTEGER :: iaux + REAL(KIND=dp) :: dedw_rpa, dedw_sigma + + dedw_sigma = 0.0_dp + dedw_rpa = 0.0_dp + + ! Loop which take each eigenvalue to: i) get E_RPA & ii) integrates the spline to get E_c correction + DO iaux = 1, rpa_sigma%dimen_RI_red + IF (rpa_sigma%sigma_eigenvalue(iaux) > 0.0_dp) THEN + IF (MODULO(iaux, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE + dedw_rpa = dedw_rpa + LOG(1.0_dp + rpa_sigma%sigma_eigenvalue(iaux)) - rpa_sigma%sigma_eigenvalue(iaux) + IF (MODULO(iaux, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE + dedw_sigma = dedw_sigma - cubic_spline_integr(rpa_sigma%sigma_eigenvalue(iaux), rpa_sigma%sigma_param) + END IF + END DO + + ! (use 2.0_dp its better for compilers) + rpa_sigma%e_sigma_corr = rpa_sigma%e_sigma_corr + (wj*dedw_sigma/(2.0_dp*pi*2.0_dp)) + rpa_sigma%e_rpa_by_eig_val = rpa_sigma%e_rpa_by_eig_val + (wj*dedw_rpa/(2.0_dp*pi*2.0_dp)) + + END SUBROUTINE compute_e_sigma_corr_by_freq_int + +! ************************************************************************************************** +!> \brief ... Save the calculated value of E_c correction to the global variable and memory clean. +!> \param rpa_sigma ... +!> \param unit_nr ... +!> \param e_sigma_corr ... +!> \param para_env ... +!> \param do_minimax_quad ... +! ************************************************************************************************** + SUBROUTINE finalize_rpa_sigma(rpa_sigma, unit_nr, e_sigma_corr, para_env, do_minimax_quad) + TYPE(rpa_sigma_type), INTENT(INOUT) :: rpa_sigma + INTEGER :: unit_nr + REAL(KIND=dp), INTENT(OUT) :: e_sigma_corr + TYPE(mp_para_env_type), INTENT(IN) :: para_env + LOGICAL, INTENT(IN) :: do_minimax_quad + + IF (do_minimax_quad) rpa_sigma%e_rpa_by_eig_val = rpa_sigma%e_rpa_by_eig_val/2.0_dp + CALL para_env%sum(rpa_sigma%e_rpa_by_eig_val) + e_sigma_corr = rpa_sigma%e_sigma_corr + IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') & + 'RI-RPA energy from eigenvalues of Q(w) = ', & + rpa_sigma%e_rpa_by_eig_val + + CALL rpa_sigma_cleanup(rpa_sigma) + + END SUBROUTINE finalize_rpa_sigma + +! ************************************************************************************************** +!> \brief ... integrates cubic spline to get eigenvalue sigma based E_c correction. +!> \param sigma ... +!> \param sigma_param ... +!> \return ... Integration value H(sigma) wrt to coupling constant eq 35 in Trushin et al JCP 2021 +! ************************************************************************************************** + FUNCTION cubic_spline_integr(sigma, sigma_param) RESULT(integral) + REAL(KIND=dp), INTENT(in) :: sigma + INTEGER :: sigma_param + REAL(KIND=dp) :: integral + + INTEGER :: i, m, n + REAL(KIND=dp) :: h + REAL(KIND=dp), ALLOCATABLE :: coeff(:, :), x_coord(:) + + SELECT CASE (sigma_param) + CASE (sigma_PBE0_S1) + n = 21 + ALLOCATE (x_coord(n)) + ALLOCATE (coeff(4, n)) + + coeff(1, 1) = 0.000000000000D+00 + coeff(1, 2) = 0.000000000000D+00 + coeff(1, 3) = -0.149500660756D-03 + coeff(1, 4) = -0.353017276233D-02 + coeff(1, 5) = -0.109810247734D-01 + coeff(1, 6) = -0.231246943777D-01 + coeff(1, 7) = -0.268999962858D-01 + coeff(1, 8) = -0.634751994007D-03 + coeff(1, 9) = 0.118792892470D-01 + coeff(1, 10) = -0.473431931326D-01 + coeff(1, 11) = -0.817589390539D-01 + coeff(1, 12) = 0.125726011069D-01 + coeff(1, 13) = 0.108028492092D+00 + coeff(1, 14) = 0.193548206759D+00 + coeff(1, 15) = 0.358395561305D-01 + coeff(1, 16) = -0.497714974829D-01 + coeff(1, 17) = 0.341059348835D-01 + coeff(1, 18) = 0.341050720155D-01 + coeff(1, 19) = 0.785549033229D-01 + coeff(1, 20) = 0.000000000000D+00 + coeff(1, 21) = 0.000000000000D+00 + coeff(2, 1) = 0.000000000000D+00 + coeff(2, 2) = 0.000000000000D+00 + coeff(2, 3) = -0.208376539581D+01 + coeff(2, 4) = -0.469755869285D+01 + coeff(2, 5) = -0.565503803415D+01 + coeff(2, 6) = -0.135502867642D+01 + coeff(2, 7) = 0.000000000000D+00 + coeff(2, 8) = 0.284340701746D+01 + coeff(2, 9) = 0.000000000000D+00 + coeff(2, 10) = -0.342695931351D+01 + coeff(2, 11) = 0.000000000000D+00 + coeff(2, 12) = 0.358739081268D+01 + coeff(2, 13) = 0.203368806130D+01 + coeff(2, 14) = 0.000000000000D+00 + coeff(2, 15) = -0.901387663218D+00 + coeff(2, 16) = 0.000000000000D+00 + coeff(2, 17) = 0.000000000000D+00 + coeff(2, 18) = 0.000000000000D+00 + coeff(2, 19) = 0.000000000000D+00 + coeff(2, 20) = 0.000000000000D+00 + coeff(2, 21) = 0.000000000000D+00 + coeff(3, 1) = -0.000000000000D+00 + coeff(3, 2) = -0.322176662524D+05 + coeff(3, 3) = -0.267090835643D+04 + coeff(3, 4) = -0.373532067350D+04 + coeff(3, 5) = -0.797121299000D+03 + coeff(3, 6) = 0.111299540119D+03 + coeff(3, 7) = 0.299284621116D+04 + coeff(3, 8) = -0.319333485618D+02 + coeff(3, 9) = -0.140910103454D+04 + coeff(3, 10) = -0.848330431187D+01 + coeff(3, 11) = 0.435025012278D+03 + coeff(3, 12) = -0.700327539634D+01 + coeff(3, 13) = 0.545486142353D+01 + coeff(3, 14) = -0.453346282407D+02 + coeff(3, 15) = 0.371921027910D+00 + coeff(3, 16) = 0.464101795796D+01 + coeff(3, 17) = -0.190069531714D-04 + coeff(3, 18) = 0.514345336660D-01 + coeff(3, 19) = -0.431543078188D-02 + coeff(3, 20) = -0.000000000000D+00 + coeff(3, 21) = 0.000000000000D+00 + coeff(4, 1) = 0.000000000000D+00 + coeff(4, 2) = 0.152897717268D+09 + coeff(4, 3) = 0.902815532735D+06 + coeff(4, 4) = 0.191760493084D+07 + coeff(4, 5) = 0.445372471512D+06 + coeff(4, 6) = 0.188362654331D+04 + coeff(4, 7) = -0.383203258784D+06 + coeff(4, 8) = -0.170027418959D+05 + coeff(4, 9) = 0.819629330224D+05 + coeff(4, 10) = 0.560228610945D+04 + coeff(4, 11) = -0.108203002413D+05 + coeff(4, 12) = -0.363378668069D+03 + coeff(4, 13) = -0.260332257619D+03 + coeff(4, 14) = 0.291068208088D+03 + coeff(4, 15) = 0.122322834276D+02 + coeff(4, 16) = -0.132875656470D+02 + coeff(4, 17) = 0.343356030115D-04 + coeff(4, 18) = -0.212958640167D-01 + coeff(4, 19) = 0.389311916174D-03 + coeff(4, 20) = 0.000000000000D+00 + coeff(4, 21) = 0.000000000000D+00 + x_coord(1) = 0.000000000000D+00 + x_coord(2) = 0.100000000000D-04 + x_coord(3) = 0.100000000000D-03 + x_coord(4) = 0.100000000000D-02 + x_coord(5) = 0.215443469000D-02 + x_coord(6) = 0.464158883000D-02 + x_coord(7) = 0.100000000000D-01 + x_coord(8) = 0.146779926762D-01 + x_coord(9) = 0.215443469003D-01 + x_coord(10) = 0.316227766017D-01 + x_coord(11) = 0.464158883361D-01 + x_coord(12) = 0.681292069058D-01 + x_coord(13) = 0.100000000000D+00 + x_coord(14) = 0.158489319246D+00 + x_coord(15) = 0.251188643151D+00 + x_coord(16) = 0.398107170553D+00 + x_coord(17) = 0.630957344480D+00 + x_coord(18) = 0.100000000000D+01 + x_coord(19) = 0.261015721568D+01 + x_coord(20) = 0.100000000000D+02 + x_coord(21) = 0.215443469000D+02 + + CASE (sigma_PBE0_S2) + n = 21 + ALLOCATE (x_coord(n)) + ALLOCATE (coeff(4, n)) + + coeff(1, 1) = 0.000000000000D+00 + coeff(1, 2) = 0.000000000000D+00 + coeff(1, 3) = -0.431405252048D-04 + coeff(1, 4) = -0.182874853131D-02 + coeff(1, 5) = -0.852003132762D-02 + coeff(1, 6) = -0.218177403992D-01 + coeff(1, 7) = -0.305777654735D-01 + coeff(1, 8) = -0.870882903969D-02 + coeff(1, 9) = 0.137878988102D-01 + coeff(1, 10) = -0.284352007440D-01 + coeff(1, 11) = -0.798812002431D-01 + coeff(1, 12) = -0.334010771574D-02 + coeff(1, 13) = 0.934182748715D-01 + coeff(1, 14) = 0.204960802253D+00 + coeff(1, 15) = 0.213204380281D-01 + coeff(1, 16) = -0.401220283037D-01 + coeff(1, 17) = 0.321629738336D-01 + coeff(1, 18) = 0.321618301891D-01 + coeff(1, 19) = 0.808763912948D-01 + coeff(1, 20) = 0.000000000000D+00 + coeff(1, 21) = 0.000000000000D+00 + coeff(2, 1) = 0.000000000000D+00 + coeff(2, 2) = 0.000000000000D+00 + coeff(2, 3) = -0.661870777583D+00 + coeff(2, 4) = -0.289752912590D+01 + coeff(2, 5) = -0.558979946652D+01 + coeff(2, 6) = -0.267765704540D+01 + coeff(2, 7) = 0.000000000000D+00 + coeff(2, 8) = 0.389592612611D+01 + coeff(2, 9) = 0.000000000000D+00 + coeff(2, 10) = -0.382296397421D+01 + coeff(2, 11) = 0.000000000000D+00 + coeff(2, 12) = 0.327772498106D+01 + coeff(2, 13) = 0.239633724310D+01 + coeff(2, 14) = 0.000000000000D+00 + coeff(2, 15) = -0.726304793204D+00 + coeff(2, 16) = 0.000000000000D+00 + coeff(2, 17) = 0.000000000000D+00 + coeff(2, 18) = 0.000000000000D+00 + coeff(2, 19) = 0.000000000000D+00 + coeff(2, 20) = 0.000000000000D+00 + coeff(2, 21) = 0.000000000000D+00 + coeff(3, 1) = -0.000000000000D+00 + coeff(3, 2) = -0.862385254713D+04 + coeff(3, 3) = -0.192306222883D+04 + coeff(3, 4) = -0.520047462362D+04 + coeff(3, 5) = -0.877473657666D+03 + coeff(3, 6) = 0.841408344046D+02 + coeff(3, 7) = 0.216516760964D+04 + coeff(3, 8) = 0.296702212913D+03 + coeff(3, 9) = -0.867733655494D+03 + coeff(3, 10) = -0.188410055380D+03 + coeff(3, 11) = 0.336084151111D+03 + coeff(3, 12) = 0.489746728744D+01 + coeff(3, 13) = 0.158746877181D+02 + coeff(3, 14) = -0.562764882273D+02 + coeff(3, 15) = 0.134759277149D+01 + coeff(3, 16) = 0.399959778866D+01 + coeff(3, 17) = -0.251917983154D-04 + coeff(3, 18) = 0.563694092760D-01 + coeff(3, 19) = -0.444296223097D-02 + coeff(3, 20) = -0.000000000000D+00 + coeff(3, 21) = 0.000000000000D+00 + coeff(4, 1) = 0.000000000000D+00 + coeff(4, 2) = 0.366429086790D+08 + coeff(4, 3) = 0.504466528222D+06 + coeff(4, 4) = 0.232980923705D+07 + coeff(4, 5) = 0.392124287301D+06 + coeff(4, 6) = 0.206173887726D+05 + coeff(4, 7) = -0.249217659838D+06 + coeff(4, 8) = -0.563519876566D+05 + coeff(4, 9) = 0.448530826095D+05 + coeff(4, 10) = 0.143140667434D+05 + coeff(4, 11) = -0.800144415404D+04 + coeff(4, 12) = -0.391685311241D+03 + coeff(4, 13) = -0.414433988077D+03 + coeff(4, 14) = 0.376550449117D+03 + coeff(4, 15) = 0.510124747789D+01 + coeff(4, 16) = -0.114511339236D+02 + coeff(4, 17) = 0.455083767664D-04 + coeff(4, 18) = -0.233390912502D-01 + coeff(4, 19) = 0.400817027790D-03 + coeff(4, 20) = 0.000000000000D+00 + coeff(4, 21) = 0.000000000000D+00 + x_coord(1) = 0.000000000000D+00 + x_coord(2) = 0.100000000000D-04 + x_coord(3) = 0.100000000000D-03 + x_coord(4) = 0.100000000000D-02 + x_coord(5) = 0.215443469000D-02 + x_coord(6) = 0.464158883000D-02 + x_coord(7) = 0.100000000000D-01 + x_coord(8) = 0.146779926762D-01 + x_coord(9) = 0.215443469003D-01 + x_coord(10) = 0.316227766017D-01 + x_coord(11) = 0.464158883361D-01 + x_coord(12) = 0.681292069058D-01 + x_coord(13) = 0.100000000000D+00 + x_coord(14) = 0.158489319246D+00 + x_coord(15) = 0.251188643151D+00 + x_coord(16) = 0.398107170553D+00 + x_coord(17) = 0.630957344480D+00 + x_coord(18) = 0.100000000000D+01 + x_coord(19) = 0.261015721568D+01 + x_coord(20) = 0.100000000000D+02 + x_coord(21) = 0.215443469000D+02 + + CASE (sigma_PBE_S1) + n = 22 + ALLOCATE (x_coord(n)) + ALLOCATE (coeff(4, n)) + + coeff(1, 1) = 0.000000000000D+00 + coeff(1, 2) = 0.000000000000D+00 + coeff(1, 3) = -0.493740326815D-04 + coeff(1, 4) = -0.136110637329D-02 + coeff(1, 5) = -0.506905111755D-02 + coeff(1, 6) = -0.127411222930D-01 + coeff(1, 7) = -0.220144968504D-01 + coeff(1, 8) = -0.239939034695D-01 + coeff(1, 9) = -0.436386416290D-01 + coeff(1, 10) = -0.117890214262D+00 + coeff(1, 11) = -0.141123921668D+00 + coeff(1, 12) = 0.865524876740D-01 + coeff(1, 13) = 0.179390274565D+00 + coeff(1, 14) = 0.269368658116D+00 + coeff(1, 15) = 0.785040456996D-01 + coeff(1, 16) = 0.490248637276D-01 + coeff(1, 17) = -0.111571911794D+00 + coeff(1, 18) = -0.197712184164D-01 + coeff(1, 19) = -0.197716870218D-01 + coeff(1, 20) = -0.372253617253D-01 + coeff(1, 21) = 0.000000000000D+00 + coeff(1, 22) = 0.000000000000D+00 + coeff(2, 1) = 0.000000000000D+00 + coeff(2, 2) = 0.000000000000D+00 + coeff(2, 3) = -0.709484897949D+00 + coeff(2, 4) = -0.197447407686D+01 + coeff(2, 5) = -0.315478745349D+01 + coeff(2, 6) = -0.229603163128D+01 + coeff(2, 7) = -0.670801534786D+00 + coeff(2, 8) = -0.704199644986D+00 + coeff(2, 9) = -0.400987325224D+01 + coeff(2, 10) = -0.269982990241D+01 + coeff(2, 11) = 0.000000000000D+00 + coeff(2, 12) = 0.472814414167D+01 + coeff(2, 13) = 0.207638470052D+01 + coeff(2, 14) = 0.000000000000D+00 + coeff(2, 15) = -0.389846972557D+00 + coeff(2, 16) = -0.298496119087D+00 + coeff(2, 17) = 0.000000000000D+00 + coeff(2, 18) = 0.000000000000D+00 + coeff(2, 19) = -0.601781536636D-06 + coeff(2, 20) = 0.000000000000D+00 + coeff(2, 21) = 0.000000000000D+00 + coeff(2, 22) = 0.000000000000D+00 + coeff(3, 1) = -0.000000000000D+00 + coeff(3, 2) = -0.104035132381D+05 + coeff(3, 3) = -0.108777473624D+04 + coeff(3, 4) = -0.219328637518D+04 + coeff(3, 5) = -0.260711341283D+03 + coeff(3, 6) = 0.132509852177D+02 + coeff(3, 7) = 0.165970301474D+03 + coeff(3, 8) = -0.460909893146D+03 + coeff(3, 9) = -0.112939707971D+04 + coeff(3, 10) = 0.465035067500D+02 + coeff(3, 11) = 0.123097490767D+04 + coeff(3, 12) = -0.876616265219D+02 + coeff(3, 13) = 0.790484996078D+01 + coeff(3, 14) = -0.624281400584D+02 + coeff(3, 15) = 0.324152775194D+01 + coeff(3, 16) = -0.632212496608D+01 + coeff(3, 17) = 0.202215332970D+01 + coeff(3, 18) = -0.308693235932D-06 + coeff(3, 19) = -0.495067060383D-02 + coeff(3, 20) = 0.116855980641D-02 + coeff(3, 21) = -0.000000000000D+00 + coeff(3, 22) = 0.000000000000D+00 + coeff(4, 1) = 0.000000000000D+00 + coeff(4, 2) = 0.478661516427D+08 + coeff(4, 3) = 0.285187385316D+06 + coeff(4, 4) = 0.971371823345D+06 + coeff(4, 5) = 0.116156741398D+06 + coeff(4, 6) = 0.172191903906D+05 + coeff(4, 7) = -0.241613612898D+05 + coeff(4, 8) = 0.213790845631D+05 + coeff(4, 9) = 0.790063233314D+05 + coeff(4, 10) = 0.201667888760D+04 + coeff(4, 11) = -0.344519214370D+05 + coeff(4, 12) = 0.963471669433D+03 + coeff(4, 13) = -0.292417702205D+03 + coeff(4, 14) = 0.433842720035D+03 + coeff(4, 15) = -0.132982468090D+02 + coeff(4, 16) = 0.199358142858D+02 + coeff(4, 17) = -0.365297127483D+01 + coeff(4, 18) = 0.434041376596D-07 + coeff(4, 19) = 0.101490424907D-02 + coeff(4, 20) = -0.796902275213D-04 + coeff(4, 21) = 0.000000000000D+00 + coeff(4, 22) = 0.000000000000D+00 + x_coord(1) = 0.000000000000D+00 + x_coord(2) = 0.100000000000D-04 + x_coord(3) = 0.100000000000D-03 + x_coord(4) = 0.100000000000D-02 + x_coord(5) = 0.215443469000D-02 + x_coord(6) = 0.464158883000D-02 + x_coord(7) = 0.100000000000D-01 + x_coord(8) = 0.146779926762D-01 + x_coord(9) = 0.215443469003D-01 + x_coord(10) = 0.316227766017D-01 + x_coord(11) = 0.464158883361D-01 + x_coord(12) = 0.681292069058D-01 + x_coord(13) = 0.100000000000D+00 + x_coord(14) = 0.158489319246D+00 + x_coord(15) = 0.251188643151D+00 + x_coord(16) = 0.398107170553D+00 + x_coord(17) = 0.630957344480D+00 + x_coord(18) = 0.100000000000D+01 + x_coord(19) = 0.237137370566D+01 + x_coord(20) = 0.562341325000D+01 + x_coord(21) = 0.153992652606D+02 + x_coord(22) = 0.316227766000D+02 + + CASE (sigma_PBE_S2) + n = 22 + ALLOCATE (x_coord(n)) + ALLOCATE (coeff(4, n)) + + coeff(1, 1) = 0.000000000000D+00 + coeff(1, 2) = 0.000000000000D+00 + coeff(1, 3) = -0.156157535801D-03 + coeff(1, 4) = -0.365199003270D-02 + coeff(1, 5) = -0.108302033233D-01 + coeff(1, 6) = -0.203436953346D-01 + coeff(1, 7) = -0.214330355346D-01 + coeff(1, 8) = 0.109617244934D-03 + coeff(1, 9) = 0.813969827075D-02 + coeff(1, 10) = -0.701367130014D-01 + coeff(1, 11) = -0.162002361715D+00 + coeff(1, 12) = 0.337288711362D-01 + coeff(1, 13) = 0.140348429629D+00 + coeff(1, 14) = 0.271234417677D+00 + coeff(1, 15) = 0.780732751240D-01 + coeff(1, 16) = 0.436066976238D-01 + coeff(1, 17) = -0.106097689688D+00 + coeff(1, 18) = -0.133141637069D-01 + coeff(1, 19) = -0.133143525246D-01 + coeff(1, 20) = -0.430994711278D-01 + coeff(1, 21) = 0.000000000000D+00 + coeff(1, 22) = 0.000000000000D+00 + coeff(2, 1) = 0.000000000000D+00 + coeff(2, 2) = 0.000000000000D+00 + coeff(2, 3) = -0.217211651544D+01 + coeff(2, 4) = -0.473638379726D+01 + coeff(2, 5) = -0.487821808504D+01 + coeff(2, 6) = -0.433631413905D+00 + coeff(2, 7) = 0.000000000000D+00 + coeff(2, 8) = 0.193813387881D+01 + coeff(2, 9) = 0.000000000000D+00 + coeff(2, 10) = -0.695060290528D+01 + coeff(2, 11) = 0.000000000000D+00 + coeff(2, 12) = 0.502541925806D+01 + coeff(2, 13) = 0.273498669354D+01 + coeff(2, 14) = 0.000000000000D+00 + coeff(2, 15) = -0.448708826169D+00 + coeff(2, 16) = -0.332102918195D+00 + coeff(2, 17) = 0.000000000000D+00 + coeff(2, 18) = 0.000000000000D+00 + coeff(2, 19) = -0.242488141082D-06 + coeff(2, 20) = 0.000000000000D+00 + coeff(2, 21) = 0.000000000000D+00 + coeff(2, 22) = 0.000000000000D+00 + coeff(3, 1) = -0.000000000000D+00 + coeff(3, 2) = -0.337014964214D+05 + coeff(3, 3) = -0.285795351280D+04 + coeff(3, 4) = -0.372723918347D+04 + coeff(3, 5) = -0.516689374427D+03 + coeff(3, 6) = 0.480322803175D+02 + coeff(3, 7) = 0.253894893657D+04 + coeff(3, 8) = -0.535684993409D+02 + coeff(3, 9) = -0.162223464755D+04 + coeff(3, 10) = -0.319667723139D+03 + coeff(3, 11) = 0.101401359817D+04 + coeff(3, 12) = -0.862770702569D+02 + coeff(3, 13) = 0.212578002151D+02 + coeff(3, 14) = -0.625949163782D+02 + coeff(3, 15) = 0.357838438707D+01 + coeff(3, 16) = -0.543078279308D+01 + coeff(3, 17) = 0.204380282001D+01 + coeff(3, 18) = -0.124376927880D-06 + coeff(3, 19) = -0.844892173480D-02 + coeff(3, 20) = 0.135295689023D-02 + coeff(3, 21) = -0.000000000000D+00 + coeff(3, 22) = 0.000000000000D+00 + coeff(4, 1) = 0.000000000000D+00 + coeff(4, 2) = 0.160253203309D+09 + coeff(4, 3) = 0.106174857663D+07 + coeff(4, 4) = 0.211694319531D+07 + coeff(4, 5) = 0.377995031873D+06 + coeff(4, 6) = -0.941771032564D+03 + coeff(4, 7) = -0.332306990184D+06 + coeff(4, 8) = -0.850176312292D+04 + coeff(4, 9) = 0.844978829514D+05 + coeff(4, 10) = 0.249933770676D+05 + coeff(4, 11) = -0.275803550484D+05 + coeff(4, 12) = 0.105308484492D+04 + coeff(4, 13) = -0.508788318128D+03 + coeff(4, 14) = 0.432758845019D+03 + coeff(4, 15) = -0.144367801061D+02 + coeff(4, 16) = 0.175904487062D+02 + coeff(4, 17) = -0.369208055752D+01 + coeff(4, 18) = 0.174842962107D-07 + coeff(4, 19) = 0.173203285755D-02 + coeff(4, 20) = -0.922652326542D-04 + coeff(4, 21) = 0.000000000000D+00 + coeff(4, 22) = 0.000000000000D+00 + x_coord(1) = 0.000000000000D+00 + x_coord(2) = 0.100000000000D-04 + x_coord(3) = 0.100000000000D-03 + x_coord(4) = 0.100000000000D-02 + x_coord(5) = 0.215443469000D-02 + x_coord(6) = 0.464158883000D-02 + x_coord(7) = 0.100000000000D-01 + x_coord(8) = 0.146779926762D-01 + x_coord(9) = 0.215443469003D-01 + x_coord(10) = 0.316227766017D-01 + x_coord(11) = 0.464158883361D-01 + x_coord(12) = 0.681292069058D-01 + x_coord(13) = 0.100000000000D+00 + x_coord(14) = 0.158489319246D+00 + x_coord(15) = 0.251188643151D+00 + x_coord(16) = 0.398107170553D+00 + x_coord(17) = 0.630957344480D+00 + x_coord(18) = 0.100000000000D+01 + x_coord(19) = 0.237137370566D+01 + x_coord(20) = 0.562341325000D+01 + x_coord(21) = 0.153992652606D+02 + x_coord(22) = 0.316227766000D+02 + + END SELECT + + ! determine to which interval sigma eigenvalue belongs + m = intervalnum(x_coord, n, sigma) + + ! Numerically evaluate integral + integral = 0.0_dp + IF (m == 1) THEN + integral = 0.5_dp*coeff(2, 1)*sigma + END IF + IF ((m > 1) .AND. (m < n)) THEN + h = sigma - x_coord(m) + integral = 0.5_dp*coeff(2, 1)*x_coord(2)**2/sigma & + + (coeff(1, m)*h + coeff(2, m)/2.0_dp*h**2 + & + coeff(3, m)/3.0_dp*h**3 + coeff(4, m)/4.0_dp*h**4)/sigma + DO i = 2, m - 1 + h = x_coord(i + 1) - x_coord(i) + integral = integral & + + (coeff(1, i)*h + coeff(2, i)/2.0_dp*h**2 + & + coeff(3, i)/3.0_dp*h**3 + coeff(4, i)/4.0_dp*h**4)/sigma + END DO + END IF + IF (m == n) THEN + integral = 0.5_dp*coeff(2, 1)*x_coord(2)**2/sigma + DO i = 2, m - 1 + h = x_coord(i + 1) - x_coord(i) + integral = integral & + + (coeff(1, i)*h + coeff(2, i)/2.0_dp*h**2 & + + coeff(3, i)/3.0_dp*h**3 + coeff(4, i)/4.0_dp*h**4)/sigma + END DO + integral = integral + coeff(1, n)*(1.0_dp - x_coord(n)/sigma) + END IF + integral = integral*sigma + + DEALLOCATE (x_coord, coeff) + + END FUNCTION cubic_spline_integr + +! ************************************************************************************************** +!> \brief ... Determine the interval which contains the sigma eigenvalue. +!> \param x_coord ... +!> \param n ... +!> \param sigma ... +!> \return ... an integer m +! ************************************************************************************************** + INTEGER FUNCTION intervalnum(x_coord, n, sigma) RESULT(inum) + + INTEGER :: n + REAL(KIND=dp), INTENT(in) :: x_coord(n), sigma + + INTEGER :: i + + IF (sigma <= 0.0_dp) STOP 'intervalnum: sigma should be positive' + + inum = -1 + IF (sigma > x_coord(n)) inum = n + DO i = 1, n - 1 + IF ((sigma > x_coord(i)) .AND. (sigma <= x_coord(i + 1))) inum = i + END DO + + IF (inum == -1) STOP 'interval: something was wrong' + + END FUNCTION intervalnum + +END MODULE rpa_sigma_functional diff --git a/src/rpa_util.F b/src/rpa_util.F index 5b2dd605c0..99f3895bb0 100644 --- a/src/rpa_util.F +++ b/src/rpa_util.F @@ -783,8 +783,8 @@ SUBROUTINE calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, homo, & IF (first_cycle) THEN ! In this case just update the matrix (symmetric form) with ! SQRT((epsi_a-epsi_i)/((epsi_a-epsi_i)**2+omega**2)) -!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) & -!$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega) + !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) & + !$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega) DO iiB = 1, ncol_local i_global = col_indices(iiB) @@ -799,8 +799,8 @@ SUBROUTINE calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, homo, & ELSE ! In this case the update has to remove the old omega component thus ! SQRT(((epsi_a-epsi_i)**2+omega_old**2)/((epsi_a-epsi_i)**2+omega**2)) -!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) & -!$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega,omega_old) + !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) & + !$OMP SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega,omega_old) DO iiB = 1, ncol_local i_global = col_indices(iiB) @@ -846,7 +846,7 @@ SUBROUTINE contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat CALL dgemm_counter_start(dgemm_counter) SELECT CASE (mm_style) CASE (wfc_mm_style_gemm) - ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm + ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm !!! CALL parallel_gemm(transa="N", transb="T", m=dimen_RI, n=dimen_RI, k=dimen_ia, alpha=alpha, & matrix_a=fm_mat_S, matrix_b=fm_mat_S, beta=0.0_dp, & matrix_c=fm_mat_Q_gemm) diff --git a/tests/QS/regtest-rpa-sigma/H2O_gas.xyz b/tests/QS/regtest-rpa-sigma/H2O_gas.xyz new file mode 100644 index 0000000000..d1b6342153 --- /dev/null +++ b/tests/QS/regtest-rpa-sigma/H2O_gas.xyz @@ -0,0 +1,5 @@ + 3 + + O 0.000000 0.000000 -0.111000 + H 0.000000 -0.744000 0.495000 + H 0.000000 0.744000 0.495000 diff --git a/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_clenshaw.inp b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_clenshaw.inp new file mode 100644 index 0000000000..bbe21a5b86 --- /dev/null +++ b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_clenshaw.inp @@ -0,0 +1,85 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT RI_RPA_H2O_minimax + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END TIMINGS +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-10 + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-7 + MAX_SCF 100 + SCF_GUESS ATOMIC + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &WF_CORRELATION + MEMORY 200. + NUMBER_PROC 1 + &INTEGRALS + ERI_METHOD GPW + &WFC_GPW + CUTOFF 100 + REL_CUTOFF 20 + &END WFC_GPW + &END INTEGRALS + &RI_RPA + RPA_NUM_QUAD_POINTS 6 + SIGMA_FUNCTIONAL PBE0_S2 + &HF + FRACTION 1.0000000 + &SCREENING + EPS_SCHWARZ 1.0E-8 + SCREEN_ON_INITIAL_P FALSE + &END SCREENING + &END HF + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.000 6.000 6.000 + PERIODIC NONE + &END CELL + &KIND H + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &TOPOLOGY + COORD_FILE_FORMAT xyz + COORD_FILE_NAME H2O_gas.xyz + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax.inp b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax.inp new file mode 100644 index 0000000000..23b53138cd --- /dev/null +++ b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax.inp @@ -0,0 +1,86 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT RI_RPA_H2O_minimax + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END TIMINGS +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-10 + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-7 + MAX_SCF 100 + SCF_GUESS ATOMIC + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &WF_CORRELATION + MEMORY 200. + NUMBER_PROC 1 + &INTEGRALS + ERI_METHOD GPW + &WFC_GPW + CUTOFF 100 + REL_CUTOFF 20 + &END WFC_GPW + &END INTEGRALS + &RI_RPA + MINIMAX + RPA_NUM_QUAD_POINTS 6 + SIGMA_FUNCTIONAL PBE0_S2 + &HF + FRACTION 1.0000000 + &SCREENING + EPS_SCHWARZ 1.0E-8 + SCREEN_ON_INITIAL_P FALSE + &END SCREENING + &END HF + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.000 6.000 6.000 + PERIODIC NONE + &END CELL + &KIND H + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &TOPOLOGY + COORD_FILE_FORMAT xyz + COORD_FILE_NAME H2O_gas.xyz + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax_NUM_INTEG_GROUPS.inp b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax_NUM_INTEG_GROUPS.inp new file mode 100644 index 0000000000..81f450cca6 --- /dev/null +++ b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax_NUM_INTEG_GROUPS.inp @@ -0,0 +1,87 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT RI_RPA_H2O_minimax + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END TIMINGS +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-10 + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-7 + MAX_SCF 100 + SCF_GUESS ATOMIC + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &WF_CORRELATION + MEMORY 200. + NUMBER_PROC 1 + &INTEGRALS + ERI_METHOD GPW + &WFC_GPW + CUTOFF 100 + REL_CUTOFF 20 + &END WFC_GPW + &END INTEGRALS + &RI_RPA + MINIMAX + NUM_INTEG_GROUPS 4 + RPA_NUM_QUAD_POINTS 6 + SIGMA_FUNCTIONAL PBE0_S2 + &HF + FRACTION 1.0000000 + &SCREENING + EPS_SCHWARZ 1.0E-8 + SCREEN_ON_INITIAL_P FALSE + &END SCREENING + &END HF + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.000 6.000 6.000 + PERIODIC NONE + &END CELL + &KIND H + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &TOPOLOGY + COORD_FILE_FORMAT xyz + COORD_FILE_NAME H2O_gas.xyz + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_clenshaw.inp b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_clenshaw.inp new file mode 100644 index 0000000000..5c5cac1837 --- /dev/null +++ b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_clenshaw.inp @@ -0,0 +1,68 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT H_atom + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + MULTIPLICITY 2 + POTENTIAL_FILE_NAME GTH_POTENTIALS + UKS + &MGRID + CUTOFF 300 + REL_CUTOFF 50 + &END MGRID + &QS + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-6 + MAX_SCF 100 + SCF_GUESS RESTART + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &WF_CORRELATION + MEMORY 1000.0 + NUMBER_PROC 1 + &INTEGRALS + &WFC_GPW + CUTOFF 200 + REL_CUTOFF 30 + &END WFC_GPW + &END INTEGRALS + &RI_RPA + RPA_NUM_QUAD_POINTS 6 + SIGMA_FUNCTIONAL PBE0_S2 + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.0000000 5.0000000 5.0000000 + &END CELL + &COORD + H 0.0000000 0.0000000 0.0000000 + &END COORD + &KIND H + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_minimax.inp b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_minimax.inp new file mode 100644 index 0000000000..b06995769c --- /dev/null +++ b/tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_minimax.inp @@ -0,0 +1,69 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT H_atom + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + MULTIPLICITY 2 + POTENTIAL_FILE_NAME GTH_POTENTIALS + UKS + &MGRID + CUTOFF 300 + REL_CUTOFF 50 + &END MGRID + &QS + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-6 + MAX_SCF 100 + SCF_GUESS RESTART + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &WF_CORRELATION + MEMORY 1000.0 + NUMBER_PROC 1 + &INTEGRALS + &WFC_GPW + CUTOFF 200 + REL_CUTOFF 30 + &END WFC_GPW + &END INTEGRALS + &RI_RPA + MINIMAX + RPA_NUM_QUAD_POINTS 6 + SIGMA_FUNCTIONAL PBE0_S2 + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 5.0000000 5.0000000 5.0000000 + &END CELL + &COORD + H 0.0000000 0.0000000 0.0000000 + &END COORD + &KIND H + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &TOPOLOGY + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-rpa-sigma/TEST_FILES b/tests/QS/regtest-rpa-sigma/TEST_FILES new file mode 100644 index 0000000000..27a125882e --- /dev/null +++ b/tests/QS/regtest-rpa-sigma/TEST_FILES @@ -0,0 +1,7 @@ +RPA_SIGMA_H2O_clenshaw.inp 11 1e-09 -17.192268141449258 +RPA_SIGMA_H2O_minimax.inp 11 1e-09 -17.189840389922761 +RPA_SIGMA_H_minimax.inp 11 1e-09 -0.515037791729464 +RPA_SIGMA_H_clenshaw.inp 11 1e-09 -0.515090906944438 +RPA_SIGMA_H2O_minimax_NUM_INTEG_GROUPS.inp 11 1e-09 -17.189840389922765 +#EOF + diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 8dbd0fa6ec..4ead07d16a 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -3,6 +3,7 @@ # Directories have been reordered according the execution time needed for a gfortran pdbg run using 2 MPI tasks # in case a new directory is added just add it at the top of the list.. # the order will be regularly checked and modified... +QS/regtest-rpa-sigma libint TMC/regtest_ana_on_the_fly parallel mpiranks>2 QS/regtest-cusolver cusolvermp QS/regtest-dlaf dlaf