Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation

• James W. Furness*
• ,
• Aaron D. Kaplan
• ,
• Jinliang Ning
• ,
• John P. Perdew
• , and
• Jianwei Sun*

Cite this: J. Phys. Chem. Lett. 2020, 11, 19, 8208–8215

Publication Date:September 2, 2020

https://doi.org/10.1021/acs.jpclett.0c02405

Copyright © 2020 American Chemical Society

RIGHTS & PERMISSIONS

📷

Article Views

9313

Altmetric

21

Citations

210

LEARN ABOUT THESE METRICS

Share

Add to

ExportRIS

PDF (1 MB)

Supporting Info (1)»Supporting Information

SUBJECTS:

• Approximation,
• Energy,
• Show More

📷

The Journal of Physical Chemistry Letters

Abstract

📷

The recently proposed rSCAN functional [ J. Chem. Phys. 2019 150, 161101] is a regularized form of the SCAN functional [ Phys. Rev. Lett. 2015 115, 036402] that improves SCAN’s numerical performance at the expense of breaking constraints known from the exact exchange–correlation functional. We construct a new meta-generalized gradient approximation by restoring exact constraint adherence to rSCAN. The resulting functional maintains rSCAN’s numerical performance while restoring the transferable accuracy of SCAN.

There is a fundamental trade-off at the heart of all large-scale chemical and material computational studies between prediction accuracy and computational efficiency. The level of theory used must simultaneously make accurate and efficient material property predictions. For many projects, Kohn–Sham density functional theory (KS-DFT) presents an appealing compromise, delivering useful accuracy and favorable algorithmic complexity.

The Materials Project database presents a case study of finding such a balance, (1) stating an ambitious mission of “removing the guesswork from materials design by computing properties of all known materials”. (2) At the time of writing, the database lists 125 000 inorganic structures calculated from KS-DFT using the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) exchange–correlation (XC) functional. (3) While GGA functionals can be impressively accurate for many properties, they cannot be systemically accurate for all properties, (4−6) and the last 10 years have shown that meta-GGA functionals can improve predictions for similar computational cost.

Meta-GGAs are commonly designed around constraints known for the exact XC functional while minimizing the number of free parameters that must be fit. Functionals derived in this way are termed “non-empirical,” and we refer the reader to the Supporting Information of ref (7) for precise definitions of all the exact constraints known for meta-GGAs. Alternatively, the functional can be built from a more flexible form that allows some exact constraints to be broken, so that free parameters can be tuned for accuracy to reference data sets. Functionals taking the latter route, termed “empirical” functionals, tend to be less reliable for systems outside their fitting sets, making a non-empirical functional desirable for large-scale applications.

The strongly constrained and appropriately normed (SCAN) functional (7) recovers all 17 exact constraints presently known for meta-GGA functionals and has shown good transferable accuracy, even for systems challenging for DFT methods. Examples include predicting accurate geometries and energetics for diverse ice and silicon phases (8) and for polymorphs of MnO2. (9) SCAN accurately reproduces the complex doping-driven metal–insulator transition, magnetic structure, and charge-spin stripe phases of cuprate (10−12) high-temperature superconductors and iridates. (13) It is one of the few functionals that predicts ice as less dense than liquid water under standard conditions, (14) and its description of intermediate range van der Waals interactions has been used to study the dynamics of liquid water. (14,15) Combination of SCAN with beyond DFT techniques such as van der Waals functionals and the Hubbard U self-interaction correction have proven effective for modeling the ionic and electronic structures of transition metal oxides. (16−18)

Despite these successes, SCAN’s utility for large-scale projects is limited by its sensitivity to the density of the numerical integration grid used during calculation. This poor numerical performance in many codes mandates the use of dense integration grids which reduces SCAN’s computational efficiency, (19,20) and divergence in the associated XC potential has hindered the generation of SCAN pseudopotentials. (21,22) Neither limitation is inherent to the meta-GGA level or SCAN-like functionals, as we will show.

Some modifications to SCAN have been proposed to improve its accuracy for specific systems. The revSCAN functional is a simple modification to the slowly varying limit of SCAN’s correlation energy to eliminate the fourth-order term in SCAN’s correlation energy density-gradient expansion. (23) The TASK functional is a complete revision of SCAN designed to accurately predict band gaps while retaining the exact constraints placed on the exchange energy, (24) though TASK uses a local spin-density approximation (LSDA) to model correlation. It is not expected that these modifications address the numerical inefficiencies of the parent functional.

In recent work, Bartók and Yates propose a regularized SCAN termed “rSCAN” that aims to control SCAN’s numerical challenges while changing as little as possible from the parent functional. (22) The resulting functional shows greatly improved numerical stability and enables pseudopotential generation. While initial testing suggested that rSCAN maintained the accuracy and transferability of SCAN, expanded testing by Mejı́a-Rodrı́guez and Trickey (25,26) shows that some transferability is lost, with accuracy for atomization energies (27) particularly degraded.

The need for a computationally efficient revision of SCAN is made plain in Figure 1. This figure shows three meta-GGAs: SCAN and rSCAN, which have already been introduced, and a novel meta-GGA, r2SCAN, that is the topic of this Letter. It illustrates a grid problem that arises for SCAN in codes with localized basis functions. The horizontal axis shows increasing integration grid density, and the vertical axis shows the mean absolute error (MAE) of the G3 test set (28) of 226 atomization energies. It would be difficult to assert that any of the grid settings present a converged SCAN energy, with SCAN errors varying unpredictably by a factor of 2. While rSCAN stabilizes SCAN numerically, its error offers little improvement over GGAs (e.g., PBE has a MAE of 22.2 kcal/mol (7) for the G3 set). The need for a meta-GGA that retains the accuracy of SCAN, with the grid efficiency of rSCAN, is evident. No such grid problem is found for SCAN in the plane-wave code VASP, as shown in the Supporting Information. However, in VASP, r2SCAN seems to converge with fewer iterations than SCAN does.

Figure 1 📷Figure 1. Mean absolute error (MAE) of atomization energies (kcal/mol) for the G3 set of 226 molecules (28) as a function of increasing numerical integration grid density expressed relative to the smallest grid. The grids were chosen from Turbomole (29,30) grid levels 1–7.

The SCAN functional is constructed using a dimensionless kinetic energy variable

𝛼(𝐫)=𝜏(𝐫)−𝜏W(𝐫)𝜏unif(𝐫)

(1)

where τ = ∑i|∇ϕi|2Θ(μ – εi)/2 is the positive kinetic energy density; ϕi(r) are the Kohn–Sham orbitals; Θ(μ – εi) is the orbital occupation; τW = |∇n|2/(8n) is the von Weizsäcker kinetic energy density; τunif = 3(3π2)2/3n5/3/10 is the kinetic energy density of a uniform electron gas; μ is the chemical potential; and εi are the orbital energies. SCAN uses α to tune functional performance for the local chemical environment. (31) While α allows SCAN to satisfy exact constraints that would be contradictory at the GGA level, (32) α can introduce numerical sensitivity and divergences in the XC potential. (33,34)

The design of rSCAN prioritizes numerical efficiency over satisfaction of exact constraints and instead uses a regularized α′

𝛼∼(𝐫)=𝜏(𝐫)−𝜏W(𝐫)𝜏unif(𝐫)+𝜏r

(2)

𝛼′(𝐫)=𝛼∼(𝐫)3𝛼∼(𝐫)2+𝛼r

(3)

where τr = 10–4 and αr = 10–3 are regularization constants. While the choice of a constant τr eliminates numerical instability as α → 0, α′ does not retain the correct uniform and nonuniform scaling properties of α, nor the correct uniform density limit.

For a uniform electron gas, α → 1, which SCAN uses to recover the LSDA exactly. In rSCAN, α̃ → 1/(1 + τr/τunif) which varies with the density, losing the correct uniform electron gas description. It has been shown that recovery of the uniform gas limit is critical for an accurate description of solids, atoms, and molecules. (35,36)

For a slowly varying electron gas, the exchange and correlation energies have well-known expansions in powers of the gradient of the density. Let s = |∇n|/(2kFn), a dimensionless density-gradient on the scale of the Fermi wavevector kF = (3π2n)1/3, and q = ∇2n/(4kF2n) a dimensionless density-Laplacian. The gradient expansion for the exchange energy per particle εx(r) is (37)

𝜀x=𝜀LDAx[1+𝜇AK𝑝+1462025(𝑞2−52𝑝𝑞)]+𝒪[(∇𝑛)6]

(4)

where εxLDA = −3kF/(4π), p = s2, and μAK = 10/81. For the correlation energy, following ref (3), we define an additional dimensionless density-gradient t = |∇n|/[2ksϕ(ζ)n] on the scale of the Thomas-Fermi screening wavevector 𝑘s=4𝑘F/𝜋⎯⎯⎯⎯⎯⎯⎯⎯⎯√. Here ϕ(ζ) = [(1 + ζ)2/3 + (1 – ζ)2/3]/2 is a spin-scaling function of the spin-polarization ζ = (n↑– n↓)/n. Then the density-gradient expansion of the correlation energy per particle εc(r) is (3,7,38)

𝜀c=𝜀LSDAc+𝜙(𝜁)3𝛽(𝑟s)𝑡2

(5)

where β(rs) is a weakly varying function of the Wigner–Seitz radius rs = (4πn/3)−1/3, with a maximum β(0) ≈ 0.066725. The kinetic energy density τ has an analogous but unwieldy density-gradient expansion. (39) It is generally understood that recovering the exact density-gradient expansion is relevant for solids. (40) These terms also affect the asymptotic behavior of Exc for atoms, (41) as the asymptotic limit for atoms of large-Z is a semiclassical limit that is described exactly by the LSDA at lowest-order, with the density-gradient terms modulating the higher-order terms (known accurately). (41) Thus, the uniform and slowly varying density limits are relevant to both solid-state and atomic systems.

SCAN eliminates erroneous contributions from α to the second- and fourth-order slowly varying density-gradient expansion (GE2 and GE4, respectively) of Exc by using a nonanalytic switching function whose value and derivatives of all orders are zero at α = 1. While theoretically convenient, constraining the interpolation function to have zero derivatives at α = 1 results in a twisted function that harms numerical performance. The SCAN interpolation function was replaced with a smooth polynomial in rSCAN (see Figure 2) to remove this source of numeric instability, at the expense of introducing second- and fourth-order contributions from α to the density-gradient expansion of Exc.

Figure 2 📷Figure 2. SCAN (blue, solid) and rSCAN (red, dashed) interpolation functions plotted for a generic stand-in iso-orbital indicator “Y” (α for SCAN, α′ for rSCAN, α̅ for r2SCAN). The functions mix Y = 0 (single orbital) and Y = 1 (uniform density limit for α and α̅ ) energy densities. The derivatives of the SCAN interpolation functions vanish to all orders in Y at Y → 1, allowing SCAN to recover the appropriate density-gradient expansions exactly in the slowly varying limit. The rSCAN interpolation functions are used with Y = α̅ in r2SCAN, and their smooth, nonvanishing first derivatives at Y = 1 necessitate changes from SCAN to r2SCAN in the Y = 1 energy densities.

It is clear then that rSCAN makes wide-ranging sacrifices in exact constraint adherence in order to make a numerically efficient meta-GGA. Here, we will show definitively that such sacrifices are needless and derive revisions to the rSCAN functional to restore exact constraint adherence without harming numerical efficiency. We apply these revisions to build a regularized–restored SCAN functional, r2SCAN, which recovers the most important exact constraints of SCAN. Table 1 summarizes the constraint satisfaction of the functionals concerned, and we stress that because only appropriate norm systems (7) were used to set the free parameters, all three functionals (SCAN, rSCAN, and r2SCAN) may be considered non-empirical. For brevity, we show only parts of the functional that are modified in this work and direct the reader to Section S2 of the Supporting Information for a full definition of the relevant equations.

Table 1. Summary of Exact Constraint Adherence for a Subset of the 17 Known Exact Constraints Applicable to Meta-GGA Functionalsa

SCANrSCANr2SCANuniform density√−√coordinate scaling√−√GE2√−√GE4X√−−

a

Here, GE2 denotes the second-order slowly-varying density-gradient expansion, and GE4X denotes the fourth-order GE for exchange.

There are many situations where the exact exchange–correlation potential and energy density can be expected to be reasonably smooth (see, e.g., the plots of highly accurate exchange–correlation potentials and energy densities of simple hydrides in ref (42)). In general, the exact Kohn–Sham exchange–correlation potential need not be smooth, as demonstrated by the Perdew–Parr–Levy–Balduz theorem: (43) the exchange–correlation potential, as a function of the number of electrons N, exhibits discontinuities across integer values of N, with steps and peaks in the low-density region between two separated dissimilar systems. However, a semilocal functional cannot recover the precise behaviors of the exact exchange–correlation energy and potential and instead averages over them. Therefore, we consider smoothness of the energy density and potential to be a necessary construction principle of semilocal approximate density functionals. A construction principle is any physically or mathematically motivated principle that can supplement the design of a first-principles density functional approximation.

The correct uniform- and nonuniform-scaling properties of α, as well as the correct uniform density limit of Exc, are recovered in r2SCAN by regularizing α as

𝛼¯=𝜏−𝜏W𝜏unif+𝜂𝜏W

(6)

where η = 10–3 is a simple regularization parameter. Note that because τ ≥ τW, α̅ has the same range as α, 0 ≤ α̅ < ∞. This is distinct from the dimensionless kinetic energy variable suggested by ref (34)

𝛽=𝜏−𝜏W𝜏+𝜏unif

(7)

which ranges between 0 ≤ β < 1 and has less rapidly varying derivatives than α. As this work seeks revisions to SCAN, we will not consider β or related iso-orbital indicators here and adopt α̅ as the iso-orbital indicator used throughout r2SCAN.

SCAN uses the iso-orbital indicator variable α to drive interpolation functions fx/c(α) for exchange and correlation. The rSCAN functional replaces the original SCAN interpolation functions with a polynomial function of α′ when α′ < 2.5 that smooths out the plateau-like behavior of the original near α = 1. The r2SCAN functional adopts the rSCAN interpolation function but uses α̅ as the indicator variable. Both the SCAN and rSCAN interpolation functions are shown in Figure 2 as functions of a generic indicator.

The SCAN interpolation function was designed to have vanishing derivatives at α = 1, but the rSCAN replacements go linearly through zero at this point. As a result of these nonvanishing derivatives, the interpolation function makes spurious contributions to the slowly varying density-gradient expans