9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-07 05:53:37 +01:00

properly added basis_correction

This commit is contained in:
Emmanuel Giner 2020-04-07 11:03:19 +02:00
parent 51cf96a506
commit 2fc294cd56
15 changed files with 1006 additions and 0 deletions

View File

@ -0,0 +1,3 @@
mu_of_r
ecmd_utils
dft_one_e

View File

@ -0,0 +1,27 @@
================
basis_correction
================
This module proposes the various flavours of the DFT-based basis set correction originally proposed in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714.
This basis set correction relies mainy on :
+) The definition of a range-separation function \mu(r) varying in space to mimic the incompleteness of the basis set used to represent the coulomb interaction. This procedure needs a two-body rdm representing qualitatively the spacial distribution of the opposite spin electron pairs.
Two types of \mu(r) are proposed, according to the strength of correlation, through the keyword "mu_of_r_potential" in the module "mu_of_r":
a) "mu_of_r_potential = hf" uses the two-body rdm of a HF-like wave function (i.e. a single Slater determinant developped with the MOs stored in the EZFIO folder).
When HF is a qualitative representation of the electron pairs (i.e. weakly correlated systems), such an approach for \mu(r) is OK.
See for instance JPCL, 10, 2931-2937 (2019) for typical flavours of the results.
Thanks to the trivial nature of such a two-body rdm, the equation (22) of J. Chem. Phys. 149, 194301 (2018) can be rewritten in a very efficient way, and therefore the limiting factor of such an approach is the AO->MO four-index transformation of the two-electron integrals.
b) "mu_of_r_potential = cas_ful" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder).
If the CAS is properly chosen (i.e. the CAS-like wave function qualitatively represents the wave function of the systems), then such an approach is OK for \mu(r) even in the case of strong correlation.
+) The use of DFT correlation functionals with multi-determinant reference (Ecmd). These functionals are originally defined in the RS-DFT framework (see for instance Theor. Chem. Acc.114, 305(2005)) and design to capture short-range correlation effects. A important quantity arising in the Ecmd is the exact on-top pair density of the system, and the main differences of approximated Ecmd relies on different approximations for the exact on-top pair density.
The two main flavours of Ecmd depends on the strength of correlation in the system:
a) for weakly correlated systems, the ECMD PBE-UEG functional based on the seminal work of in RSDFT (see JCP, 150, 084103 1-10 (2019)) and adapted for the basis set correction in JPCL, 10, 2931-2937 (2019) uses the exact on-top pair density of the UEG at large mu and the PBE correlation functional at mu = 0. As shown in JPCL, 10, 2931-2937 (2019), such a functional is more accurate than the ECMD LDA for weakly correlated systems.
b) for strongly correlated systems, the ECMD PBE-OT, which uses the extrapolated on-top pair density of the CAS wave function thanks to the large \mu behaviour of the on-top pair density, is accurate, but suffers from S_z dependence (i.e. is not invariant with respect to S_z) because of the spin-polarization dependence of the PBE correlation functional entering at mu=0.
An alternative is ECMD SU-PBE-OT which uses the same on-top pair density that ECMD PBE-OT but a ZERO spin-polarization to remove the S_z dependence. As shown in ???????????, this strategy is one of the more accurate and respects S_z invariance and size consistency if the CAS wave function is correctly chosen.

View File

@ -0,0 +1 @@
change all correlation functionals with the pbe_on_top general

View File

@ -0,0 +1,27 @@
program basis_correction
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
read_wf = .True.
touch read_wf
no_core_density = .True.
touch no_core_density
provide mo_two_e_integrals_in_map
call print_basis_correction
! call print_e_b
end
subroutine print_e_b
implicit none
print *, 'Hello world'
print*,'ecmd_lda_mu_of_r = ',ecmd_lda_mu_of_r
print*,'ecmd_pbe_ueg_mu_of_r = ',ecmd_pbe_ueg_mu_of_r
print*,'ecmd_pbe_ueg_eff_xi_mu_of_r = ',ecmd_pbe_ueg_eff_xi_mu_of_r
print*,''
print*,'psi_energy + E^B_LDA = ',psi_energy + ecmd_lda_mu_of_r
print*,'psi_energy + E^B_PBE_UEG = ',psi_energy + ecmd_pbe_ueg_mu_of_r
print*,'psi_energy + E^B_PBE_UEG_Xi = ',psi_energy + ecmd_pbe_ueg_eff_xi_mu_of_r
print*,''
print*,'mu_average_prov = ',mu_average_prov
end

View File

@ -0,0 +1,92 @@
BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_eff_xi_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_pbe_ueg_eff_xi_mu_of_r = multi-determinantal Ecmd within the PBE-UEG and effective spin polarization approximation with mu(r),
!
! see Eqs. 30 in ???????????
!
! Based on the PBE-on-top functional (see Eqs. 26, 27 of J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.5082638)
!
! and replaces the approximation of the exact on-top pair density by the exact on-top of the UEG
!
! !!!! BUT !!!! with an EFFECTIVE SPIN POLARIZATION DEPENDING ON THE ON-TOP PAIR DENSITY
!
! See P. Perdew, A. Savin, and K. Burke, Phys. Rev. A 51, 4531 (1995). for original Ref., and Eq. 29 in ???????????
END_DOC
implicit none
double precision :: weight,density
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),g0_UEG_mu_inf,on_top
ecmd_pbe_ueg_eff_xi_mu_of_r = 0.d0
print*,'Providing ecmd_pbe_ueg_eff_xi_mu_of_r ...'
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight=final_weight_at_r_vector(ipoint)
mu = mu_of_r_prov(ipoint,istate)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! We use the effective spin density to define rho_a/rho_b
rho_a = 0.5d0 * (density + effective_spin_dm(ipoint,istate))
rho_b = 0.5d0 * (density - effective_spin_dm(ipoint,istate))
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
! We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
! with the effective rho_a and rho_b
on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_PBE)
ecmd_pbe_ueg_eff_xi_mu_of_r(istate) += eps_c_md_PBE * weight
enddo
enddo
double precision :: wall1, wall0
call wall_time(wall1)
print*,'Time for the ecmd_pbe_ueg_eff_xi_mu_of_r:',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_lda_eff_xi_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_lda_eff_xi_mu_of_r = multi-determinantal Ecmd within the LDA and effective spin polarization approximation with mu(r),
!
! corresponds to equation 40 in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714
!
! !!!! BUT !!!! with an EFFECTIVE SPIN POLARIZATION DEPENDING ON THE ON-TOP PAIR DENSITY
!
! See P. Perdew, A. Savin, and K. Burke, Phys. Rev. A 51, 4531 (1995). for original Ref., and Eq. 29 in ???????????
END_DOC
implicit none
integer :: ipoint,istate
double precision :: rho_a, rho_b, ec
logical :: dospin
double precision :: wall0,wall1,weight,mu,density
dospin = .true. ! JT dospin have to be set to true for open shell
print*,'Providing ecmd_lda_eff_xi_mu_of_r ...'
ecmd_lda_eff_xi_mu_of_r = 0.d0
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
mu = mu_of_r_prov(ipoint,istate)
weight = final_weight_at_r_vector(ipoint)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
rho_a = 0.5d0 * (density + effective_spin_dm(ipoint,istate))
rho_b = 0.5d0 * (density - effective_spin_dm(ipoint,istate))
call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
if(isnan(ec))then
print*,'ec is nan'
stop
endif
ecmd_lda_eff_xi_mu_of_r(istate) += weight * ec
enddo
enddo
call wall_time(wall1)
print*,'Time for ecmd_lda_eff_xi_mu_of_r :',wall1-wall0
END_PROVIDER

View File

@ -0,0 +1,121 @@
BEGIN_PROVIDER [double precision, ecmd_pbe_on_top_mu_of_r, (N_states)]
BEGIN_DOC
!
! Ecmd functional evaluated with mu(r) and depending on
! +) the on-top pair density
!
! +) the total density, density gradients
!
! +) the spin density
!
! Defined originally in Eq. (25) of JCP, 150, 084103 1-10 (2019) for RS-DFT calculations, but evaluated with mu(r).
!
! Such a functional is built by interpolating between two regimes :
!
! +) the large mu behaviour in cst/(\mu^3) \int dr on-top(r) where on-top(r) is supposed to be the exact on-top of the system
!
! +) mu= 0 with the usal ec_pbe(rho_a,rho_b,grad_rho_a,grad_rho_b)
!
! Here the approximation to the exact on-top is done through the assymptotic expansion (in \mu) of the exact on-top pair density (see Eq. 29) but with a mu(r) instead of a constant mu
!
! Such an asymptotic expansion was introduced in P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
!
END_DOC
implicit none
double precision :: weight
double precision :: eps_c_md_on_top_PBE,on_top_extrap,mu_correction_of_on_top
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top
ecmd_pbe_on_top_mu_of_r = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
mu = mu_of_r_prov(ipoint,istate)
! depends on (rho_a, rho_b) <==> (rho_tot,spin_pol)
rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
ecmd_pbe_on_top_mu_of_r(istate) += eps_c_md_on_top_PBE * weight
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_pbe_on_top_su_mu_of_r, (N_states)]
BEGIN_DOC
!
! Ecmd functional evaluated with mu(r) and depending on
! +) the on-top pair density
!
! +) the total density, density gradients
!
! +) !!!!! NO SPIN POLAIRIZATION !!!!!
!
! Defined originally in Eq. (25) of JCP, 150, 084103 1-10 (2019) for RS-DFT calculations, but evaluated with mu(r).
!
! Such a functional is built by interpolating between two regimes :
!
! +) the large mu behaviour in cst/(\mu^3) \int dr on-top(r) where on-top(r) is supposed to be the exact on-top of the system
!
! +) mu= 0 with the usal ec_pbe(rho_a,rho_b,grad_rho_a,grad_rho_b)
!
! Here the approximation to the exact on-top is done through the assymptotic expansion (in \mu) of the exact on-top pair density (see Eq. 29) but with a mu(r) instead of a constant mu
!
! Such an asymptotic expansion was introduced in P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
!
END_DOC
implicit none
double precision :: weight
double precision :: eps_c_md_on_top_PBE,on_top_extrap,mu_correction_of_on_top
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top,density
ecmd_pbe_on_top_su_mu_of_r = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
mu = mu_of_r_prov(ipoint,istate)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! rho_a = rho_b = rho_tot/2 ==> NO SPIN POLARIZATION
rho_a = 0.5d0 * density
rho_b = 0.5d0 * density
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
ecmd_pbe_on_top_su_mu_of_r(istate) += eps_c_md_on_top_PBE * weight
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,82 @@
subroutine print_basis_correction
implicit none
integer :: istate
provide mu_average_prov
if(mu_of_r_potential.EQ."hf")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
provide ecmd_pbe_on_top_mu_of_r ecmd_pbe_on_top_su_mu_of_r
endif
print*, ''
print*, ''
print*, '****************************************'
print*, '****************************************'
print*, 'Basis set correction for WFT using DFT Ecmd functionals'
print*, 'These functionals are accurate for short-range correlation'
print*, ''
print*, 'For more details look at Journal of Chemical Physics 149, 194301 1-15 (2018) '
print*, ' Journal of Physical Chemistry Letters 10, 2931-2937 (2019) '
print*, ' ???REF SC?'
print*, '****************************************'
print*, '****************************************'
print*, 'mu_of_r_potential = ',mu_of_r_potential
if(mu_of_r_potential.EQ."hf")then
print*, ''
print*,'Using a HF-like two-body density to define mu(r)'
print*,'This assumes that HF is a qualitative representation of the wave function '
print*,'********************************************'
print*,'Functionals more suited for weak correlation'
print*,'********************************************'
print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
enddo
print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo
else if(mu_of_r_potential.EQ."cas_ful")then
print*, ''
print*,'Using a CAS-like two-body density to define mu(r)'
print*,'This assumes that the CAS is a qualitative representation of the wave function '
print*,'********************************************'
print*,'Functionals more suited for weak correlation'
print*,'********************************************'
print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
enddo
print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo
print*,''
print*,'********************************************'
print*,'********************************************'
print*,'+) PBE-on-top Ecmd functional : (??????? REF-SCF ??????????)'
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
enddo
print*,''
print*,'********************************************'
print*,'+) PBE-on-top no spin polarization Ecmd functional : (??????? REF-SCF ??????????)'
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)
enddo
print*,''
endif
print*,''
print*,'**************'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) , state ',istate,' = ',mu_average_prov(istate)
enddo
end

View File

@ -0,0 +1,83 @@
BEGIN_PROVIDER [double precision, ecmd_lda_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_lda_mu_of_r = multi-determinantal Ecmd within the LDA approximation with mu(r) ,
!
! see equation 40 in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714
END_DOC
implicit none
integer :: ipoint,istate
double precision :: rho_a, rho_b, ec
logical :: dospin
double precision :: wall0,wall1,weight,mu
dospin = .true. ! JT dospin have to be set to true for open shell
print*,'Providing ecmd_lda_mu_of_r ...'
ecmd_lda_mu_of_r = 0.d0
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
mu = mu_of_r_prov(ipoint,istate)
weight = final_weight_at_r_vector(ipoint)
rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! Ecmd within the LDA approximation of PRB 73, 155111 (2006)
call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
if(isnan(ec))then
print*,'ec is nan'
stop
endif
ecmd_lda_mu_of_r(istate) += weight * ec
enddo
enddo
call wall_time(wall1)
print*,'Time for ecmd_lda_mu_of_r :',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_pbe_ueg_mu_of_r = multi-determinantal Ecmd within the PBE-UEG approximation with mu(r) ,
!
! see Eqs. 13-14b in Phys.Chem.Lett.2019, 10, 2931 2937; https://pubs.acs.org/doi/10.1021/acs.jpclett.9b01176
!
! Based on the PBE-on-top functional (see Eqs. 26, 27 of J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.5082638)
!
! but it the on-top pair density of the UEG as an approximation of the exact on-top pair density
END_DOC
implicit none
double precision :: weight
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top
double precision :: g0_UEG_mu_inf
ecmd_pbe_ueg_mu_of_r = 0.d0
print*,'Providing ecmd_pbe_ueg_mu_of_r ...'
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight=final_weight_at_r_vector(ipoint)
! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
mu = mu_of_r_prov(ipoint,istate)
rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
! We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)
! The form of interpolated (mu=0 ---> mu=infinity) functional originally introduced in JCP, 150, 084103 1-10 (2019)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_PBE)
ecmd_pbe_ueg_mu_of_r(istate) += eps_c_md_PBE * weight
enddo
enddo
double precision :: wall1, wall0
call wall_time(wall1)
print*,'Time for the ecmd_pbe_ueg_mu_of_r:',wall1-wall0
END_PROVIDER

1
src/ecmd_utils/NEED Normal file
View File

@ -0,0 +1 @@

View File

@ -0,0 +1,4 @@
==========
ecmd_utils
==========

View File

@ -0,0 +1,152 @@
!****************************************************************************
subroutine ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,e)
!*****************************************************************************
! Short-range spin-dependent LDA correlation functional with multideterminant reference
! for OEP calculations from Section V of
! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006)
!
! Input: rhot : total density
! rhos : spin density
! mu : Interation parameter
! dospin : use spin density
!
! Ouput: e : energy
!
! Created: 26-08-11, J. Toulouse
!*****************************************************************************
implicit none
double precision, intent(in) :: rho_a,rho_b,mu
logical, intent(in) :: dospin
double precision, intent(out):: e
double precision :: e1
double precision :: rhoa,rhob
double precision :: rhot, rhos
rhoa=max(rho_a,1.0d-15)
rhob=max(rho_b,1.0d-15)
rhot = rhoa + rhob
rhos = rhoa - rhob
call ec_only_lda_sr(mu,rho_a,rho_b,e1)
if(isnan(e1))then
print*,'e1 is NaN'
print*,mu,rho_a,rho_b
stop
endif
call DELTA_LRSR_LDAERF (rhot,rhos,mu,dospin,e)
if(isnan(e))then
print*,'e is NaN'
print*,mu,rhot,rhos
stop
endif
e = e1 + e
end
!****************************************************************************
subroutine DELTA_LRSR_LDAERF (rhot,rhos,mu,dospin,e)
!*****************************************************************************
! LDA approximation to term Delta_(LR-SR) from Eq. 42 of
! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006)
!
! Input: rhot : total density
! rhos : spin density
! mu : Interation parameter
! dospin : use spin density
!
! Ouput: e : energy
!
! Warning: not tested for z != 0
!
! Created: 26-08-11, J. Toulouse
!*****************************************************************************
implicit none
double precision rhot, rhos, mu
logical dospin
double precision e
double precision f13, f83, pi, rsfac, alpha2
double precision rs, rs2, rs3
double precision rhoa, rhob, z, z2, onepz, onemz, zp, zm, phi8
double precision g0f, g0s
double precision bd2, bd3
double precision c45, c4, c5
double precision bc2, bc4, bc3t, bc5t, d0
double precision delta2,delta3,delta4,delta5,delta6
double precision delta
parameter(f13 = 0.333333333333333d0)
parameter(f83 = 2.6666666666666665d0)
parameter(pi = 3.141592653589793d0)
parameter(rsfac = 0.620350490899400d0)
parameter(alpha2 = 0.2715053589826032d0)
rs = rsfac/(rhot**f13)
rs2 = rs*rs
rs3 = rs2*rs
! Spin-unpolarized case
if (.not.dospin) then
z = 0.d0
! Spin-polarized case
else
rhoa=max((rhot+rhos)*.5d0,1.0d-15)
rhob=max((rhot-rhos)*.5d0,1.0d-15)
z=min((rhoa-rhob)/(rhoa+rhob),0.9999999999d0)
endif
z2=z*z
bd2=dexp(-0.547d0*rs)*(-0.388d0*rs+0.676*rs2)/rs2
bd3=dexp(-0.31d0*rs)*(-4.95d0*rs+rs2)/rs3
onepz=1.d0+z
onemz=1.d0-z
phi8=0.5d0*(onepz**f83+onemz**f83)
zp=onepz/2.d0
zm=onemz/2.d0
c45=(zp**2)*g0s(rs*zp**(-f13))+(zm**2)*g0s(rs*zm**(-f13))
c4=c45+(1.d0-z2)*bd2-phi8/(5.d0*alpha2*rs2)
c5=c45+(1.d0-z2)*bd3
bc2=-3.d0*(1-z2)*(g0f(rs)-0.5d0)/(8.d0*rs3)
bc4=-9.d0*c4/(64.d0*rs3)
bc3t=-(1-z2)*g0f(rs)*(2.d0*dsqrt(2.d0)-1)/(2.d0*dsqrt(pi)*rs3)
bc5t = -3.d0*c5*(3.d0-dsqrt(2.d0))/(20.d0*dsqrt(2.d0*pi)*rs3)
d0=(0.70605d0+0.12927d0*z2)*rs
delta2=0.073867d0*(rs**(1.5d0))
delta3=4*(d0**6)*bc3t+(d0**8)*bc5t;
delta4=4*(d0**6)*bc2+(d0**8)*bc4;
delta5=(d0**8)*bc3t;
delta6=(d0**8)*bc2;
delta=(delta2*(mu**2)+delta3*(mu**3)+delta4*(mu**4)+delta5*(mu**5)+delta6*(mu**6))/((1+(d0**2)*(mu**2))**4)
! multiply by rhot to get energy density
e=delta*rhot
end
!*****************************************************************************
double precision function g0s(rs)
!*****************************************************************************
! g"(0,rs,z=1) from Eq. 32 of
! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006)
!
! Created: 26-08-11, J. Toulouse
!*****************************************************************************
implicit none
double precision rs
double precision rs2, f53, alpha2
parameter(f53 = 1.6666666666666667d0)
parameter(alpha2 = 0.2715053589826032d0)
rs2=rs*rs
g0s=(2.d0**f53)*(1.d0-0.02267d0*rs)/((5.d0*alpha2*rs2)*(1.d0+0.4319d0*rs+0.04d0*rs2))
end

View File

@ -0,0 +1,53 @@
subroutine ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_on_top_PBE)
implicit none
BEGIN_DOC
!
! General e_cmd functional interpolating between :
!
! +) the large mu behaviour in cst/(\mu^3) on-top
!
! +) mu= 0 with the usal ec_pbe at
!
! Depends on : mu, the density (rho_a,rho_b), the square of the density gradient (grad_rho_a,grad_rho_b)
!
! the flavour of on-top densiyt (on_top) you fill in: in principle it should be the exact on-top
!
! The form of the functional was originally introduced in JCP, 150, 084103 1-10 (2019)
!
END_DOC
double precision, intent(in) :: mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top
double precision, intent(out) :: eps_c_md_on_top_PBE
double precision :: pi, e_pbe,beta,denom
double precision :: grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
integer :: m
pi = 4.d0 * datan(1.d0)
eps_c_md_on_top_PBE = 0.d0
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do m = 1, 3
grad_rho_a_2 += grad_rho_a(m)*grad_rho_a(m)
grad_rho_b_2 += grad_rho_b(m)*grad_rho_b(m)
grad_rho_a_b += grad_rho_a(m)*grad_rho_b(m)
enddo
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,sigmaoo,sigmacc,sigmaco)
! usual PBE correlation energy using the density, spin polarization and density gradients for alpha/beta electrons
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
denom = (-2.d0+sqrt(2d0))*sqrt(2.d0*pi)* on_top
if (dabs(denom) > 1.d-12) then
! quantity of Eq. (26)
beta = (3.d0*e_PBE)/denom
eps_c_md_on_top_PBE = e_PBE/(1.d0+beta*mu**3)
else
eps_c_md_on_top_PBE =0.d0
endif
end

View File

@ -0,0 +1,72 @@
subroutine ec_md_on_top_PBE_mu_corrected(mu,r,two_dm,eps_c_md_on_top_PBE)
implicit none
BEGIN_DOC
! enter with "r(3)", and "two_dm(N_states)" which is the on-top pair density at "r" for each states
!
! you get out with the energy density defined in J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.508263
!
! by Eq. (26), which includes the correction of the on-top pair density of Eq. (29).
END_DOC
double precision, intent(in) :: mu , r(3), two_dm
double precision, intent(out) :: eps_c_md_on_top_PBE(N_states)
double precision :: two_dm_in_r, pi, e_pbe(N_states),beta(N_states),mu_correction_of_on_top
double precision :: aos_array(ao_num), grad_aos_array(3,ao_num)
double precision :: rho_a(N_states),rho_b(N_states)
double precision :: grad_rho_a(3,N_states),grad_rho_b(3,N_states)
double precision :: grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo,on_top_corrected
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_on_top_PBE = 0.d0
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,rho_a,rho_b, grad_rho_a, grad_rho_b, aos_array, grad_aos_array)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do istate = 1, N_states
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate)*grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate)*grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate)*grad_rho_b(m,istate)
enddo
enddo
do istate = 1, N_states
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
! usual PBE correlation energy using the density, spin polarization and density gradients for alpha/beta electrons
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE(istate))
! correction of the on-top pair density according to Eq. (29)
on_top_corrected = mu_correction_of_on_top(mu,two_dm)
! quantity of Eq. (27) with a factor two according to the difference of normalization
! between the on-top of the JCP paper and that of QP2
beta(istate) = (3.d0*e_PBE(istate))/( (-2.d0+sqrt(2d0))*sqrt(2.d0*pi)*2.d0* on_top_corrected)
! quantity of Eq. (26)
eps_c_md_on_top_PBE(istate)=e_PBE(istate)/(1.d0+beta(istate)*mu**3)
enddo
end
double precision function mu_correction_of_on_top(mu,on_top)
implicit none
BEGIN_DOC
! mu-based correction to the on-top pair density provided by the assymptotic expansion of
!
! P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
!
! This is used in J. Chem. Phys.150, 084103 (2019); Eq. (29).
END_DOC
double precision, intent(in) :: mu,on_top
double precision :: pi
pi = 4.d0 * datan(1.d0)
mu_correction_of_on_top = on_top / ( 1.d0 + 2.d0/(dsqrt(pi)*mu) )
mu_correction_of_on_top = max(mu_correction_of_on_top ,1.d-15)
end

View File

@ -0,0 +1,194 @@
subroutine ecmd_pbe_ueg_at_r(mu,r,eps_c_md_PBE)
implicit none
BEGIN_DOC
! provides the integrand of Eq. (13) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! !!! WARNING !!! This is the total integrand of Eq. (13), which is e_cmd * n
!
! such a function is based on the exact behaviour of the Ecmd at large mu
!
! but with the exact on-top estimated with that of the UEG
!
! You enter with r(3), you get out with eps_c_md_PBE(1:N_states)
END_DOC
double precision, intent(in) :: mu , r(3)
double precision, intent(out) :: eps_c_md_PBE(N_states)
double precision :: pi, e_PBE, beta
double precision :: aos_array(ao_num), grad_aos_array(3,ao_num)
double precision :: rho_a(N_states),rho_b(N_states)
double precision :: grad_rho_a(3,N_states),grad_rho_b(3,N_states)
double precision :: grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
double precision :: g0_UEG_mu_inf, denom
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_PBE = 0.d0
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,rho_a,rho_b, grad_rho_a, grad_rho_b, aos_array, grad_aos_array)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do istate = 1, N_states
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate)*grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate)*grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate)*grad_rho_b(m,istate)
enddo
enddo
do istate = 1, N_states
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
if(mu == 0.d0) then
eps_c_md_PBE(istate)=e_PBE
else
! note: the on-top pair density is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
denom = (-2.d0+sqrt(2d0))*sqrt(2.d0*pi) * 4.d0*rho_a(istate)*rho_b(istate)*g0_UEG_mu_inf(rho_a(istate),rho_b(istate))
if (dabs(denom) > 1.d-12) then
beta = (3.d0*e_PBE)/denom
eps_c_md_PBE(istate)=e_PBE/(1.d0+beta*mu**3)
else
eps_c_md_PBE(istate)=0.d0
endif
endif
enddo
end
subroutine eps_c_md_PBE_from_density(mu,rho_a,rho_b, grad_rho_a, grad_rho_b,eps_c_md_PBE) ! EG
implicit none
BEGIN_DOC
! provides the integrand of Eq. (13) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! !!! WARNING !!! This is the total integrand of Eq. (13), which is e_cmd * n
!
! such a function is based on the exact behaviour of the Ecmd at large mu
!
! but with the exact on-top estimated with that of the UEG
!
! You enter with the alpha/beta density and density gradients
!
! You get out with eps_c_md_PBE(1:N_states)
END_DOC
double precision, intent(in) :: mu(N_states) , rho_a(N_states),rho_b(N_states), grad_rho_a(3,N_states),grad_rho_b(3,N_states)
double precision, intent(out) :: eps_c_md_PBE(N_states)
double precision :: pi, e_PBE, beta
double precision :: aos_array(ao_num), grad_aos_array(3,ao_num)
double precision :: grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
double precision :: g0_UEG_mu_inf, denom
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_PBE = 0.d0
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do istate = 1, N_states
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate)*grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate)*grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate)*grad_rho_b(m,istate)
enddo
enddo
do istate = 1, N_states
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
if(mu(istate) == 0.d0) then
eps_c_md_PBE(istate)=e_PBE
else
! note: the on-top pair density is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
denom = (-2.d0+sqrt(2d0))*sqrt(2.d0*pi) * 4.d0*rho_a(istate)*rho_b(istate)*g0_UEG_mu_inf(rho_a(istate),rho_b(istate))
if (dabs(denom) > 1.d-12) then
beta = (3.d0*e_PBE)/denom
eps_c_md_PBE(istate)=e_PBE/(1.d0+beta*mu(istate)**3)
else
eps_c_md_PBE(istate)=0.d0
endif
endif
enddo
end
subroutine eps_c_md_PBE_at_grid_pt(mu,i_point,eps_c_md_PBE)
implicit none
BEGIN_DOC
! provides the integrand of Eq. (13) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! !!! WARNING !!! This is the total integrand of Eq. (13), which is e_cmd * n
!
! such a function is based on the exact behaviour of the Ecmd at large mu
!
! but with the exact on-top estimated with that of the UEG
!
! You enter with the alpha/beta density and density gradients
!
! You get out with eps_c_md_PBE(1:N_states)
END_DOC
double precision, intent(in) :: mu
double precision, intent(out) :: eps_c_md_PBE(N_states)
integer, intent(in) :: i_point
double precision :: two_dm, pi, e_pbe,beta,mu_correction_of_on_top
double precision :: grad_rho_a(3),grad_rho_b(3)
double precision :: grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
double precision :: rhoc,rhoo,ec_pbe_88
double precision :: delta,two_dm_corr,rho_a,rho_b
double precision :: grad_rho_2,denom,g0_UEG_mu_inf
double precision :: sigmacc,sigmaco,sigmaoo
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_PBE = 0.d0
do istate = 1, N_states
! total and spin density
rhoc = one_e_dm_and_grad_alpha_in_r(4,i_point,istate) + one_e_dm_and_grad_beta_in_r(4,i_point,istate)
rhoo = one_e_dm_and_grad_alpha_in_r(4,i_point,istate) - one_e_dm_and_grad_beta_in_r(4,i_point,istate)
! gradients of the effective spin density
grad_rho_a_2 = 0.D0
grad_rho_b_2 = 0.D0
grad_rho_a_b = 0.D0
do m = 1, 3
grad_rho_a_2 += one_e_dm_and_grad_alpha_in_r(m,i_point,istate)**2.d0
grad_rho_b_2 += one_e_dm_and_grad_beta_in_r(m,i_point,istate) **2.d0
grad_rho_a_b += one_e_dm_and_grad_alpha_in_r(m,i_point,istate) * one_e_dm_and_grad_beta_in_r(m,i_point,istate)
enddo
sigmacc = grad_rho_a_2 + grad_rho_b_2 + 2.d0 * grad_rho_a_b
sigmaco = 0.d0
sigmaoo = 0.d0
rho_a = one_e_dm_and_grad_alpha_in_r(4,i_point,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,i_point,istate)
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
if(e_PBE.gt.0.d0)then
print*,'PBE gt 0 with regular dens'
endif
if(mu == 0.d0) then
eps_c_md_PBE(istate)=e_PBE
else
! note: the on-top pair density is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
denom = (-2.d0+dsqrt(2d0))*sqrt(2.d0*pi) * 4.d0*rho_a*rho_b*g0_UEG_mu_inf(rho_a,rho_b)
if (dabs(denom) > 1.d-12) then
beta = (3.d0*e_PBE)/denom
! Ecmd functional with the UEG ontop pair density when mu -> infty
! and the usual PBE correlation energy when mu = 0
eps_c_md_PBE(istate)=e_PBE/(1.d0+beta*mu**3)
else
eps_c_md_PBE(istate)=0.d0
endif
endif
enddo
end

View File

@ -0,0 +1,94 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
double precision function correction_to_on_top_from_UEG(mu,r,istate)
implicit none
integer, intent(in) :: istate
double precision, intent(in) :: mu,r(3)
double precision :: rho_a(N_states),rho_b(N_states)
double precision :: g0_UEG_mu_inf, g0_UEG_mu
call dm_dft_alpha_beta_at_r(r,rho_a,rho_b)
correction_to_on_top_from_UEG = g0_UEG_mu_inf(rho_a(istate),rho_b(istate)) / g0_UEG_mu(mu,rho_a(istate),rho_b(istate))
end
double precision function g0_UEG_mu_inf(rho_a,rho_b)
BEGIN_DOC
! Pair distribution function g0(n_alpha,n_beta) of the Colombic UEG
!
! Taken from Eq. (46) P. Gori-Giorgi and A. Savin, Phys. Rev. A 73, 032506 (2006).
END_DOC
implicit none
double precision, intent(in) :: rho_a,rho_b
double precision :: rho,pi,x
double precision :: B, C, D, E, d2, rs, ahd
rho = rho_a+rho_b
pi = 4d0 * datan(1d0)
ahd = -0.36583d0
d2 = 0.7524d0
B = -2d0 * ahd - d2
C = 0.08193d0
D = -0.01277d0
E = 0.001859d0
if (dabs(rho) > 1.d-12) then
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
x = -d2*rs
g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*exp(x)
else
g0_UEG_mu_inf= 0.d0
endif
end
double precision function g0_UEG_mu(mu,rho_a,rho_b)
implicit none
BEGIN_DOC
! Pair distribution function g0(n_alpha,n_beta) of the UEG interacting with the long range interaction erf(mu r12)/r12
!
! Taken from P. Gori-Giorgi and A. Savin, Phys. Rev. A 73, 032506 (2006).
END_DOC
double precision, intent(in) :: rho_a,rho_b,mu
double precision :: zeta,pi,rho,x,alpha
double precision :: B, C, D, E, d2, rs, ahd, h_func, kf
pi = 4d0 * datan(1d0)
rho = rho_a+rho_b
alpha = (4d0/(9d0*pi))**(1d0/3d0)
ahd = -0.36583d0
d2 = 0.7524d0
B = -2d0 * ahd - d2
C = 0.08193d0
D = -0.01277d0
E = 0.001859d0
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
kf = (alpha*rs)**(-1d0)
zeta = mu / kf
x = -d2*rs*h_func(zeta)/ahd
g0_UEG_mu = (exp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) )
end
double precision function h_func(zeta)
implicit none
double precision, intent(in) :: zeta
double precision :: pi
double precision :: a1, a2, b1, b2, b3, ahd, alpha
pi = 4d0 * datan(1d0)
ahd = -0.36583d0
alpha = (4d0/(9d0*pi))**(1d0/3d0)
a1 = -(6d0*alpha/pi)*(1d0-log(2d0))
b1 = 1.4919d0
b3 = 1.91528d0
a2 = ahd * b3
b2 = (a1 - (b3*alpha/sqrt(pi)))/ahd
h_func = (a1*zeta**2d0 + a2*zeta**3d0) / (1d0 + b1*zeta + b2*zeta**2d0 + b3*zeta**3d0)
end