For pedagogical reasons, I started adding libxc to the atomic DFT dftatom. In this program only spherical symmetric charges. That means that the final potential is just radial.
To my understanding, if the LDA xc energy is defined as
$$ E_{xc}[n(r)] = \int n(\mathbf{r}) \varepsilon_{xc}[n(\mathbf{r})]\ d\mathbf{r}, $$ the radial xc potential contribution will be $$ V_{xc}[n(r)] = \varepsilon_{xc}(r) + n(r)\,\frac{d \varepsilon_{xc}(r)}{d n(r)} $$
According to the manual of libxc from the function with the interface
void xc_lda_exc_vxc(const xc_lda_type *p, inp np, double *rho,
double *exc, double *vrho);
the $\varepsilon_{xc}$ from exc
and $\frac{d \varepsilon_{xc}}{d n}$ from vrho
. Judging from comparison to simple analytical expressions (e.g. Slater exchange), I would say the vrho
actually returns the whole radial expression $\varepsilon_{xc} + n\frac{d \varepsilon_{xc}}{d n}$ already.
I am starting to wonder is the gradient correction term for GGA returned from the interface. Generally, GGA potential is written as
$$ V_{xc}[n] = \varepsilon_{xc} + n\frac{\partial \varepsilon_{xc}}{\partial n} - \nabla\cdot\left[n\frac{\partial \varepsilon_{xc}}{\partial \nabla n}\right]. $$
and with the contracted gradient $\sigma$
$$ V_{xc}[n] = \varepsilon_{xc} + n\frac{\partial \varepsilon_{xc}}{\partial n} - 2 \, \nabla\cdot\left[n\frac{\partial \varepsilon_{xc}}{\partial \sigma} \nabla n\right]. $$
The interface for GGA is
void xc_gga_exc_vxc(const xc_func_type *p, int np, double *rho, double *sigma,
double *exc, double *vrho, double *vsigma)
Also here I would say that vrho
returns the sum of the first two terms and vsigma
is $n\frac{\partial \varepsilon_{xc}}{\partial \sigma}$.
Can somebody verify the correct implementation for GGA and LDA from the libxc functions?