r²SCAN-D4: Dispersion corrected meta-generalized gradient approximation for general chemical applications Free

The quantum mechanical description of physical and chemical materials at the electronic resolution is an increasingly important task for in silico simulations. Here, density functional theory (DFT) has emerged in the past decades as one of the most versatile methodological frameworks.^1,2 This leading position in both materials and chemical applications is largely due to the excellent accuracy over the computational cost ratio as well as the broad applicability across system classes of today’s density functional approximations (DFAs).^3–5 

The Jacob’s ladder hierarchy⁶ is commonly used to classify DFAs. In this hierarchy, DFAs are systematically improved by ascending rungs of different approximations: the local density approximation (LDA), generalized gradient approximations (GGAs), meta-GGAs, hybrid functionals (including a fraction of nonlocal exact exchange), and double-hybrid functionals (including nonlocal correlation). In terms of efficiency, meta-GGAs are in a favorable spot, as they have the same cubic scaling with the system size as LDA. Yet, many of the meta-GGAs proposed so far cannot truly leverage the full potential of their rung. Some shortcomings of the existing functionals are increased sensitivity to the numeric integration grid as observed in the strongly constrained and appropriately normed (SCAN) functional⁷ or several Minnesota type functionals,^8–10 purely empirical parameters as present in the B97M functional,¹¹ and sensitivity to the kinetic energy density.^12,13 Recent developments of semi-local DFAs combine exact constraints with various degrees of parametrization to improve descriptions of short- to medium-range electron correlation.^7,14,15

The SCAN functional⁷ is constructed to rigorously satisfy all known exact constraints suitable for a meta-GGA. While the functional itself has shown excellent performance in previous studies, the severe numerical instabilities inherent to the functional impeded its adoption for many computational studies. With the recently proposed regularized SCAN (rSCAN)¹⁶ and the subsequent restoration of exact constraints in r²SCAN,¹⁷ the main drawback of the SCAN functional seems to be resolved. Shortcomings of the SCAN functional might still be present in its successors, rSCAN and r²SCAN. Notably, the description of water clusters (H₂O)_n (2 ≤ n ≤ 8) show the overbinding tendency of SCAN.¹⁸ A recent work by Sharkas et al.,¹⁸ however, demonstrates that this can be mended by the Perdew–Zunger self-interaction correction (PZ-SIC). SCAN and r²SCAN show similar behaviors for self-interaction error prone systems, which may make r²SCAN amenable to a PZ-SIC correction as well. However, r²SCAN is often more accurate than SCAN, as in the extensive benchmarking here, in the atomization energies of molecules¹⁷ and in the spin-crossover energies of molecules.¹⁹ 

Nevertheless, semilocal functionals cannot include long-range correlation effects such as London dispersion interactions.²⁰ To truly judge its applicability, we extensively tested r²SCAN combined with the state-of-the-art D4 dispersion correction,²¹ which shows unprecedented performance for a range of diverse chemical and physical properties. To investigate the development of the SCAN-type functionals, we include both SCAN-D4 and rSCAN-D4 in the comparison to r²SCAN-D4 and can attribute improvements in non-covalent interactions (NCIs) mainly to the regularization and improvements in thermochemistry and barrier heights to the restoration of the exact constraints.

We give a concise methodological overview (Sec. II) on r²SCAN and D4 before testing the full method against established DFAs over a wide range of benchmarks (Sec. III), with particular focus on molecular geometries, thermochemistry, kinetics, and non-covalent interactions in small and large complexes.

The rSCAN¹⁶ functional regularizes the severe numerical instability or inefficiency of the otherwise successful SCAN⁷ functional at the expense of breaking exact constraints SCAN was constructed to obey. This problem arises in many codes that employ localized basis sets and is less problematic in many codes that employ plane wave basis sets. While numerical challenges are indeed resolved, a rigorous adherence to exact constraints is core to the design of the SCAN functional and likely important for transferable accuracy across domains of applicability.¹⁹ This seems to be reflected in rSCAN’s relatively poor performance for molecular atomization energies compared to other tests.^22,23 The r²SCAN functional¹⁷ combines the good accuracy of SCAN with the numerical efficiency of rSCAN by directly restoring exact constraint satisfaction to the rSCAN regularizations.

The SCAN functional is constructed as an interpolation between single orbital and slowly varying energy densities designed to maximize exact constraint satisfaction.⁷ The interpolation is controlled by an iso-orbital indicator

where τ_W = |∇ρ|²/(8ρ) and $τU=3(3π2)2/3ρ5/3/10$ are the von Weizsäcker and uniform electron gas kinetic energy densities, respectively.²⁴ In subsequent studies, α has been shown to contribute to numerical instability.^25,26 To remove these effects, a regularized α′ was used in rSCAN that removes single orbital divergences at the expense of breaking exact coordinate scaling conditions^27–29 and the uniform density limit. These conditions are restored in r²SCAN by adopting a different regularization,

where η = 10⁻³ is a regularization parameter.

The second regularization made in the rSCAN functional is to substitute the twisted piece-wise exponential interpolation of the original SCAN with a smooth polynomial function. This removes problematic oscillations in the exchange-correlation potential but introduces spurious terms in the slowly varying density gradient expansion that deviate from the exact expansion^30,31 recovered by SCAN. A corrected gradient expansion term is used in r²SCAN that cancels these spurious terms, so the functional recovers the slowly varying density gradient expansion to the second order. A recent modification of SCAN for improved bandgap accuracy from Aschebrock and Kümmel named “TASK”³² is able to enforce the fourth-order gradient expansion for the exchange energy without apparent numerical problems,³³ resolving the dominant source of numerical inefficiency. The importance of the fourth-order exchange terms is not established, however, and we are thus satisfied using one less exact constraint compared to SCAN. TASK uses an LSDA for correlation, however, and consequently violates many important exact constraints for correlation, e.g., the second order gradient expansion, that are obeyed by SCAN and r²SCAN.

Numerical instabilities are revealed by SCAN’s sensitivity to the choice of the numerical integration grid, often requiring dense computationally costly grids.^16,25,34 This issue has been addressed with the rSCAN and r²SCAN functionals. Figure 1 shows that the regularization indeed leads to two orders of magnitude error reduction when comparing r²SCAN with SCAN. This holds for both total energy and nuclear gradients for all chosen numerical settings. In practice, this allows for more computationally favorable settings. To give a rough estimate of the computational cost of r²SCAN compared to SCAN, we consider system 10 of the S30L³⁵ with 158 atoms and 8250 atomic orbitals in a def2-QZVP basis set. A SCAN calculation using Turbomole’s grid 4 and a rad size of 50 (8.5 × 10⁶ grid points) would take ∼10 h, while an r²SCAN calculation with the Turbomole grid m4 and a rad size of 6 (1.6 × 10⁶ grid points) takes only three and a half hours for the same numerical accuracy, resulting in a computational saving of a factor of three to five.³⁶ We recommend using r²SCAN with a rad size of 6 and potentially increasing it to 10 for problematic geometry optimizations.³⁷ We also compared SCAN and r²SCAN with different energy cutoffs in a PAW expansion and found that r²SCAN is not as sensitive as SCAN, i.e., the total energy converges significantly faster.³⁸ 

As London dispersion interactions arise from nonlocal electron correlations, they cannot be captured by any meta-GGA. In the past years, a range of schemes have been developed to capture these interactions in the DFT framework.^20,39–43 Here, we combine r²SCAN with the semi-classical D4 dispersion correction.²¹ Its energy contribution is calculated by

where R_AB is the atomic distance, $CnAB$ is the n-th order dispersion coefficient, and $fnBJ(RAB)$ is the Becke–Johnson damping function.^44,45R_ABC and $C9ABC$ denote the geometrically averaged distance and dispersion coefficient, respectively, and θ_ABC is the angle dependent term of the triple-dipole contribution.^46,47 The s₈ parameter for the two-body dispersion and the a₁ and a₂ parameters entering the critical radius in the damping function are adjusted to match the local description of a specific DFA. Damping parameters are fitted using a Levenberg–Marquardt least-squares minimization to reference interaction energies, as described in Ref. 21. Optimized parameters are given in Table I. The b parameters for rSCAN-VV10 and r²SCAN-VV10 were determined to be 10.8 and 12.3, respectively, on the same set.^42,51

All ground state molecular DFT calculations were performed with a development version of Turbomole.^52,53 The resolution of identity (RI) approximation^54,55 was applied in all calculations for the electronic Coulomb energy contributions. For all functionals except SCAN, Turbomole’s modified grids of type m4 were used. For all SCAN calculations, grid 4 with an increased radial integration size of 50 was used instead. Self-consistent field convergence criteria of 10⁻⁷ Hartree were applied. Ahlrichs’ type quadruple-zeta basis sets, def2-QZVP,⁵⁶ were used throughout if not stated otherwise.

The periodic electronic structure calculations were conducted with vasp 6.1^57,58 with projector-augmented plane waves (PAWs) with an energy cutoff of 800 or 1000 eV (hard PAWs^59,60). Tight self-consistent field settings and large integration (and fine FFT) grids are used. The Brillouin zone sampling has been increased to converge the interaction energy to 0.1 kcal/mol. The non-periodic directions use a vacuum spacing of 12 Å.

To evaluate the description of covalent bond distances, we compare experimental and calculated ground–state equilibrium distances R_e (in pm) for 35 light main group bonds (LMGB35⁶¹), 11 heavy main group bonds (HMGB11⁶¹), and 50 bonds in 32 3d transition metal complexes (TMC32⁶²). Additionally, we investigate the bond distances and angles for a set of simple organic molecules against accurate semi-experimental Refs. 63 and 64. Extended statistics and optimized geometries are made freely available.⁶⁵ 

We include r²SCAN-D4, PBE0-D4,⁶⁶ TPSS-D4,⁶⁷ and PBE-D4⁶⁸ in the comparison shown in Fig. 2. For organic molecules, we find exceptional performance for all functionals, with errors smaller than 1 pm in the bond distances and half a degree in the bond angles. While all methods reproduce the reference values closely, we observe the best agreement from r²SCAN-D4 with a mean absolute deviation (MAD) of 0.4 pm and 0.3° for the bond distances and angles, respectively. For light main group elements, all methods give a mean absolute deviation of less than 1 pm as well, which was also observed in previous studies.^21,25 In comparison with the other methods tested here, r²SCAN-D4 also yields the lowest MAD of only 0.7 pm. Finally, for transition metal complexes, r²SCAN-D4 performs reasonably well with an MAD of 1.9 pm. Overall, the performance of r²SCAN is similar to, and sometimes even better than, the hybrid PBE0-D4, which in turn is one of the best performing hybrid functionals for molecular geometries.⁶¹ 

To investigate the performance of r²SCAN-D4 for general main group chemistry, we use the main group thermochemistry, kinetics, and non-covalent interactions (GMTKN55) database.⁶⁹ The GMTKN55 database is a compilation of 55 benchmark sets to assess the performance of DFAs and allows a comprehensive comparison of DFAs. It contains five categories: namely, basic properties, barrier heights, isomerizations and reactions, and intermolecular and intramolecular non-covalent interactions (NCIs), see Fig. 3 for a performance overview of selected methods. Usual weighted total MADs (WTMAD2s) range from 2 kcal/mol to 3 kcal/mol for double hybrid functionals, over 3 kcal/mol–4 kcal/mol for hybrid functionals to 8 kcal/mol–9 kcal/mol for (meta-)GGAs, while the lowest rung functionals such as PWLDA yield WTMADs of 17 kcal/mol on the GMTKN55. With the exception of the semi-empirical B97M-V¹¹ (and its B97M-D4 variant⁷⁰), r²SCAN-D4 is the best non-hybrid functional on the GMTKN55 so far with a WTMAD2 of 7.5 kcal/mol, compared to other meta-GGAs such as SCAN-D4 (8.61 kcal/mol) or TPSS-D4 (9.36 kcal/mol). For the isomerization and reaction category as well, as for the intramolecular NCIs, r²SCAN-D4 can even compete with the performance of the hybrid PBE0-D4 (WTMAD2 6.66 kcal/mol).

We additionally evaluated rSCAN-D4 and SCAN-D4⁷¹ on the GMTKN55 set to monitor the development in the SCAN-family of functionals. The main difference between SCAN-D4 and rSCAN-D4 is the general improvement in the description of non-covalent interactions, while both functionals perform similarly well in all other categories. Here, rSCAN-D4 improves for both NCI categories with a weighted MAD of 6.8 kcal/mol over SCAN-D4, which yields a weighted MAD of 7.6 kcal/mol. This improvement in rSCAN-D4 is mainly responsible for the smaller WTMAD2 of 8.3 kcal/mol compared to the WTMAD2 of 8.6 kcal/mol for SCAN-D4. For r²SCAN-D4, the improved description of NCI in rSCAN-D4 is preserved (weighted MAD of 6.6 kcal/mol), but r²SCAN-D4 bests its predecessor in all three remaining categories, resulting in its exceptional WTMAD2 of 7.5 kcal/mol. The mindless benchmark (MB16-43 subset of GMTKN55) is specifically useful for testing a method’s robustness to deal with the unusual chemistry in artificial molecules. Here, we see that enforcing exact constraints in non-empirical DFAs yields generally lower errors for artificial molecules than their empirical counterparts (see Table II).

To stress the importance of including a dispersion correction, we test the plain dispersion-uncorrected r²SCAN, which yields a significantly worse WTMAD2 of 8.8 kcal/mol, a difference similar in magnitude to the improvement from SCAN-D4 to r²SCAN-D4. In summary, r²SCAN-D4 shows a systematic improvement over its predecessor SCAN-D4 in all categories of GMTKN55 and can preserve improvements present in rSCAN-D4. This makes r²SCAN-D4 one of the best non-empirical meta-GGAs that have been broadly benchmarked so far.

Metal organic chemistry is one of the major application areas of non-hybrid DFAs. Here, we use the MOR41 benchmark set that contains 41 closed-shell metal–organic reactions representing common chemical reactions relevant in transition-metal chemistry and catalysis.⁷² We compare the statistical deviations from high-level references of r²SCAN-D4 to PBE0-D4, TPSS-D4, and PBE-D4 in Table III. The r²SCAN-D4 functional is one of the best meta-GGAs tested so far on the MOR41 benchmark set, with an MAD of 3.3 kcal/mol. While the SCAN-D4 method provides a slightly lower MAD, the analysis of other statistical quantities, such as the standard deviation (SD) and the maximum absolute error (AMAX), suggests less systematic results compared to r²SCAN-D4. This is confirmed by the Gini coefficient,⁷³ which is 0.44 for r²SCAN-D4 and 0.50 for SCAN-D4. Compare this to B97M-D4, one of the best meta-GGAs tested on the GMKTN55 set, which yields a larger MAD of 3.8 kcal/mol.⁷⁰ 

With the improved description of non-covalent interactions (NCIs), while retaining the computational efficiency of a meta-GGA, r²SCAN-D4 is a promising choice for interaction and association energies of large complexes. The results for the S30L,³⁵ L7,⁷⁴ and X40 × 10⁷⁵ benchmark set are shown in Fig. 4.

We choose the recently revised L7 benchmark⁷⁴ set to assess the performance of r²SCAN-D4 against local natural orbital based coupled cluster including singles, doubles, and perturbative triple excitations in the complete basis set limit [LNO-CCSD(T)/CBS].⁷⁶ Close agreement with an MAD of 0.9 kcal/mol is reached for r²SCAN-D4. This is a significant improvement over other meta-GGAs such as SCAN-D4 and TPSS-D4 with MADs of 1.3 kcal/mol and 1.4 kcal/mol, respectively.

We also investigated the description of association energies for large supramolecular complexes using the S30L benchmark set.³⁵ SCAN-D4 proved to be one of most accurate meta-GGAs in the previous benchmarks,²¹ giving a remarkable MAD of 2.0 kcal/mol, close to the uncertainty of the provided reference interactions; r²SCAN-D4 further improves upon this.

In particular, the association energies of the halogen-bonded complexes (15 and 16) are improved with r²SCAN-D4. The same trend can be observed in the HAL59 benchmark set of the GMTKN55, which shows an MAD of 1.0 kcal/mol with SCAN-D4 and improves with r²SCAN-D4 to an MAD of 0.8 kcal/mol. To confirm this trend, we additionally evaluated the X40 × 10 benchmark⁷⁵ containing 40 halogen bond dissociation curves with SCAN-D4 and r²SCAN-D4. Again, r²SCAN-D4 gives the lowest MAD of 0.36 kcal/mol, showcasing on overall improved description of this kind of NCIs.

To evaluate if the good performance for non-covalent interactions transfers from the gas phase to solids, molecular crystals and their polymorphic forms provide useful test cases.^77–79 Here, we investigate the lattice energy benchmark DMC8⁸⁰ shown in Table IV. The DMC8 benchmark contains a subset of the X23^78,81–83 and ICE10⁸⁴ benchmark sets with accurate structures and corresponding highly accurate fixed node diffusion Monte Carlo (FN-DMC) results. Due to SCAN’s tendency to overbind hydrogen bonded systems, such as ice polymorphs or hydrogen bonded molecular crystals, dispersion corrected SCAN was problematic for these systems. With the improved description of non-covalent interactions in r²SCAN-D4, this issue is mitigated and we find an overall improved MAD of 0.7 kcal/mol. This MAD is only half of the SCAN-D4 error of 1.5 kcal/mol for these systems and close to the very good performance of the hybrid PBE0-D4 of 0.5 kcal/mol.⁸⁵ Only the ice polymorphs are systematically overbound by r²SCAN-D4, which is, however, a problem of many functionals^86,87 and may be a self-interaction error.¹⁸ Both r²SCAN-D4 and SCAN-D4 yield similar results for the self-interaction subset (SIE4x4) of the GMTKN55, and therefore, a recent work investigating self-interaction corrections for SCAN might also be transferable to r²SCAN as well.⁸⁸ In contrast, the relative stability of the ice polymorph is reproduced correctly. The energy differences of ice II and ice VIII with respect to ice Ih are 0.03 kcal/mol and 0.70 kcal/mol, respectively, in good agrement with the reference values of 0.05 kcal/mol and 0.41 kcal/mol, respectively.

The benzene crystal has been frequently used for electronic benchmark purposes.^89–91 Here, we evaluated the equation of state (EOS) to compare with experimental measurements and the Murnaghan EOS fit to the FN-DMC from Ref. 80. The resulting EOS is shown in Fig. 5 and agrees excellently with the high-level method as well as the experimental estimate. A slight underestimation of the unit cell volume by 2.6% and overestimation of the bulk modulus by 6.4% can be seen. To highlight once again the importance of London dispersion on properties beyond the mere energy, we report plain r²SCAN results as well. The r²SCAN EOS has a significant offset equilibrium volume that is overestimated by 5.4% and a bulk modulus underestimated by 34.0%.

We have presented an accurate and robust combination of the non-empirical r²SCAN DFA with the state-of-the-art D4 dispersion correction. The resulting r²SCAN-D4 electronic structure method shows exceptional performance across several diverse categories of chemical problems assessed by thousands of high-level data points in a number of comprehensive benchmark sets. Included in the assessment were molecular thermochemistry for both main group and transition metal compounds, barrier heights, structure optimizations, lattice energies of molecular crystals, and both inter- and intramolecular non-covalent interactions of small to large systems, creating an extensive coverage of chemically relevant problems.

For the large GMTKN55 benchmark collection of about 1500 data points, r²SCAN-D4 is one of the most accurate meta-GGAs tested so far. Unlike the best meta-GGA on this set, the dispersion corrected B97M functional, r²SCAN-D4 can transfer this accuracy to chemically distinct systems such as metalorganic reactions. We find significant improvements in NCIs, which were one of the weak spots of SCAN based methods. More detailed analysis showed that improvements can mainly be found in the description of hydrogen and halogen bonded systems. The same trend is found for molecular crystals, where SCAN-D4’s tendency to overbind is mostly resolved in r²SCAN-D4, giving close to hybrid DFT results for lattice energies.

We found r²SCAN-D4 to be an accurate and (more importantly) consistent DFA for a large variety of problems and chemical systems. The already good performance of the original SCAN functional is kept and systematically improved in r²SCAN, while the numeric stability is almost on par with the established GGA functionals. We were able to gain some insight into the improvement from SCAN over rSCAN to r²SCAN, where we can attribute the improved description of non-covalent interactions to the regularization in the step from SCAN to rSCAN and the improved thermochemistry and barrier heights to the constraint restoration in the step from rSCAN to r²SCAN. Like SCAN, r²SCAN is not fitted to molecules, so its accuracy in extensive molecular tests demonstrates the predictive power of its exact constraints and appropriate norms.

With r²SCAN-D4, a meta-GGA method is finally available that truly leverages the advantages of its rung in Jacob’s ladder while retaining favorable numerical properties and fulfilling important exact constraints. We anticipate r²SCAN-D4 to be a valuable electronic structure method with broad applications in computational chemistry and material science.

We thank Stefan Grimme for valuable discussions. S.E. was supported by the DFG in the framework of the priority program 1807 “Control of London Dispersion Interactions in Molecular Chemistry.” J.N., J.W.F., and J.S. acknowledge the support of the U.S. DOE, Office of Science, Basic Energy Sciences (Grant No. DE-SC0019350, core research). A.D.K. acknowledges the support of the U.S. DOE, Office of Science, Basic Energy Sciences, through Grant No. DE-SC0012575 to the Energy Frontier Research Center: Center for Complex Materials from First Principles. J.P.P. was supported by the U.S. National Science Foundation under Grant No. DMR-1939528 (CMMT with a contribution from CTMC).

The data that support the findings of this study are openly available at https://github.com/awvwgk/r2scan-d4-paper and from the corresponding author upon reasonable request.