from sympy import pi, var, S, init_printing, log, Function, sqrt
init_printing(use_latex=True)
The basic relation for any LDA potential is the contribution to the total/free energy:
$$ F_{xc}[n_e] = \int n_e(\mathbf{x}) \epsilon_{xc}^{LD}(n_e) d^3 x\,. $$From this we calculate the XC potential $V_{xc}(\mathbf{x})$ as a functional derivative (this follows from minimizing the total/free energy): $$ V_{xc}(\mathbf{x}) \equiv {\delta F_{xc}[n_e] \over \delta n_e} = \epsilon_{xc}^{LD}(n_e) + n_e(\mathbf{x}) {d \epsilon_{xc}^{LD}(n_e) \over d n_e} = {d \over d n_e}\left(n_e \epsilon_{xc}^{LD}(n_e)\right) $$ The $n_e(\mathbf{x})$ is the (positive) electronic particle density, $\epsilon_{xc}^{LD}(n_e)$ is the XC energy density and $V_{xc}(\mathbf{x})$ is the XC potential.
Sometimes the potential can be decomposed into exchange and correlation parts as follows: $$\epsilon_{xc}^{LD}(n_e)= \epsilon_{x}^{LD}(n_e) + \epsilon_{c}^{LD}(n_e)$$
var("n")
ex = -3/(4*pi) * (3*pi**2*n)**(S(1)/3)
The corresponding $V_x$ is then calculated as:
Vx = (n*ex).diff(n)
Vx
Finally, we notice:
Vx/ex
So we have: $$V_{x}(\mathbf{x}) = {4\over 3}\epsilon_{x}^{LD}(n_e)$$ So this show that the exchange part can be calculated using the following Fortran statements:
ex = -3/(4*pi) * (3*pi**2*n)**(1.0_dp/3)
Vx = 4*ex/3
The following is taken from [1]. The constants are:
$$ \begin{eqnarray*} \gamma &=& -0.1423 \\ \beta_1 &=& 1.0529 \\ \beta_2 &=& 0.3334 \\ A &=& 0.0311 \\ B &=& -0.048 \\ C &=& 0.0020 \\ D &=& -0.0116 \\ \end{eqnarray*} $$[1] Perdew, J. P., & Zunger, A. (1981). Self-interaction correction to density-functional approximations for many-electron systems. Physical Review B, 23(10), 5048–5079.
rs = (3/(4*pi*n))**(S(1)/3)
rs
var("A B C D")
ec = A*log(rs) + B + C*rs*log(rs) + D*rs
ec
Vc = (n*ec).diff(n).expand()
Vc
Which is a little complicated, so we try the following approach:
rs.diff(n) / rs
So ${d r_s \over d n} = - {1\over 3n} r_s$.
rs = Function("rs")(n)
ec = A*log(rs) + B + C*rs*log(rs) + D*rs
ec
Vc = (n*ec).diff(n).expand()
Vc
Now we substitute for the derivative:
Vc = Vc.subs(rs.diff(n), (-1/(3*n))*rs)
Vc
So this part of $V_c$ can be calculated using the following Fortran statements:
rs = (3/(4*pi*n))**(1.0_dp/3)
ec = A*log(rs) + B + C*rs*log(rs) + D*rs
Vc = A*log(rs) + (B-A/3) + 2*C*rs*log(rs)/3 + (2*D-C)*rs/3
var("gamma beta1 beta2")
ec = gamma / (1+beta1*sqrt(rs)+beta2*rs)
ec
Vc = (n*ec).diff(n)
Vc
Vc = Vc.subs(rs.diff(n), (-1/(3*n))*rs)
Vc
We express $V_c$ as $V_c = e_c \cdot C$. Arguably the simplest expression for $C$ is:
C = (Vc / ec).simplify()
C
But in literature the following is used (both numerator and denumerator must be divided by 6 below):
C.normal()
So the following Fortran statements implement that:
ec = gam / (1+beta1*sqrt(rs)+beta2*rs)
Vc = ec * (1+7*beta1*sqrt(rs)/6 + 4*beta2*rs/3) / (1+beta1*sqrt(rs) + beta2*rs)
The final Fortran subroutine is:
subroutine xc_pz(n, exc, Vxc)
! Calculates XC LDA density and potential from the charge density "n".
! Uses the Perdew Zunger [1] parametrization.
!
! [1] Perdew, J. P., & Zunger, A. (1981). Self-interaction correction to
! density-functional approximations for many-electron systems. Physical Review
! B, 23(10), 5048–5079.
real(dp), intent(in) :: n ! charge density
real(dp), intent(out) :: exc ! XC density
real(dp), intent(out) :: Vxc ! XC potential
real(dp), parameter :: gam = -0.1423_dp
real(dp), parameter :: beta1 = 1.0529_dp
real(dp), parameter :: beta2 = 0.3334_dp
real(dp), parameter :: A = 0.0311_dp
real(dp), parameter :: B = -0.048_dp
real(dp), parameter :: C = 0.0020_dp
real(dp), parameter :: D = -0.0116_dp
real(dp) :: ex, ec, Vx, Vc, rs
if (n == 0) then
exc = 0
Vxc = 0
return
end if
ex = -3/(4*pi) * (3*pi**2*n)**(1.0_dp/3)
Vx = 4*ex/3
rs = (3/(4*pi*n))**(1.0_dp/3)
if (rs >= 1) then
ec = gam / (1+beta1*sqrt(rs)+beta2*rs)
Vc = ec * (1+7*beta1*sqrt(rs)/6 + 4*beta2*rs/3) / &
(1+beta1*sqrt(rs) + beta2*rs)
else
ec = A*log(rs) + B + C*rs*log(rs) + D*rs
Vc = A*log(rs) + (B-A/3) + 2*C*rs*log(rs)/3 + (2*D-C)*rs/3
end if
exc = ex + ec
Vxc = Vx + Vc
end subroutine