Regularized SCAN functional

First principles modeling of electronic structure has become a standard tool in studying the structure, stability, and dynamics of matter on the atomistic scale, with Density Functional Theory (DFT) being particularly popular, due to the balance of computational accuracy and cost.¹ The major source of inaccuracy in Kohn-Sham DFT calculations² is the necessity of using the exchange-correlation functional, which for general systems, only exists in approximate forms. Semilocal functionals based on the Generalized Gradient Approximation (GGA), for example, the Perdew-Burke-Ernzerhof (PBE) functional,³ model the electronic structure at a reasonable accuracy for a wide range of problems. However, there is a need for functionals with yet greater accuracy. Compared to GGAs, the meta-generalized Gradient Approximations (mGGA) provide more flexibility in the approximate functional form by introducing another local property on which the exchange-correlation functional depends, the orbital kinetic energy density, in addition to the electron density and its gradients. The recently proposed Strongly Constrained and Appropriately Normed (SCAN) functional is the first mGGA constructed such that it satisfies all known constraints that a semi-local functional can satisfy, and the remaining free parameters are fitted to reproduce exact or accurate reference values, or norms, of exchange and correlation energies. The resulting functional has proved broadly transferable⁴ and improved the DFT description of a wide range of systems, such as liquid water and ice,⁵ semiconductor materials,⁶ or metal oxides.⁷ 

Despite the tremendous success of SCAN, its implementation in DFT packages intended for condensed matter simulations is, at the time of writing this manuscript, somewhat limited. For example, most recent versions of the all-electron general potential linearized augmented plane wave (LAPW) codes elk⁸ and WIEN2K⁹ only allow non-selfconsistent calculations with mGGA functionals. To date, in plane-wave pseudopotential DFT implementations, the availability of SCAN-based pseudopotentials has also been limited, and to our knowledge, only a norm-conserving library exists,¹⁰ which is lacking kinetic energy density augmentation terms and nonlinear core corrections. For this reason, many calculations published on condensed phase simulations use PBE pseudopotentials,^5,11–15 which at best is an uncontrolled approximation. This type of inconsistency in using pseudopotentials has been studied by Fuchs et al.,¹⁶ and they have shown that using local density approximation (LDA) pseudopotentials in GGA calculations leads to significant errors in the calculated structural properties. We have also found earlier¹⁷ that all-electron properties are much more accurately reproduced when consistent pseudopotentials are used.

Our motivation for this current work was to generate a library of SCAN ultrasoft pseudopotentials for the entire periodic table, based on our previous work.¹⁷ However, we found severe numerical instabilities in both the solution of the atomic all-electron generalized Kohn-Sham equation, which is normally the first step in the pseudopotential generation workflow and again during the pseudopotential construction itself. Indeed, it has been previously observed that SCAN is numerically less stable than GGA exchange-correlation functionals,¹⁸ and recent work has identified shortcomings of the isoorbital indicator component of some mGGA functionals.¹⁹ To remedy this situation, we propose a regularized form of the original SCAN functional (rSCAN), which retains the accuracy of the original form, while improving its stability.

In this paper, we analyze the properties of the isoorbital indicator of SCAN used to connect different approximations of the exchange-correlation energy based on the local bonding environment. We describe a modification that eliminates the unphysical divergence of the exchange-correlation potential, which occurs in some free atoms, while keeping the isoorbital indicator close to the original expression for most regions. We also identify a feature of the switching function in SCAN, which introduces rapidly oscillating regions in the exchange-correlation potential. This causes instabilities in the pseudopotential generation procedure and also affects the discrete representation of the potential on a Fourier grid, which is pivotal in DFT programs using a plane-wave basis set. We propose a small modification, which provides smoother switching, while retaining the superb description of the exchange-correlation energy of the original SCAN functional. We test the rSCAN to establish its closeness to the original form and provide benchmark calculations of fully consistent plane-wave pseudopotential DFT with rSCAN.

A crucial ingredient in SCAN²⁰ and some other mGGA functionals²¹ is the isoorbital indicator function defined as $α=τ−τWτU$⁠, with the definitions of the used quantities listed in Table I. α detects the different local bonding environments, such as covalent single, metallic, or weak bonds, and is used to switch between various local approximations of the exchange and correlation energies, derived for the appropriate bonding type.

However, Furness and Sun¹⁹ have found that derivatives of α with respect to the electron density display divergent behavior at rapidly decreasing electron densities, such as in some free atoms at large distance from the nucleus. The consequence is that the potential itself becomes divergent in these cases, as the derivatives of α appear in the expressions for the potential and are not dampened sufficiently by other terms. Therefore, the resulting exchange or correlation potentials can diverge, for example, in the case of the hydrogen 1s type orbital. Figure 1 shows the SCAN exchange-correlation potential corresponding to the density and kinetic energy density of $ψ(r)=e−rπ$⁠, the H 1s orbital, exhibiting the unphysically divergent behavior.

Furness and Sun¹⁹ suggested an alternative isoorbital indicator function $β=τ−τWτ+τU$⁠, which still displays some divergence, but at a significantly smaller rate, therefore in total resulting in a physically well-behaved potential. In this paper, however, we intend to propose the least amount of modification in the SCAN functional form; hence, we resorted to regularizing the original isoorbital indication function.

The worst divergence occurs in the low-density, single-orbital region, where α ≈ 0, or in case of the 1s orbital example, α = 0 exactly. It is partially due to the rapidly decreasing τ_U in the denominator of α, which leads to numerical instabilities at low-density regions in α. We propose our first regularization in the kinetic energy density of the uniform electron gas as $τU′=τU+τr$⁠, where τ_r = 1 × 10⁻⁴ is a small constant, which only affects α at very low densities.

The second proposed regularization is described as $α′=α3α2+αr$⁠, where α_r = 1 × 10⁻³ is a small constant, and the regularized isoorbital indicator function α′ only differs from the original α function at small values. However, in the single-orbital region, this construction allows vanishing derivatives of α′ with respect to n, ∇n, and τ, therefore minimizing the interference of the switching construction with the physically motivated parts of the exchange and correlation functional expressions. It should be noted, however, that upon introducing τ′ and α′, the exchange energy no longer scales exactly under uniform scaling of the density, although in a practical calculation this effect is expected to remain negligible.

In low-density regions, rSCAN corrects the divergence of derivatives, as well as adjusting the physical interpretation that the isoorbital indicator provides. For example, in the case of isolated noble gas atoms, with the exception of helium, the tail of the valence p orbitals tends to dominate far from the nucleus. According to the original definition, this results in α ≫ 1 at greater distances from the nucleus, corresponding to weak bonds,²⁰ whereas the regularized α indicator returns to zero. This is more similar to helium, where α = 0 everywhere, by construction. Figure 2 compares the original and the regularized isoorbital indicator functions for the isolated Kr atom, also indicating the proportion of the highest contributing orbital type.

The original SCAN functional form includes a switching function, based on the isoorbital indicator. The switching function facilitates a smooth transition between limiting cases, which are constructed observing the constraints based on exact density functional. The functional form of the switching function had been carefully selected, and its parameters were fitted such that the resulting exchange and correlation energies reproduce those of accurate model systems. Even so, the actual form is arbitrary, and we identified the region corresponding to α ≈ 1 as another source of numerical instability. Figure 3 shows the switching function and its first and second derivatives, both contributing to the resulting exchange and correlation potentials. The region around α ≈ 1 is constructed so flat that f⁽ⁿ⁾(1) = 0 for every n, in order to preserve the gradient expansion for the exchange energy in the slowly varying limit. However, this in turn introduces severe oscillations in the derivatives at the surrounding region. These oscillations also manifest in the exchange-correlation potential, as shown for GGA part of the potential in Fig. 4 and mentioned in Ref. 18.

Our intent is to make minimal changes to the switching function, and we found that replacing the region 0 < α < 2.5 by a 7-th degree polynomial removes the oscillatory behavior, while keeping the performance of SCAN similar to the original functional form, although recognizing that we lose the gradient expansion in the slowly varying limit. We fitted the coefficients (supplementary material) of the polynomials such that the derivatives f^(0,1,2)(0) and f^(0,1,2,3)(2.5) are retained and the additional constraint f(1) = 0 is satisfied. Figure 3 compares the original and modified switching functions and their derivatives, demonstrating the improved smoothness, as also evidenced in the practical case of two isolated atoms in Fig. 4.

We implemented rSCAN in the CASTEP²² plane-wave pseudopotential DFT program and the PySCF quantum chemistry package.²³ Self-consistent calculations were performed by solving the generalized Kohn-Sham equations iteratively.^18,24 In CASTEP, ultrasoft pseudopotentials were generated on-the-fly, using the methodology we described elsewhere.¹⁷ We have also pseudized the τ-dependent part, V_τ of the exchange-correlation potential. In our solid-state calculations, we used Monkhorst-Pack k-point grids²⁵ with a 0.02 Å⁻¹ (0.014 Å⁻¹ in case of metals) spacing to sample the Brillouin zone, and the basis_precision: extreme setting in CASTEP for the energy cutoff of the plane-wave basis. PySCF was used to compute the Ar dimer dissociation energies, using the aug-cc-PVQZ basis set²⁶ at the standard grid settings. We used CASTEP to optimize the geometry of the water monomer and hexamer configurations, using a cubic box with 15 Å sides, 750 eV plane-wave cutoff, and the Γ point in the Brillouin zone.

The parameters in the exchange and correlation switching function of the original SCAN were fitted to reproduce the exchange and correlation energies of isolated Ne, Ar, Kr, and Xe atoms, the interaction energies of compressed Ar dimers, and the jellium surface exchange-correlation energy. We compared the accuracy of these quantities, with the exception of the jellium surface exchange-correlation energy, and summarized the results in Table II. For the relative binding energy curve of the Ar dimer at 1.6 Å, 1.8 Å, and 2.0 Å, the mean absolute error of rSCAN is 1.1 kcal/mol, while the figure for the original SCAN was below 1 kcal/mol.

We also benchmarked the rSCAN on some model systems in the literature where results with the original SCAN are available. The set is far from complete, and we note that the literature figures are not consistent: they were obtained by a broad range of codes using different basis sets, in some cases with inconsistent projector augmented-wave (PAW) pseudopotentials. However, our results demonstrate that rSCAN has a performance comparable to the original SCAN functional.

Table III lists the lattice constants of a set of simple solids as calculated with rSCAN and compares them to experiment as well as the original SCAN figures reproduced from the supplementary material of Ref. 20, showing good agreement. A recently published shortcoming of SCAN is the overestimation of magnetic energies of ferromagnetic systems.¹¹ We have found that rSCAN performs similarly, obtaining m = 2.62 μ_B of the spin moments for bcc iron at an optimized lattice constant of 2.84 Å, in good agreement of the SCAN values presented in Ref. 11 m = 2.60 μ_B at the optimized 2.85 Å lattice constant.

Interaction energies of water systems are a very strict test of density functionals, and the original SCAN functional performs remarkably well, predicting the correct energetic ordering of ice polymorphs and water hexamer conformations. With the rSCAN, the water monomer geometry is very close to that of the original SCAN and the dipole moments of the isolated molecule are also in close agreement. We have also calculated the dissociation energies of four low-energy water hexamers, as shown in Table IV, recovering the same energetic ordering as predicted by CCSD(T)³¹ and SCAN, and somewhat improving the absolute values of the energies.

Exchange-correlation functionals based on the mGGA have become increasingly successful, but their implementation in solid-state DFT packages lags behind the theoretical developments. We have implemented the SCAN mGGA functional in a plane-wave DFT program, using ultrasoft pseudopotentials generated with the same functional and solving the electronic problem self-consistently via the generalized Kohn-Sham scheme. To achieve this it was necessary to introduce a regularized form of the SCAN functional that has an improved numerical stability while retaining the accuracy of the original form. We note that the few adjustable parameters, which we imported from SCAN may be reoptimized to further improve the performance, but that is outside of the scope of our current work. Our benchmark calculations illustrate that the proposed rSCAN functional remains transferable and accurate for a broad range of solid state and molecular systems. rSCAN will make the generation of pseudopotential and PAW datasets more straightforward in other packages, while its improved smoothness properties should improve the stability of any DFT implementation where the exchange-correlation functionals need to be represented on a grid.

See supplementary material for the numerical values of the polynomial coefficients of the modified exchange and correlation switching functions.

We would like to thank Chris Pickard, Philip Hasnip, and Dominik Jochym for useful discussions. Both authors acknowledge support from the Collaborative Computational Project for NMR Crystallography (CCP-NC) and UKCP Consortium, both funded by the Engineering and Physical Sciences Research Council (EPSRC) under Grant Nos. EP/M022501/1 and EP/P022561/1, respectively. Computing resources were provided by the STFC Scientific Computing Department’s SCARF cluster.