Dispersion-corrected PBEsol exchange-correlation functional

We study the compatibility between the PBEsol exchange-correlation energy functional of Phys. Rev. Lett. 100 , 136406 (2008) and the rVV10 van der Waals nonlocal correlation functional of Phys. Rev. B 87 , 041108 (2013). By applying a density-gradient dependence in the expression of rVV10, we develop a new functional: The PBEsol + rVV10s functional, which is remarkably accurate for layered solids, rare-gas crystals, densely packed bulk solids, and the adsorption of noble-gas atoms on metal surfaces, and is also competitive for noncovalent interactions between molecular complexes as well as for potential-energy curves of noble-gas dimers. The capacity of this dispersion-corrected functional to describe various interactions in solids makes it useful for electronic-structure calculations in solid-state physics and materials science.

One of the most popular GGA XC functionals is the PBEsol functional approximation [15,23], which is used in many electronic-structure calculations in condensed-matter physics and materials science. This approximated functional describes well the short-range correlation [25,26], showing important binding in layered materials [27]; but, as in the case of other semilocal functionals, the PBEsol functional completely misses the long-range van der Waals (vdW) interaction, which is of a nonlocal nature. On the other hand, the PBEsol functional is, in general, not compatible with available dispersion corrections. The main reason for this limitation is the behavior of the PBEsol exchange functional at large values of the reduced gradient s (s = |∇n|/[2(3π 3 ) 1/3 n 4/3 ]) where it is locally constrained by the Lieb-Oxford bound [28,29]. In order to solve this issue, various GGA exchange functionals have been constructed to be compatible with vdW correlation functionals. Here we mention the C09 and cx13 exchange functionals [30,31], which are close to the PBEsol exchange functional at small s but recover the large-s limit of the revPBE and PW86r exchange functionals [32], respectively. A comparison of their exchange enhancement factors [defined by F x = x / unif x with unif x = −3(3π 2 n) 1/3 /(4π )] is shown in Fig. 1. The corresponding van der Waals functionals C09-vdW-DF [30] and vdW-DF-cx [31,[33][34][35] are accurate for various interactions in molecules and solids [36][37][38].
Among most popular methods for adding dispersion effects to a semilocal functional, we find: (i) The DFT-D semiempirical approach proposed by Grimme [39], which is based on the introduction of a forcefield correction of the form Here, R represents the interatomic distance, s 6 is a scaling factor that only depends on the semilocal functional, f dmp (R) is a dumping function, and E LR vdW represents the long-range vdW interaction, C 6 being the dispersion coefficient for each atom pair. This method has been further improved [40][41][42][43] and applied to molecules [44] and solids [45,46]. The densitydependent model of Johnson and Becke [47] incorporates C 8 and C 10 as well. The PBEsol-D2 functional has been used in various applications [48][49][50][51][52] but not as often as the PBE-D2 functional.
(ii) The Rutgers-Chalmers approach [53], which is based on the introduction of a nonlocal correlation functional of the form where ε nl c [n] represents a correlation energy per particle, Here, the kernel (r, r ) = [n(r), n(r ), ∇n(r), ∇n(r ), |r − r |] is obtained in the framework of an adiabatic-coupling fluctuation-dissipation formalism by using a simplified treatment of the electronic polarizability. Several simplifications and improvements [54][55][56][57] have made this approach very accurate and efficient. In particular, the VV10 and rVV10 kernels [54] depend on the so-called b and c parameters. Based on these kernels, various approaches have been developed recently for the study of the properties of solid systems, e.g., the AM05-VV10sol [27], SCAN + rVV10 [58], PBE + rVV10L [59], and SG4 + rVV10m [60] functionals.
In this paper, we study the compatibility between PBEsol and rVV10 functionals. Upon using a simple gradient dependence of the parameter c, we construct a functional, the PBEsol + rVV10s functional, which is found to be accurate for noncovalent interactions in molecules and solids, being always comparable with and often better than the PBE + rVV10L functional. The paper is organized as follows: In Sec. II, we construct the PBEsol + rVV10 and PBEsol + rVV10s functionals. In Sec. III, we present the computational details of our calculations. In Sec. IV, we report the results we have obtained for layered materials, noble-gas solids, strongly bound solids, the adsorption of xenon atoms on (111) metal surfaces, and various molecular tests. In Sec. V, we summarize our conclusions.
In order to investigate the compatibility between PBEsol (semilocal XC functional) and rVV10 (nonlocal correlation functional) functionals, we consider two combined functionals: (i) The PBEsol + rVV10 functional, which consists of simply adding to the PBEsol functional the rVV10 functional E rV V 10 c,nl [n] dictated by Eq. (5) with c = 0.0093 (as before) but with a midrange parameter b changed from the original value (fitted to the S22 molecular test) to the new value of b = 20 that we obtain by fitting to various layered materials. This construction is similar to the PBE + rVV10L functional case, which was constructed with the parameters c = 0.0093 and b = 10.
(ii) The PBEsol + rVV10s functional, which consists of adding to the PBEsol functional the rVV10 functional E rV V 10 c,nl (1) c(s) → 0.0093 when s 2. We note that most of the systems where the vdW interaction is important have rapid fluctuations of the electronic density and, consequently, the reduced gradient can be large (s 2).
(2) The short-range correlation of the PBEsol functional is already accurate, and to avoid the double-counting problem we require c(s) 0.0093 at 0 < s 1. This condition attempts to remove part of the rVV10 nonlocal correlation in order to be compatible with the PBEsol functional.
(3) The case where 0.0093 < c(s) < max c(s) at s = 0. This case is present in the middle of any bonding region. We recall that, when s = 0, the PBEsol functional behaves as the local density approximation (LDA) [1].
These three conditions are fulfilled by the simple expression, where c 0 = 0.0093, c 1 = 0.5, and c 2 = 300. As for the parameter b, here we take the value of b = 10, which has been found by fitting to layered materials. In order to visualize these functionals, we consider two interacting parallel jellium slabs of bulk parameter r s = 2.07, each with a thickness of L = 5 a.u. The slabs are separated by the distance D (between the nearest edges). This model system was used by Dobson and Wang to show the van der Waals accuracy of the random-phase approximation (RPA) [66] and for the study of various exchange-correlation kernels in the framework of time-dependent DFT [67]. We use this system (with slabs that are narrow enough) in order to model layered materials.  (7) versus the distance z perpendicular to two interacting parallel jellium slabs of bulk parameter r s = 2.07. Each slab is L = 5-a.u. thick. The distance between the slabs is D = 7 a.u. Middle panel: Correlation energy per particle ε c for the PBEsol, PBEsol + rVV10, and PBEsol + rVV10s functionals versus the distance z for the same system. Lower panel: Interaction correlation energy (E system c − 2E slab c ) as a function of the distance between the slabs D (in atomic units). The RPA values are taken from Table I of Ref. [67]. The inset shows absolute relative errors (in percentages) of the correlation interaction energy with respect to the RPA versus D.
In the upper panel of Fig. 2, we show the reduced gradient s and the function c(s) of Eq. (7) versus the distance z when the slabs are separated by the distance D = 7 a.u. Inside the slabs, c ≈ 0.015 when s ≈ 0, being slightly bigger than c 0 .
Close to the surface of the slabs where the Friedel oscillations are present, c → c 1 = 0.5 such that the rVV10 term is dumped. This is a good feature because in this region there is no vdW interaction. Between the slabs, c ≈ c 0 with the exception of the middle distance where the PBEsol functional becomes a LDA.
In the middle panel of Fig. 2, we show the correlation energy per particle c (z) for the same jellium system. Nonlocal correlation yields a small positive correction [54], which is slightly more pronounced in the case of the PBEsol + rVV10s functional.
In the lower panel of Fig. 2, we report the interaction correlation energy (E system c − 2E slab c ) as a function of the slab distance D. The PBEsol functional correlation energy decays too fast, as expected, being at D = 10 a.u. several orders of magnitude smaller than the RPA (which is correct at D → ∞) with a relative error of 100% (see the inset). Our PBEsol + rVV10s functional approach yields, however, a remarkable agreement with the RPA (see also the inset) at distances D 3 a.u., whereas the PBEsol + rVV10 functional exhibits a less accurate behavior at long distances D. This is due to the fitted b value (b = 20), which makes this functional dump various vdW effects. We recall that the SCAN + rVV10 functional [58] has b = 15.7 (and the same c = 0.0093); but the SCAN meta-GGA is significantly more accurate than the PBEsol functional GGA for the study of the dispersion interaction in molecular systems.
We observe that at D = 0 there is an important difference between the RPA and the other shown functionals. Even if the RPA misses part of the short-range correlation, its result at D = 0 is in agreement with an energy-optimized exchangecorrelation kernel (see Table I of Ref. [67]). The RPA is also more accurate than the PBEsol functional for the atomization energy of molecules [68] but is slightly worse than the PBEsol functional for jellium surface correlation energies [69]. Due to this last point, we may expect that the PBEsol functional can be better than the RPA at D = 0. Note that the rVV10 correction worsens the PBEsol functional result; in this case there is no van der Waals interaction between the jellium slabs.
Finally, we mention that rVV10 is a two-body approximation that does not include multilayer screening effects [70], nontrivial long-range scalings of the vdW interaction [70][71][72], or nontrivial momentum dependence for two-dimensional materials [73,74]. On the other hand, the RPA can naturally account for many-body effects up to infinite order. In this sense, for a multilayer jellium system, the PBEsol + rVV10s functional may be less accurate than in the case of two interacting jellium slabs.

III. COMPUTATIONAL DETAILS
All calculations have been performed with QUANTUM ESPRESSO (QE) [75] with the use of ultrasoft pseudopotentials from the standard QE library. An energy cutoff of 120 Ry and a charge-density cutoff of 960 Ry were used for all materials under study in order to achieve well-converged results. Our Monkhorst-Pack k-point meshes have the following values: for SET27 and noble-gas solids, 16 × 16 × 16; for layered materials graphite and hexagonal boron nitride (h-BN), 16 × 16 × 8; for 2H where E bulk represents the energy of the bulk unit cell, E unit is the energy of the unit element (molecule, atom, or monolayer), and n is the number of unit elements in the bulk unit cell. In order to calculate the interaction energy for SET27 and noble-gas solids, atoms were placed in a box with dimensions 13 × 14 × 15 Å; for layered materials, the monolayer (with fixed relaxed bulk geometry) calculations were performed with a vacuum space of about 33 Å for graphite and h-BN, 50 Å for NbS 2 , and 37 Å for all others; for the molecular S22x5 set, a simple cubic cell of a 20-Å side was used for all the configurations; and for the noble-gas dimers, a box was used with the dimensions 16 × 17 × 18 Å.
To study the adsorption of Xe gas atoms on (111) metal surfaces (Pd, Pt, Cu, and Ag), we used a ( √ 3 × √ 3)R30 • supercell periodic structure with a vacuum space of about 18 Å along the stacking direction to avoid the interaction between the periodic images. A surface was modeled using repeated slabs containing six atomic layers. The Brillouin zone was sampled using a 8 × 8 × 1 k-point mesh. We fix the metal atoms in their relaxed bulk positions considering that the topmost surface layer for the Xe/metal system does not play a significant role in the interaction mechanism between the Xe adatoms and the metal surface as shown in previous theoretical investigations [76]. The adsorption energy was calculated by where E Xe-metal , E metal , and E Xe represent the total energy of the Xe + metal surface, the bare metal surface, and an isolated Xe atom, respectively.

A. Layered materials
Layered materials are nowadays of utmost importance in science and technology. Recently, more than 5000 layered materials have been identified. Among them, about 1000 materials are easily exfoliable with noncovalent layer-layer binding energies of less than 100 meV/Å 2 [77].
In Table I, we show the relative errors of layer-layer binding energy E b , interlayer lattice constant c, and intralayer lattice constant a for 28 layered materials. This test set has been also used in Refs. [58,59], being slightly larger than the benchmark of Ref. [27]. Overall, the PBEsol + rVV10s functional is the most accurate functional of Table I, being competitive with the best-known methods.

Layer-layer binding energy E b
The RPA correlation energy is known to be accurate for long-range interactions, including van der Waals interactions in metal slabs [66], layered materials [78,83,92], and noblegas solids [93]. The RPA can also correctly describe anion-π interactions [94] and surface asymptotics [95], but it misses short-range correlations in both three and two dimensions [96,97]. Moreover, the RPA is not, in general, accurate for the description of total, ionization, and atomization energies of molecules [68,98,99]. Hence, one needs to be cautious when judging functionals against the RPA layer-layer binding energy E b [100] especially when the noncovalent layer-layer interaction is not only of the van der Waals type.
Let us consider the well-known material graphite. As shown in Table I, dispersion corrected PBEsol functionals in this material underestimate the RPA layer-layer binding energy E RP A b by more than 30%. In Table II, we present various predictions of the graphite binding energy together with experimental results. From a comparison of quantum Monte Carlo calculations with experimental measurements, one may think that E b = 20 meV/Å 2 is the best estimate for the graphite binding energy so that all dispersion-corrected functionals are reasonably accurate for graphite as they all predict binding energies in the range of 11 E b 27 meV/Å 2 .
As for the remaining layered materials, the overall MAREs indicate that SCAN + rVV10, PBE + rVV10L, AM05 + VV10sol, and PBEsol + rVV10s functional calculations are the closest to the RPA reference values, whereas rVV10 and C09-vdW-DF overestimate the RPA reference values by more than 40%. Finally, we stress again that the RPA is not a gold standard for layered systems. Then a comparison for the equilibrium structure that is experimentally based may better reveal the accuracy of functionals.

Lattice constants a and c
In the upper panel of Fig. 3, the relative error in the calculation of the interlayer lattice constant c is presented versus the relative error in the calculation of the intralayer lattice constant a for several functionals: PBEsol + rVV10s, PBE + rVV10L, SCAN + rVV10, AM05 + VV10sol, and vdW-DF-cx functionals. The best performance is given by the vdW-DF-cx functional, with a MARE of 0.9% for both lattice constants (see Ref. [36]). With the exception of PtS 2 and PtSe 2 , the second best performance is achieved with the use of the new PBEsol + rVV10s functional, being even superior to SCAN + rVV10. This figure shows that errors in the calculation of the lattice constants a and c are clearly correlated [105]. In order to understand the relatively large errors that are made in the case of PtS 2 and PtSe 2 , in the lower panel of Fig. 3 we show the results one obtains from various functionals. The LDA strongly underestimates the interlayer lattice constant c for both PtS 2 and PtSe 2 . The PBEsol functional improves essentially with respect to the LDA, underestimating c by about 5% and 4% for PtS 2 and PtSe 2 , respectively. Adding a rVV10-like dispersion correction to the PBEsol functional has the effect of a further underestimation of c, i.e., PBEsol + rVV10 and PBEsol + rVV10s functionals are both worse than the PBEsol functional, although the PBEsol + rVV10s functional is better than the PBEsol + rVV10 functional; indeed, the PBEsol + rVV10s functional is close to AM05 + VV10sol. PtS 2 and PtSe 2 both represent difficult materials for all the GGA functionals that are constructed for bulk solids, such as the PBEsol, AM05, even C09-vdW-DF [30], and vdW-DF-cx functionals. Figure 3 shows that PtS 2 and PtSe 2 are both outliners for the PBEsol + rVV10s functional data; indeed, by eliminating them we obtain a MARE of 1.0% for the PBEsol + rVV10s functional interlayer lattice constant c, lower than, e.g., SCAN + rVV10 which gives a MARE of 1.34% [58].

B. Noble-gas solids
Noble-gas solids are pure van der Waals materials, and they are (unlike layered materials) completely isotropic and highly symmetric (face-centered-cubic structure). Hence, they represent an important benchmark test for dispersion-corrected functionals [107,108].
In Table III, we show relative errors in the calculation of equilibrium lattice constants and cohesive energies with the use of a variety of functionals. Remarkably, the PBEsol + rVV10s functional is found to be among the most accurate functionals, competing with the RPA@PBE functional and performing significantly better than the PBE + rVV10L and SCAN + rVV10 functionals.
In Fig. 4, the cohesive energy is exhibited as a function of the lattice constant for Ne, Ar, and Kr, as obtained with the use of the PBEsol, PBEsol + rVV10, and PBEsol + rVV10s functionals (see also Fig. 2 of Ref. [107]). The CCSD(T) results [109] are also represented for comparison. In the region near equilibrium, the PBEsol + rVV10s functional is found to be more realistic than the PBEsol and PBEsol + rVV10 functionals. When the lattice constant is stretched, the performances of the PBEsol + rVV10s and PBEsol + rVV10 functionals are comparable. The same behavior is found for molecules as shown in Sec. IV E below.

C. Strongly bound solids
Here we use the SET27 benchmark [60,110], which contains simple metals, transition metals, semiconductors, ionic crystals, and insulators. In Table IV, we show that for strongly bound solids both PBEsol + rVV10 and PBEsol + rVV10s functionals preserve the high accuracy of the PBEsol functional in the calculation of structural properties (lattice constants and bulk moduli), whereas their performance is slightly worse in the calculation of cohesive energies (with the exception of ionic crystals for which the rVV10 correction is beneficial). We note that the vdW-DF-cx performs better than the PBEsol functional for the class of nonmagnetic transition metals for thermophysical properties [38]. Figure 5 exhibits the correlation between the relative errors (in percentages) of lattice constants and bulk moduli. The use of the PBEsol + rVV10 and PBEsol + rVV10s functionals yield results that are always close to the results one obtains by using the PBEsol functional. Nonetheless, the inclusion of dispersion corrections (in the PBEsol + rVV10 and PBEsol + rVV10s functionals) brings an underestimation of the lattice constant and an overestimation of the bulk modulus. We note that there are no points (for the PBEsol, PBEsol + rVV10, and PBEsol + rVV10s functionals) in the first quadrant of Fig. 5, which means that within the benchmark set under study there is no material for which these functionals overestimate both the lattice constant and the bulk modulus.

Lattice
In the case of the equilibrium distance d, the use of the PBEsol, PBE + vdW, and PBE + vdW surf functionals yields accurate results in close agreement with experiment. The use of vdW-DF1 and vdW-DF2, however, yields an overestimation of the equilibrium distance d by ∼22% for Xe/Pd(111), ∼17% for Xe/Pt (111), and ∼11% for Xe/Cu (111). The present PBEsol + rVV10 and PBEsol + rVV10s functional calculations yield equilibrium distances that are only slightly compressed with respect to the more accurate PBEsol functional distances by less than 5%.
In the case of the adsorption energy E ad , the PBEsol + rVV10s functional yields the best results, describing accurately all systems under study (see Table V). The worst results for this quantity are obtained with the use of the PBEsol and DFT-D2 functionals. Overall, the present PBEsol + rVV10s functional is one of the most promising functionals shown in Table V.

S22x5 set
The S22x5 set [62] of intermolecular interaction energies of noncovalently bonded complexes represents a generalization of the popular S22 set [63,125]. The S22x5 set consists of five subsets (denoted S22x5-d) with d being the relative interaction distance as compared to S22: In S22x5-0.9 the complexes are compressed with respect to S22, S22x5-1.0 is exactly S22, and in S22x5-2.0 the complexes are considerably stretched with respect to S22.
In the case of S22x5-1.0 (the regular S22), the PBEsol + rVV10s functional is found to be the most accurate dispersioncorrected GGA for solids, being better than AM05-VV10sol, SG4 + rVV10m, PBE + rVV10L, BEEF-vdW, and PBEsol + rVV10 functionals. This encouraging result suggests that the PBEsol + rVV10s functional can be used in calculations of hybrid interfaces (e.g., interactions between molecular complexes and metal clusters). Finally, in the case of S22x5-1.2, S22x5-1.5, and S22x5-2.0 all functionals have MAEs below 1.5 kcal/mol, so they do not fail badly even when the MARE is relatively large. For example, vdW-DF2 has a MARE of 34% for the S22x5-2.0 set. For the S22x5-2.0 test, the interaction energy is extremely small, and a correct description must include RPA-like manybody terms in order to obtain the long-range scaling of vdW interactions [70,128].

Noble-gas dimers
In Table VII, we report the accuracy of several XC functionals for the calculation of equilibrium bond lengths of the noble-gas dimers Ne 2 , Ar 2 , and Kr 2 . The PBEsol + rVV10s functional is found to yield the best results, systematically improving over the results obtained with the use of the PBEsol + rVV10 functional. The worst results are obtained with the use of C09-vdW-DF, which yields a large overestimation of the bond length as in the case of noble-gas solids (see Table III).
In Fig. 6, we show the potential curves of the noble-gas dimers Ne 2 , Ar 2 , and Kr 2 . In the case of Ne 2 , although the PBEsol functional is close to exact the dispersion-corrected functionals, the PBEsol + rVV10 and PBEsol + rVV10s functionals clearly overestimate the magnitude of the interaction energy. In the case of Ar 2 and Kr 2 , however, the best result is obtained with the use of the PBEsol + rVV10s functional, although still underestimating (as with the use of the other functionals under study here) the magnitude of the interaction energy. Note that a similar underbinding was obtained for noble-gas crystals (see Fig. 4).

V. CONCLUSIONS
We have investigated the compatibility between the PBEsol XC energy functional and the rVV10 van der Waals nonlocal correlation functional. By applying a density-gradient dependence in the expression of rVV10, we have developed a new functional: the PBEsol + rVV10s functional with the parameters b = 10 and c = c(s) given by Eq. (7). We have found that, although both PBEsol + rVV10 (with b = 20 and c = 0.0093) and PBEsol + rVV10s functionals perform at the level of the PBEsol functional for strongly bound solids, they improve considerably over the PBEsol functional in the case of noncovalent interactions. Our test calculations for layered materials, noble-gas solids, the adsorption of a Xe atom on a metal surface, molecular complexes, and rare-gas dimers show that the density-gradient-dependent PBEsol + rVV10s functional is superior to the PBEsol + rVV10 functional and, therefore, our construction for c = c(s) is realistic.
Due to its performance for the structural and energetical properties of bulk solids and hybrid interfaces (i.e., Xe/metal surface), we conclude that the present PBEsol + rVV10s functional can be used in many applications especially for electronic-structure calculations in the framework of solidstate physics and materials science.
Further improvements of our functional may be feasible by considering both parameters of the rVV10 nonlocal correlation functional to be dependent on the reduced density gradient: b = b(s) and c = c(s). We will address this idea in future work. Nevertheless, although the PBEsol + rVV10s functional is easily implementable in any code that contains the original rVV10, the functional derivative of a van der Waals functional with b = b(s) will be significantly more complicated.