# 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

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

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

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'\).

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

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

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

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

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

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

The energy is evaluated efficiency as a convolution in Fourier space

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

with