Regularized SCAN functional

TL;DR

This paper introduces a regularized version of the SCAN density functional to improve numerical stability, enabling reliable pseudopotential generation while maintaining similar performance to the original functional.

Contribution

The authors develop a modified SCAN functional that eliminates numerical instabilities without significantly altering its accuracy.

Findings

Regularized SCAN matches original performance.

Enhanced numerical stability for pseudopotential generation.

Maintains accuracy while improving reliability.

Abstract

We propose modifications to the functional form of the SCAN density functional to eliminate numerical instabilities. This is necessary to allow reliable, automatic generation of pseudopotentials (including PAW potentials). The regularized SCAN is designed to match the original form very closely, and we show that its performance remains comparable.

Tables

Table 1. Table 1: Definition of the quantities based on the Kohn-Sham orbitals used in this work.

Kohn-Sham orbitals	${\psi }_i$
orbital kinetic energy density	$\tau =1/2{\sum }_i^{\text{occ}}{\left|\nabla {\psi }_i\right|}^2$
electron density	$n={\sum }_i^{\text{occ}}{\left|{\psi }_i\right|}^2$
Weizsäcker kinetic energy density	${\tau }_W=\frac{{\left|\nabla n\right|}^2}{8n}$
kinetic energy density of the uniform electron gas	${\tau }_U=\left(3/10\right){\left(3{\pi }^2\right)}^{2/3}{n}^{5/3}$

Table 2. Table 2: Exchange and correlation energies of isolated noble gas atoms, in hartrees. Original SCAN values are obtained from Ref. Sun  et al. , 2015 , reference values from Refs. Becke, 1988 ; Chakravorty  et al. , 1993 ; McCarthy and Thakkar, 2011 .

		Ne	Ar	Kr	Xe
$E_{\text{x}}$	SCAN	-12.164	-30.263	-94.068	-179.325
	rSCAN	-12.163	-30.298	-94.199	-179.632
	ref.	-12.108	-30.188	-93.890	-179.200
$E_{\text{c}}$	SCAN	-0.345	-0.691	-1.756	-2.899
	rSCAN	-0.345	-0.695	-1.768	-2.914
	ref.	-0.391	-0.723	-1.850	-3.000
$E_{\text{xc}}$	SCAN	-12.508	-30.954	-95.826	-182.218
	rSCAN	-12.508	-30.993	-95.966	-182.546
	ref.	-12.499	-30.911	-95.740	-182.200

Table 3. Table 3: Equilibrium lattice constants (Å) of a selection of metallic and semiconductor solids (a subset of “LC20” in Ref. Sun  et al. , 2015 ), computed using the rSCAN functional. Experimental values, corrected for zero point anharmonic expansion, were taken from Ref. Hao  et al. , 2012 , and reference SCAN values from Ref. Sun  et al. , 2015 .

	Li	Na	Ag	C	Si	SiC	LiF	MgO
Expt.	3.451	4.207	4.063	3.555	5.422	4.348	3.974	4.188
SCAN	3.460	4.190	4.079	3.550	5.424	4.349	3.980	4.206
rSCAN	3.453	4.197	4.039	3.555	5.441	4.353	3.964	4.200

Table 4. Table 4: Dissociation energies (meV/monomer) of a few low-energy water hexamers conformations, the equilibrium bond length (Å), bond angle and dipole moment (Debye) of the water molecule. Reference hexamer dissociation values are computed by CCSD(T) Santra  et al. ( 2008 ) , while the geometry of the water molecule is from Ref.  Benedict  et al. , 1956 and its dipole moment from Ref.  Dyke and Muenter, 1973 . SCAN values were obtained from Ref.  Sun  et al. , 2016 .

	prism	cage	book	chair	$r_{\text{OH}}$	${\theta }_{\text{HOH}}$	$\mu$
ref.	348	346	339	332	0.957	104.5°	1.855
SCAN	377	376	370	360	0.961	104.5°	1.847
rSCAN	359	358	356	348	0.959	104.4°	1.847

Full text

I Introduction

First principles modelling 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 costJones (2015). The major source of inaccuracy in Kohn-Sham DFT calculationsKohn and Sham (1965) 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) functionalPerdew et al. (1996), 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 transferableSun et al. (2016), and improved the DFT description of a wide range of systems, such as liquid water and ice,Chen et al. (2017) semiconductor materialsRemsing et al. (2017) or metal oxidesGautam and Carter (2018).

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 elkelk and WIEN2KBlaha et al. (2018) 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 existsYao and Kanai (2017), which is lacking kinetic energy density augmentation terms and non-linear core corrections. For this reason, many calculations published on condensed phase simulations use PBE pseudopotentialsFu and Singh (2018); Chen et al. (2017); Janthon et al. (2014); Dorner et al. (2018); Bonati and Parrinello (2018); Isaacs and Wolverton (2018), which at best is an uncontrolled approximation. This type of inconsistency in using pseudopotentials has been studied by Fuchs et al Fuchs et al. (1998), and they have shown that using LDA pseudopotentials in GGA calculations leads to significant errors in the calculated structural properties. We have also found earlierBartók and Yates (2019) 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 Bartók and Yates (2019). 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 functionalsYang et al. (2016), and recent work has identified shortcomings of the iso-orbital indicator component of some mGGA functionalsFurness and Sun (2019). 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 iso-orbital indicator of SCAN, used to connect different approximations of the exchange-correlation energy based on the local bonding environment. We describe a modification which eliminates the unphysical divergence of the exchange-correlation potential which occurs in some free atoms, while keeping the iso-orbital 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.

II Regularized SCAN

II.1 The iso-orbital indicator function

A crucial ingredient in SCANSun et al. (2015) and some other mGGA functionalsSun et al. (2013) is the iso-orbital indicator function, defined as $\alpha=\frac{\tau-\tau_{W}}{\tau_{U}}$, with the definitions of the used quantities listed in Table 1.

$\alpha$ detects 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 SunFurness and Sun (2019) have found that derivatives of $\alpha$ with respect to the electron density display divergent behaviour 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 $\alpha$ 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 1$s$ type orbital. Figure 1 shows the SCAN exchange-correlation potential corresponding to the density and kinetic energy density of $\psi(r)=\nicefrac{{e^{-r}}}{{\sqrt{\pi}}}$, the H $1s$ orbital, exhibiting the unphysically divergent behavior.

Furness and Sun suggested an alternative iso-orbital indicator function $\beta=\frac{\tau-\tau_{W}}{\tau+\tau_{U}}$, that still displays some divergence, but at a significantly smaller rate, therefore in total resulting in a physically well-behaved potential.Furness and Sun (2019) 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 iso-orbital indication function.

The worst divergence occurs in the low-density, single-orbital region, where $\alpha\approx 0$, or in case of the 1s orbital example, $\alpha=0$ exactly. It is partially due to the rapidly decreasing $\tau_{U}$ in the denominator of $\alpha$, which leads to numerical instabilities at low-density regions in $\alpha$. We propose our first regularization in the kinetic energy density of the uniform electron gas, as $\tau_{U}^{\prime}=\tau_{U}+\tau_{r}$, where $\tau_{r}=1\times 10^{-4}$ is a small constant, which only affects $\alpha$ at very low densities.

The second proposed regularization is described as $\alpha^{\prime}=\frac{\alpha^{3}}{\alpha^{2}+\alpha_{r}}$, where $\alpha_{r}=1\times 10^{-3}$ is a small constant, and the regularized iso-orbital indicator function $\alpha^{\prime}$ only differs from the original $\alpha$ function at small values. However, in the single-orbital region this construction allows vanishing derivatives of $\alpha^{\prime}$ with respect to $n$, $\nabla n$ and $\tau$, 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 $\tau^{\prime}$ and $\alpha^{\prime}$, 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 iso-orbital indicator provides. For example, in the case of isolated noble gas atoms, with the exception of helium, the tail of the valence $p$ orbitals tend to dominate far from the nucleus. According to the original definition, this results in $\alpha\gg 1$ at greater distances from the nucleus, corresponding to weak bondsSun et al. (2015), whereas the regularized $\alpha$ indicator returns to zero. This is more similar to helium, where $\alpha=0$ everywhere, by construction. Figure 2 compares the original and the regularized iso-orbital indicator functions for the isolated Kr atom, also indicating the proportion of the highest contributing orbital type.

II.2 The switching function

The original SCAN functional form includes a switching function, based on the iso-orbital 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 $\alpha\approx 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 $\alpha\approx 1$ is constructed so flat that $f^{(n)}(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 Figure 4 and mentioned in Ref. Yang et al., 2016.

Our intent is to make minimal changes to the switching function, and we found that replacing the region $0<\alpha<2.5$ by an 7-th degree polynomial removes the oscillatory behaviour, 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 coefficientsSI of the polynomials such that the derivatives $f^{(0,1,2)}(0)$, $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 Figure 4.

III Results

We implemented rSCAN in the CASTEPClark et al. (2005) planewave-pseudopotential DFT program and the PySCF quantum chemistry packageSun et al. (2017). Self-consistent calculations were performed by solving the generalized Kohn-Sham equations iterativelyYang et al. (2016); Perdew et al. (2017). In CASTEP, ultrasoft pseudopotentials were generated on-the-fly, using the methodology we described elsewhereBartók and Yates (2019). We have also pseudized the $\tau$-dependent part, $V_{\tau}$ of the exchange-correlation potential. In our solid-state calculations, we used Monkhorst-Pack $k$-point gridsMonkhorst and Pack (1976) with a 0.02 Å*-1* (0.014 Å*-1* in case of metals) spacing to sample the Brillouin zone, and the basis_precision : extreme setting in CASTEP for the energy cutoff of the planewave basis. PySCF was used to compute the Ar dimer dissociation energies, using the aug-cc-PVQZ basis setWoon and Dunning Jr. (1993) 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 planewave cutoff and the $\Gamma$ 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 2. 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 PAW pseudopotentials. However, our results demonstrate that rSCAN has a performance comparable to the original SCAN functional.

Table 3 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. Sun et al., 2015, showing good agreement. A recently published shortcoming of SCAN is the overestimation of magnetic energies of ferromagnetic systems.Fu and Singh (2018) We have found that rSCAN performs similarly, obtaining $m=2.62\,\mu_{B}$ of the spin moments for bcc iron at the optimized lattice constant of 2.84 Å, in good agreement of the SCAN values presented in Ref. Fu and Singh, 2018 $m=2.60\,\mu_{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 4, recovering the same energetic ordering as predicted by CCSD(T)Santra et al. (2008) and SCAN, and somewhat improving the absolute values of the energies.

IV Conclusions

Exchange-correlation functionals based on the meta-Generalized Gradient Approximation 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 re-optimized 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.

Supplementary Material

See Supplementary MaterialSI for the numerical values of the polynomial coefficients of the modified exchange and correlation switching functions.

Acknowledgements.

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 numbers EP/M022501/1 and EP/P022561/1, respectively. Computing resources were provided by the STFC Scientific Computing Department’s SCARF cluster.