How does libxc calculate the potential of GGA functionals?

Question

If the GGA xc energy is defined as $$ E_{xc}[n] = \int n(\mathbf{r}) \varepsilon_{xc}[n(\mathbf{r}), \nabla n(\mathbf{r})]\ d\mathbf{r}, $$ the potential can be 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]. $$ The third term (if I understand it correctly) requires the evaluation of a divergence term $\left(\nabla\cdot\right)$ which (I think) needs the spatial information of the points in the grid.

However, from LibXC manual, to calculate the GGA potential the user only needs to supply the density $n$ and the density gradients $\nabla n$, without any further information on grid points.

My question is: how to calculate the divergence term above if only $n$ and $\nabla n$ are given? Or do I miss something?

Answer 1

Libxc does not compute the full potential. Instead, it provides necessary ingredients for the code that calls Libxc to compute the potential. This strategy allows the library to work with a variety of codes employing different grid types.

For example, a code for single atoms might compute the divergence term on a radial grid with finite differences. A plane wave code could compute the divergence term with Fourier Transforms. There are many possibilities, but the design of Libxc allows it to work in all of these circumstances. The only downside is that the client code has to do a bit of the work itself.

---

For practical purposes—assuming the simplest case without spin polarization—it helps to recognize that Libxc recasts GGAs into the form $$ E_{xc}[n] = \int n \, \epsilon(n, \sigma) \, d\mathbf{r}, $$ where $\sigma=|\nabla n|^2$. The associated potential (valid for most boundary conditions) is then $$ v_{xc} = \epsilon + n \frac{\partial \epsilon}{\partial n} - 2 \nabla \cdot \left(n \frac{\partial \epsilon}{\partial \sigma} \nabla n\right) . $$ To compute the real space potential using Libxc:

1. your code provides $n$ and $\sigma$ at points in space;
2. Libxc returns $\epsilon$ and the derivatives $\partial \epsilon / \partial n$ and $\partial \epsilon / \partial \sigma$;
3. your code assembles the potential, computing the divergence term in some suitable manner.

---

There are a few other strategies worth knowing about too. For example, another form for the real space potential is $$ v_{xc} = \epsilon + n \frac{\partial \epsilon}{\partial n} - 2 \left(\frac{\partial \epsilon}{\partial \sigma} + n \frac{\partial^2 \epsilon}{\partial n \partial \sigma}\right) \sigma - 2 n \frac{\partial^2 \epsilon}{\partial \sigma^2} \nabla n \cdot \nabla \sigma - 2n \frac{\partial \epsilon}{\partial \sigma} \nabla^2 n, $$ which is obtained by expanding the divergence analytically. I have on rare occasions found this form useful (in a plane wave code for systems with vacuum), getting $\epsilon$ and its first and second derivatives from Libxc.

Finally, as Susi Lehtola (who would know) points out, it's worth asking if you actually need the potential in real space. You might only need integrals involving the real space potential, in which case you can avoid the divergence term with integration by parts.