61 lines
2.1 KiB
Fortran
61 lines
2.1 KiB
Fortran
|
|
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
|
|
! convertion from (alpha,beta) formalism to (closed, open) formalism for the density
|
|
call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
|
|
if(rhoc.lt.1.d-10)then
|
|
return
|
|
else if(on_top/(rhoc**2) .lt. 1.d-6)then
|
|
return
|
|
endif
|
|
|
|
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
|
|
! same same for gradients : convertion from (alpha,beta) formalism to (closed, open) formalism
|
|
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
|
|
|
|
|