comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

RPA based SIGMA Functionals (#3695)

Browse source 
Co-authored-by: Raviraj Chandreshbhai Mandalia <mandal53@localhost.localdomain>
This commit is contained in:
Raviraj26 2024-09-26 11:15:54 +02:00 • committed by GitHub
parent bce21c9c6e
commit a4477e416d
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
18 changed files with 1233 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

@ -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

6
src/input_constants.F
View file

@ -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

20
src/input_cp2k_mp2.F
View file

@ -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 @@ CONTAINS
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.", &

12
src/mp2.F
View file

@ -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 @@ CONTAINS
! 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 @@ CONTAINS
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 @@ CONTAINS
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

1
src/mp2_gpw.F
View file

@ -265,7 +265,6 @@ CONTAINS
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

5
src/mp2_setup.F
View file

@ -100,6 +100,10 @@ CONTAINS
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 @@ CONTAINS
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)

6
src/mp2_types.F
View file

@ -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

38
src/rpa_main.F
View file

@ -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 @@ CONTAINS
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 @@ CONTAINS
! 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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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

746
src/rpa_sigma_functional.F Normal file
View file

@ -0,0 +1,746 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! 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

10
src/rpa_util.F
View file

@ -783,8 +783,8 @@ CONTAINS
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 @@ CONTAINS
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 @@ CONTAINS
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)

5
tests/QS/regtest-rpa-sigma/H2O_gas.xyz Normal file
View file

@ -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

85
tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_clenshaw.inp Normal file
View file

@ -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

86
tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax.inp Normal file
View file

@ -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

87
tests/QS/regtest-rpa-sigma/RPA_SIGMA_H2O_minimax_NUM_INTEG_GROUPS.inp Normal file
View file

@ -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

68
tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_clenshaw.inp Normal file
View file

@ -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

69
tests/QS/regtest-rpa-sigma/RPA_SIGMA_H_minimax.inp Normal file
View file

@ -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

7
tests/QS/regtest-rpa-sigma/TEST_FILES Normal file
View file

@ -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

1
tests/TEST_DIRS
View file

@ -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