mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-31 05:03:40 +02:00
54 lines
1.9 KiB
Fortran
54 lines
1.9 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
|
||
|
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
|
||
|
|
||
|
|