Why Density Functionals Should Not Be Judged Primarily by Atomization Energies

While most molecules and solids are spin-unpolarized, most chemically-active atoms are partly spin-polarized. As a result, the errors of the spin-dependence of a density functional are much more troublesome for atomization energies than they are for typical reaction or formation energies. This observation explains why the atomization energy errors of approximate functionals do not correlate with their other errors, and why the errors of atomization energies for a given functional can be radically reduced by fitting the energies of the atoms. We present an illustrative example from the recent nonempirical construction of the SCAN meta-generalized gradient approximation.


Introduction
This year is the 50th anniversary of Kohn-Sham density functional theory [1], which from a humble beginning has grown up to be the most widely-used method of electronic structure calculation in condensed matter physics, quantum chemistry, and materials engineering.The popularity of this theory arises first because it is computationally efficient, with an orbital structure like that of Hartree or Hartree-Fock theory, and second because it is usefully accurate.There is an underlying exact theory for the ground-state energy and electron density, which inspires the continuing and successful (if slow) search for more accurate but still computationally tractable approximations to the needed density functional for the exchange-correlation energy.This is also the tenth anniversary of collaboration among the authors of this paper.In 2005, three of us (GIC, AR, JPP) published the first [2] of 35 collaborative research papers in density functional theory.(One year later, another of us (JS) started research in this theory as a graduate student with JPP.)That first collaborative paper raised questions about the atomization energies of molecules, which up to that time had been widely regarded as a "gold standard" of accuracy for density functional approximations.In 2015, when three of us (JS, AR, and JPP) have proposed what we regard as an optimally accurate functional that preserves computational efficiency [3], we return to the issue in the title of this article and present an explanation which is supported by numerical evidence at this high level of approximation.

Density and Spin-Density Functional Theories
The original Kohn-Sham theory [1] assumed that the electrons (with mutual Coulomb repulsion) see a static, spin-independent, and multiplicative external potential v r ( )  (typically but not necessarily the Coulomb attraction to the nuclei).Under this condition, the exact exchange-correlation energy E xc is a functional E xc [n] of the electron density distribution n r ( )  ; i.e., there is a rule (in practice an uncomputable one) that assigns one energy E xc to each possible electron density distribution n r ( )  . Apart from E xc , the rest of the total energy is computable exactly.
Later this theorem was generalized [4] to electrons that see spin-dependent external potentials v r ↑ ( ) is a functional of the electron spin-density distributions n r ↑ ( )  and n r ↓ ( )  .In nature, non-relativistic elec- trons in the absence of a magnetic field always see an exter- , so these two exact theories should yield the same exact total density n n n = + ↑ ↓ and the same exact total energy (which includes E xc as one term).
But approximations to these two exact theories need not agree.We can see this already at the simplest level, the local spin-density approximation [5] , is the exchange-correlation energy per particle of an electron gas with uniform or  r -independent spin densities n ↑ and n ↓ .Eq. ( 1) is exact for spin densities that are uniform or that vary slowly-enough over space.The local density approximation will give the same total energy as Eq. ( 1) only when the system is fully spin-unpolarized, with n r n r n r  2. This is the case for the He atom ground state, but not for the fully spin-polarized H atom ground-state with n r n r Yet both atoms see a spin-independent external po- tential.
Tong [5] made an early LDA calculation of the cohesive energy of metallic sodium, i.e., the energy per atom to convert the solid to a vapor of atoms.He found a value 23% larger than the experimental result.Kohn [6] argued that LDA should be accurate for solid sodium (in which the spin-unpolarized valence electron density resembles a uniform gas), and that the error must be in the energy of the atom.Gunnarsson, Lundqvist, and Wilkins [7] applied LSDA to the partly spinpolarized sodium atom, and found a cohesive-energy error of only 4%.It is understandable that approximate spin-density functionals can be more accurate than the corresponding approximate density functionals.The spin-density functionals get more input information about the electrons, and don't have to be as "smart" as the density functionals would have to be to achieve the same level of accuracy.
The cohesive energy of a solid was never of great interest to condensed matter physicists, who were typically more interested in crystal structures, lattice constants, bulk moduli, phonon frequencies, band structures, etc.However, around 1993 chemists also started to be interested in density functional theory.To a chemist, the atom is rightly the building block of all ordinary matter.As a result, an era began in which the atomization energies of molecules were a kind of "gold standard" for density functionals.(In theoretical chemistry, atomization energies were traditionally used for the calculation of the standard enthalpies of formation, because early wavefunction theory was unable to yield reliable energies for elements in the solid state.In physical chemistry, atomization energies are not required, as the enthalpies of formations require only the standard-state enthalpies of the reagents and the product.)In those earlier days, empirical functionals were fitted to large data sets of atomization energies, and all functionals were judged by their error statistics for atomization energies.That era is now ending, with the availability of much larger and more diverse molecular data sets and direct calculations of the energies of solid state elements.
Perdew and Schmidt [8] have proposed a "Jacob's Ladder" of spin-density functional approximations: The inputs at each point  r are the local spin densities (the only inputs in LSDA), their gradients (added in a generalized gradient approximation or GGA), the orbital kinetic energy densities , (further added in a meta-GGA), etc.The inputs displayed explicitly in Eq. ( 3) are all available at each point  r in a Kohn- Sham calculation, making the first three rungs (LSDA, GGA, and meta-GGA) semilocal and computationally efficient in this hierarchy.The ingredients not displayed explicitly in Eq. ( 3) are fully nonlocal, so the higher rungs (hybrid-like and random phase approximation (RPA)-like approximations) are computationally more demanding.
Given that the exact total energies for atoms, molecules and solids in the absence of magnetic fields must be the same, we can ask whether for example MAE(DFT)-MAE(SDFT) gets smaller as we climb up the ladder and approach the exact theory.
Here MAE is the mean absolute error for some property over a data set, applying the same functional in density functional (DFT) and spin-density functional theory (SDFT).We hope to answer this question numerically in future work.The answer is currently unknown, since nearly all studies are based on SDFT.

Atomization Energy Puzzle
As we climb the Jacob's Ladder of approximate spin-density functionals from LSDA (which is more or less unique) to various flavors of GGA and meta-GGA, there is a clear decrease in the MAE of the atomization energies of molecules and solids, as one would expect.
Following suggestions from Dewar [9] and Delley [10], as further explored by Delley [11], and Grimme [12], Csonka, Ruzsinszky, Tao and Perdew [2] found that the MAE's of the atomization energies for functionals on all rungs of the ladder could be reduced, sometimes drastically, by using atomic energies shifted by the addition of empirical constants, with one set of atomic constants for each functional.{Very recently, optimized atomic correction energies were successfully applied to direct RPA (dRPA) and dual-hybrid dRPA75 [13].}Large atomic shifts indicated an inconsistency between a functional's description of atoms and its description of molecules built up from those atoms.(The shifts were almost negligibly small for the TPSS [14] meta-GGA.)The fact that the empirically-corrected MAE's were so small indicated that most functionals were more consistent for molecules than for their constituent atoms.Thus the functionals were expected to be (and are) much more accurate for reaction energies (e.g., 2H 2 O → 2H 2 + O 2 ) and formation energies involving no free atom than for atomization energies.
Goerigk and Grimme [15] framed the problem in a different way.They looked at many large data sets of molecular properties, and at many functionals, and found that the errors that a functional makes for atomization energies are almost uncorrelated with the errors that it makes for other molecular properties.
Those are the facts, but what is the explanation?Atoms obviously differ from molecules and solids in several ways: (1) The density of an atom is more confined, and further from the slowly-varying limit, than is the density of a molecule or solid.(2) Open-subshell atoms like carbon might be expected to have smaller gaps between the highest-occupied and lowest unoccupied orbital energies than molecules or semiconducting solids do, leading to near-degeneracy correlation, although SCAN calculations (Table 1) do not support this expectation.(3) Open-subshell atoms like have groundstate densities that are non-spherical, with shapes very different from the shape of the external potential.(4) Open-subshell atoms are partly spin-polarized, while most molecules and solids are spin-unpolarized, so spin-polarization errors can be more troublesome for atomization energies than for reaction or formation energies that involve no free atom.So which of these differences explains the atomization energy puzzle?The early work [7] on the cohesive energy of sodium metal, described in section 2 of this paper, suggests that the spindependence of the approximate functionals provides the right explanation.But that work (largely unknown to chemists) was done on the lowest and least-accurate rung of the ladder.Complicating the answer is the fact that the exchange-correlation energy can be written as = + , the sum of a large exchange energy and a small correlation energy.Nearly every approximate density functional for exchange satisfies the exact spin-scaling relation [16] E Thus, given a sufficiently accurate E x [n], errors in the spin dependence must arise from the relatively small correlation en- . To show that these are responsible for the atomization energy puzzle on the higher rungs of the ladder, we need a very accurate spin-density functional.As discussed in the next section, we have one now.

Exact Constraints and the SCAN Meta-GGA
Although the exact is incomputable, there are exact formal expressions for it which can be used to derive its mathematical features.These features can be regarded as exact constraints, and used to construct approximations that are fully or partly nonempirical.Since the exact constraints are universal (independent of the external potential), constraint-based approximations can be reliable even for systems unlike those on which they have been tested and benchmarked.
The early successes of LSDA, even for metal surfaces, were surprising and demanded an explanation.The explanation [17,18] was not only that LSDA was appropriately normed on the only density for which it could be exact (the electron gas of uniform spin densities), but also that LSDA satisfied a certain subset of exact constraints that it inherited from that appropriate norm.Generalized gradient approximations (GGA's) like the PBE of Ref. [19] and meta-GGA's like the TPSS of Ref. [14] were then constructed to satisfy additional exact constraints.The added ingredients on higher rungs of Jacob's Ladder can be used to satisfy exact constraints that cannot be satisfied on lower rungs.
Recently three of us have proposed SCAN [2], a strongly constrained and appropriately normed meta-GGA.We believe that SCAN is nearly as accurate as a computationally-efficient semilocal exchange-correlation functional can be.SCAN is the first functional to satisfy all 17 known exact constraints that a semilocal functional can, including a tight new lower bound [20] on the exchange energy.But there are still infinitely many ways to satisfy these 17 exact constraints.Thus, to the previous meta-GGA appropriate norms (the electron gas of uniform (5) (6) or slowly-varying spin densities, and the one-electron atom), SCAN adds the exchange energies and correlation energies of rare-gas atoms and other nonbonded systems in which the exact exchange-correlation hole remains close to its electron (as it must if a semilocal functional is to be very accurate for a system).It is important that, unlike empirical functionals, SCAN is not fitted to any bonded system.Then the SCAN predictions for bonded systems are genuine predictions and not fits.
SCAN is constructed as an interpolation/extrapolation on the dimensionless meta-GGA ingredient where τ τ τ = + ↑ ↓ is the exact positive orbital kinetic energy density, τ W = |Ñn| 2 /(8n) is the von Weizsäcker kinetic energy density (exact for one-and spin-unpolarized twoelectron ground states), and τ unif is the kinetic energy density of a uniform electron gas.We have argued [21] that α is the right meta-GGA ingredient because it can recognize and give an appropriate description to three kinds of bonds: covalent single (α = 0), metallic (α ≈ 1), and weak (α >> 1).Since many of the exact constraints are for α = 0 or α ≈ 1, we need an interpolation/ extrapolation guided by the remaining constraints and the appropriate norms.
We have tested SCAN and other semilocal functionals on the G3 set (223 atomization energies of molecules), the BH76 set (76 barrier heights to chemical reactions), the S22 set (22 binding energies for weakly-bound molecular complexes, with hydrogen bonds or van der Waals bonds), and the LC20 set (20 lattice constants of solids).Table 2 shows the MAE's in kcal/mol (first three sets) or Å (fourth set).LSDA is the local spin-density approximation, PBE [19] is a GGA, while TPSS [14] is our earlier nonempirical meta-GGA and SCAN [2] is our current nonempirical meta-GGA.pSCAN (proto-SCAN) is identical to SCAN except in the way it interpolates the α = 0 correlation energy between fully spin-unpolarized and fully spin-polarized limits.pSCAN is in fact the version of SCAN that we first submitted for publication, while SCAN is the version finally accepted for publication.A first observation is that most calculated properties (the atomization energies of G3, the barrier heights of BH76, and the lattice constants of LC20) improve as we climb from the first rung (LSDA) to the second (GGA PBE) and then to the third (meta-GGA's TPSS, pSCAN, and SCAN).An exception is the weak-interaction binding energies of S22, which actually worsen from LSDA to PBE to TPSS but then improve dramatically in pSCAN and SCAN.The improvement in the pSCAN and SCAN lattice constants is also dramatic.The added constraints and norms beyond those in TPSS are doing what they should do.
A second observation, and one more relevant to the conclusions of this article, is that proto-SCAN and SCAN are identical for S22 and LC20 (where all systems are fully spin-unpolarized).They make very similar errors for BH76, even though the transition states of chemical reaction are typically partly spin-polarized.But SCAN is much better than proto-SCAN for the G3 atomization energies, where the open-shell free atoms are typically partly spin-polarized.
The difference between proto-SCAN and SCAN arises only from the way in which we interpolate the α = 0 correlation energy between fixed fully spin-unpolarized and fully spinpolarized limits.Let us define the relative spin polarization The function that exactly interpolates the exchange energy per particle of a uniform electron gas between its ς = 0 and | ς | = 1 limits is , which varies from 1 at ς = 0 to 2 1/3 = 1.26 at | ς | = 1.For α = 0 correlation, proto-SCAN uses the interpolation , which varies from 1 at ς = 0 to 0 at | ς | = 1.The latter limit is needed to make the correlation energy per particle vanish for all one-electron (α = 0 and | ς | = 1) densities.But SCAN uses a different interpolation between the same two limits, , designed to make the low-density limit of the SCAN exchangecorrelation energy nearly independent of ς over the range 0 ≤ | ς | ≤ 0.7, as in TPSS.In this way, SCAN and TPSS satisfy, as best they can, the exact constraint that the exchange-correlation energy per particle in the low-density limit should be independent of ς.Thus both SCAN and TPSS yield accurate atomization energies, while proto-SCAN is less accurate.
The two interpolation functions are compared in Table 3.For a given | ς | between 0 and 1, G G c c p > , so SCAN provides a little more negative correlation energy in a partly spin-polarized atom than proto-SCAN does.This slightly lowers the energy of the atom and reduces the too-high atomization energies found in proto-SCAN.While some exact constraints can only be satisfied on a sufficiently high rung of Jacob's Ladder, others (including the absence of spin-polarization dependence in the low-density limit) are already satisfied on the first or LSDA rung and get harder to satisfy on higher rungs.
The effect on the atomization energies of stricter enforcement of the low-density-limit constraint is rather subtle.It is magnified in the G3 set by the presence of numerous large molecules that contain many C atoms that "compound the error".While the cohesive energy of a solid is an intensive (per atom) quantity, the atomization energy of a molecule is an extensive one; it and its density functional errors grow with the size of the molecule.

Conclusions
The atomization energy errors of a density functional can strongly magnify minor errors in its spin-polarization dependence which are of little importance for other properties of many molecules and solids.Thus density functionals should not be judged primarily by their atomization energy errors, but by a wider spectrum of tests.Although SCAN atomization energies are not significantly better than TPSS atomization energies, we believe that SCAN is significantly better than TPSS.We already have evidence [2] that SCAN can reach a new level of accuracy for the energy differences between molecules or between solids at fixed atomic composition.In future work, we will test SCAN and other functionals for reaction energies, isomerization energies, structural energy differences in solids, and formation energies (and thus for relative stabilities).

Table 1
Energy gap (eV) between the lowest-unoccupied and highest-occupied atomic orbital energies, from SCAN

Table 2
[2]n absolute errors of semilocal spin-density functionals for three molecular (G3, BH76, S22, in kcal/mol) and one solid-state (LC20, in Å) test sets, as described in the text.References for the test sets are given in Ref.[2].

Table 3
Functions that interpolate between the fully spinunpolarized and fully spin-polarized limits for exchange, and for α = 0 correlation in proto-SCAN and SCAN spin-density functionals for any system.