How can the GGA functionals in libxc be evaluated?

Question

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?

Answer 1

The xc energy is given by

$$ E_\text{xc} = \int n({\bf r}) \epsilon_\text{xc} ({\bf r}) {\rm d}^3r \equiv \int f_\text{xc} ({\bf r}) {\rm d}^3r $$

where $\epsilon_\text{xc}$ is the energy density per particle and $f_\text{xc}$ is the energy density. What Libxc gives you is $\epsilon_\text{xc}$ and derivatives of $f_\text{xc}$, such as $\partial f_\text{xc} / \partial n_\sigma$, $\partial f_\text{xc} / \partial \gamma_{\sigma \sigma'}$, $\partial f_\text{xc} / \partial \nabla^2 n_\sigma$, and $\partial f_\text{xc} / \partial \tau_\sigma$ that you need for e.g. self-consistent field calculations.

How you actually use these variables depends on the numerical approach. The approach you have described above is the one based on a real-space potential $v({\bf r})$, which is simple for the LDA but gets a bit hairy already for GGAs since you need to evaluate the divergence of a term that contains both the density and the density functional.

If you have a basis set, you can use integration by parts to get rid of the divergences and end up with a straightforward quadrature problem as discussed e.g. in our recent open access review article: Molecules 25, 1218 (2020). This also works nicely with numerical basis sets, as I've demonstrated e.g. in Int. J. Quantum Chem. 119, e25945 (2019).

If you have a basis set, you can also evaluate the real-space potential $v({\bf r})$ in an analytical fashion. This leads to higher-order derivatives of $f_\text{xc}$ which are also provided by Libxc as well as derivatives of the density; I have demonstrated this for atoms in Phys. Rev. A 101, 012516 (2020).

In contrast, as far as I know, programs employing the finite difference approach like dftatom evaluate the GGA contribution to the potential by finite differences.