Background

Motivation

In a previous effort Marques et al. created a c-library called libxc. This library consists of 150 different local and semi-local functionals, and it is linked by 18 different codes. In other words, when using this library, scientists have over 1500 different functional and code combinations to choose from. The publication process of a functional becomes easier, because it needs just to be added into one place. In a similar spirit, libvdwxc is intended as a library that enables calculations of a particular family of functionals, which cannot be readily added to libxc.

Currently, vdW-DF is implemented at least in Abinit, Quantum Espresso [GBB+09], VASP [KJ99] and GPAW [MHJ05] [ERM+10] codes.

The implementations vary and so do the results. Typical van der Waals energy landscapes have weak curvature, and thus bonding distances are very sensitive to even small differences in implementations. Such differences can be caused by various means. These include,

Van der Waals forces

It is widely acknowledged that van der Waals (vdW) forces are important with regard to the ground state geometry in a wide range of systems. Accordingly, there exists a variety of empirical, semi-empirical and ab initio approaches in describing effects of van der Waals forces on atomic geometry. One of the most widely applied one is the density functional theory [KS65] (DFT), which is an efficient approach for computer simulations. There exists a huge variety of different DFT-codes and different xc-functional approximations. The reason for large amount of these functionals is that in practice, the exact xc-functional needs to be approximated. Most of the xc-functional approximations are local or semi-local and thus do not have any van der Waals interactions. Within the context of DFT, these van der Waals forces are part of the exchange corelation (xc) functional.

The apparent problem in denity functional theory community is that the approximations in xc-functionals are necessary, but their errors are mathematically not well defined. Therefore, a general approach is to choise a suitable functional by experience depending on the suitability and desired accuracy. In this work, we will focus on the vdW-DF family of van der Waals xc-functionals. After the first vdW-DF density functional for general geometries [DRSchroder+04], several approximations developed. These include vdW-DF2 [LMK+10], vdW-DF-CX [BH14], optPBE and optB88 [KlimevsBM10].

Following the original paper, there are now several variants of this functional available. For physisorption and weak binding, local and semi-local density functionals are not well suited for these problems. However, there a general trend in the vdW-DF fuctional family is becoming a generally applicable functional [BH14] [KlimevsBM11] (bulk properties, surfaces, clusters, dispersive interactions) instead being profiled as the functional for weak dispersive interactions. Below, in this article, we only focus on the technical details of the algorithm, but we recommend other sources for a general review [BCL+15] and interpretation [HBSchroder14] of the vdW-DF functional.

For completeness, even if we here only focus on the vdW-DF functional family, there exist other ab initio vdW-DF approaches including Wannier function approach [Sil08], and semi-empirical such as the TS-method [TS09] and Grimme D-correction series [Gri06] [GEG11].

One of the benefits of open source development is the use of specialized libraries. For example, it is unnecessary for a density functional theory code to implement fast Fourier transform when it can be linked as an external library. The downside of external libraries are the extra effort often needed in compiling and dynamic/static linking of the external libraries. A known success story is, LibXC [MOB12], which is a C-library which allows easy inclusion of hundreds of different xc-functionals to tens of density functional theory codes.

To be precise, with vdW-DF functional family, we mean all functionals which an be written as

\[E_{\rm xc}[n({\bf r})] = E_{\rm x}^{\rm GGA}[n({\bf r})] + E_{\rm c}^{\rm LDA}[n({\bf r})] + E_{\rm c}^{\rm nl}[n({\bf r})],\]

where the vdW-interaction is in the non-local part of the functional

\[E_{\rm c}^{\rm nl}[n({\bf r})] = \frac{1}{2} \int {\rm d}{\bf r}_1 \int {\rm d}{\bf r}_2 n({\bf r}_1) n({\bf r}_2) \phi(q_1, q_2, r_{12}), \label{enl}\]

where \(r_{12} = |{\bf r}_1 - {\bf r}_2|\), and \(q_1\) and \(q_2\) are functions of local density and gradient at \({\bf r}_1\) and \({\bf r}_2\) respectively. This convolution, involving 6D-integration can be solved by an algorithm by Roman-Perez and Soler [RomanPerezS09] with the 2D-spline decomposition of the kernel \(\phi\) and Fourier transforms. Currently, there are many codes which implement a part of the vdW-DF functional family by separate implementations. However, our goal is to write a robust library, which would allow an easy inclusions of vdW-DF family of functionals to any DFT code.

Theory

The van der Waals energy functional consists of a semilocal exchange and correlation functional plus a non-local contribution to the correlation. This latter takes the form of a double integral

\[E_\text{C}^{\text{nl}} = \frac12 \iint n(\boldsymbol r) \phi(q_0(\boldsymbol r), q_0(\boldsymbol r'), \Vert\boldsymbol r - \boldsymbol r'\Vert) n(\boldsymbol r') d r d \boldsymbol r' \label{eq:realspace_vdw_integral}\]

featuring the van der Waals kernel \(\phi(q_0(\boldsymbol r), q_0(\boldsymbol r'), \Vert\boldsymbol r - \boldsymbol r'\Vert)\). Here, \(q_0(\boldsymbol r)\) is a function \(q_0(n(\boldsymbol r),\Vert\boldsymbol \nabla n(\boldsymbol r)\Vert)\) of both the density and its gradient in \(\boldsymbol r\). Thus, each pair of points \(\boldsymbol r\) and \(\boldsymbol r'\) contributes an energy determined by the density semilocally at the two points \(\boldsymbol r\) and \(\boldsymbol r'\) and the distance \(\Vert\boldsymbol r - \boldsymbol r'\Vert\) between them.

Eq.~eqref{eq:realspace_vdw_integral} is prohibitively expensive to evaluate by direct numerical integration, but it can be efficiently approximated by the method of Roman-Perez and Soler.

In this method, the kernel is parametrized in terms of splines that interpolate its radial dependence of both \(\boldsymbol r\) and \(\boldsymbol r'\).

\[E_\text{C}^{\text{nl}} = \frac12 \sum_{\alpha\beta} \iint \theta_\alpha(\boldsymbol r) \phi_{\alpha\beta}(\Vert\boldsymbol r - \boldsymbol r'\Vert) \theta_\beta(\boldsymbol r) d \boldsymbol r d \boldsymbol r'\]

In the method by Roman-Perez and Soler, the kernel is expanded as

\[\phi(q_1, q_2, r_{12}) = \sum_{\alpha \beta} \phi(q^{\rm m}_\alpha, q^{\rm m}_{\beta}, r_{12}) p_{\alpha}(q_1) p_{\beta}(q_2). %\label{phiint}\]

Here, \(q^{\rm m}_{\alpha}\) defines a discrete mesh, typically of 20 interpolation points points, \(\alpha \in (0,1,2, \ldots,19)\). In this fashion, the 3-dimensional kernel is only stored in \(20^2\) one dimensional slices, since interpolation is performed over \(q_1\) and \(q_2\). We shorten the notation to

\[\phi(q^{\rm m}_\alpha, q^{\rm m}_{\beta}, r_{12}) = \phi_{\alpha\beta}(r_{12}).\]

In Eq.~ref{phiint}, \(p_{\alpha}(q)\)‘s are a cubic interpolation splines, i.e. which are defined as piecewise continuous polynomials

\[ \begin{align}\begin{aligned}\begin{split}p_{\alpha}(q^{\rm m}_\beta + {\rm d}q) = \sum_{c=0}^{3} a^c_{\alpha\beta} ({\rm d}q)^c, 0 \leq {\rm dq} \leq q^{\rm m}_{\beta+1}-q^{\rm m}_{\beta}, \\ p_{\alpha}(q^{\rm m}_\beta) = \delta_{\alpha \beta},\end{split}\\\begin{split}%\label{inter} \lim_{\epsilon \rightarrow 0} p_{\alpha}(q_\beta^{\rm m}+\epsilon) - p_{\alpha}(q_\beta^{\rm m}-\epsilon) = 0, \label{a1} \\ \lim_{\epsilon \rightarrow 0} p'_{\alpha}(q_\beta^{\rm m}+\epsilon) - p'_{\alpha}(q_\beta^{\rm m}-\epsilon) = 0, \label{a2} \\ \lim_{\epsilon \rightarrow 0} p''_{\alpha}(q_\beta^{\rm m}+\epsilon) - p''_{\alpha}(q_\beta^{\rm m}-\epsilon) = 0. \label{a3}\end{split}\end{aligned}\end{align} \]

It is easy to show that due to Eq.~ref{inter}, the interpolation polynomial will be exact on the mesh points, but approximate in between. Equations ref{a1}, ref{a2} and ref{a3} just ensure continuity of the interpolation polynomial up to second derivative.

We define an auxillary quantity, \(\theta_\alpha\) which is the key quantity in actual computation

\[\theta_\alpha({\bf r}) = n({\bf r}) p_{\alpha}(q({\bf r})).\]

By inserting Eq.~ref{phiint} into Eq.~ref{enl}, yields

\[E_{\rm c}^{\rm nl}[n({\bf r})] = \frac{1}{2} \sum_{\alpha \beta} \int {\rm d}{\bf r}_1 \int {\rm d}{\bf r}_2 \theta_{\alpha}({\bf r}_1) \theta_{\beta}({\bf r}_2) \phi_{\alpha\beta}(r_{12}),\]

which is just a convolution, and can be written written using convolution theorem

\[E_{\rm c}^{\rm nl}[n({\bf r})] = \frac{1}{2} \int {\rm d}{\bf k} \theta^*_\alpha({\bf k}) \theta_{\beta}({\bf k}) \phi_{\alpha \beta}(k), %\label{enlk}\]

The energy is evaluated efficiency as a convolution in Fourier space

\[E_\text{C}^{\text{nl}} = \frac12 [\Omega?] \sum_{\alpha\beta} \theta_\alpha^*(\boldsymbol k) \phi_{\alpha\beta}(k) \theta_\beta(\boldsymbol k) d \boldsymbol k,\]

where \(\theta_\alpha(\boldsymbol k)\) is the Fourier transform of

\[\theta_\alpha(\boldsymbol r) = n(\boldsymbol r) p_\alpha(q_0(\boldsymbol r)),\]

with

\[q_0(\boldsymbol r) = \textrm{[definition of q0]}\]