Comparative study of the PBE and SCAN functionals: The particular case of alkali metals

The SCAN meta-generalized gradient approximation (GGA) functional is known to describe multiple properties of various materials with different types of bonds with greater accuracy, compared to the widely used PBE GGA functional. Yet, for alkali metals, SCAN shows worse agreement with experimental results than PBE despite using more information about the system. In the current study, this behavior for alkali metals is explained by identifying an inner semicore region which, within SCAN, contributes to an underbinding. The inner semicore push toward larger lattice constants is a general feature but is particularly important for very soft materials, such as the alkali metals, while for harder materials the valence region dominates. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5092748


I. INTRODUCTION
Presently, the most common theoretical approach to calculating the properties of solids and molecules is Kohn-Sham density functional theory (KS-DFT). 1 The main accuracy restricting factor of this method is the functional form of the exchange-correlation energy, Exc. For local and semilocal functionals, Exc is given by where exc is the exchange-correlation energy per volume unit and is a function of local electronic properties, such as the electron density, electron density gradient, or kinetic-energy density (KED). The simplest approximation of exc is the local density approximation (LDA). 1 The next step of functional development was to add a functional dependence on the gradient of the density. This led to the generalized gradient approximations (GGAs), 2,3 which have better accuracy in multiple cases. In the meta-GGA (MGGA) functionals, the KED and/or the Laplacian of the density are also used in the parameterization of exc. Several MGGAs with different constraints and goals have been developed (see Ref. 4 for a review), and benchmarks of these different functionals have shown how MGGAs can improve the overall accuracy compared to GGAs. [5][6][7][8] The improved performance can, depending on the point of view, be related to the MGGAs being able to distinguish more bonding situations, [9][10][11] better fit reference data, 12 or satisfy more exact constraints. 13 One of the recent MGGAs, which has gained considerable attention, is the SCAN functional. 13 For instance, it has shown successes in calculating the formation enthalpy of various solids 14 or the structural and energetic properties of ferroelectric materials. 15 On the other hand, SCAN performs poorly for the magnetic properties of transition metals. [16][17][18][19] It is natural to compare SCAN with the PBE GGA. 3 First of all, they are constructed following a similar philosophy of constraint satisfaction. Furthermore, PBE can be considered as a good functional for solids since it gives reasonable equilibrium lattice constants, a 0 , and cohesive energies, E coh (see, e.g., Ref. 5). While it is possible to construct GGA functionals which give better results than PBE for the lattice constants, 20-23 these will tend to overestimate the cohesive energies of solids. 5 Thus, for SCAN to be a systematic improvement on PBE for solids, one requirement would be that it simultaneously improves on the lattice constants and the cohesive energy. Numerical tests have shown that on average the SCAN does exactly this for a wide range of solids. 5,24 The improvement of SCAN over PBE is, however, not universal. A close look at the results in Ref. 5 also reveals how SCAN performs disappointingly for most alkali metals. This is illustrated in Fig. 1 where SCAN is compared to LDA and PBE for a few selected materials. Considering first Si and Ge, which we use for illustrative purposes as representatives for a large group of systems (covalent semiconductors), the well known tendencies of LDA to overestimate the cohesive energies, thereby underestimating the lattice constant, and PBE to overestimate the lattice constants, thereby underestimating the cohesive energies, can be seen. For Si and Ge, SCAN systematically improves both the lattice constant and cohesive energies over PBE. For the alkali metals on the other hand, we can see that SCAN underbinds even more, i.e., gives even larger lattice constants (except for Na) and even smaller cohesive energies, than PBE ( Fig. 1).
In the present study, we aim at a detailed understanding of how the poor performance of SCAN for the alkali metals is related to its functional form. Apart from the obvious interest in the alkali metals themselves, understanding the disappointing performance of SCAN for this class of materials is also important for developing more accurate density functionals in general. While SCAN may perform well in statistical studies, where the focus is on average errors for databases containing a large number of strongly ionic and covalently bonded materials, these studies may somewhat hide a systematic problem for the more weakly bonded alkali metals because these systems only make up a small subset of the database. Actually, the density distribution in the alkali metals is rather particular. The bonding region is characterized by both the density and reduced gradients being low. This means that the correlation energy becomes comparable to the exchange and that regions of exc that are otherwise not sampled are probed.

II. METHODOLOGY
As will be discussed below, we will focus on the exchange energy in the present analysis. To describe the analytical form of an exchange functional, it is common to define the enhancement factor where e LDA x = −Cxn 4 3 [Cx = (3/4)(3/π) 1/3 , atomic units are used throughout this work] is the LDA exchange-energy density for the homogeneous electron gas 1 and n = ∑ N i=1 ψi 2 is the electron density. For GGA functionals, Fx depends on the gradient of the density, ∇n, while for MGGA functionals it also depends on the noninteracting KS KED τ KS = (1 2) ∑ N i=1 ∇ψi 2 (in the present work, we are not concerned with Laplacian-dependent MGGAs). In the following, we will use dimensionless expressions to characterize the density, namely, the reduced density gradient and reduced KED where τ TF = (3 10)(3π 2 ) 2 3 n 5 3 is the Thomas-Fermi (TF) KED 25,26 which is exact for the homogeneous electron gas. Here, we note that in our previous studies 8, 10 and others, 11,12 τ KS /τ TF was instead labeled as t −1 . In iso-orbital regions where the density is dominated by one orbital, the KED is given exactly by the von Weizsäcker form 27 which makes a convenient measure of how much the density n at a point of space is dominated by a single orbital. 9 Since one can write τ vW /τ TF = 5p/3, then Note that τ vW is a strict lower bound to the KED 9,28-30 so that 5p/3 is a lower bound to t. As mentioned in Sec. I, the goal of the present work is to rationalize the SCAN results on the alkali metals and to understand the worsening in the performance compared to PBE. For this purpose, potassium is the case study that will be considered in Sec. III. The analysis will consist of a careful comparison of the PBE and SCAN enhancement factors Fx. The calculations were carried out with the WIEN2k code. 31 The SCAN calculations were done nonself-consistently using the PBE densities and orbitals 5,32 so that the only difference in the total energy stems from the functional form of exc. Note that in our previous work, 5 the self-consistent effects were estimated to be quite small, below 0.01 Å in most cases, except for the van der Waals systems where they could be larger. For the spatial distribution analysis of the exchange-correlation energy, the sampling of the Voronoi cell of one atom was done on an equidistant radial mesh for 400 different directions from the atom. The distance between the sample points is the same for both volumes, resulting in a larger number of samples for the expanded structure.

III. ANALYSIS
We start by showing in Fig. 2 the difference between PBE and SCAN of the exchange-correlation energy Exc as a function of lattice constant a for the alkali metal potassium. Since all terms in the total energy except Exc are the same for the PBE and SCAN calculations and the SCAN energy is evaluated with PBE density, the slope of E SCAN xc − E PBE xc is directly related to the difference between the equilibrium lattice constants a PBE 0 and a SCAN 0 . As seen in Fig. 2, the slope of the exchange-energy difference E SCAN As a direct consequence, the SCAN equilibrium lattice constant a 0 is "pushed" toward a larger value than the one obtained with PBE. As PBE already overestimates a 0 of potassium (as well as Rb and Cs, see Fig. 1), then SCAN worsens the agreement. Thus, it is the exchange component of E SCAN xc that is responsible for the overestimated lattice constant of K. An interesting feature about the alkali metals is that the low density in the bonding region means that the correlation energy density is comparable to the exchange energy density. Figure 2 shows that the correlation energy exhibits the opposite trend and somewhat compensates for the "push" toward larger volume by the SCAN exchange. However, the compensation is only partial and the slope of the total exchange-correlation energy curve remains negative. In the following, we will thus focus on a detailed analysis of the exchange energy, which is the driving force behind the overestimated lattice constant.
The variation of Ex with respect to the lattice constant a can be explained in terms of changes in the density n and enhancement factor Fx [Eq. (2)]. We separate these two effects by expanding the exchange energy shifts and keeping only terms that are first order in the perturbation δex ≈ −Cx(n + δn) 4 3 ( The first part, δe enha respectively. This choice of volumes for calculating δex is somehow arbitrary; however, the linear monotonic behavior of Ex seen in Fig. 2 shows that it is unimportant for the conclusion. For the sampling of δex, we have chosen grids of equidistant points starting at the atomic positions. Since the grid contains more points for the larger volume, an additional contribution, δe new x = −Cxn 4 3 Fx, representing the new sample points has to be taken into account and added to δe enha x and δe dens x to get the full δex. Figure 3(a) illustrates these components to the difference δe SCAN x − δe PBE x as functions of distance from the potassium atom. δe new x is small and only appears for distances larger than 4 bohrs because of the shape of the Voronoi cell.
In the valence region that we define as the distance beyond 2 bohrs, the density and enhancement terms tend to cancel each other. Actually, δe enha,SCAN x −δe enha,PBE x is positive, which indicates that the SCAN exchange enhancement factor is less sensitive to a change in the volume. On the other hand, δe dens,SCAN In Fig. 3(a), one can also identify a region, between 1 and 2 bohrs, where the total exchange energy density difference δe SCAN x − δe PBE x is negative, thus forcing SCAN lattice constant to be larger than the PBE one. This region of negative values is clearly due to the component δe enha x , i.e., a faster increase in the magnitude of F SCAN x in a region with a high density when the volume gets bigger, which is particularly important for the lattice constant [see Fig. 3 Actually, a strong influence on the equilibrium lattice constant coming from the region between 1 and 2 bohrs is at first sight somewhat surprising as one would associate it with an inner semicore region. To understand its origin, we first show 5p/3 and t, Eqs. (3) and (4), and the normalized orbital densities of a free potassium atom as functions of the distance to the atom in Fig. 4. It is seen that the 1-2 bohr region is indeed dominated by the 3s and 3p semicore orbitals. In this region, the electron density n is very large compared to the valence region such that even small changes δFx in the enhancement factor have a large impact on the exchange energy [since δFx is multiplied by n 4/3 , Eq. (8)] and thereby lead to large values of δe enha x . From Fig. 4, it is also possible to understand why the 4s shell also contributes to δe SCAN x − δe PBE x below 2 bohrs. Indeed, at a distance around 1.6 bohrs from the atom, the outer lobe of the 4s shell starts to become important. Since the 4s shell is strongly perturbed by the chemical bonding, then its influence on δFx in the 1.6-2.0 bohr region should be important.
To obtain insight into the individual contributions 33 to δFx in δe enha x , we proceed by expanding it as We first consider the variations δp and δt due to an expansion of the volume. These are depicted for the inner semicore region in Fig. 5(a). As expected, the reduced density gradient p gets larger when the volume increases, i.e., δp > 0, especially for a distance larger than 1.
that is not available in the gradient of the electron density, something that is a premise for the development of MGGA. The behavior of t can be understood from Fig. 4. As observed earlier, 33,34 one can clearly identify peaks in the reduced density gradient p that are located in transition regions where the dominating shell is switching from one to another. In these inter-shell regions, t is substantially larger than 5p/3 (see Fig. 4) so that α is, as expected from Eq. (7), larger in such regions with contributions coming from different shells. 9 In the inner semicore (1-2 bohrs) region, α is small, reflecting how it is dominated by orbitals of similar shape (n = 3). In the solid, the inner semicore region becomes increasingly dominated by the 3s and 3p orbitals as the unit cell expands, thereby becoming more "atomiclike." As p hardly changes [δp ≈ 0 for r < 1.6 bohrs, Fig. 5(a)], the smaller values of α are the result of smaller values of t (δt < 0). The partial derivatives in Eq. (9), depicted in Fig. 5(b), reflect the dependence of the functional on changes in p and t around their values at the equilibrium lattice constant a 0 . Figure 5(b) also shows that ∂F SCAN x ∂p is approximately twice larger than ∂F PBE x ∂p. Actually, the large derivatives of SCAN are somewhat surprising because earlier illustrations (see, e.g., Fig. 1 of Ref. 13) give more the impression of a smooth and subdued functional form. However, in Fig. 6(a), we show that the smooth behavior is mainly along lines of constant values of α. Perpendicular to these lines, F SCAN x shows a somewhat more "bumpy" behavior. Such bumps lead to an erratic behavior of the derivatives, as shown in Figs. 6(b) and 6(c).
The large positive ∂F SCAN x ∂p will, when multiplied by the negative −Cxn 4/3 , Eq. (8), and the positive δp in the inner semicore region, Fig. 5(a), contribute to the negative slope of E SCAN observed in Fig. 2. We have already discussed how δt shows a different behavior than δp in the inner semicore region and, in principle, a MGGA could compensate for this contribution in its dependence on the KED. However, as both δp and the corresponding partial derivative ∂F SCAN Equation (9) underlines how the partial derivatives are an important factor in determining energy differences and thereby the performance of a MGGA. This would suggest that they should be routinely shown when reporting a new functional. We should also point out that the partial derivatives in Eq. (9) are part of the analytical expression of the MGGA potential for self-consistent calculations, 35 and the behavior observed in Figs. 6(b) and 6(c) could thus be responsible for SCAN resulting in a large overestimation of the magnetic moment in itinerant transition metals. [16][17][18][19] In this context, it is interesting to note that a fixed-spin moment calculation, which involved only the SCAN energy and used the PBE potential, resulted in the same overestimation of the magnetic moment as a self-consistent calculation. 19 The expansion with respect to volume, Eq. (8), highlights how features of the analytic form of the energy functional are directly related to the potential. Similar to the expansion with respect to volume, the exchange- correlation energy can also be expanded with respect to the magnetic moment.
The observation of an erratic behavior of the functional also falls in line with recent observations of relative strong grid dependence of SCAN results. 36,37 Such a grid dependence has also been analyzed previously in Ref. 38 for other MGGA functionals, where it was also pointed out that poor convergence with grid size can lead to unintended contributions to the energy differences.
It is noteworthy that among the 44 solids tested in Ref. 5, the alkali and alkaline earth metals are the only ones for which lattice for Al is shown in Fig. 7. The inner semicore region can be identified between 0.5 and 1.0 bohrs and does indeed have a negative total δe SCAN x − δe PBE x due to the δFx contribution. It will therefore push SCAN to have a larger lattice constant compared to PBE, similar to what was observed for potassium (Fig. 3). However, contrary to potassium, the influence of the valence region is much larger than the inner semicore region. The valence contribution is mainly positive which results overall in a smaller SCAN equilibrium lattice constant than with PBE. We have performed the same analysis for FCC-Ca (not shown) Si in the diamond lattice. Also here, an inner semicore push toward larger lattice constants due to δe enha x can be identified. However, this is compensated by the valence region, which means that SCAN and PBE lead to very similar lattice constants for Ca.
Finally, one could also argue that the SCAN underbinding of the alkali metals should be cured by explicitly adding contributions for the long-range dispersion interactions. 36 Such corrections will however universally strengthen the bonding and thus lead to a worse performance in cases such as Ca and Al where SCAN already tends to overbind. Thus, two strategies could be followed to cure the problem of SCAN for the alkali metals: either by modifying the functional form such that the results for the alkali metals are improved or by adding a term that explicitly accounts for the dispersion term. However, in the latter case, the functional form of SCAN should also be modified in order to avoid an overbinding for other systems like Ca or Al.

IV. SUMMARY
In the current study, we have analyzed in detail the results obtained with the MGGA functional SCAN for the alkali metals. For these systems, SCAN is less accurate than the more simple GGA functional PBE. SCAN has a clear tendency to underbind the alkali metals; i.e., the equilibrium lattice constants are too large and the cohesive energies are too small. We have shown that this behavior of SCAN is attributed to an inner semicore push toward larger lattice constants, which was revealed by a careful comparison of the PBE and SCAN enhancement factors. Such an inner semicore push toward larger lattice constants can probably be identified for many materials; however, it is the most important mechanism for soft materials such as alkali metals, while for harder materials (e.g., semiconductors and ionic insulators) the valence region dominates (as shown for Al).
A detailed analysis, such as the one that we have presented, leads to a clear understanding of the failures or unexpected results that a functional produces. A functional may have an analytical form that is inappropriate within a particular regime, e.g., for low densities or high density gradients, and the precise identification of the problem in the functional form may give a clue of how to modify the functional form to cure the problem. Our study furthermore highlights the importance of the partial derivatives in determining energy differences and suggests that these should be routinely shown when reporting a new functional.