# How does libxc calculate the potential of GGA functionals?

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?

• Well, I think you miss that they actually don't have the density gradients as input but something they call contracted gradients $\sigma$. Oct 27 '20 at 19:34

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.

• Good answer. However, if one has a properly defined basis set, one doesn't even need to calculate the potentials in real space, as integrating by parts leads to equations that don't involve the nasty divergence terms, see e.g. doi.org/10.1016/0009-2614(92)85009-Y. This is the way e.g. my fully numerical HelFEM program works, see e.g. doi.org/10.1002/qua.25945 for calculations on atoms. Oct 28 '20 at 17:01
• Thanks, edited in response.
– wcw
Oct 29 '20 at 17:17