9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00

added pbe_ueg_self_contained.irp.f

This commit is contained in:
eginer 2022-06-15 19:39:01 +02:00
parent 900e010395
commit e99e976c22
4 changed files with 115 additions and 1 deletions

View File

@ -0,0 +1,62 @@
double precision function ecmd_pbe_ueg_self_cont(dens,spin_pol,mu,e_PBE)
implicit none
! dens = total density
! spin_pol = spin_polarization (n_a - n_b)/dens
! e_PBE = PBE correlation (mu=0) energy evaluated at (dens,spin_pol,grad_rho)
! e_PBE = epsilon_PBE * dens
! dens = a + b
! spin_pol = (a - b)/(a+b)
! spin_pol * dens = a - b
! a - b + a+b = 2 a
! a - b - a - b = - 2b
double precision, intent(in) :: dens,spin_pol,mu,e_PBE
double precision :: rho_a,rho_b,pi,g0_UEG_func,denom,beta
pi = dacos(-1.d0)
rho_a = (dens * spin_pol + dens)*0.5d0
rho_b = (dens - dens * spin_pol)*0.5d0
if(mu == 0.d0) then
ecmd_pbe_ueg_self_cont = 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*rho_b*g0_UEG_func(rho_a,rho_b)
if (dabs(denom) > 1.d-12) then
beta = (3.d0*e_PBE)/denom
ecmd_pbe_ueg_self_cont=e_PBE/(1.d0+beta*mu**3)
else
ecmd_pbe_ueg_self_cont=0.d0
endif
endif
end
double precision function g0_UEG_func(rho_a,rho_b)
! 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).
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
x = -d2*rs
if (dabs(rho) > 1.d-20) then
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
x = -d2*rs
if(dabs(x).lt.50.d0)then
g0_UEG_func= 0.5d0 * (1d0+ rs* (-B + rs*(C + rs*(D + rs*E))))*dexp(x)
else
g0_UEG_func= 0.d0
endif
else
g0_UEG_func= 0.d0
endif
g0_UEG_func = max(g0_UEG_func,1.d-14)
end

View File

@ -20,6 +20,7 @@ subroutine print_su_pbe_ot
integer :: istate
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)
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ecmd_pbe_ueg_test , state ',istate,' = ',ecmd_pbe_ueg_test(istate)
enddo
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)

View File

@ -81,3 +81,54 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
print*,'Time for the ecmd_pbe_ueg_mu_of_r:',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_test, (N_states)]
BEGIN_DOC
! test of the routines contained in pbe_ueg_self_contained.irp.f
END_DOC
implicit none
double precision :: weight
integer :: ipoint,istate,m
double precision :: mu,rho_a,rho_b
double precision :: dens,spin_pol,grad_rho,e_PBE,delta_rho
double precision :: ecmd_pbe_ueg_self_cont,eps_c_md_PBE
ecmd_pbe_ueg_test = 0.d0
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)
! conversion from rho_a,rho_b --> dens,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)
dens = rho_a + rho_b
spin_pol = (rho_a - rho_b)/(max(dens,1.d-12))
delta_rho = rho_a - rho_b
! conversion from grad_rho_a ... to sigma
double precision :: grad_rho_a(3),grad_rho_b(3),grad_rho_a_2(3),grad_rho_b_2(3),grad_rho_a_b(3)
double precision :: sigmacc,sigmaco,sigmaoo
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
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
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,sigmaoo,sigmacc,sigmaco)
! call the PBE energy
call ec_pbe_only(0.d0,dens,delta_rho,sigmacc,sigmaco,sigmaoo,e_PBE)
eps_c_md_PBE = ecmd_pbe_ueg_self_cont(dens,spin_pol,mu,e_PBE)
ecmd_pbe_ueg_test(istate) += eps_c_md_PBE * weight
enddo
enddo
!
END_PROVIDER