```
from sympy import pi, var, S, init_printing, log, Function, sqrt
init_printing(use_latex=True)
```

# Exchange and Correlation Potential

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)$$

## Exchange part

$$\epsilon_{x}^{LD}(n_e) = -{3\over 4\pi} (3 \pi^2 n_e)^{1\over 3}$$


```
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

## Correlation part of Perdew and Zunger (PZ)

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
```

### Case $r_s < 1$


```
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

### Case $r_s \ge 1$


```
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)

## Fortran Subroutine

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

502072cf9818119602df3838b1115c8a244a7c5b
ipynb
Jupyter Notebook
src/tests/fem/plots/Vxc derivation.ipynb
certik/hfsolver
b4c50c1979fb7e468b1852b144ba756f5a51788d
BSD-2-Clause
2015-03-24T13:06:39.000Z
2022-03-29T00:14:02.000Z
src/tests/fem/plots/Vxc derivation.ipynb
certik/hfsolver
b4c50c1979fb7e468b1852b144ba756f5a51788d
BSD-2-Clause
2015-03-25T04:59:43.000Z
2017-06-06T23:00:09.000Z
src/tests/fem/plots/Vxc derivation.ipynb
certik/hfsolver
b4c50c1979fb7e468b1852b144ba756f5a51788d
BSD-2-Clause
2016-01-20T13:38:22.000Z
2020-11-24T15:35:43.000Z
Qwen/Qwen-72B
1. YES
2. YES