.. _introduction:
.. index:: Introduction
Background
************************************
Motivation
======================
In a previous effort Marques *et al.* created a c-library called
:program:`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, :program:`libvdwxc` is
intended as a library that enables calculations of a particular family
of functionals, which cannot be readily added to :program:`libxc`.
Currently, vdW-DF is implemented at least in `Abinit
`_, `Quantum Espresso
`_ :cite:`GiaBarBon09`, `VASP
`_ :cite:`KreJou99` and `GPAW
`_ :cite:`MorHanJac05`
:cite:`EnkRosMor10` 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 :cite:`KohSha65` (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
:cite:`DioRydSch04`, several approximations developed. These
include vdW-DF2 :cite:`LeeMurKon10`, vdW-DF-CX
:cite:`BerHyl14`, optPBE and optB88
:cite:`KliBowMic10`.
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
:cite:`BerHyl14` :cite:`KliBowMic11` (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
:cite:`BerCooLee15` and interpretation
:cite:`HylBerSch14` 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 :cite:`Sil08`, and semi-empirical such as
the TS-method :cite:`TkaSch09` and Grimme D-correction series
:cite:`Gri06` :cite:`GriEhrGoe11`.
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
:cite:`MarOliBur12`, 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
.. math::
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
.. math::
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 :math:`r_{12} = |{\bf r}_1 - {\bf r}_2|`, and :math:`q_1` and
:math:`q_2` are functions of local density and gradient at :math:`{\bf
r}_1` and :math:`{\bf r}_2` respectively. This convolution, involving
6D-integration can be solved by an algorithm by Roman-Perez and Soler
:cite:`RomSol09` with the 2D-spline decomposition of the
kernel :math:`\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
.. math::
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 :math:`\phi(q_0(\boldsymbol r),
q_0(\boldsymbol r'), \Vert\boldsymbol r - \boldsymbol r'\Vert)`.
Here, :math:`q_0(\boldsymbol r)` is a function
:math:`q_0(n(\boldsymbol r),\Vert\boldsymbol \nabla n(\boldsymbol
r)\Vert)` of both the density and its gradient in :math:`\boldsymbol
r`. Thus, each pair of points :math:`\boldsymbol r` and
:math:`\boldsymbol r'` contributes an energy determined by the density
semilocally at the two points :math:`\boldsymbol r` and
:math:`\boldsymbol r'` *and* the distance :math:`\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 :math:`\boldsymbol r` and
:math:`\boldsymbol r'`.
.. math::
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
.. math::
\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, :math:`q^{\rm m}_{\alpha}` defines a discrete mesh, typically of
20 interpolation points points, :math:`\alpha \in (0,1,2,
\ldots,19)`. In this fashion, the 3-dimensional kernel is only stored
in :math:`20^2` one dimensional slices, since interpolation is
performed over :math:`q_1` and :math:`q_2`. We shorten the notation to
.. math::
\phi(q^{\rm m}_\alpha, q^{\rm m}_{\beta}, r_{12}) = \phi_{\alpha\beta}(r_{12}).
In Eq.~\ref{phiint}, :math:`p_{\alpha}(q)`'s are a cubic interpolation
splines, i.e. which are defined as piecewise continuous polynomials
.. math::
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},
%\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}
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, :math:`\theta_\alpha` which is the
key quantity in actual computation
.. math::
\theta_\alpha({\bf r}) = n({\bf r}) p_{\alpha}(q({\bf r})).
By inserting Eq.~\ref{phiint} into Eq.~\ref{enl}, yields
.. math::
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
.. math::
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
.. math::
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 :math:`\theta_\alpha(\boldsymbol k)` is the Fourier transform of
.. math::
\theta_\alpha(\boldsymbol r) = n(\boldsymbol r) p_\alpha(q_0(\boldsymbol r)),
with
.. math::
q_0(\boldsymbol r) = \textrm{[definition of q0]}