Versatile van der Waals Density Functional Based on a Meta-Generalized Gradient Approximation

A “best-of-both-worlds” van der Waals (vdW) density functional is constructed, seamlessly supplementing the strongly constrained and appropriately normed (SCAN) meta-generalized gradient approximation for shortand intermediate-range interactions with the long-range vdW interaction from rVV10, the revised Vydrov–van Voorhis nonlocal correlation functional. The resultant SCANþrVV10 is the only vdW density functional to date that yields excellent interlayer binding energies and spacings, as well as intralayer lattice constants in 28 layered materials. Its versatility for various kinds of bonding is further demonstrated by its good performance for 22 interactions between molecules; the cohesive energies and lattice constants of 50 solids; the adsorption energy and distance of a benzene molecule on coinage-metal surfaces; the binding energy curves for graphene on Cu(111), Ni(111), and Co(0001) surfaces; and the rare-gas solids. We argue that a good semilocal approximation should (as SCAN does) capture the intermediate-range vdW through its exchange term. We have found an effective range of the vdW interaction between 8 and 16 Å for systems considered here, suggesting that this interaction is negligibly small at the larger distances where it reaches its asymptotic power-law decay.


I. INTRODUCTION
In 2004, graphene, the first two-dimensional (2D) material, was experimentally realized [1], triggering the renaissance of layered materials in both condensed matter physics and materials science.In the same year, the first general-geometry van der Waals (vdW) density functional was devised [2], followed by a series of efforts [3][4][5][6][7][8][9][10] to include the vdW interaction within the framework of Kohn-Sham density functional theory (DFT) [11], the current mainstream first-principles approach.Layered materials including graphene, h-BN, transition-metal dichalcogenides, black phosphorus, etc., have demonstrated various new concepts and potential technical applications [12][13][14][15].By contrast, despite much progress, the effort continues to find a vdW density functional with an acceptable accuracy for layered materials [16][17][18][19].For example, none of the density functionals tested in Refs.[17] and [19] can simultaneously reproduce the high-level theoretical interlayer binding energy, as well as the experimental interlayer and intralayer lattice constants with a satisfying accuracy for 28 layered materials.An accurate description of the layerlayer vdW interaction is self-evidently important for studying the physical properties of layered materials, especially the evolution of property as a function of configuration.Exploring the interactions between the layered materials and the environment, such as substrates and molecular absorbates, further requires accurate treatment of the vdW and other kinds of chemical bonding on the same footing.In short, a better versatile vdW density functional has been long-awaited and is urgently needed.
The difficulty in describing the layered materials arises from the nonlocal and long-range nature of the vdW interaction, as well as the coexistence of the weak vdW bonding and much stronger chemical bonding therein.The vdW interaction originates from dynamic electron correlations, causing a net attraction between fragments [20].It has minor effects in most bulk systems, becomes noticeable in some bulk solids like soft alkali metals [21], and is significant in sparse matter including molecular complexes, molecular crystals, layered materials, and surface-adsorbate systems, as well as the so-called "soft matter."Fully accounting for the vdW interaction is achievable by highlevel methods such as quantum Monte Carlo (QMC) [22], coupled-cluster singles and doubles with perturbative triples [CCSD(T)] [23], and adiabatic-connection fluctuationdissipation theorem within the random-phase approximation (RPA) [24], which are, however, only feasible for limitedsize systems because of high computational cost.
Aiming for a better efficiency, various approaches have been proposed within the framework of DFT.The popular ones include the DFTþD series [25,26], Tkatchenko-Scheffler (TS) methods [27][28][29], the Rutgers-Chalmers vdW-DF family [2][3][4][5][6][7][8], as well as the VV10 [10] and the rVV10 [30] methods.DFTþD2 and its predecessors correct the DFT total energy by adding an empirical atom-pairwise interaction correction, parametrized by atomic C 6 coefficients.The DFTþD3 and TS methods improve over them by taking account of the dependence of the vdW coefficients on the chemical environment.The vdW density functionals (vdW-DFs, VV10, and rVV10) require only the electron density and its first derivative as inputs and hence are conceptually applicable to any chemical environment with general geometry.In these methods, the total exchange-correlation (xc) energy consists of the semilocal xc and nonlocal correlation components: Comparing with the Rutgers-Chalmers vdW-DFs, the VV10 and rVV10 methods are more flexible in form thanks to two empirical parameters, C and b.The original VV10 [10] takes [31,32], C ¼ 0.0093, and b ¼ 5.9.The rVV10 method differs from VV10 only by a slightly revised form, together with a different b ¼ 6.3 [30].With such settings, VV10 and rVV10 describe the tested molecular systems very well [10,30] and reproduce remarkably well the experimental intralayer and interlayer lattice constants for 28 layered materials [17].However, VV10 and rVV10 overestimate the interlayer binding energies with respect to the RPA results by around 40%, although in a preferred systematic way [17].The interlayer binding energies from VV10 and rVV10 can be improved by a refitted b parameter (b ¼ 9.15 for both), however, at the undesired loss of accuracy for the lattice constants of layered materials and for other systems like molecular complexes [18].
Here, we propose starting with an accurate meta-GGA [33] that includes even intermediate-range vdW interactions, and then adding a long-range vdW correction [30] in which the intermediate-range contribution is tuned out.Note that a meta-GGA, without an intermediate-range contribution [34], has been corrected at intermediate and long ranges in Ref. [34] (as fitted and tested only for molecules).
Recently, the "strongly constrained and appropriately normed" (SCAN) meta-GGA [33] was constructed based on all 17 known exact constraints appropriate for a semilocal functional, and a set of "appropriate norms" for which a semilocal function can be exact or nearly exact, like slowly varying electron densities and compressed Ar 2 .Presenting remarkable accuracy and versatility for the structural and energetic properties in strongly bonded systems, SCAN also gets the intermediate-range vdW interaction about right, which was otherwise severely overestimated by the local density approximation (LDA) and severely underestimated by other nonempirical semilocal functionals.(The intermediate range is roughly the distance between nearest-neighbor atoms at equilibrium.)This has been demonstrated by the improved interaction energies for the S22 (cf.Table IV) data set [33] and by better geometric and energetic properties of ice in different phases [35].With SCAN for E 0 xc and appropriate parametrization, the VV10 or rVV10 E nl c will take effect only for the long-range vdW interaction and keep the advantages of SCAN intact.In this study, we demonstrate this concept with the rVV10 method, which allows for a fast implementation [36].The resultant SCANþrVV10, with a single b parameter and fixed C, turns out to be a successful versatile vdW density functional, with a high accuracy for different systems, including molecular complexes, bulk solids, benzene adsorbed on metal surfaces, and the more challenging layered materials.(Computational details are given in Appendix A.) While some of the SCAN exact constraints could be violated by SCANþrVV10, the violation would necessarily be small.
It is now common to pair a nonlocal vdW functional with a semilocal exchange-correlation functional but in a way that to us seems to make "the tail wag the dog": Start from the nonlocal vdW functional and then add a GGA chosen to minimize its "contamination" from intermediate-range vdW.Our approach instead is to start from a meta-GGA that can recognize and describe diverse bonds, including intermediate-range vdW, and then add an exclusively longranged nonlocal vdW functional.Underlying the conventional approach is the belief that a semilocal density functional should not capture an intermediate-range vdW interaction, or at least that it should not capture this contribution to the correlation energy through its exchange term.In Sec.IVA, we argue that this belief is incorrect, and we show how satisfaction of a tight lower bound on the exchange energy, as in SCAN, can lead to a correct and understandable description of intermediate-range vdW.

II. THEORY: PARAMETERS IN SCANþrVV10
The rVV10 nonlocal correlation functional takes a similar form as the Rutgers-Chalmers vdW-DFs, where nðrÞ is the electron density, and Φðr; r 0 Þ is the kernel describing the density-density interactions.To ensure zero E nl c for the uniform electron gas, 3 4 in Hartree is required for VV10 and rVV10.(β is not needed for Rutgers-Chalmers vdW-DFs.)The two empirical parameters C and b appear in the kernel: C chosen for accurate −C 6 =R 6 vdW interactions between molecules at large separation R, and b controlling the damping of E nl c at short range.For a semilocal E 0 xc , C ¼ 0.0093 was recommended [10], and the b parameter was determined as 5.9 and 6.3 by fitting to the interaction energies of the S22 set [37,38] for the original VV10 and rVV10 [10,30].Increasing C or b generally results in smaller vdW correction.Björkman et al. further proposed b ¼ 9.15 for VV10 by fitting to the binding energies of 26 layered materials [18], and we obtained the same number for rVV10.Reference [39] found b ¼ 9.3 for the structural properties of water.These results highlight the flexibility of the VV10 and rVV10 methods but also point out the need for a new approach with a single set of vdW empirical parameters for all materials.(In the following, we use "rVV10" to specifically denote the original rVV10 density functional with 0093, and b ¼ 6.3 [30].)Here, we determine the b parameter for SCANþrVV10 by fitting to eleven data points around the equilibrium bond length from the CCSD(T) binding curve of the Ar dimer [40].(One point with a binding energy close to zero was excluded.)Note that three data points on the repulsive wall have been used as an "appropriate norm" for constructing SCAN [33].With this choice, the more diverse S22 data set is reserved for benchmarking in the following, and the computational cost for the fitting is largely reduced.In the end, we obtain a mean absolute relative error of 2.6% for these selected data points with b ¼ 15.7.The binding curve for Kr 2 was also computed with the same parametrization.In Fig. 1, the two binding curves obtained from PBE, rVV10, SCAN, and SCANþrVV10 are compared with the CCSD(T) results [40,41].Varying the value of C does not noticeably improve the binding curve, and hence we stick to the recommended value of 0.0093.
Such a large b parameter in SCANþrVV10 is not unexpected.From Fig. 1, as well as previous SCAN calculations on S22 [33] (also cf.Fig. 2), noticeable improvement over PBE by SCAN itself for molecular systems is clear.With the kinetic energy density as an extra input, meta-GGA can recognize covalent, metallic, and weak bonds [42], so it is conceptually not surprising that some meta-GGAs can include-to some extent-the intermediate-range vdW interactions.For example, the M06-L [43] has a performance comparable to that of SCAN for the S22 data set.However, the M06-L has molecular systems in its construction, while SCAN is not fitted to any bonded system.Constrained by the 17 exact constraints and guided by a set of appropriate norms, the resulting mathematical form of SCAN achieves excellent control of error cancellation between semilocal exchange and correlation, not only for covalent bonds as LDA and GGA do, but also for the intermediate-range vdW interaction, although the resulting binding still comes from the exchange part instead of the correlation [44,45].A similar situation occurs in covalent molecules, where semilocal exchange rather reliably imitates exact exchange plus leftright correlation [46,47].(A more detailed explanation is given in Sec.IVA.)With such improvement for the intermediate-range vdW interaction, SCAN requires a smaller correction from E nl c and therefore a larger b parameter in the rVV10 kernel.The large b parameter, in turn, keeps intact the excellent capability of SCAN to describe non-vdW systems, leaving SCANþrVV10 as a "best-of-both-worlds" method.

III. RESULTS: PERFORMANCE OF SCANþrVV10
A. SCANþrVV10 for molecular complexes We first evaluate the performance of SCANþrVV10 for molecular complexes by calculating the interaction energies for the S22 set, which includes seven hydrogenbonded, eight dispersion-bound, and seven mixed complexes.In Fig. 2, we compare the absolute relative errors and relative errors for the interaction energies from PBE, rVV10, vdW-DF2, SCAN, and SCANþrVV10, with the CCSD(T) results [37,38] as the reference (raw data provided in Appendix B).With the large b ¼ 15.7 parameter from Ar 2 binding curve fitting, SCANþrVV10 approaches the accuracy of the original rVV10 (for which b ¼ 6.3 was fitted to these data) for these molecular systems, better than the vdW-DF2 (numerical results from Ref. [10]), which underestimates the interaction energies for S22 by about 15% on average.

B. SCANþrVV10 for bulk solids
To benchmark the performance of SCANþrVV10 for systems where vdW has little or only a slight effect, we compiled a set of 50 solid systems studied previously in Refs.[48,49], which includes (i) 13 group-IV and III-V semiconductors, (ii) 5 insulators, (iii) 8 main-group metals, (iv) 3 ferromagnetic transition metals, Fe, Co, and Ni, and (v) 21 other transition metals for which non-spin-polarized calculations were performed.For this set of solids, atomization energies and lattice volumes from RPA, and the corresponding experimental values after the zero-point correction, are available [48,49].In Fig. 3, the comparison between the RPA, PBE, SCAN, and SCANþrVV10 is given (raw data provided in Appendix B).For the atomization energy, both SCANþrVV10 and SCAN are only slightly better than PBE and RPA.(However, atomization energy may not be a good choice to assess a semilocal functional [50].)For the lattice volume, we found that SCANþrVV10 is essentially as good as RPA, it behaves similarly to SCAN for most solids, and it is much better than the PBE functional that overestimates, with the mean absolute relative error about 3%.SCAN overestimates the volume by 6%-10% for K, Rb, and Cs (the three outliers) because of the long-range vdW effect in these soft alkali metals.This vdW effect originates from the long-range attraction between semicore p orbitals [21], which SCAN fails to describe.Not unexpectedly, SCANþrVV10 properly improves the volumes for these three and does not skew the distribution for the other solids as desired.

C. SCANþrVV10 for benzene adsorbed on metal surfaces
The outstanding performance of SCANþrVV10 for the molecular and solid systems above is encouraging.To be a successful general-geometry density functional, SCANþrVV10 should pass more stringent tests, where both strong local atomic bonds and weak vdW interactions need to be well described.We first consider the widely studied benzene ring adsorbed on the (111) surface of coinage metals [51][52][53][54][55][56].This is also one of the most common systems for benchmarking a vdW functional's capability to simultaneously describe the vdW bonding and the metallic bonding.At first, we want to emphasize that the lattice constants from SCANþrVV10 are 3.544, 4.058, and 4.073 Å for Cu, Ag, and Au, in excellent agreement with the experimental values of 3.595, 4.062, and 4.062 Å.These results are much better than the vdW-DF1 methods with different exchange density functionals [54].In Table I, we compare the calculated adsorption energy E ad and distance Δz between the benzene and the surface from SCAN and SCANþrVV10, with available experimental values [29,[57][58][59][60][61].SCANþrVV10 agrees well with experiment for E ad , while SCAN systematically underestimates it by about 0.4 eV.The experimental values for the distance between benzene and the surfaces were measured for Ag [29] and estimated for Au [61], both of which are close to our SCANþrVV10 results.FIG. 2. Box plots for the absolute relative errors and relative errors of the interaction energies from PBE, rVV10, vdW-DF2 (numerical results from Ref. [10]), SCAN, and SCAN þ rVV10 with respect to the CCSD(T) results [37,38], for the molecular dimers in the S22 data set.The PBE errors, too big to display on this scale, have been scaled down by a factor of 3 in this figure.
The Tukey box plot used here summarizes the overall distribution of a set of data points: The bottom and the top of the box are the first (Q1) and third (Q3) quartiles (25% of data points lies below Q1, and another 25% above Q3); the band inside the box denotes the median; the circles, if any, denote outliers that lie further than 1.5 Ã jQ3 − Q1j away from the box; the vertical line extends from the minimum to the maximum, except for the outliers.We also denote the mean value with a shape inside the box.[53], which still overestimates the adsorption for the case of Cu and apparently cannot apply to general geometries.The recent rev-vdW-DF2 [7] is close to SCANþrVV10 based on its results for benzene on Cu(111), but more thorough calculations are needed for confirmation.Since the adsorption of organic molecules on metallic surfaces plays an important role in catalysis, molecular sensors, and switches, etc., SCANþrVV10 will be a very helpful tool in these fields.Benzene absorption on Pt(111) surfaces is another widely studied system, both experimentally [62][63][64] and theoretically [53,54,65].Unlike the case for the coinage metal surfaces, strong chemisorption occurs between carbon (hydrogen) atoms and Pt(111) surface because of the covalent contribution.The range of theoretical adsorption energies is quite large.When a molecule adsorbs in this way on a metal surface, there could be a spurious charge transfer creating a "density-driven error" [66] in the adsorption energy from any semilocal functional, which could be removed, for example, by a self-interaction correction.We aim to investigate this possibility in future work.

D. SCANþrVV10 for graphene on Ni, Co, and Cu
The binding energy curve for a graphene sheet in registry with a Ni(111) [67][68][69], Co(0001) [68,69], or Cu(111) [68,69] surface provides possibilities for both physisorption and weak chemisorption.(The assumed registry is reasonable for chemisorption but not at large separations; see Table I of Ref. [69].) Figure 4 compares the binding energy curves calculated with PBE, SCAN, SCANþrVV10, and RPA@PBE (the random-phase approximation with PBE orbitals and orbital energies).Reference [67] shows that PBE yields no net charge transfer between the graphene sheet and the Ni substrate, suggesting that the density-driven errors discussed in the previous paragraph may be absent for these systems.
Experimental information for these systems is scarce.LEED [70] suggests that the equilibrium position of the graphene sheet on Ni is 2.1 Å above the outermost plane of Ni nuclei, as found by all the methods shown in Fig. 4. PBE certainly underbinds, and the binding increases from PBE to RPA to SCAN to SCANþrVV10.At large separations, RPA and SCANþrVV10 show a long-range vdW, while PBE and SCAN do not.For graphene on Cu, RPA and SCANþrVV10 agree well in their description of a single physisorption minimum.
For graphene on Ni and Co, SCANþrVV10 excellently reproduces the long-range part of the binding curve from RPA, while the chemisorption is stronger; hence, the physisorption looks more like a shoulder than a minimum.It is often assumed that the RPA binding energy curve is a reliable reference.Certainly, RPA is the highest-level approximation in Fig. 4. But, as Ref. [68] points out, "RPA tends to underestimate atomization energies for molecules."Indeed, the atomization energies of molecules are much more accurate in SCAN and SCANþrVV10 than in RPA."One could thus expect that the physisorption minimum for graphene on metals is rather well described by RPA, whereas the chemisorption minimum is most likely underestimated" [68].Without further information, we cannot say whether RPA or SCANþrVV10 is more correct.

E. SCANþrVV10 for rare-gas solids
All binding in the rare-gas solids comes from the vdW interaction.Table II shows the equilibrium lattice constants for Ne, Ar, and Kr in the LDA, PBE, SCAN, and SCANþrVV10, in comparison with experiment and CCSD(T) [71].LDA consistently overbinds, and PBE consistently underbinds.SCANþrVV10 overbinds Ne (because SCAN does), but not as much as LDA, and is excellent for Ar and Kr (just as it is for the dimers Ar 2 and Kr 2 , but not Ne 2 [72]).The rVV10 correction to SCAN is important.For example, it doubles the binding energy of solid Ar.

F. SCANþrVV10 for layered materials
Layered materials have become one of the main arenas for condensed matter physics as well as materials science during the last decade [12][13][14][15].When embarking on firstprinciples exploration of known or unknown 2D materials [17,73], three of the most fundamental quantities have to be correctly predicted: the interlayer binding energy (E b ), the interlayer lattice constant (c), and the intralayer lattice constant (a).The first two are mainly determined by vdW interactions, and the last one by the stronger covalent or ionic bonding.Björkman et al. [17,19] have tested the performance for a series of vdW density functionals, using E b from RPA calculation and experimental lattice constants as the references for 28 hexagonal or tetragonal layered materials.One conclusion is that no single functional yet can yield a mean absolute relative error <10% for E b , <2% for c, and <1% for a.Among the tested functionals, the recently proposed rev-vdW-DF2 [7] was found to be the best for the lattice constants, together with a relatively small mean absolute relative error of 16% for E b .In Fig. 5, we compare the performance of this functional, SCAN, and SCANþrVV10 for 28 layered materials for which the RPA E b are available [17,19] (raw data provided in Appendix B). (Here, we found the mean absolute relative error of E b from rev-vdW-DF2 to be close to 26%, which may arise from the different implementation.)As expected, SCAN predicts an accurate intralayer lattice constant for which vdW has little effect (there is one outlier from VTe 2 with a underestimated by 3%, which may be due to the self-interaction error).However, SCAN does strongly underbind along the c direction, as illustrated by the systematically nearly 60% underestimated E b , and the largely scattered absolute relative errors of c.With the b parameter determined by the Ar 2 binding curve, the mean absolute relative error for c from SCANþrVV10 drastically reduces to 1.4%, comparable with that from rev-vdW-DF2, and the mean absolute relative error for E b dramatically reduces to 8%, only about 1=3 of that from rev-vdW-DF2.As expected, the long-range rVV10 correction is much more important for the interlayer binding TABLE II.Calculated lattice constants and binding energies for the fcc rare-gas solids from LDA, PBE, SCAN, and SCAN þ rVV10, compared with the experimental and CCSD(T) values [71].

Ne
Ar Kr  energy than it is for the S22 binding energies of small molecular complexes.Meanwhile, the mean absolute relative error for a is as small as 0.6%, the same as those from SCAN and rev-vdW-DF2.Hence, SCANþrVV10 becomes the first functional reducing the mean absolute relative error <10% for E b , <2% for c, and <1% for a.This significant success suggests that SCANþrVV10 can be the method of choice for layered materials, especially in the case of large-scale computation.RPA calculations are currently the best-available benchmark for layered materials, although RPA is not an exact method [19].Very recently, the binding energy of black phosphorus was calculated with the more reliable QMC method [74].Black phosphorus has an orthorhombic layered structure consisting of puckered layers, with experimental lattice constants of a ¼ 3.313, b ¼ 4.374, c ¼ 10.473 Å (the interlayer one), and the layer-layer binding energy is 81 meV=atom from QMC [74].Our SCANþrVV10 agrees with these values very well, reproducing the intralayer lattice constants with a ¼ 3.29, b ¼ 4.45 Å, the interlayer lattice constant with c ¼ 10.599 Å, and interlayer binding energy with 82.9 meV=atom (all calculated with full geometric relaxation).

How can SCAN describe the intermediate-range attraction between overlapped but nonbonded densities?
The exact exchange-correlation energy arises from a Coulomb interaction between the electron density at each r and the density at each r 0 of the exchange-correlation hole around an electron at r [75,76].In many-electron systems, the exact exchange-correlation hole is more localized around its electron than is the exact exchange hole and thus more amenable to semilocal approximation.The cancellation of error between semilocal exchange and semilocal correlation is what makes semilocal approximations work for molecules [46] and for vdW-bound systems [45].It has been known for a long time that the exchange term of a semilocal functional can imitate intermediaterange correlation [46].
The LDA is the simplest nonempirical semilocal density functional for the exchange-correlation energy (taking a local functional as a limit of a semilocal one).LDA overbinds atoms in molecules and also overbinds small molecular complexes at equilibrium, exaggerating the intermediaterange attraction between closed-shell molecules.Note that the vdW binding of all semilocal density functionals essentially comes from the exchange part, and of course, the LDA attraction vanishes exponentially as the molecules move far apart, and their density overlap tends to zero.
When we pass to the GGA, this intermediate-range attraction is largely or completely lost because GGA makes the exchange energy density much more negative than in LDA for regions of large reduced density gradient s ¼ j∇nj=½2ð3π 2 Þ 1=3 n 4=3 , such as the nonoverlapped density tails of molecules [45].For example, the PBE [32] GGA exchange energy density e PBE x ðrÞ approaches 1.804e LDA x ðrÞ as s → ∞ for a spin-unpolarized density.Then, the molecules can lower their energy by moving further apart, increasing the region of nonbonded tail density.
While this behavior of GGA is suggested by the gradient expansion and greatly improves the atomization energies of molecules, there are reasons to believe that it is not correct.There is a rigorous tight lower bound, E x ½n ≥ 1.174E LDA x ½n, on the exchange energy of any spinunpolarized two-electron density [77], which is strongly suspected to hold for all spin-unpolarized densities [77].One argument [78] is as follows: (1) Start with a fixed spin-independent external potential that confines electrons to a given volume.(2) Add N electrons to this volume, one at a time.When N equals 2, the exact exchange energy E x ½n satisfies (4) Unless there is an unanticipated nonmonotonicity, the tight bound E x ½n ≥ 1.174E LDA x ½n will hold for all intermediate densities.This tight bound is strongly violated by GGA, as well as by many pre-SCAN meta-GGAs, which also capture little or no intermediate-range attraction.
This tight bound can only be satisfied for all spinunpolarized densities within a semilocal approximation by imposing a similar bound on the exchange energy density, e x ðrÞ ≥ 1.174e LDA x ðrÞ [77].It was imposed on SCAN for two-electron densities in this way, and it then turned out to be satisfied by SCAN for all spin-unpolarized densities.Thus, SCAN, while improving atomization energies of molecules even more than GGA, keeps part of the LDA intermediate-range attraction between molecules.

B. vdW corrections to meta-GGA
The benchmarking calculations have established SCANþrVV10 as a successful versatile vdW density functional, especially in the near-equilibrium cases.The success builds on three facts: the high quality of SCAN as a semilocal density functional even for the intermediaterange vdW interactions, the high quality of the nonlocal correlation functional in rVV10 for long-range vdW interactions, and the flexibility in the rVV10 method to adjust the short-range damping via the b parameter.Specifically, the high quality of SCAN requires a downscaled long-range correction such that the error introduced by the correction is also downscaled.Based on this argument, we expect similar success from a combination of SCAN with other long-range vdW corrections that have such flexibility: the DFTþD [25,26] and Tkatchenko-Scheffler (TS) methods [27][28][29].This idea has been confirmed by the very recent development of SCAN-D3 [72].In these methods, pairwise interatomic terms are added to the DFT potentials, with explicit usage of C 6 parameters and a damping function.These methods are, in principle, ion or atom based, instead of electron density based, and hence may not be applied as generally as the rVV10, VV10, and Rutgers-Chalmers vdW-DF methods [2,3,20].Furthermore, the DFTþD methods only introduce corrections to the energy and forces but not to the Kohn-Sham effective potential, which is important to describe metal surfaces, as shown [79] by the self-consistent implementation of a TS [27] method.However, incorporating SCAN within these methods can be of special interest.For example, incorporating SCAN within the DFTþD methods, one can readily include high-order C 8 , C 10 terms.
Our version of rVV10 is fitted to systems with fundamental energy gaps (Ar 2 and pairs of well-separated molecules), but it performs just as well for the gapless systems studied here (the coinage metals interacting with benzene, and the zero-gap layered materials).We would not, however, expect rVV10 to be accurate for the longrange vdW asymptotics in systems of small or zero gap with more challenging geometries, such as pairs of large fullerene molecules [80,81].While taking the electrondensity distribution into account [27], the more-recent TS methods [28,29] further include screening effects and many-body dispersion, which can be important for vdW asymptotics [81].Even when SCANþrVV10 is not accurate for the long-range vdW coefficient C 6 , it could still be accurate for near-equilibrium binding.This appears to be so for the zero-gap layered material graphite, although the interaction between two isolated distant graphene planes is derived to be D −3 [82], but ∼D −4 in pairwise-interaction models including SCANþrVV10.
The Rutgers-Chalmers vdW-DF family [20], which inspired the VV10 and rVV10 methods, currently sticks to the LDA correlation in E 0 xc .Much effort has been made to design a better accompanying semilocal exchange functional [2][3][4][5][6][7][8], which should ideally include no vdW effects.We cannot expect improvement of the Rutgers-Chalmers vdW-DF1 and vdW-DF2 methods by directly incorporating SCAN, either only in the exchange part or the whole semilocal exchange correlation energy.Other meta-GGAs [83,84] may provide better options for the semilocal part.
SCANþrVV10 is the first method to deliberately pair a nonlocal vdW correlation functional with a semilocal functional that includes intermediate-range vdW.Its success suggests new routes to include vdW within the framework of DFT: The semilocal part of the functional can be a meta-GGA, which properly describes all bonding types, including intermediate-range vdW, leaving the nonlocal part to account only for the long-range vdW contribution.Besides, SCANþrVV10 inherits the improvements of SCAN meta-GGA over other nonempirical or semiempirical GGAs, including better energetics and geometries (as shown here as well as in Ref. [35]), better band gap for semiconductors and insulators [85], etc.

C. Effective range of the vdW interaction
Note that rVV10 is a sophisticated pairwise-interaction model for the long-range vdW interaction, in which the interaction is not between pairs of atoms but between pairs of volume elements.No pairwise-interaction model will always give the correct asymptotics for the vdW interaction, as discussed in Sec.IV B. Yet SCANþrVV10 seems to be very accurate for the binding energies and geometries studied here.How can this be?
Recently, Ambrosetti et al. [86] calculated the power-law exponent of the vdW interaction as a function of the distance between layers in bilayers of materials such as graphene and MoS 2 .They used the many-body dispersion approach [28] and found that the asymptotic exponent does indeed differ from that predicted by a pairwise-interaction model, but the asymptote is only approached at very large separations on the order of 100-200 Å.
Concerns about the pairwise-interaction nature of SCANþrVV10 can be addressed to some extent by estimating an effective range for the vdW interaction, which can be found by introducing a truncated kernel Φðr; r 0 Þ with value zero for jr − r 0 j > r cut .As shown in Table III, r cut ¼ 8 Å agrees with r cut ¼ 53 Å for small molecules, including (stretched) Ar 2 , and the S22 set of molecular complexes.For solids and layered materials, r cut ¼ 16 Å agrees with r cut ¼ 53 Å for the cohesive energies and interlayer binding energies within 3 meV=atom, and the geometric properties remain stable with an even shorter r cut .Such a short r cut indicates that the incorrect vdW asymptotic behavior of pairwise-interaction TABLE III.The mean relative error (MRE) and mean absolute relative error (MARE) of the truncated SCANþrVV10 predictions for the Ar 2 binding curve, S22 interaction energy, atomization energy and lattice volume for 50 solids, and interlayer binding energy and lattice constants of 28 layered materials calculated with the default r cut ¼ 53 Å, compared with those from a much smaller r cut .The MRE and MARE are in percentage.models is not sampled significantly by the energetic and geometric properties studied here.

V. CONCLUSION
In conclusion, we have devised a promising vdW density functional SCANþrVV10 based on the SCAN meta-GGA.It works for general geometries as exemplified by its excellent performances for the molecular complexes, solids, benzene adsorption on metal surfaces, and layered materials.SCANþrVV10 achieves an accuracy comparable to that of higher-level methods like RPA and CCSD(T) for various situations benchmarked here, with a computational cost reduced by several orders of magnitude.SCANþrVV10 outperforms other currently available methods with comparable computational efficiency for the energetics and structures of layered materials, making it the right choice for the computational study of 2D layered materials.The success of SCANþrVV10 suggests new routes to include vdW within the framework of density functional theory, as summarized in Sec.IV B.

APPENDIX A: COMPUTATIONAL DETAILS
All calculations in this study were performed with the projector augmented wave (PAW) method [87] as implemented in the VASP code (version 5.3.5)[88][89][90].The PAW Perdew-Burke-Ernzerhof (PBE) [32] pseudopotentials (version .52)recommended in the VASP manual [91] were employed.The rVV10 nonlocal correlation functional has been implemented within the VASP code based on the vdW-DF implementation by Klimeš et al. [5] and has been benchmarked by comparing with the original rVV10 implementation [30] within the Quantum-Espresso code.
To calculate the binding curve of Ar 2 , the dimer and the atom were put in a cubic supercell with a length of 25 Å, an energy cutoff for the plane-wave basis as high as 1200 eV was chosen to ensure high quality for the fitting, and the single Γ point was used for the Brillouin-zone sampling.
To calculate the interaction energies of the S22 data set [37,38], the molecular dimer or complex was put into an orthorhombic supercell, where the periodic images were separated with a 20-Å-thick vacuum slab along all three directions to avoid interaction between them.An energy cutoff of 900 eV was chosen.The same simulation cell and the same computational parameters were used for the components of the molecular dimer or complex to minimize errors.For all molecular calculations, the single Γ point was used for the Brillouin-zone sampling.
For the 50 solid systems, the initial atomic structures for the solids and the spin configuration for the isolated atoms were chosen to be the same as those in Refs.[48,49], and 800 eV was chosen for the basis cutoff energy.The isolated atom was simulated with a 14 × 15 × 16 Å 3 orthogonal cell.For a solid, a Γ-centered Monkhorst-Pack [92] k-mesh corresponding to at least 4000 per reciprocal atom was used for the Brillouin-zone integration.Both the lattice vectors and internal atomic coordinates have been relaxed until the residual atomic force was less than 0.01 eV=Å.
The initial atomic models for the benzene molecule on Cu, Ag, and Au (111) surfaces were constructed based on those in Ref. [51].The metallic surfaces were modeled with five 3 × 3 atomic layers, with the bottom three layers fixed to their bulk coordinates.A 20-Å-thick vacuum slab was inserted between the metal slab and its periodic images.During the relaxation, the lattice vectors were fixed, with only the internal atomic coordinates of the top two metal atomic layers and the benzene molecule optimized.The energy cutoff was 600 eV, and the Γ-centered 4 × 4 × 1 Monkhorst-Pack [92] k-meshes were used.For SCAN, since the binding curve is very shallow, we further calculated the adsorption energies at distances Δz from 2.8 to 3.6 Å with a step of 0.1 Å.This is the reason why we only keep one decimal digit for the SCAN Δz in Table I.The same was done for SCANþrVV10, but we found SCANþrVV10 was able to find the minima with a standard atomic relaxation.
To calculate the binding curves for graphene in registry with the Ni(111), Co(0001), and Cu(111) surfaces, we put the graphene sheet on top of the surface (as shown in Fig. 1 of Ref. [69]).The metal surface was modeled with a four-layer slab generated with the experimental lattice constants and a 20-Å-thick vacuum layer.The energy cutoff of 600 eV was used, and the Γ-centered 16 × 16 × 1 Monkhorst-Pack [92] k-meshes were used.As in the RPA calculations [68,69] the PBE orbitals were used to evaluate the total energies without self-consistent iterations for the other functionals.
The binding energy and the equilibrium lattice constants of the face-centered cubic (fcc) Ne, Ar, and Kr solids were extracted from the binding-energy curves near the equilibrium lattice constants.These curves have been fitted to an eighth-order polynomial.An energy cutoff of 750 eV was used, and the Γ-centered 12 × 12 × 12 Monkhorst-Pack [92] k-meshes were used.The reference free atom was modeled with a single rare-gas atom in a 6 × 6 × 6 fcc supercell.
For the layered-structure materials, the initial atomic structures were taken from Refs.[17,19].The energy cutoff for the plane-wave basis was 800 eV, and the Γ-centered 12 × 12 × 6 and 12 × 12 × 1 Monkhorst-Pack [92] k-meshes were used for the bulk form and monolayer, respectively.A 30-Å-thick vacuum slab was utilized to model the monolayers.The internal coordinates were always fully relaxed.To determine the lattice vectors of the bulk structures, we further relaxed all three lattice vectors.To determine the binding energy, we fixed the intralayer lattice constants to their experimental values and relaxed the interlayer lattice constants only for the bulk structures, as in the RPA calculations [17,19].Table VI shows that the theoretical and experimental intralayer lattice constants agree very closely.The relaxation criteria were the same as those for the 50 solids.

APPENDIX B: RAW DATA FOR THE BOX PLOTS
The raw data for generating the box plots in Figs.2,3, and 5 are given in Tables IV-VI.TABLE IV.Interaction energies, in kcal=mol, for the molecular dimers in the S22 data set from PBE, rVV10, vdW-DF2 (numerical results from Ref. [10]), SCAN, and SCANþrVV10 with respect to the CCSD(T) results [37,38].[48,49].

FIG. 5 .
FIG.5.Box plots for the absolute relative errors and relative errors of the interlayer binding energies, and interlayer and intralayer lattice constants (c and a) from rev-vdW-DF2, SCAN, and SCAN þ rVV10, for 28 layered materials.The reference values are from RPA for the binding energy and from experiment for the lattice constants[17,19].
This work was supported as part of the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Grant No. DE-SC0012575.This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy.
[53]chemical trend of the ionic radii, we are confident in our SCANþrVV10 prediction for Δz in the case of Cu.The results from other methods, including PBE, PBE þ vdW, PBE þ vdW surf , MP2, vdW-DF1, DFTþD, and M06-L, etc., have been compiled in Ref.[53].Unfortunately, none of them simultaneously predicts the correct E ad and Δ z , with the best results coming from PBE þ vdW surf Abs. Rel.Err.for E a Rel.Err.for E a Abs.Rel.Err.for V 0 FIG. 3. Box plots for the absolute relative errors and relative errors of the atomization energies E a , and the equilibrium lattice volumes V 0 , from RPA, PBE, SCAN, and SCAN þ rVV10 for 50 solids, with respect to the experimental values.The RPA, PBE, and experimental values are from Refs.[48,49].Considering

TABLE V .
Atomization energies E a and equilibrium lattice volumes V 0 from RPA, PBE, SCAN, and SCANþrVV10 for 50 solids, compared to the experimental values.The RPA, PBE, and experimental values are from Refs.

TABLE VI .
[17,19]ayer binding energy E b in meV=Å 2 , interlayer lattice constant c in Å, and intralayer lattice constant a in Å for 28 layered materials from rev-vdW-DF2, SCAN, and SCANþrVV10.The reference values for E b are from RPA calculations and from experiments for c and a[17,19].