Simple way to apply nonlocal van der Waals functionals within all-electron methods

TL;DR

This paper introduces a simple, basis-set independent smoothing method to efficiently apply nonlocal van der Waals functionals within all-electron methods, enabling accurate solid-state property calculations without extensive plane-wave expansions.

Contribution

The authors develop a straightforward smoothing procedure that extends the Román-Pérez and Soler method to all-electron densities, facilitating nonlocal vdW functional calculations in all-electron frameworks.

Findings

The smoothing method yields lattice constants and binding energies in good agreement with existing literature.

Results demonstrate the method's simplicity, efficiency, and compatibility with all-electron calculations.

The approach is basis-set independent and easily implementable.

Abstract

The method based on fast Fourier transforms proposed by G. Rom\'an-P\'erez and J. M. Soler [Phys. Rev. Lett. 103, 096102 (2009)], which allows for a computationally fast implementation of the nonlocal van der Waals (vdW) functionals, has significantly contributed to making the vdW functionals popular in solid-state physics. However, the Rom\'an-P\'erez-Soler method relies on a plane-wave expansion of the electron density; therefore it can not be applied readily to all-electron densities for which an unaffordable number of plane waves would be required for an accurate expansion. In this work, we present the results for the lattice constant and binding energy of solids that were obtained by applying a smoothing procedure to the all-electron density calculated with the linearized augmented plane-wave method. The smoothing procedure has the advantages of being very simple to implement, basis-set independent, and allowing the calculation of the potential. It is also shown that the results agree very well with those from the literature that were obtained with the projector augmented wave method.

Tables

Table 1. Table 1: Equilibrium lattice constant a 0 subscript 𝑎 0 a_{0} (in Å) and cohesive energy E coh subscript 𝐸 coh E_{\text{coh}} (in eV/atom) of 35 cubic solids. The space group number is indicated in parenthesis. An estimate of the error bar for the WIEN2k results is 0.003 Å for a 0 subscript 𝑎 0 a_{0} and 0.03 eV/atom for E coh subscript 𝐸 coh E_{\text{coh}} .

	$a_0$				$E_{\text{coh}}$
	PBE		optB88-vdW		PBE		optB88-vdW
Solid	WIEN2k	VASP	WIEN2k	VASP	WIEN2k	VASP	WIEN2k	VASP
Li (229)	3.437	3.437,111Ref. Klimeš et al., 2011. 3.439333Ref. Park et al., 2015.	3.433	3.432,111Ref. Klimeš et al., 2011. 3.435333Ref. Park et al., 2015.	1.61	1.60,111Ref. Klimeš et al., 2011. 1.61333Ref. Park et al., 2015.	1.58	1.57,111Ref. Klimeš et al., 2011.1.59333Ref. Park et al., 2015.
Na (229)	4.200	4.200,111Ref. Klimeš et al., 2011. 4.193333Ref. Park et al., 2015.	4.153	4.169,111Ref. Klimeš et al., 2011. 4.152333Ref. Park et al., 2015.	1.08	1.08,111Ref. Klimeš et al., 2011. 1.09333Ref. Park et al., 2015.	1.05	1.04,111Ref. Klimeš et al., 2011.1.07333Ref. Park et al., 2015.
K (229)	5.283	5.284,111Ref. Klimeš et al., 2011.5.278,222Ref. Schimka et al., 2013.5.284333Ref. Park et al., 2015.	5.170	5.168,111Ref. Klimeš et al., 2011.5.159,222Ref. Schimka et al., 2013.5.162333Ref. Park et al., 2015.	0.87	0.86,111Ref. Klimeš et al., 2011.0.87,222Ref. Schimka et al., 2013.0.87333Ref. Park et al., 2015.	0.88	0.88,111Ref. Klimeš et al., 2011.0.89333Ref. Park et al., 2015.
Rb (229)	5.671	5.671,111Ref. Klimeš et al., 2011.5.667,222Ref. Schimka et al., 2013.5.668333Ref. Park et al., 2015.	5.506	5.506,111Ref. Klimeš et al., 2011.5.535,222Ref. Schimka et al., 2013.5.501333Ref. Park et al., 2015.	0.78	0.77,111Ref. Klimeš et al., 2011.0.77,222Ref. Schimka et al., 2013.0.77333Ref. Park et al., 2015.	0.82	0.81,111Ref. Klimeš et al., 2011.0.83333Ref. Park et al., 2015.
Cs (229)	6.162	6.160,111Ref. Klimeš et al., 2011.6.156,222Ref. Schimka et al., 2013.6.161333Ref. Park et al., 2015.	5.915	5.899,111Ref. Klimeš et al., 2011.5.900,222Ref. Schimka et al., 2013.5.899333Ref. Park et al., 2015.	0.72	0.70,111Ref. Klimeš et al., 2011.0.72,222Ref. Schimka et al., 2013.0.72333Ref. Park et al., 2015.	0.79	0.79,111Ref. Klimeš et al., 2011.0.66333Ref. Park et al., 2015.
Ca (225)	5.528	5.533,111Ref. Klimeš et al., 2011.5.524,222Ref. Schimka et al., 2013.5.532333Ref. Park et al., 2015.	5.446	5.450,111Ref. Klimeš et al., 2011.5.443,222Ref. Schimka et al., 2013.5.450333Ref. Park et al., 2015.	1.91	1.90,111Ref. Klimeš et al., 2011.1.91,222Ref. Schimka et al., 2013.1.90333Ref. Park et al., 2015.	1.88	1.88,111Ref. Klimeš et al., 2011.1.88333Ref. Park et al., 2015.
Sr (225)	6.025	6.019,111Ref. Klimeš et al., 2011.6.020,222Ref. Schimka et al., 2013.6.026333Ref. Park et al., 2015.	5.919	5.917,111Ref. Klimeš et al., 2011.5.917,222Ref. Schimka et al., 2013.5.911333Ref. Park et al., 2015.	1.61	1.61,111Ref. Klimeš et al., 2011.1.61,222Ref. Schimka et al., 2013.1.61333Ref. Park et al., 2015.	1.60	1.61,111Ref. Klimeš et al., 2011.1.61333Ref. Park et al., 2015.
Ba (229)	5.020	5.028,111Ref. Klimeš et al., 2011.5.018,222Ref. Schimka et al., 2013.5.030333Ref. Park et al., 2015.	4.905	4.917,111Ref. Klimeš et al., 2011.4.903,222Ref. Schimka et al., 2013.4.915333Ref. Park et al., 2015.	1.89	1.88,111Ref. Klimeš et al., 2011.1.88,222Ref. Schimka et al., 2013.1.88333Ref. Park et al., 2015.	1.99	1.99,111Ref. Klimeš et al., 2011.2.00333Ref. Park et al., 2015.
Al (225)	4.041	4.041,111Ref. Klimeš et al., 2011. 4.040333Ref. Park et al., 2015.	4.054	4.054,111Ref. Klimeš et al., 2011. 4.054333Ref. Park et al., 2015.	3.44	3.50,111Ref. Klimeš et al., 2011. 3.54333Ref. Park et al., 2015.	3.24	3.34,111Ref. Klimeš et al., 2011.3.38333Ref. Park et al., 2015.
C (227)	3.575	3.574,111Ref. Klimeš et al., 2011. 3.573333Ref. Park et al., 2015.	3.577	3.577,111Ref. Klimeš et al., 2011. 3.575333Ref. Park et al., 2015.	7.71	7.70,111Ref. Klimeš et al., 2011. 7.85333Ref. Park et al., 2015.	7.72	7.70,111Ref. Klimeš et al., 2011.7.89333Ref. Park et al., 2015.
Si (227)	5.471	5.465,111Ref. Klimeš et al., 2011. 5.469333Ref. Park et al., 2015.	5.464	5.460,111Ref. Klimeš et al., 2011. 5.469333Ref. Park et al., 2015.	4.57	4.62,111Ref. Klimeš et al., 2011. 4.61333Ref. Park et al., 2015.	4.67	4.74,111Ref. Klimeš et al., 2011.4.74333Ref. Park et al., 2015.
Ge (227)	5.764	5.766,111Ref. Klimeš et al., 2011. 5.783333Ref. Park et al., 2015.	5.762	5.762,111Ref. Klimeš et al., 2011. 5.798333Ref. Park et al., 2015.	3.73	3.72,111Ref. Klimeš et al., 2011. 3.74333Ref. Park et al., 2015.	3.93	3.90,111Ref. Klimeš et al., 2011.3.92333Ref. Park et al., 2015.
Sn (227)	6.657	6.652333Ref. Park et al., 2015.	6.640	6.637333Ref. Park et al., 2015.	3.17	3.20333Ref. Park et al., 2015.	3.42	3.44333Ref. Park et al., 2015.
Pb (225)	5.034	5.031333Ref. Park et al., 2015.	5.026	5.018333Ref. Park et al., 2015.	2.98	2.98333Ref. Park et al., 2015.	3.32	3.32333Ref. Park et al., 2015.
V (229)	2.998	2.996,222Ref. Schimka et al., 2013.2.978333Ref. Park et al., 2015.	2.986	2.984,222Ref. Schimka et al., 2013.2.969333Ref. Park et al., 2015.	5.37	5.37,222Ref. Schimka et al., 2013.5.41333Ref. Park et al., 2015.	5.63	5.74333Ref. Park et al., 2015.
Fe (229)	2.832	2.832,222Ref. Schimka et al., 2013.2.832333Ref. Park et al., 2015.	2.820	2.821,222Ref. Schimka et al., 2013.2.821333Ref. Park et al., 2015.	4.92	4.89,222Ref. Schimka et al., 2013.5.16333Ref. Park et al., 2015.	5.03	5.11333Ref. Park et al., 2015.
Ni (225)	3.518	3.507,222Ref. Schimka et al., 2013.3.518333Ref. Park et al., 2015.	3.513	3.512,222Ref. Schimka et al., 2013.3.511333Ref. Park et al., 2015.	4.76	4.75,222Ref. Schimka et al., 2013.4.80333Ref. Park et al., 2015.	4.85	4.98333Ref. Park et al., 2015.
Cu (225)	3.632	3.635,111Ref. Klimeš et al., 2011.3.631,222Ref. Schimka et al., 2013.3.635333Ref. Park et al., 2015.	3.630	3.632,111Ref. Klimeš et al., 2011.3.622,222Ref. Schimka et al., 2013.3.629333Ref. Park et al., 2015.	3.52	3.49,111Ref. Klimeš et al., 2011.3.50,222Ref. Schimka et al., 2013.3.49333Ref. Park et al., 2015.	3.59	3.52,111Ref. Klimeš et al., 2011.3.57333Ref. Park et al., 2015.
Nb (229)	3.312	3.310,222Ref. Schimka et al., 2013.3.308333Ref. Park et al., 2015.	3.306	3.307,222Ref. Schimka et al., 2013.3.303333Ref. Park et al., 2015.	6.98	6.96,222Ref. Schimka et al., 2013.7.06333Ref. Park et al., 2015.	7.42	7.49333Ref. Park et al., 2015.
Mo (229)	3.162	3.161,222Ref. Schimka et al., 2013.3.151333Ref. Park et al., 2015.	3.162	3.162,222Ref. Schimka et al., 2013.3.154333Ref. Park et al., 2015.	6.28	6.28,222Ref. Schimka et al., 2013.7.76333Ref. Park et al., 2015.	6.85	6.95333Ref. Park et al., 2015.
Rh (225)	3.832	3.830,111Ref. Klimeš et al., 2011.3.831,222Ref. Schimka et al., 2013.3.824333Ref. Park et al., 2015.	3.835	3.831,111Ref. Klimeš et al., 2011.3.834,222Ref. Schimka et al., 2013.3.829333Ref. Park et al., 2015.	5.74	5.82,111Ref. Klimeš et al., 2011.5.70,222Ref. Schimka et al., 2013.6.02333Ref. Park et al., 2015.	6.06	6.10,111Ref. Klimeš et al., 2011.6.34333Ref. Park et al., 2015.
Pd (225)	3.943	3.943,111Ref. Klimeš et al., 2011.3.933,222Ref. Schimka et al., 2013.3.942333Ref. Park et al., 2015.	3.940	3.941,111Ref. Klimeš et al., 2011.3.930,222Ref. Schimka et al., 2013.3.938333Ref. Park et al., 2015.	3.71	3.71,111Ref. Klimeš et al., 2011.3.76,222Ref. Schimka et al., 2013.3.74333Ref. Park et al., 2015.	4.03	3.96,111Ref. Klimeš et al., 2011.4.04333Ref. Park et al., 2015.
Ag (225)	4.148	4.154,111Ref. Klimeš et al., 2011.4.145,222Ref. Schimka et al., 2013.4.147333Ref. Park et al., 2015.	4.136	4.141,111Ref. Klimeš et al., 2011.4.127,222Ref. Schimka et al., 2013.4.130333Ref. Park et al., 2015.	2.53	2.50,111Ref. Klimeš et al., 2011.2.52,222Ref. Schimka et al., 2013.2.52333Ref. Park et al., 2015.	2.81	2.76,111Ref. Klimeš et al., 2011.2.82333Ref. Park et al., 2015.
Ta (229)	3.320	3.317,222Ref. Schimka et al., 2013.3.309333Ref. Park et al., 2015.	3.311	3.307,222Ref. Schimka et al., 2013.3.306333Ref. Park et al., 2015.	8.22	8.11,222Ref. Schimka et al., 2013.8.41333Ref. Park et al., 2015.	8.33	8.50333Ref. Park et al., 2015.
W (229)	3.185	3.182,222Ref. Schimka et al., 2013.3.172333Ref. Park et al., 2015.	3.184	3.181,222Ref. Schimka et al., 2013.3.178333Ref. Park et al., 2015.	8.34	8.39,222Ref. Schimka et al., 2013.8.48333Ref. Park et al., 2015.	8.84	9.01333Ref. Park et al., 2015.
Ir (225)	3.874	3.868,222Ref. Schimka et al., 2013.3.873333Ref. Park et al., 2015.	3.882	3.879,222Ref. Schimka et al., 2013.3.886333Ref. Park et al., 2015.	7.35	7.31,222Ref. Schimka et al., 2013.7.28333Ref. Park et al., 2015.	7.63	7.60333Ref. Park et al., 2015.
Pt (225)	3.971	3.969,222Ref. Schimka et al., 2013.3.968333Ref. Park et al., 2015.	3.979	3.978,222Ref. Schimka et al., 2013.3.980333Ref. Park et al., 2015.	5.55	5.51,222Ref. Schimka et al., 2013.5.58333Ref. Park et al., 2015.	5.85	5.90333Ref. Park et al., 2015.
Au (225)	4.161	4.154,222Ref. Schimka et al., 2013.4.157333Ref. Park et al., 2015.	4.166	4.156,222Ref. Schimka et al., 2013.4.161333Ref. Park et al., 2015.	3.03	3.05,222Ref. Schimka et al., 2013.3.04333Ref. Park et al., 2015.	3.82	3.40333Ref. Park et al., 2015.
LiF (225)	4.070	4.068111Ref. Klimeš et al., 2011.	4.032	4.033111Ref. Klimeš et al., 2011.	4.33	4.32111Ref. Klimeš et al., 2011.	4.54	4.53111Ref. Klimeš et al., 2011.
LiCl (225)	5.152	5.152111Ref. Klimeš et al., 2011.	5.113	5.114111Ref. Klimeš et al., 2011.	3.37	3.42111Ref. Klimeš et al., 2011.	3.58	3.61111Ref. Klimeš et al., 2011.
NaF (225)	4.706	4.708111Ref. Klimeš et al., 2011.	4.642	4.647111Ref. Klimeš et al., 2011.	3.84	3.82111Ref. Klimeš et al., 2011.	4.05	4.02111Ref. Klimeš et al., 2011.
NaCl (225)	5.700	5.701111Ref. Klimeš et al., 2011.	5.617	5.622111Ref. Klimeš et al., 2011.	3.10	3.15111Ref. Klimeš et al., 2011.	3.31	3.33111Ref. Klimeš et al., 2011.
MgO (225)	4.259	4.257111Ref. Klimeš et al., 2011.	4.234	4.231111Ref. Klimeš et al., 2011.	4.99	4.97111Ref. Klimeš et al., 2011.	5.24	5.21111Ref. Klimeš et al., 2011.
SiC (216)	4.385	4.377111Ref. Klimeš et al., 2011.	4.380	4.375111Ref. Klimeš et al., 2011.	6.40	6.44111Ref. Klimeš et al., 2011.	6.49	6.52111Ref. Klimeš et al., 2011.
GaAs (216)	5.749	5.752111Ref. Klimeš et al., 2011.	5.742	5.751111Ref. Klimeš et al., 2011.	3.15	3.15111Ref. Klimeš et al., 2011.	3.37	3.36111Ref. Klimeš et al., 2011.

Table 2. Table 2: Equilibrium lattice constant a 0 subscript 𝑎 0 a_{0} (in Å) and cohesive energy E coh subscript 𝐸 coh E_{\text{coh}} (in meV/atom) of rare-gas solids. The space group number is indicated in parenthesis. An estimate of the error bar for the WIEN2k results is 0.01 Å for a 0 subscript 𝑎 0 a_{0} and 5 meV/atom for E coh subscript 𝐸 coh E_{\text{coh}} .

	$a_0$					$E_{\text{coh}}$
	LDA		optB88-vdW			LDA		optB88-vdW
Solid	WIEN2k	CP2K111Ref. Tran and Hutter, 2013.	WIEN2k	CP2K111Ref. Tran and Hutter, 2013.	VASP222Ref. Callsen and Hamada, 2015.	WIEN2k	CP2K111Ref. Tran and Hutter, 2013.	WIEN2k	CP2K111Ref. Tran and Hutter, 2013.	VASP222Ref. Callsen and Hamada, 2015.
Ne (225)	3.86	3.86	4.27	4.24	4.25	87	92	49	59	45
Ar (225)	4.94	4.94	5.24	5.24	5.22	138	136	136	143	143
Kr (225)	5.33	5.36	5.63	5.63	5.61	169	164	179	181	180

Table 3. Table 3: Equilibrium lattice constant c 0 subscript 𝑐 0 c_{0} (in Å) and interlayer binding energy E b subscript 𝐸 b E_{\text{b}} (in meV/atom) of hexagonal layered solids. The intralayer lattice constant a 𝑎 a was kept fixed at the experimental value of 2.462 and 2.503 Å for graphite and h-BN, respectively. The space group number is indicated in parenthesis. An estimate of the error bar for the WIEN2k results is 0.03 Å for c 0 subscript 𝑐 0 c_{0} and 5 meV/atom for E b subscript 𝐸 b E_{\text{b}} .

	$c_0$				$E_{\text{b}}$
	LDA		optB88-vdW		LDA		optB88-vdW
Solid	WIEN2k	VASP	WIEN2k	VASP	WIEN2k	VASP	WIEN2k	VASP
Graphite (194)	6.68	6.62,111Ref. Graziano et al., 2012.6.75222Ref. Björkman et al., 2012.	6.72	6.72,111Ref. Graziano et al., 2012.6.76333Ref. Björkman, 2014.	24	24,111Ref. Graziano et al., 2012.25222Ref. Björkman et al., 2012.	68	65,111Ref. Graziano et al., 2012.66333Ref. Björkman, 2014.
h-BN (194)	6.48	6.42,111Ref. Graziano et al., 2012.6.58222Ref. Björkman et al., 2012.	6.59	6.60,111Ref. Graziano et al., 2012.6.64333Ref. Björkman, 2014.	28	28,111Ref. Graziano et al., 2012.28222Ref. Björkman et al., 2012.	69	65,111Ref. Graziano et al., 2012.67333Ref. Björkman, 2014.

Equations

$$E_{\text{xc}}=E_{\text{xc}}^{\text{SL/hybrid}}+E_{\text{c,disp}},$$

$$E_{\text{c,disp}}^{\text{NL}}=\frac{1}{2}\int \int \rho ({\mathbf{r}}_1)\Phi \left({\mathbf{r}}_1,{\mathbf{r}}_2\right)\rho ({\mathbf{r}}_2)d^3{r}_1{d}^3{r}_2,$$

$${\psi }_{n\mathbf{k}}(\mathbf{r})=\sum _{\mathbf{K}}{c}_{n\mathbf{k}\mathbf{K}}{\varphi }_{\mathbf{k}\mathbf{K}}(\mathbf{r})+\sum _i{c}_{n\mathbf{k}i}^{\text{LO}}{\varphi }_i^{\text{LO}}(\mathbf{r}),$$

$$\phi_{\bm{k}\bm{K}}(\bm{r})=\left\{\begin{array}[]{l@{\quad}l}\sum\limits_{\ell,m}\sum\limits_{f}d_{f\bm{k}\bm{K}}^{\alpha\ell m}D_{f}^{\alpha\ell}(r_{\alpha})Y_{\ell m}(\hat{\bm{r}}_{\alpha}),&\bm{r}\in\text{S}_{\alpha}\\ \frac{1}{\sqrt{\Omega}}e^{i\left(\bm{k}+\bm{K}\right)\cdot\bm{r}},&\bm{r}\in\text{I}\end{array}\right.$$

$$\phi_{i}^{\text{LO}}(\bm{r})=\left\{\begin{array}[]{l@{\quad}l}\sum\limits_{f}d_{\text{LO},if}^{\alpha_{i}\ell_{i}m_{i}}D_{if}^{\alpha_{i}\ell_{i}}(r_{\alpha_{i}})Y_{\ell_{i}m_{i}}(\hat{\bm{r}}_{\alpha_{i}}),&\bm{r}\in\text{S}_{\alpha_{i}}\\ 0,&\bm{r}\in\text{I}\end{array}\right.,$$

$$\rho_{\text{s}}(\bm{r})=\left\{\begin{array}[]{l@{\quad}l}\rho(\bm{r}),&\rho(\bm{r})\leqslant\rho_{\text{c}}\\ \frac{\rho(\bm{r})+A\rho_{\text{c}}\left(\rho(\bm{r})-\rho_{\text{c}}\right)^{n}}{1+A\left(\rho(\bm{r})-\rho_{\text{c}}\right)^{n}},&\rho(\bm{r})>\rho_{\text{c}}\end{array},\right.$$

$$\Phi \left(q_1,q_2,r_{12}\right)=\sum _{\alpha ,\beta }\Phi \left(q_{\alpha },q_{\beta },r_{12}\right)p_{\alpha }(q_1)p_{\beta }(q_2),$$

$$\nabla {\rho }_{\text{s}}=\frac{d{\rho }_{\text{s}}}{d\rho }\nabla \rho$$

$${\nabla }^2{\rho }_{\text{s}}=\frac{d^2{\rho }_{\text{s}}}{d{\rho }^2}{\left|\nabla \rho \right|}^2+\frac{d{\rho }_{\text{s}}}{d\rho }{\nabla }^2\rho ,$$

$$\frac{d{\rho }_{\text{s}}}{d\rho }=\frac{1+\left(1-n\right)A{\left(\rho -{\rho }_{\text{c}}\right)}^n}{{\left(1+A{\left(\rho -{\rho }_{\text{c}}\right)}^n\right)}^2}$$

$$\frac{d^2{\rho }_{\text{s}}}{d{\rho }^2}=-\frac{n(1+n)A{\left(\rho -{\rho }_{\text{c}}\right)}^{n-1}+n(1-n)A^2{\left(\rho -{\rho }_{\text{c}}\right)}^{2n-1}}{{\left(1+A{\left(\rho -{\rho }_{\text{c}}\right)}^n\right)}^3}.$$

$$v_{\text{c,disp}}^{\text{NL}}(\mathbf{r})=\frac{\delta E_{\text{c,disp}}^{\text{NL}}}{\delta \rho (\mathbf{r})}=\frac{\delta E_{\text{c,disp}}^{\text{NL}}}{\delta {\rho }_{\text{s}}(\mathbf{r})}\frac{d{\rho }_{\text{s}}(\mathbf{r})}{d\rho (\mathbf{r})},$$

$$v_{\text{c,disp}}^{\text{NL}}\left(\mathbf{r}\right)=\sum _{\mathbf{G}}{v}_{\text{c,disp}}^{\text{NL},\mathbf{G}}{e}^{i\mathbf{G}\cdot \mathbf{r}},$$

$$v_{\text{c,disp}}^{\text{NL}}\left(\mathbf{r}\right)=\sum _{\alpha }\sum _L\sum _{M=0}^L\sum _p{v}_{\text{c,disp}}^{\text{NL},\alpha LMp}\left(r_{\alpha }\right)Z_{LMp}\left({\hat{\mathbf{r}}}_{\alpha }\right),$$

$$Z_{L0}=Y_{L0}$$

$$Z_{LM+}=\frac{(-1{)}^M}{\sqrt{2}}\left(Y_{LM}+Y_{LM}^{\ast }\right),$$

$$Z_{LM-}=\frac{(-1{)}^M}{i\sqrt{2}}\left(Y_{LM}-Y_{LM}^{\ast }\right)$$

$$e^{i\mathbf{G}\cdot \mathbf{r}}$$

$$v_{\text{c,disp}}^{\text{NL},\alpha LMp}\left(r_{\alpha }\right)=4\pi i^L\sum _{\mathbf{G}}{v}_{\text{c,disp}}^{\text{NL},\mathbf{G}}{e}^{i\mathbf{G}\cdot {\tau }_{\alpha }}{j}_L\left(G{r}_{\alpha }\right)Z_{LMp}^{\ast }\left(\hat{\mathbf{G}}\right),$$

Full text

I Introduction

It is well known that in density functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) the exchange-correlation (xc) functionals $E_{\text{xc}}$ of the semilocal and hybrid approximations do not account properly for the van der Waals (vdW) interactions since the attractive part, the London dispersion forces, is missing and, as a consequence, lead to erratic results when applied to systems where the vdW interactions play an important role.Grimme et al. (2016) In this respect, two extreme and opposite functionals are the local density approximation (LDA)Kohn and Sham (1965) and the generalized gradient approximation (GGA) of Becke, Lee, Yang, and Parr Becke (1988); Lee et al. (1988) that lead to severe overbinding and underbinding (or even no binding at all), respectively.Kristyán and Pulay (1994); Pérez-Jordá and Becke (1995) Therefore, dispersion correction terms to be added to a semilocal (SL) or hybrid xc functional,

[TABLE]

have been proposed, such that much more reliable results can be obtained for vdW systems with DFT (see Refs. Klimeš and Michaelides, 2012; Dobson, 2014; Berland et al., 2015; Grimme et al., 2016 for reviews).

The most simple and popular type of correction for $E_{\text{c,disp}}$ consists of an atom-pairwise termWu et al. (2001); Wu and Yang (2002); Hasegawa and Nishidate (2004); Grimme (2004) for each atom pair $A$-$B$ (separated by $R_{AB}$): $-f^{\text{damp}}\left(R_{AB}\right)C_{6}^{AB}/R_{AB}^{6}$, where $C_{6}^{AB}$ are the dispersion coefficients and $f^{\text{damp}}$ is a damping function. Such corrections lead to little increase in the computational time and have proven to be very useful for improving the reliability of DFT calculations.Grimme et al. (2016) Depending on the particular scheme, the coefficients $C_{6}^{AB}$ are precomputed Wu et al. (2001); Wu and Yang (2002); Hasegawa and Nishidate (2004); Grimme (2004) or evaluated on the fly by using properties of the system like the atomic positions or the electron density.Becke and Johnson (2005); Tkatchenko and Scheffler (2009); Grimme et al. (2010) Nowadays, the most used atom-pairwise methods are DFT-D3 from Grimme and co-workers Grimme et al. (2010, 2011) and PBE+TS from Tkatchenko and Scheffler.Tkatchenko and Scheffler (2009); Tkatchenko et al. (2012) The atom-pairwise methods can be applied to molecules and extended systems.

Also popular, are the so-called nonlocal (NL) vdW dispersion terms of the formDion et al. (2004); Langreth et al. (2009)

[TABLE]

where the kernel $\Phi$ depends on the electron density $\rho$ and its derivative $\nabla\rho$ as well as on $\left|\bm{r}_{1}-\bm{r}_{2}\right|$. The first functional of the form of Eq. (2) that could be applied to all kinds of system was proposed by Dion et al.Dion et al. (2004) (DRSLL) and was derived as a simplification of the adiabatic connection formula.Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977) The DRSLL term was originally added to the semilocal functional consisting of the GGA revPBEZhang and Yang (1998) for exchange and LDA Vosko et al. (1980); Perdew and Wang (1992) for correlation. Since then, other kernels $\Phi$ in Eq. (2) or associated semilocal functionals have been proposed by various authors. Lee et al. (2010); Klimeš et al. (2010); Cooper (2010); Vydrov and Van Voorhis (2009, 2010); Wellendorff et al. (2010); Sabatini et al. (2013); Hamada (2014); Berland and Hyldgaard (2014); Peng et al. (2016) Overall, the most recent versions of nonlocal vdW functionals and pair-wise methods are rather similar in terms of accuracy (see, e.g., Refs. Grimme, 2011; Goerigk, 2014; Rêgo et al., 2015; Tran et al., 2016; Lozano et al., 2017), nevertheless, the double integration in Eq. (2) makes the calculations with the nonlocal vdW functionals more expensive.

Brute-force methods to carry out the double integration in Eq. (2) have been proposed,Lazić et al. (2010); Nabok et al. (2011) however an important step in the development of nonlocal vdW functionals for periodic systems was made by Román-Pérez and SolerRomán-Pérez and Soler (2009) (RPS) who proposed a very efficient method for evaluating Eq. (2). Their method, which is based on fast Fourier transforms (FFT) and the convolution theorem, is now the standard method and is implemented in various solid-state codes. Sabatini et al. (2012); Wellendorff et al. (2010); Klimeš et al. (2011); Tran and Hutter (2013); Larsen et al. (2017) Furthermore, Román-Pérez and Soler also showed that their method allows for a straightforward calculation of the functional derivative of Eq. (2) that is much simpler than the complicated ways from Refs. Thonhauser et al., 2007; Gulans et al., 2009. Also, in Ref. Sabatini et al., 2012, the formula for the contribution to the stress tensor was derived. For these reasons, the RPS method is the preferred one compared to the others that have been proposed.Gulans et al. (2009); Lazić et al. (2010); Nabok et al. (2011)

However, since the RPS method relies on a plane-wave expansion of $\rho$, it is not obvious how to apply it to all-electron densities, since their large variations close to the nuclei would require an unrealistically large plane-wave expansion. This is the reason why, to our knowledge, the RPS method has been implemented only into codes relying exclusively on a plane-wave expansion of $\rho$, i.e., pseudopotential codes. In order to use the RPS method within an all-electron code, a simple solution would be to smooth the electron density close to the nuclei, such that a pure plane-wave expansion of $\rho$ is possible.

The goal of this work is to present a smoothing procedure to all-electron densities and to study in detail which degree of smoothing should be applied in order to make the calculations affordable without sacrificing too much accuracy in the calculation of Eq. (2). This is done in the framework of the linearized-augmented plane-waveAndersen (1975); Singh and Nordström (2006) (LAPW) method, which provides very accurate all-electron densities. However, our smoothing method is basis-set independent and, therefore, can be implemented in any all-electron code. Thus, we expect our work to contribute to a much more widespread use of the nonlocal vdW functionals in the all-electron solid-state community.

The properties that will be considered for our tests are the lattice constant and binding energy of various types of solids, and the results will be compared to results that were obtained with the projector augmented-wave (PAW) method.Blöchl (1994)

The paper is organized as follows. Section II gives details about the used methods. Then, the results are presented and discussed in Sec. III, while Sec. IV summarizes the main conclusions of this work.

II Methodology

We start with a brief introduction to the LAPW method that is used in the present work to solve the Kohn-ShamKohn and Sham (1965) (KS) equations of DFT. Adopting a notation that is common to all flavors of the LAPW method, Singh and Nordström (2006); Sjöstedt et al. (2000); Michalicek et al. (2013) the Bloch orbitals ($n$ is the band index and $\bm{k}$ is a vector in the first Brillouin zone) are expanded as

[TABLE]

where

[TABLE]

are basis functions (indexed with the reciprocal lattice vector $\bm{K}$) that consist of a linear combination of products between a radial function $D_{f}^{\alpha\ell}$ and a spherical harmonics $Y_{\ell m}$ inside the atomic sphere $\text{S}_{\alpha}$ centered at nucleus $\alpha$, and of a plane wave in the interstitial region I. The coefficients $d_{f\bm{k}\bm{K}}^{\alpha\ell m}$ are determined from the requirement of continuity of $\phi_{\bm{k}\bm{K}}$ [and its derivative(s) depending on the LAPW flavor] at the sphere boundary. The basis functions

[TABLE]

are the so-called local orbitals (LO) that are defined only inside the spheres and set to zero in the interstitial region.Singh and Nordström (2006) The number of basis functions in Eq. (3) is determined by the cutoff $K_{\text{max}}$ such that $\left|\bm{k}+\bm{K}\right|\leqslant K_{\text{max}}$ and the number of LO in the second term. Note that the basis-set expansion Eq. (3) is used for the valence and unoccupied orbitals, but not for the deep-lying core orbitals which are calculated by solving the radial KS equations inside the MT spheres numerically without basis-set expansion.Singh and Nordström (2006)

Once the orbitals are obtained, the electron density $\rho$ is calculated and expanded in (real) spherical harmonics and plane waves inside the atomic spheres and interstitial regions, respectively. The great advantage of the LAPW method is to provide a virtually exact all-electron solution of the KS equations for a given xc functional.

On the other hand, however, the LAPW method does not allow for a straightforward use of the RPS method to calculate Eq. (2) since a pure plane-wave expansion of the all-electron density $\rho$ is out of reach. In order to use the RPS method with the LAPW method, the most obvious choice is to smooth $\rho$ in the core region such that a plane-wave expansion is possible. A simple procedure to construct a smooth density $\rho_{\text{s}}$ is given by

[TABLE]

where $\rho_{\text{c}}$ and $n$ are parameters, and $A$ is set to 1 Bohr3n and is introduced only for the sake of consistency of the units [e.g., $A$ would be $(0.529177)^{3n}$ Å3n if $\rho$ was expressed in Å*-3*]. $\rho_{\text{c}}$ is the density cutoff that has to be chosen carefully such that the plane-wave expansion of $\rho_{\text{s}}$ everywhere in the unit cell is small enough to make the calculation with the RPS method affordable, but without loosing too much accuracy with respect to the reference calculation with the original density. $n$ determines until which order the derivatives of $\rho_{\text{s}}$ are continuous, and in Appendix A it is shown explicitly that $\nabla\rho_{\text{s}}$ and $\nabla^{2}\rho_{\text{s}}$ are continuous if $n\geqslant 2$.

Figures 1 and 2 show the effect of the smoothing procedure for different values of $n$ and $\rho_{\text{c}}$ in solid Li and Au, which represent cases of a very light and a very heavy atom, respectively. It is clear that a pure plane-wave expansion of such densities $\rho_{\text{s}}$ without the peaks close to the nuclei should converge much faster than the original density $\rho$. From Fig. 1, we can see that the smoothest curve is obtained with $n=1$ (continuity of $\nabla\rho_{\text{s}}$), while requiring higher continuous derivatives of $\rho_{\text{s}}$ leads to a maximum and overall to a less smooth density $\rho_{\text{s}}$. Convergence tests of Eq. (2) with respect to the cutoff for the reciprocal lattice vector $G_{\text{max}}$ of the plane-wave expansion of $\rho_{\text{s}}=\sum_{G<G_{\text{max}}}\rho_{\text{s}}^{\bm{G}}e^{i\bm{G}\cdot\bm{r}}$ showed that, indeed, $n=1$ requires a smaller $G_{\text{max}}$ than larger $n$. Therefore, for efficiency $n=1$ is a better choice, but the price to pay is to have second and higher derivatives of $\rho_{\text{s}}$ that are discontinuous, which is not elegant but of no consequence from the practical point of view according to our tests. $n=1$ in Eq. (6) has been chosen for all calculations presented in Sec. III.

Figure 2 shows $\rho_{\text{s}}$ for different values of $\rho_{\text{c}}$ (with $n=1$). For Li for instance, the cutoffs $\rho_{\text{c}}=0.5$, 1.0, and 1.5 Bohr*-3* lead to a change in the density in the region corresponding to a distance $d$ from the nucleus that is smaller than approximately 0.55, 0.45, and 0.35 Bohr, respectively. The crucial point is to choose $\rho_{\text{c}}$ such that a good balance between computational efficiency and accuracy is achieved. This will be studied in detail in Sec. III.

Like the all-electron density $\rho$, $\rho_{\text{s}}$ has a kink at the nucleus, however, if $\rho$ at the nucleus is much larger than $\rho_{\text{c}}$, then this kink is strongly attenuated [see Eqs. (8) and (10)]. On a scale like in Figs. 1 and 2, a kink is clearly visible only for the very lightest atoms.

The RPS method is described in detail in Ref. Román-Pérez and Soler, 2009, thus only the basic idea is now briefly summarized. The kernel $\Phi$ in Eq. (2) depends on $\rho$ and $\nabla\rho$ at $\bm{r}_{1}$ and $\bm{r}_{2}$ via a function $q_{0}\left(\rho,\left|\nabla\rho\right|\right)$ that is evaluated at $\bm{r}_{1}$ and $\bm{r}_{2}$. Such a dependency of $\Phi$ on $\bm{r}_{1}$ and $\bm{r}_{2}$ individually, and not only on $r_{12}=\left|\bm{r}_{1}-\bm{r}_{2}\right|$, prevents Eq. (2) to be evaluated by a single convolution. Therefore, Román-Pérez and Soler proposed to expand $\Phi$ with an interpolation formula consisting of factorized terms,

[TABLE]

allowing the use of convolution that is performed efficiently with FFT. The sum on the right-hand side of Eq. (7) runs over a two-dimensional mesh $(q_{\alpha},q_{\beta})$ of values of $q_{0}$ at which $\Phi$ is pre-calculated and multiplied by cubic polynomials $p_{\alpha}$ evaluated at $q_{i}=q_{0}\left(\rho(\bm{r}_{i}),\left|\nabla\rho(\bm{r}_{i})\right|\right)$ ($i=1,2$). The accuracy of the interpolation is determined by the number of points $N_{q}$ on the $q$-mesh, and the chosen cutoff value $q_{0}^{\text{c}}$ should be larger than all values of $q_{0}$ in the unit cell. Not explicitly shown in Eq. (7), an interpolation over $r_{12}$ [or equivalently over $k$ for the Fourier transform $\Phi\left(q_{\alpha},q_{\beta},k\right)]$ is also done. A careful study of the influence of $N_{q}$ and $q_{0}^{\text{c}}$ on the lattice constant and cohesive energy of solids was reported in Ref. Klimeš et al., 2011. It was concluded that $N_{q}=30$ and $q_{0}^{\text{c}}=10$ lead to results that are well converged, therefore these parameters were chosen for our calculations. However, in two cases, namely K and Cs, we observed that a rather large change in the lattice constant (of the order of 0.02 Å) was obtained by increasing $N_{q}$ to 40. Thus, the results in Sec. III for K and Cs were obtained with $N_{q}=40$ and $q_{0}^{\text{c}}=10$, while for all other solids $N_{q}=30$ was used. We mention that Wu and GygiWu and Gygi (2012) pointed out that $q_{1}q_{2}\Phi$ is smoother than $\Phi$, such that it is computationally more advantageous to expand the former instead of the latter, since the results converge faster with respect to $N_{q}$. In Ref. Corsetti et al., 2013, a similar idea was applied to the kernel VV10 of Vydrov and Van Voorhis.Vydrov and Van Voorhis (2010)

The kernel $\Phi$ that we have considered for the present work is DRSLL.Dion et al. (2004) Extensions of Eq. (2) for spin-polarized systems were proposed in Refs. Obata et al., 2013; Thonhauser et al., 2015. However, since our results will be compared to the results from Refs. Klimeš et al., 2011; Schimka et al., 2013; Park et al., 2015, we followed the procedure in these studies which consists of simply using the sum of the up and down electron densities to evaluate Eq. (2).Klimeš et al. (2011) This concerns the calculations for solid Fe and Ni, and most free atoms.

Concerning the potential $v_{\text{c,disp}}^{\text{NL}}=\delta E_{\text{c,disp}}^{\text{NL}}/\delta\rho$, it was shown in Ref. Thonhauser et al., 2007 that adding $v_{\text{c,disp}}^{\text{NL}}$ to the semilocal component of the exchange-correlation potential for self-consistent calculations leads to no visible change in the total energy curve. In other words, the effect of $v_{\text{c,disp}}^{\text{NL}}$ on the orbitals and electron density is too small to affect properties calculated with the total energy. However, for the geometry optimization with the forces and stress tensor,Sabatini et al. (2012) the complete potential is required, therefore it is still desirable to have access to $v_{\text{c,disp}}^{\text{NL}}$. As explained in more detail in Appendix B, the proper calculation of $v_{\text{c,disp}}^{\text{NL}}$ requires the calculation of $d\rho_{\text{s}}/d\rho$, which is trivially done with Eq. (6) of the present scheme.

The LAPW calculations were done with the WIEN2k codeBlaha et al. (2001) and the parameters of the calculations like the basis-set size or number of $\bm{k}$-points for Brillouin zone integrations were chosen to be very well converged. For the RPS method, the subroutines available in the QUANTUM ESPRESSO codeGiannozzi et al. (2009); Sabatini et al. (2012) were used and modified. The FFT were done using version 3.3.5 of the FFTW software package,FFT which is efficiently parallelized with MPI.

III Results and discussion

The objective of this section is twofold. First, we will discuss the convergence behavior of the lattice constant and binding energy with respect to the density cutoff $\rho_{\text{c}}$ and plane-wave expansion cutoff $G_{\text{max}}$. Details will be shown for a few representative solids. Then, for a large set of solids (listed in Tables 1-3), we will compare our results with the values that were obtained with other implementations of the nonlocal vdW functionals. All results presented in this section were obtained with the functional optB88-vdW,Klimeš et al. (2010, 2011) whose semilocal component consists of a modification of the GGA B88Becke (1988) for exchange (optB88) and LDAPerdew and Wang (1992) for correlation, while the nonlocal term Eq. (2) uses the DRSLLDion et al. (2004) kernel $\Phi$. optB88-vdW is one of the functionals proposed in Refs. Klimeš et al., 2010, 2011 that have been shown to be more accurate for both molecules and solids than the older functionals vdW-DF1Dion et al. (2004) and vdW-DF2.Lee et al. (2010)

III.1 Convergence tests

As described in Sec. II, the cutoff $\rho_{\text{c}}$ in Eq. (6), which determines the degree of smoothness applied to $\rho$, has to be chosen. This is a tradeoff between computational efficiency and accuracy, since a smaller $\rho_{\text{c}}$ allows for a smaller $G_{\text{max}}$ (see below), but increases the deviation of the smooth density $\rho_{\text{s}}$ from the true density $\rho$ and, therefore, also the value of the vdW energy [Eq. (2)].

Starting the discussion with the convergence with respect to the density cutoff $\rho_{\text{c}}$, we have observed that for all cases in Tables 1-3 except Ne, a cutoff $\rho_{\text{c}}=0.3~{}$Bohr*-3* provides results for the geometry and binding energy that are already quite well converged. Compared to fully converged results (that are obtained with $\rho_{\text{c}}$ between 1 Bohr*-3* and 3 Bohr*-3* depending on the solid), the difference is below $\sim 0.005$ Å for the lattice constant and 0.04 eV/atom for the binding energy. The only exception is the very weakly bound rare gas Ne, since the converged lattice constant is about 0.02 Å larger than the value obtained with $\rho_{\text{c}}=0.3$ Bohr*-3*. Nevertheless, for vdW systems, an uncertainty in the lattice constant or bond length of the order of a few 0.01 Å is usually acceptable.

In the examples shown in Figs. 3-5 for the cubic solids Ca, Au, and GaAs, the lattice constants $a_{0}$ calculated with $\rho_{\text{c}}=0.3$ Bohr*-3* differ from the converged values by 0.003, 0.002, and 0.001 Å, respectively. For the cohesive energy $E_{\text{coh}}$, the error is the largest for Au (0.03 eV/atom). Also in the case of hexagonal BN (Fig. 6), the results for the lattice constant $c_{0}$ (the interlayer distance is $c_{0}/2$) and interlayer binding energy $E_{\text{b}}$ are very well converged with $\rho_{\text{c}}=0.3$ Bohr*-3*.

Thus, these results show that $\rho_{\text{c}}=0.3$ Bohr*-3* is a pretty good choice in terms of accuracy for all solids, including vdW systems. This shows that all-electron benchmark results can be obtained even with a smooth density $\rho_{\text{s}}$ which differs considerably from the true density in the region close to the nuclei (see Figs. 1-2).

At that point we should mention that Ref. Klimeš et al., 2011 reports a comparison between various schemes for the calculation of nonlocal vdW functionals. The goal was to estimate the accuracy of the PAW method (the density consists of a sum of pseudovalence and soft-core components) compared to all-electron results. A smoothing procedure is also used, however, no details are given, and also very little is mentioned about the procedure to get the all-electron results. The conclusion from this study is that the PAW results for the lattice constant and cohesive energy show good agreement with the all-electron results, similar to our conclusion, as shown below. However, note that no van der Waals systems were considered in Ref. Klimeš et al., 2011.

Concerning the convergence with the plane-wave cutoff $G_{\text{max}}$, we have observed that for most solids, $G_{\text{max}}=15$ Bohr*-1* is enough for the density cutoff $\rho_{\text{c}}=0.3$ Bohr*-3*. This is what is shown in Figs. 3-5 for Ca, Au, and GaAs. With $G_{\text{max}}=15$ Bohr*-1* the lattice constant and cohesive energy are usually converged within a few 0.001 Å and 0.01 eV/atom, respectively. However, in the case of weakly bound systems like h-BN (Fig. 6), a larger $G_{\text{max}}$ of 20-25 Bohr*-1* should be used for well converged results. In order to be on the safe side, a $G_{\text{max}}$ of about 20 Bohr*-1* is also recommended for the alkali metals, which are rather soft systems with a core-core vdW attraction that should not be neglected.Tao et al. (2010) As clearly shown in Figs. 3-6, the larger $\rho_{\text{c}}$ is, the larger $G_{\text{max}}$ should be in order to properly expand the density $\rho_{\text{s}}$.

From these convergence tests, it can be concluded that smoothing the all-electron density with a cutoff $\rho_{\text{c}}=0.3$ Bohr*-3* is quite safe in terms of accuracy. Furthermore, this value of $\rho_{\text{c}}$ is appropriate for the very light as well as for the very heavy atoms. With $\rho_{\text{c}}=0.3$ Bohr*-3*, a density plane-wave expansion in the range 15$-$20 Bohr*-1* has to be used (25 Bohr*-1* for the very weakly bound Ne), which seems to be similar to what may be needed to get converged results for weakly bound systems with pseudopotentials methods (see Ref. Sorescu et al., 2014).

III.2 Comparison with other codes

Turning now to the comparison with results from the literature, Klimeš et al. (2011); Schimka et al. (2013); Park et al. (2015); Tran and Hutter (2013); Callsen and Hamada (2015); Graziano et al. (2012); Björkman et al. (2012); Björkman (2014) Tables 1-3 show our converged (with respect to $\rho_{\text{c}}$ and $G_{\text{max}}$) optB88-vdW results for the lattice constant and binding energy along with results obtained with the PAW and Gaussian augmented plane wave (GAPW) methods as implemented into the VASPKresse and Furthmüller (1996) and CP2K codes,VandeVondele et al. (2005) respectively. Results obtained with the PBEPerdew et al. (1996) or LDA functionals are also shown in order to provide an idea how much (dis)agreement should be expected between the different codes.

The results in Table 1 for the strongly bound cubic solids show that the agreement between the WIEN2k and VASP codes is in general excellent. The differences in the lattice constant $a_{0}$ between the two codes are of the order of only a few 0.001 Å for most solids. We just note that in some cases, the various VASP results do not agree with each others, as for instance for Rb and Ge, where the differences are above 0.03 Å with optB88-vdW. In such cases, our WIEN2k benchmark results may help to indicate which one of the VASP results should be more correct. The WIEN2k and VASP results for the cohesive energy $E_{\text{coh}}$ are overall in quite good agreement as well, since the differences are typically of the order of a few $0.01$ eV/atom, which is very small. An exception is Au for which a large discrepancy is obtained with optB88-vdW (3.82 eV/atom for WIEN2k and 3.40 eV/atom for VASPPark et al. (2015)), while the PBE results differ by only $\sim 0.01$ eV/atom. Let us also mention that the two VASP results for Mo with PBE differ considerably (6.28 eV/atom from Ref. Schimka et al., 2013 and 7.76 eV/atom from Ref. Park et al., 2015). Such large disagreements in $E_{\text{coh}}$, as obtained for Mo and Au, may be due to different electronic configurations of the $d$-electrons in the isolated atoms.

As a side note, we mention that for the solids in Table 1, the omission of the vdW term [Eq. (2)] in the optB88-vdW functional leads to lengthening of the lattice constant and reduction of the cohesive energy that are rather substantial. Our calculations with the functional $E_{\text{xc}}=E_{\text{x}}^{\text{optB88}}+E_{\text{c}}^{\text{LDA}}$ (results not shown) lead to lattice constants that are larger by 0.05-0.1 Å for the transition metals, while for all systems containing alkali or alkaline earth atoms (except Li and MgO) $a_{0}$ is larger by 0.1-0.25 Å. For the cohesive energy, the values without vdW term are smaller by 0.3-0.6 eV/atom for the systems with alkali and alkaline earth atoms and 1-1.5 eV/atom for the transition metals. Therefore, not only vdW systems, but also those with supposedly unimportant vdW interactions are useful to test the implementation of functionals specifically designed for vdW systems. As shown in Ref. Klimeš et al., 2011, optB88-vdW is, compared to PBE, of the same accuracy for the lattice constant and slightly more accurate for the cohesive energy, such that it can be considered as rather good for solids. Without the vdW term, the lattice constant and cohesive energy are compared to experiment largely overestimated and underestimated, respectively.

Tables 2 and 3 show the results for the weakly bound rare-gas and hexagonal layered solids, respectively. For vdW systems, it is reasonable to tolerate uncertainties of a few 0.01 Å for the lattice constant and up to 10 or 20 meV/atom for the binding energy, depending on the system. From the results it can be inferred that the agreement between the various codes is rather good and, actually, on average the discrepancies do not seem to be larger for optB88-vdW than for LDA. Nevertheless, we note that for Ne, the lattice constant is noticeably larger (by 0.02-0.03 Å) with WIEN2k. As noticed above, among all systems that we have considered, Ne is the only one for which a density cutoff $\rho_{\text{c}}=0.3$ Bohr*-3* is not large enough to get a lattice constant that is within 0.01 Å of the converged value. With $\rho_{\text{c}}=0.3$ Bohr*-3*, $a_{0}=4.25$ Å which is closer to the VASP and CP2K values and would possibly indicate that the PAW and GAPW smooth densities plugged into Eq. (2) correspond more to our density with $\rho_{\text{c}}=0.3$ Bohr*-3*. In the case of the layered systems graphite and h-BN, Björkman (2012, 2014); Rêgo et al. (2015) the largest disagreements are for the lattice constant $c_{0}$ calculated with LDA, since the various VASP values differ by more than 0.15 Å which is quite large. Our LDA results for $c_{0}$ are in between the values from Refs. Graziano et al., 2012; Björkman et al., 2012. The optB88-vdW values for $c_{0}$ do not differ that much ($\sim 0.05$ Å), but the agreement between the codes is also not perfect. We mention that in our calculations, the intralayer lattice constant $a$ was kept fixed at the experimental value of 2.462 and 2.503 Å for graphite and h-BN, respectively. The same procedure was used in Ref. Björkman et al., 2012, while no details are given in Ref. Graziano et al., 2012. On the other hand, the interlayer binding energies $E_{\text{b}}$ calculated with WIEN2k and VASP agree quite well.

III.3 Further discussion

It was mentioned in Sec. II that the potential $v_{\text{c,disp}}^{\text{NL}}=\delta E_{\text{c,disp}}^{\text{NL}}/\delta\rho$ can be calculated (see Appendix B for details), which is necessary for the proper calculation of the forces acting on the nuclei. The unit cells of the systems that we have considered for our benchmark calculations do not contain free internal parameters. For such cases it does not matter how accurately the potential is calculated since it has basically no effect on the total energy curve.Thonhauser et al. (2007) A way to check the correctness of the implementation of $v_{\text{c,disp}}^{\text{NL}}$ is to consider the forces in systems with free internal parameters. We could check that the forces are correct, thus validating our implementation of the potential. Figure 7 shows the example of the Ne dimer, where we can see that the calculated forces agree very well with the numerical derivative of the total energy. The agreement is excellent in both cases, with or without the vdW term [Eq. (2)] in the functional optB88-vdW. According to some tests, consistency between the total energy and forces is also achieved with small plane-wave cutoff $G_{\text{max}}$ that maybe not be large enough for very well converged lattice constant for instance.

As an alternative to the smooth density $\rho_{\text{s}}$ given by Eq. (6), we may use the pseudocharge density generated by the method of WeinertWeinert (1981) for the calculation of the Coulomb potential in the LAPW method. In this method, the total (i.e., electronic plus nuclear) charge density inside the atomic spheres Sα is replaced by a pseudodensity $\rho_{\text{s}}^{\text{pseudo}}$ having the same multipole moments (and therefore producing the same Coulomb potential in the interstitial), but that is much smoother such that it can be expanded in plane waves. The only modification that has to be done for the present purpose is to remove the nuclear contribution. For comparison, the optB88-vdW lattice constants were also calculated using the pseudodensity, and the results (not shown) for most strongly bound solids are very similar to the reference results in Table 1 obtained with Eq. (6). The largest discrepancies, which are of the order of 0.02 Å, are found for the alkali metals (except Li) and Ag, and for these cases the lattice constant is usually larger with the pseudodensity $\rho_{\text{s}}^{\text{pseudo}}$ than Eq. (6). For the rare-gas solids, the lattice constants using $\rho_{\text{s}}^{\text{pseudo}}$ are noticeably larger by 0.03 (Ne), 0.05 (Ar), and 0.05 Å (Kr), while for graphite and h-BN, $c_{0}$ is larger with $\rho_{\text{s}}^{\text{pseudo}}$ by 0.01 and 0.03 Å, respectively, which can be considered as reasonably small. So, overall the use of the pseudodensity $\rho_{\text{s}}^{\text{pseudo}}$ leads to results which are quite close to the reference results, except for the rare gases and a couple of other systems.

Now, let us enumerate the pros and cons of both schemes, Eq. (6) and the pseudocharge method, to generate a smooth density. Equation (6) is trivial to implement and basis-set independent, while the pseudocharge method is a rather complicated method to implement and was designed specifically for the LAPW method.Weinert (1981) An unambiguous way to reach the converged all-electron results with Eq. (6) is to increase the density cutoff $\rho_{\text{c}}$, while it is not clear how to do it with the pseudocharge method (however, we admit that we have not checked the influence of the parameter $n$ in this methodWeinert (1981)). Finally, the third advantage with Eq. (6) is the possibility to calculate $d\rho_{\text{s}}/d\rho$ for the potential [see Eq. (12)], while it is at present not clear for us how to calculate $d\rho_{\text{s}}^{\text{pseudo}}/d\rho$ for the pseudodensity. Actually, if $\rho_{\text{s}}^{\text{pseudo}}\neq\rho_{\text{s}}^{\text{pseudo}}(\rho$), then this would be impossible even numerically. On the other hand, the advantage of the pseudocharge method is to lead to faster calculations, since due to the way $\rho_{\text{s}}^{\text{pseudo}}$ is constructed, it is sufficient to expand it with the plane-wave cutoff $G_{\text{max}}$ for the expansion of the Coulomb potential in the interstitial region ($G_{\text{max}}=12$ Bohr*-1* is the default in the WIEN2k code). The ideal method should have all advantages of Eq. (6) and produce converged results with a density plane-wave cutoff that is below 15 Bohr*-1* for all kinds of systems including weakly bound systems.

IV Summary and conclusion

In summary, a procedure for combining the efficient FFT-based method of RPS for nonlocal vdW functionals with all-electron methods has been presented and implemented into an LAPW code. It is based on the simple idea which consists of smoothing the all-electron density in the core region of the atoms, and then to use the resulting smooth density in the RPS method. There are obviously many ways to smooth an all-electron density, but the one that we have proposed [Eq. (6)] has several advantages: it is trivial to implement and contains a parameter (the density cutoff $\rho_{\text{c}}$) that gives us control for approaching all-electron benchmark results. Furthermore, the smoothing procedure is quite general in the sense that it is basis-set independent. However, the conversion of the potential from plane waves (obtained as output of the RPS method) to the required format will depend on the basis set. Nevertheless, most basis sets other than plane waves use spherical harmonics, such that the use of the Rayleigh formula for $e^{i\bm{G}\cdot\bm{r}}$ should be the solution in most codes.

A detailed study of the convergence of the results with respect to the density and plane-wave expansion cutoffs has shown that quite-well converged results can be obtained with a plane-wave cutoff that is reasonable and tremendously smaller than if the all-electron density was used. Interestingly, a density cutoff $\rho_{\text{c}}=0.3$ Bohr*-3* seems to be universally good across the whole periodic table of elements.

Then, very well converged results obtained with the optB88-vdW functional for a large set of strongly bound and vdW solids were compared to results obtained with other codes based on the PAW or GAPW methods. Overall, an excellent agreement between the various codes for the lattice constant and binding energy was obtained.

To conclude, we have shown that it is possible to use the efficient RPS method within an all-electron framework. Our work should also pave the way for proposing other methods based on the same idea. The pseudocharge density method for the Coulomb potential would also be a possibility, however, one would need to figure out how to calculate the potential.

Acknowledgements.

This work was supported by the project F41 (SFB ViCoM) of the Austrian Science Fund (FWF).

Appendix A Derivatives of $\rho_{\text{s}}$

At $\rho=\rho_{\text{c}}$, a derivative of Eq. (6) is continuous if the corresponding derivatives for $\rho>\rho_{\text{c}}$ and $\rho\leqslant\rho_{\text{c}}$ are equal. The first and second derivatives of $\rho_{\text{s}}$ are given by

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

From Eqs. (8)-(11), we can see that if $n\geqslant 1$, $\nabla\rho_{\text{s}}$ is continuous since $d\rho_{\text{s}}/d\rho=1$ at $\rho=\rho_{\text{c}}$. For $\nabla^{2}\rho_{\text{s}}$ it is the case if $n\geqslant 2$ since $d^{2}\rho_{\text{s}}/d\rho^{2}=0$. In general, for any value of $n\geqslant 1$ in Eq. (6), $d\rho_{\text{s}}/d\rho=1$ and $d^{m}\rho_{\text{s}}/d\rho^{m}=0$ for $m=2,3,\ldots,n$ when $\rho=\rho_{\text{c}}$, such that the $n$ first derivatives of Eq. (6) are continuous.

Appendix B Calculation of the potential $v_{\text{c,disp}}^{\text{NL}}$

In addition of being computationally efficient, the RPS method also leads to a straightforward calculation of the functional derivative of Eq. (2) [see Eq. (10) in Ref. Román-Pérez and Soler, 2009]. However, since in the present case $E_{\text{c,disp}}^{\text{NL}}$ is evaluated with a density ($\rho_{\text{s}}$) that is not the one ($\rho$) used to minimize the total energy, the chain rule has to be applied to get the correct expression for the potential:

[TABLE]

where $\delta E_{\text{c,disp}}^{\text{NL}}/\delta\rho_{\text{s}}$ is the part provided by the RPS method and $d\rho_{\text{s}}/d\rho$ is given by Eq. (10). Actually, it should be stressed that it was crucial to use a smoothing scheme that makes possible the derivation of $d\rho_{\text{s}}/d\rho$.

The second point concerns the conversion of $v_{\text{c,disp}}^{\text{NL}}$ (provided as a pure plane-wave expansion by the RPS procedure) into the LAPW format, namely (see Sec. II), as plane-wave and spherical harmonics expansions inside the interstitial region and atomic spheres, respectively:

[TABLE]

[TABLE]

where $Z_{LMp}$ are real spherical harmonics ($p=\{+,-\}$ for $M\geqslant 1$ or absent for $M=0$) that are defined as followsKara and Kurki-Suonio (1981):

[TABLE]

for $M=0$ and

[TABLE]

[TABLE]

for $M\geqslant 1$. Of course, the Fourier coefficients $v_{\text{c,disp}}^{\text{NL},\bm{G}}$ in Eq. (13) are the same as those obtained from the RPS method. The radial functions in Eq. (14) can be obtained by using the Rayleigh formula ($j_{\ell}$ are spherical Bessel functions)Arfken and Weber (2005)

[TABLE]

in Eq. (13), such that

[TABLE]

where $\tau_{\alpha}$ is the position of nucleus $\alpha$.