Using van der Waals Density Functionals

Author:: Lampros Andrinopoulos, Imperial College London
Date:: November 2012

Activating vdW-DF

The van der Waals energy is calculated in ONETEP using the van der Waals Density Functional method, developed by Dion et al [Dion2004].

The only input variable needed to activate the vdW-DF functional is to set xc_functional VDWDF. If a vdW_df_kernel file is not present in the working directory, it will be automatically generated.

Theory

The form for the exchange-correlation functional proposed by Dion et al is:

(10)\[E_{xc} = E_{x}^{\rm{revPBE}} + E_c^{\rm{PW92}} + E_c^{\rm{nl}}\]

where the non-local exchange-correlation energy is given by:

(11)\[E_{c}^{\mathrm{nl}} = \frac{1}{2}\int\int d{\mathbf{r}}d{\mathbf{r}}'\rho({\mathbf{r}})\phi({\mathbf{r}},{\mathbf{r}}')\rho({\mathbf{r}}')\]

where \(\rho({\mathbf{r}})\) is the electron density at \({\mathbf{r}}\) and \(\phi({\mathbf{r}},{\mathbf{r}}')\) is the nonlocal exchange correlation kernel whose form is explained in [Dion2004].

Non-local correlation energy

The direct calculation of the integral in the form of Eq. (11) is very computationally expensive, as it involves a six-dimensional spatial integral.

The algorithm proposed later by Roman-Perez and Soler [Roman-Perez2009] improves the efficiency of the calculation. They observed that with the form used by Dion et al for \(\phi\), the above expression can be re-written as:

\[E_{c}^{\mathrm{nl}} = \frac{1}{2}\int\int d{\mathbf{r}}d{\mathbf{r}}'\rho({\mathbf{r}})\phi(q,q',r)\rho({\mathbf{r}}')\]

where \(r=|{\mathbf{r}}-{\mathbf{r}}'|\), and \(q\) and \(q'\) are the values of a universal function \(q_0[\rho({\mathbf{r}}),|\nabla \rho({\mathbf{r}})|]\) at \({\mathbf{r}}\) and \({\mathbf{r}}'\). They thus proposed a way to expand the kernel \(\phi\) using interpolating polynomials \(p_\alpha(q)\) for chosen values \(q_\alpha\) of \(q\), and tabulated functions \(\phi_{\alpha\beta}(r)\) for the kernel corresponding to each pair of interpolating polynomials. The interpolating polynomials \(p_{\alpha}\) are cubic splines that evaluate to a Kronecker delta on each respective interpolating point. A mesh of 20 interpolation points is used in Soler’s implementation. The Soler form of the nonlocal energy can be written as:

(12)\[\phi(q_1,q_2,r) = \sum_{\alpha\beta}\phi_{\alpha\beta}(r) p_{\alpha}(q_1) p_{\beta}(q_2)\]

The universal function \(q_0({\mathbf{r}})\) is in practice given by:

(13)\[q_0({\mathbf{r}}) = \Bigg(1 + \frac{\epsilon_c^{\rm{PW92}}}{\epsilon_x^{\rm{LDA}}} + \frac{0.8491}{9}\Big(\frac{\nabla\rho}{2\rho k_F}\Big)^2\Bigg) k_F\]

with \(k_F=(3\pi^2\rho)^{1/3}\). The quantity \(q_0\) is first “saturated” to limit its maximum value, according to:

\[q_0^{\text{sat}}(\rho,{|\nabla{\rho}|}) = q_c \Bigg(1-\exp\Big(-\sum_{m=1}^{m_c}\frac{(q/q_c)^m}{m}\Big)\Bigg)\]

where \(q_c\) is the maximum value of the mesh of \(q_{\alpha}\).

To evaluate this, we first define a quantity \(\theta_{\alpha}({\mathbf{r}}) = \rho({\mathbf{r}}) p_{\alpha}(q(\rho({\mathbf{r}}),\nabla\rho({\mathbf{r}}))\) in real space. In terms of this, Eq. (11) can be written as:

(14)\[E_c^{\mathrm{nl}} = \frac{1}{2} \sum_{\alpha\beta} \int \int d{\mathbf{r}}d{\mathbf{r}}' \theta_{\alpha}({\mathbf{r}}) \theta_{\beta}({\mathbf{r}}') \phi_{\alpha\beta}(r)\]

It can be shown that this can be written as a reciprocal space integral:

(15)\[ E_c^{\mathrm{nl}} = \frac{1}{2} \sum_{\alpha\beta}\int d\mathbf{k} \theta^{*}_{\alpha}(\mathbf{k})\theta_{\beta}(\mathbf{k})\phi_{\alpha\beta}(k)\]

Since the kernel is radially dependent in real space, it is only dependent on the magnitude of the G-vectors, hence the kernel need only be evaluated as a 1-dimensional function \(\phi_{\alpha\beta}(k)\) for each \(\alpha\), \(\beta\).

The kernel \(\phi\) and its second derivatives are tabulated for a specific set of radial points and transformed to reciprocal space. These values are then used to interpolate the kernel at every point \(\mathbf{k}\) in reciprocal space required to calculate Eq. (15).

Kernel

This section details the evaluation of the NLXC kernel. The kernel \(\phi({\mathbf{r}},{\mathbf{r}}')\) as specified by Dion et al [Dion2004] is given by (in atomic units):

\[\phi({\mathbf{r}},{\mathbf{r}}') = \frac{1}{\pi^2}\int_{0}^{\infty}a^2da \int_0^{\infty}b^2db W(a,b) T(\nu(a),\nu(b),\nu'(a),\nu'(b))\]

where

\[T(w,x,y,z) = \frac{1}{2}\Big[\frac{1}{w+x} + \frac{1}{y+z}\Big]\Big[\frac{1}{(w+y)(x+z)}+\frac{1}{(w+z)(y+x)}\Big],\]

and

\[\begin{split}\begin{aligned} W(a,b) = 2\Big[ & (3-a^2)b\cos b \sin a \\ + & (3-b^2)a\cos a \sin b \\ + & (a^2+b^2-3)\sin a\sin b \\ - & 3ab\cos a \cos b \Big]/(a^3b^3),\end{aligned}\end{split}\]

and

\[\nu(y) = 1- e^{-\gamma y^2/d^2}; \quad \nu'(y) = 1- e^{-\gamma y^2/d'^2};\]

where \(d=|{\mathbf{r}}-{\mathbf{r}}'|q_0({\mathbf{r}})\), \(d'=|{\mathbf{r}}-{\mathbf{r}}'|q_0(\mathbf{r'})\)

Following this chain of logic, it is clear that this the kernel can in fact be considered as a function only of \(|{\mathbf{r}}-{\mathbf{r}}'|\), \(q_0({\mathbf{r}})\) and \(q_0({\mathbf{r}}')\), since all other variables are dummy variables which are integrated over. The kernel can therefore be written as

(16)\[\phi(r,q_0({\mathbf{r}}),q_0({\mathbf{r}}'))\]

This makes it possible to evaluate the integrals above so as to tabulate the kernel values numerically for a pre-chosen set of radial points and \(q_0\) values.

Non-local potential

Starting from (15), one can evaluate the potential \(v^{\mathrm{nl}}({\mathbf{r}})\) corresponding to this energy, by evaluating all terms in \(\partial E_{\mathrm{nl}} / \partial n({\mathbf{r}})\). The non-local potential \(v_i^{\mathrm{nl}}\) at point \({\mathbf{r}}_i\) on the grid is thus written explicitly in terms of the derivatives of the \(\theta_{\alpha}\) quantities with respect to the values \(\rho_j\) at all other points on the grid:

(17)\[v_i^{\mathrm{nl}} = \sum_{\alpha}(u_{\alpha i}{\frac{\partial{\theta_{\alpha i}}}{\partial{\rho_i}}}+\sum_j u_{\alpha j} {\frac{\partial{\theta_{\alpha j}}}{\partial{\nabla\rho_j}}}{\frac{\partial{\nabla\rho_j}}{\partial{\rho_i}}})\]

This makes use of the quantities \(u_\alpha({\mathbf{r}})= \sum_{\beta}\mathcal{F}(\theta_{\beta}(\mathbf{k})\phi_{\alpha\beta}(k))\): which are already calculated in the evaluation of the energy.

Using the White and Bird [White1994] approach, Eq. (17) can be written as:

(18)\[v_{\mathrm{nl}}({\mathbf{r}}) = \sum_{\alpha} \Big( u_{\alpha}({\mathbf{r}}){\frac{\partial{\theta_{\alpha}({\mathbf{r}})}}{\partial{\rho({\mathbf{r}})}}} - \int\int i\mathbf{G}\cdot \frac{\nabla\rho({\mathbf{r}}')}{|\nabla\rho({\mathbf{r}}')|} {\frac{\partial{\theta_{\alpha}({\mathbf{r}}')}}{\partial{|\nabla\rho({\mathbf{r}}')|}}}e^{i\mathbf{G}\cdot ({\mathbf{r}}-{\mathbf{r}}')} d{\mathbf{r}}d\mathbf{G} \Big)\]

For this we need to calculate \({\frac{\partial{\theta}}{\partial{\rho}}}\) and \({\frac{\partial{\theta}}{\partial{{|\nabla{\rho}|}}}}\):

(19)\[\begin{split}\begin{aligned} {\frac{\partial{\theta_\alpha}}{\partial{\rho}}} &=p_\alpha + \rho {\frac{\partial{p_\alpha}}{\partial{\rho}}} \nonumber \\ &=p_\alpha + \rho {\frac{\partial{p_\alpha}}{\partial{q}}}{\frac{\partial{q}}{\partial{\rho}}} \nonumber \\ &=p_\alpha + \rho {\frac{\partial{p_\alpha}}{\partial{q}}} \frac{q}{k_F} {\frac{\partial{k_F}}{\partial{\rho}}} + \rho {\frac{\partial{p_\alpha}}{\partial{q}}} k_F ({\frac{\partial{{\varepsilon}_c}}{\partial{\rho}}}{\varepsilon}_x^{-1}-{\varepsilon}_c{\varepsilon}_x^{-2}{\frac{\partial{{\varepsilon}_x}}{\partial{\rho}}} - \frac{8}{3(3\pi^2)^{2/3}}\frac{Z}{4}(\nabla\rho)^2 \rho^{-11/3}) \nonumber \\ &=p_\alpha + q/3{\frac{\partial{p_\alpha}}{\partial{q}}} + k_F\rho {\frac{\partial{p_\alpha}}{\partial{q}}} ({\frac{\partial{{\varepsilon}_c}}{\partial{\rho}}}{\varepsilon}_x^{-1}-{\varepsilon}_c{\varepsilon}_x^{-2}{\frac{\partial{{\varepsilon}_x}}{\partial{\rho}}}- \frac{2Z}{3(3\pi^2)^{2/3}} {|\nabla{\rho}|}^2\rho^{-11/3})\end{aligned}\end{split}\]

(20)\[ {\frac{\partial{\theta_\alpha}}{\partial{{|\nabla{\rho}|}}}} = \rho {\frac{\partial{p_\alpha}}{\partial{q}}} {\frac{\partial{q}}{\partial{{|\nabla{\rho}|}}}} = \frac{Z}{2\rho k_F} \rho {\frac{\partial{p_\alpha}}{\partial{q}}}{|\nabla{\rho}|}\]

Combining Eqs. (18), (19) and (20) gives us the final expression for the nonlocal potential.

Overview of computational algorithm

Module `nlxc_mod`

The main module for the calculation of the non-local energy and potential is nlxc_mod. The tabulation of the kernel \(\phi\) is performed only if a kernel file is not found, by vdwdf_kernel.

The input required to calculate the non-local energy and potential is essentially just the density and its gradient on the fine grid. The calculation of \(q\) and the Fourier transformed \(\theta_\alpha\) from Eq. (15) is performed first, in the routine nlxc_q0_theta. The derivatives of the \(\theta_\alpha\)s with respect to the density and the module of its gradient are calculated on-the-fly in the real-space loop for the calculation of the non-local potential \(v_{nl}\) in Eq. (17). This is to avoid storing unnecessary arrays. From Eq. (18) two transforms are required per \(\alpha\) value, a forward FFT, followed by a backward FFT for calculating the non-local potential.

Subroutines to interpolate the polynomials as well as the kernel using cubic splines are used (spline_interp and interpolate). The interpolating polynomials \(p_\alpha\) used are Kronecker deltas, so they take the value 1 on the interpolating point and are zero at the other points.

Module `vdwdf_kernel`

The kernel \(\phi_{\alpha\beta}(k)\) is tabulated for 1024 radial reciprocal space points and 20 \(q_0\) points. Gaussian quadrature is used to calculate Eq. (16) and then the result is Fourier transformed. The second derivatives of the kernel are calculated by interpolation, and also tabulated. The default name of the file is vdw_df_kernel. The program will first check if this file exists: if it does, it will be loaded in and need not be calculated. If it does not, it will be generated from scratch (which only takes a few minutes) and then it is written out for future re-use in the current working directory.

The format of the vdw_df_kernel file is:

N_alpha  N_radial
max_radius
q_points(:)
kernel(0:N_radial,alpha,beta)
kernel2(0:N_radial,alpha,beta)

where kernel2 is the array of second derivatives of the kernel.

[Dion2004] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).

[Roman-Perez2009] G. Román-Pérez and J. M. Soler, Phys. Rev. Lett. 103, 096102 (2009).

[White1994] J. A. White and D. M. Bird, Phys. Rev. B 50, 4954 (1994).