Locality of correlation in density functional theory Free

Kohn-Sham (KS) density functional theory¹ (DFT) enjoys remarkable popularity, being used in more than 30 000 papers last year.² Essentially all these calculations use some approximation to the exchange-correlation energy, E_XC, as a functional of the (spin)-densities.³ Although many hundreds of such approximations exist and appear in standard codes, most calculations are run with one of a few approximations. These standard approximations^4–8 have been around for almost twenty years and their successes and failures are well-known. Of course there are many excellent improvements beyond these approximations for specific purposes.^9,10

Although the exact theory of DFT is well-established^11,12 and rarely questioned these days, it is sometimes claimed that there is no systematic approach to the overall construction of approximate functionals. In fact, this is not true. More than 40 years ago, Lieb and Simon (LS) proved that Thomas-Fermi (TF) theory becomes relatively exact for the total (non-relativistic) energy of any system in a limit in which both the particle number N becomes large and the coordinates are scaled. For atoms, this amounts to Z → ∞, where Z is the nuclear charge, while keeping the atom neutral (Z = N). But the limit also applies to all molecules and solids, once their bond lengths are simultaneously scaled.¹³ This limit then implies a systematic approach to the construction of density functionals, by finding the density functional corrections to the TF approximation that yield the exact leading corrections to the TF energy. Unfortunately, even in very simple situations, deriving and using these leading corrections can be very demanding, and they will not, in general, be explicit density functionals.^14–16 Nonetheless, such functionals are, in principle, well-defined, if not easily approximated.

The modern world uses the KS variant of DFT, in which only the XC energy needs approximating. Schwinger demonstrated^17,18 that, for atoms, the local density approximation (LDA) for exchange becomes relatively exact in the large-Z limit, and a proof (for non-singular potentials) was given by Conlon in the general case.¹⁹ A rigorous proof for atoms was given by Fefferman and Seco.²⁰ The next correction was estimated by Elliott and Burke²¹ and shown to be accurately reproduced by both the B88 functional⁴ and the exchange part of PBE.⁸ Thus, at the generalized gradient approximation (GGA) level, the most commonly used approximations (at least 70% of present calculations²) reproduce the leading correction to LDA in this limit. This insight was used (in an inverted piece of logic) to construct the exchange functional in PBEsol,²² which is useful for calculating the equilibrium properties of solids (but not their thermochemistry). The tension between accurate energetics and bond lengths^23,24 suggests that the GGA form is insufficient to capture all aspects of the leading correction to LDA accurately.²² The asymptotic expansion for exchange has been used to construct new functional approximations.^25–27 

The present paper explores the next logical step in this analysis: Does the correlation energy alone also become relatively exact in LDA in this limit? If so, what is the leading correction, and do standard approximations capture this? These are very fundamental questions about the nature of the XC functional and our ability to approximate it with semilocal functionals. Such questions concern properties of the exact functional that are highly relevant to approximations.

We cannot hope to answer this question in general at present, but we can try to answer this question in the one case where we have sufficient data: atoms. By very carefully analyzing accurate correlation energies²⁸ for a sequence of non-relativistic atoms up to Z = 86, we find that they can be fitted to the following form:

The first constant, A_C, is about 20.7 mhartree and is derived from the uniform electron gas and confirmed analytically in Ref. 29. The second constant, B_C, is a property of all atoms and is highly inaccurate in LDA. From theory and our calculations, we estimate B_C to be 37.2 mhartree, whereas in LDA, it is found analytically to be −4.5 mhartree. However, even this very limited information yields much insight into the nature of the correlation energy functional. At a very elementary level, we see why it has been so much more difficult to identify than the exchange form: the 1/|r − r′| behavior of the Coulomb potential produces a ZlnZ term, so that only at astronomical values of Z does the LDA correlation become relatively exact. The fact that $B C LDA$ is negative helps to explain the huge overestimate of LDA correlation energies. We also note that it is only the high-density limit of LDA that is relevant to atoms, not the entire density dependence of the uniform gas. Thus, within this analysis, the LDA is only particularly relevant for Coulomb matter in the extreme high-density limit. Of course, for a simple valence metal solid which is well-approximated by a uniform gas, the LDA is needed to produce an accurate correlation energy. Moreover, for energy differences and derivatives with respect to nuclear coordinates that are relevant to bond lengths, more terms in the LDA correlation energy might be relevant.

A first modern attempt to extract the behavior of the correlation energies for large non-relativistic atoms was made by Kunz and Rueedi.²⁹ They complemented analytical work based on a perturbative many-body approach with numerical values from Clementi and Corongiu³⁰ up to Z = 54. Their estimate for B_C is inaccurate, due to inaccuracies in the data and an insufficiently precise extrapolation.³¹ 

The layout of this article is a little unusual. Almost half of the material is introductory, explaining how the Lieb-Simon limit is relevant to approximating density functionals for any component of the energy: kinetic, exchange, or correlation. This is needed to understand the relevance of the results for correlation. Part of doing this is a simple illustration of how the large-N limit of a series can be used to accurately approximate its value for all values of N, including even N = 1 (or less). We also use this to give our definition of non-empirical parameters. We next give a substantial amount of background material which also serves to define our notation. The main body of the paper is the careful extrapolation, from the results we have for finite Z, to the large-Z limit. We test these extrapolations on cases where we know the answer and show that the results are independent of the details of our procedure. In Sec. IV, we discuss the relevance of our results for the density functional approximations.

We use atomic units throughout, so all energies are in Hartrees and all distances in Bohr radii. All calculations employ spin DFT and are entirely non-relativistic.

We beg the forbearance of the reader, as we close with a simple mathematical illustration of the nature and the usefulness of expansions as N → ∞, which is the type of expansion we will apply to density functionals. Later, we will tie this particular example to DFT, but for now, consider only the mathematics.

Consider an infinite series with terms,

i.e., a list of 1/n², but where each value is repeated n² times. If f_n is the nth term in this series, suppose we wish to find the partial sum,

This function is plotted in Fig. 1.

Clearly the sum is unbounded as N → ∞. To expand the series for large N, consider those values at which S is an integer, a condition analogous to filling an atomic shell,

where

Expanding j(M) in inverse powers of M^−1/3, we find

yielding

where b₀(N) = (3N)^1/3, b₁(N) = − 1/2, etc. We then define the sum up to a given order in N^1/3,

and compare these with the exact curve in Fig. 1. The most important point for now is that including only the first two terms yields a very reasonable approximation to S_N, i.e., the large-N expansion can be used to approximate the series for any value of N, all the way down to N = 1. Thus, even though we might only care about, e.g., N < 100, we can find good approximations by studying the approach to the large-N limit. We will return to this example several times in this article. For now, we also point out several features:

• The leading term, B₀(N) = (3N)^1/3, is always an overestimate, being 44% too high at N = 1, with the error decreasing with increasing N.

• Inclusion of the second term reduces the error to a 6% underestimate for N = 1, and now the sign of the error varies. This is almost always a better approximation than B₀(N), except near N = 0, where B₀ gives the exact answer, while B₁ = − 1/2. A physical realization of this expansion and interpretation of B₀ and B₁ is discussed in Sec. II D.

• The first two terms are smooth functions of N, but the exact partial sum S_N has kinks wherever S_N passes through an integer.

• The next term in the series, b₂(N), behaves as N^−1/3. It is not smooth but instead periodic across a shell; thus, it (and all subsequent terms) requires a more careful analysis than that for filled shells.³² We discuss b₂ further in Sec. III F.

The formal basis of DFT was established first by Hohenberg and Kohn,³³ but we use the more general formulation of Levy here.^11,12 The ground-state energy of any electronic system can be found from an exact variational principle in the density,

where v(r) is the one-body potential, the minimization is over all normalized non-negative densities with finite kinetic energy,¹² and

where $T ˆ$ is the kinetic energy operator, $V ˆ ee$ is the electron-electron repulsion, and the minimization is over all normalized antisymmetric wavefunctions.

The original DFT was that of Thomas-Fermi,^34,35 in which the kinetic energy of the electrons is approximated with a local approximation,

where A_S = 3(3π²)^2/3/10, and their interaction energy is approximated by the Hartree electrostatic self-energy of the density,

Applying such approximations to the total energy, and minimizing to find the density of a system, leads to self-consistent TF theory, which was used for several decades in material calculations.³⁶ The TF equation for an atom is now iconic,^34,35,37 and was solved numerically in the original papers. Energies are typically accurate to within about 10%, but the theory is too crude for modern chemical and material purposes.³⁶ 

The variant in use in almost all modern electronic structure calculations is Kohn-Sham (KS) DFT, which posits a fictitious set of non-interacting electrons with the same ground-state density.¹ These electrons satisfy the Pauli principle and obey the KS equations,

Writing the energy in terms of these orbitals,

where T_S is the kinetic energy of the KS orbitals, U their Hartree energy, and E_XC is defined by Eq. (14). For simplicity, we write all equations here in terms of DFT, but in practice, and in our calculations, we use spin-DFT.³⁸ 

Kohn and Sham also wrote down the most primitive approximation to use in their equations,¹ the local density approximation (LDA) for XC,

where $e XC unif ( n )$ is the XC energy density of a uniform gas of density n, which is now known very accurately.^39–41 This was used for a generation in solid-state physics,^42,43 but almost always significantly overbinds molecules, so it never became widespread in chemistry. However, most of the properties of molecules in LDA calculations (such as bond lengths and vibrational frequencies) are surprisingly accurate.⁴² 

At first, it was thought that LDA might be improved using the gradient expansion,⁴⁴ which includes the leading corrections to the LDA energy for a sufficiently slowly varying gas,

where s = |∇n|/(2k_Fn), k_F is the local Fermi wavevector, and g_XC(n) has been derived from many-body theory,^45,46 at least in the high-density limit. However, in many cases, these corrections worsened LDA results for atoms and molecules. This leads ultimately to the development of modern generalized gradient approximations,⁴⁷ which include more general functions of s (but not higher gradients) and led to the widespread use of KS-DFT in chemistry and materials today.⁴⁸ In this language, we can then understand that $T TF =T S LDA$⁠, i.e., the Thomas-Fermi theory approximates the kinetic energy as a local functional of the density for non-interacting electrons.

The standard way to think about XC is in terms of the strength of the Coulomb interaction.³ In DFT, we put a coupling constant λ in front of V_ee and imagine varying it, keeping the density fixed. Then, for finite systems,

where the correlation energy contains only λ² and higher powers of λ. Truncation at first order yields exact exchange (EXX), the DFT version of Hartree-Fock⁴⁹ (HF), which yields ground-state energies almost identical to those of HF.⁵⁰ 

One can relate this λ-dependence to coordinate scaling of the density,⁵¹ 

via

This can be used to write E_XC entirely in terms of a purely potential contribution, the adiabatic connection formula.^52–54 Furthermore, in the absence of degeneracies,

and

We begin with a definition of scaling to large Z. For any system of N non-relativistic electrons with one-body potential v(r), we define a ζ-scaled system following Lieb and Simon,¹³ 

where ζ is a positive real number. As ζ → ∞, this corresponds to simultaneously scaling the coordinates and increasing the particle number. For atoms or ions,

i.e., ζ-scaling is simply the same as changing the number of protons and electrons by the same fraction. For molecules and solids, all nuclear separations R also scale as R/ζ^1/3. The complementary density scaling is

This differs from the usual coordinate scaling of Eq. (18) as the particle number also changes. Analogous to that case, the ζ-scaled density is, in general, not the ground-state density of the ζ-scaled potential, except within the TF approximation.

Lieb and Simon⁵⁵ (LS) rigorously proved that, as ζ → ∞,

for all non-relativistic Coulombic systems. For example, the percentage error of the energy in a TF calculation for an atom vanishes for large Z. Moreover, any well-behaved integral over the density of the form¹² 

will also become relatively exact. Although n^TF(r) has many well-known deficiencies (divergence at the nucleus, non-exponential decay at large distances, missing quantum oscillations), integrals such as I over n^TF(r) have vanishing relative error in this limit.

Without taking too large a leap, we assume that the LS result is also true for non-interacting electrons in any potential,⁵⁶ yielding

again, a universal limit of all systems. Thus any approximation that fails to satisfy this limit is unlikely to be a useful starting point for higher accuracy approximations for a large variety of systems.

The simplest useful example of T_S is to consider non-interacting electrons in a hydrogen potential, −1/r. Their energies are −1/2n², where n is the principal quantum number. Via the virial theorem, T_S = − E, and accounting for double occupation of the orbitals, in fact,

where S_N is the partial sum from the Introduction. Thus our purely mathematical example is in fact a relevant example. The TF density for this problem is simply⁵⁷ 

where r_c = (18/Z)^1/3, and Θ is the Heaviside step function. This situation is similar to the interacting case, in which the TF density is singular at the origin, missing oscillations, and does not decay correctly at large r. Insertion of Eq. (29) into Eq. (11) yields precisely (3N/2)^1/3/2, the leading term B₀(N/2)/2 in the expansion for large N. This is consistent with Eq. (27). Thus, this reasoning says that, for an orbital-free theory to take advantage of the LS result, it should recover both this value and the next correction⁵⁸ (−1/2) in the large-Z limit. Of course, this is unlikely to be sufficient to determine a usefully accurate orbital-free DFT. But if the reasoning commonly applied to XC is applied to this problem, i.e., that approximate functionals should build in those conditions that they can satisfy, it is a necessary exact condition. No approximate kinetic energy functional known at present produces the correct value for the coefficients of this expansion for all electronic systems.³⁷ 

This section addresses a fundamental question about DFT in chemistry and materials science, which is as follows: How can such simple approximations, such as LDA, GGA, or even hybrids, possibly yield usefully accurate results for such a demanding many-body problem? While their failures for many specific properties and systems are myriad, their successes are legion. How can this be, since everyone knows how demanding and complicated the many-body problem is?

The seeds of the answer are in the Lieb-Simon theorem. As Z grows, local approximations become more accurate for the energy and other integrated quantities (but not pointwise for densities, energy-densities, or potentials). The many-body problem becomes simple in two extremes: N = 1 and N → ∞. All the detailed behaviors of the wavefunctions are averaged away, as quantum oscillations become more and more rapid in space, having less and less effect on system averages.

This then provides a systematic approach to constructing DFT approximations, but one which is entirely different from the traditional approach in terms of many-body theory and expansions in λ, the electron-electron interaction. The Thomas-Fermi approximation becomes relatively exact for large N. Note that we do not claim that our systems of interest have large N, and indeed, a direct large-N expansion of the energy is often not sufficiently accurate. Instead, we apply the large-N DFT approximation self-consistently to the problem at hand. We are using the density functional that is exact at large N, which is much more accurate than simply taking the large-N behavior of our specific system. In terms of our hydrogenic illustration, a GGA for T_S can be constructed that automatically recovers the first two terms in the expansion of the energy, but that is more accurate for small N than the asymptotic expansion alone.³⁷ 

We now conjecture that the following statement is true for KS-DFT:

We call this a conjecture, because the Lieb-Simon theorem for TF theory has been proven rigorously in a mathematical physics sense,⁵⁹ but no such general demonstration yet exists for E_XC. However, all evidence suggests that this conjecture is correct. For example, few would doubt that, in a limit in which the number of particles is growing, E_C/E_X → 0, so that the statement needs only be proven for E_X. This was demonstrated very carefully by Schwinger^{17,18,60–66} for atoms and detailed in the pioneering book of Englert.³² It was later proven rigorously for atoms.²⁰ It can be proven for all systems if the singularity at the nuclei is smoothed,¹⁹ but no general mathematical proof has been given for Coulomb interactions. For our purposes, this is sufficient.

Virtually all modern XC approximations reduce to $E X LDA$ for uniform densities, which means that they also do so for inhomogeneous systems as ζ → ∞. Assuming the conjecture is correct, this makes sense as $E X LDA$ is then a universal limit for all systems, and its relevance has nothing to do with the rapidity of the density variation (which is never slow for realistic finite systems). By studying the exchange energies of noble gas atoms of increasing Z, one can²¹ deduce $E X LDA [ n ]$⁠, i.e., it is a simple limit of all systems, one in which the many-body problem becomes very difficult (the number of particles is diverging), but in which the density functional for the integrated quantity simplifies enormously.

As ζ → ∞, the coupling between electrons becomes weak, just as it happens for γ → ∞, in coordinate scaling. Thus this limit is weakly correlated. But the number of particles is also growing, and the number of pair-interactions is growing even faster. In this limit,⁶⁷ 

If we compare this with the usual coordinate scaling, and speak in terms of the coupling constant, as ζ grows, the effective coupling constant is shrinking. But due to the increased number of interactions, the net effect is as given above. The analog to Eq. (20) using the Lieb-Simon theorem is

while for XC, the analog to Eq. (21) is

We can then argue that, if standard DFT approximations work reasonably well for weakly correlated systems, it is because the density-dependence of the XC energy, found self-consistently in the KS scheme, is well-approximated by the limiting density functional. Most molecular systems at or near equilibrium bond lengths are weakly correlated, so LDA yields usefully accurate values in this region. While LDA is insufficiently accurate for the energy difference between an equilibrium bond and the dissociated atoms, its error is extremely systematic, and the next correction to the energy in the form of a GGA works reasonably well.

This is entirely analogous to the starting approximation for most many-body treatments. Hartree-Fock becomes relatively exact as the electron-electron repulsion becomes weak. But, in practice, we apply it directly to any many-electron system with the full electron-electron repulsion and find usefully accurate descriptions (although not bond-energetics) from the self-consistent solution of those equations. LDA becomes relatively exact in a different limit, as discussed above. In this sense, LDA is as non-empirical as HF and is substantially more accurate for bond energies.

Here we claim that a locality principle is at work in DFT: Many of the successes and failures of DFT approximations to XC can be understood in terms of the expansion of the functional around the large-Z, i.e., local limit. The success of LDA for real materials (not slowly varying gases) has nothing to do with the uniform gas per se. It is a universal limit of all quantum systems, and that limit is most easily calculated from the uniform gas. The real question is as follows: For real molecules and materials, does LDA dominate, and if so, how large are the corrections? LDA works as well (or as badly) as it does for weakly correlated systems because, for those systems, the density-dependence of the XC energy functional is moderately accurately approximated by its limiting form (i.e., a local one). For such systems, and only such systems, inclusion of the next correction should improve its accuracy. Thus most molecules and many materials at equilibrium are weakly correlated. LDA works reasonably well, and standard GGA’s usually improve energetics substantially. On the other hand, static correlation is a signal that this expansion is failing, and hybrid functionals^6,68 are a crude attempt to account for this.

We therefore define the beyond-local (BL) XC energy functional as

and an obvious question arises: If LDA becomes relatively exact for E_XC at large ζ, what approximation for $E XC BL [ n ]$ becomes relatively exact in this limit? This is the natural expansion of which LDA is the leading order, and one can hope that accurate inclusion of the next order will lead to accurate energetics for weakly correlated systems. Presumably, standard GGA’s are crude examples of such corrections.

The answers for exchange are simple and well-established.²¹ In this case, for atoms,

The dominant term is exact within LDA, given by

To evaluate this, we write the TF atomic density in dimensionless form,³² 

where x = Z^1/3r/a and a = (3π/4)^2/3/2. Inserting this into Eq. (36) yields³⁷ 

where M is known numerically to be about 0.615 434 679.³⁷ Thus A_X ≈ 0.220 827 4. The precise value of B_X has not been derived, but estimates for B_X have been extracted numerically by careful extrapolation to the large Z limit of atoms. $B X LDA ≈0$⁠, but B_X ≈ 0.224, i.e., there is a large beyond-local contribution. The two most commonly used^4,8 generalized gradient approximations for E_X recover this beyond-local contribution rather accurately, whereas the gradient expansion for the slowly varying gas, when applied as an approximation to atoms, is smaller by about a factor of 2.^22,67

At this point, we make a small aside about what we mean by the term non-empirical. For simplicity, return to our mathematical example. Suppose we argue (or notice from examples) that the leading term must grow as N^1/3, and we use some set of tabulated values of S_N to fit the coefficient in B₀(N). If we used just the point N = 1, we would get this coefficient quite wrong (1 instead of 3^1/3 ≈ 1.4422). If we use the first 10 points, we get a more accurate answer. But if, instead, we use tabulated values of S_N to extrapolate the N → ∞ value, we get our best estimate of the non-empirical value. For example, knowing S_N for N = 1–100, and the form AN^1/3 + B, we can estimate A = 1.4458 (versus the exact value of 3^1/3 = 1.4422) and B = 0.525 (versus 1/2 = 0.5), by using only the last 50 values.

By contrast, simple fitting (such as a polynomial fit) to a range of finite data, including small N but without accounting for the asymptotic analytic form, will usually produce a numerically more accurate value within the range of the data, but will fail as N grows larger, and might be disastrous when attempting to find the next correction. Note also that it is vital to know the correct form, to guarantee relative exactness as N grows, and to allow investigation of higher-order terms.

Thus, by using tabulated values and the correct form of the expansion, we can estimate a non-empirical coefficient, even though we might not have the ability or energy to derive its value analytically. The larger the values of N that we use, and the more terms in the form, the more accurate our estimate can be. We denote the result as non-empirical, even though it has not been derived analytically, and we do not get the exact analytic value, just a best estimate of it.

We are almost ready to begin our analysis of the locality of the correlation energy. All our data will come from non-relativistic atomic calculations. Our starting point is the relatively recent publication of total correlation energies of spherical atoms²⁸ up to Z = 86 and non-spherical ones⁶⁹ up to Z = 36. These add to the well-known benchmark set⁷⁰ up to Z = 18, although they are not quite as accurate. It is the existence of these results that make the following analysis possible. Their values appear (with many others) in Tables I and II and are plotted in Fig. 2, in the form of correlation energy per electron.

The behavior of correlation energies across rows of the periodic table, and how it changes as one goes down a column, will play an important role in our analysis, but the data set becomes sparse for large Z. We have therefore performed pure random-phase-approximation (RPA) calculations for all elements up to Z = 86, and these results are included in Tables I and II and Fig. 2. We find that these RPA results can be “asymptotically corrected” using a simple formula (68) described later in the manuscript. Using only noble gases as a fitting set, this asymptotically corrected RPA (acRPA) nearly matches QC trends for all atoms for which QC data are available. We thus use the acRPA to “fill in” gaps in the reference data set.

The adiabatic connection formula^52–54 combined with the fluctuation-dissipation theorem (ACFD) and the RPA is fast becoming a de facto standard^71–73 for calculations of nearly quantum chemical accuracy energy differences. Unfortunately, it is well known that RPA gives poor estimates of absolute correlation energies. As we will show later, this erroneous contribution can, at least in the atomic systems considered here, be cancelled by a simple fit, Eq. (68), depending on the number of electrons only. This thus renders corrected RPA suitable for accurate correlation energy calculations of atoms, on a par with the most accurate quantum chemical benchmarks.

RPA energies are calculated using the ACFD correlation energy formula

Here the pair-density $n 2 C λ$ is found from the fluctuation dissipation theorem via

where stars indicate spatial convolutions. The non-interacting response χ₀ takes as input the Kohn-Sham orbitals ϕ_i and energies ϵ_i obeying ${ − 1 2 ∇ 2 + v S − ϵ i } ϕ i =0$⁠, and KS Greens functions G(r, r′; ϵ_i − ω) obeying ${ − 1 2 ∇ 2 + v S − ϵ i + i ω } G ( r , r ′ ; ϵ i − i ω ) =−δ ( r − r ′ )$⁠. We work from spherical and spin symmetric groundstates which make all calculations essentially one-dimensional. All equations can thus be carried out using a large radial grid which helps to reduce numerical errors in (39). Errors are estimated to be well under 1.1 mhartree/electron. We use the same algorithms and code as previous work.^74,75

Equation (39) is an indirect orbital function which formally maps a Kohn-Sham potential v_S and orbital occupation factors f_i to a correlation energy. We must thus start from a reasonably accurate potential if we are to expect accurate energies. For this work, we perform RPA calculations using strictly spherical Kohn-Sham potentials calculated using LEXX⁷⁴ theory extended to d and f shells. LEXX yields a Hartree and exchange energy functional of the same form as HF theory, but in which all orbitals are systematically calculated in a multiplicative spherically symmetric potential. Its associated correlation energy is thus slightly larger in magnitude. For most elements, the occupation factors are assigned according to Hund’s rules, in accordance with the theoretical and experimental evidence. For the transition metals, we use the lowest energy orbital filling (in exact exchange theory) of the s and d shells (e.g., 4s¹3d¹⁰ for Cu and 5s¹4d⁴ for Nb), which allow ready comparison with the results of McCarthy and Thakkar.⁶⁹ For the atoms in row 6, we simply fill the orbitals according to Hund’s rules to avoid the issues of (near-)degeneracy.

Ideally the orbitals would be calculated on the exact Kohn-Sham potential v_s(r), but this is available for only a very limited number of atoms, forcing us to use an appropriate approximation to the potential. The LEXX potential was chosen due to its good asymptotic form and inclusion of static correlation in open shell systems. To ensure that our RPA energies are appropriate, we performed some additional tests to study dependence on the potential. First, for C and F we evaluated RPA energies using both the exact⁷⁶ potential and the LEXX potential and found the difference to be less than 0.3 mhartree/electron. We also compared energies calculated using the PBE potential and LEXX potential across many species and found a maximum difference of 1.0 mhartree/electron, a likely worst case scenario. We thus assume our total error (numerical and error from the potential) to be under 1.5 mhartree/electron or 3%, whichever is the larger. This gives us, through RPA and the corrections discussed later, accurate benchmarks for all non-relativistic atoms with up to 86 electrons.

At the opposite end of the scale, one can perform atomic DFT calculations of correlation energies rather simply and go to much larger Z than with correlated wavefunction treatments. For the purposes of our study, we include all Z up to 86 and use both LDA and PBE results. We use the PW92 parametrization of the uniform gas,⁴¹ the relevance of which will be discussed below. We include PBE because, as shown below, its form mimics that of the exact functional in the large Z limit, so that the large-Z behavior can be most easily extracted by comparison with that approximation.

Tables I and II list all our results. In Fig. 2, we plot correlation energies per electron as a function of Z, as well as the RPA and LDA results. The shapes are well-known and unsurprising. The RPA significantly overcorrelates the atoms, and the LDA is even worse. Many quantum chemists question the value of LDA for such calculations, as the relevance of the uniform gas for such systems is far from obvious. Indeed, from this figure alone, it is unclear that the QC correlation energy per electron is diverging logarithmically — the limiting behavior of Eq. (1) derived from the uniform gas. One could reasonably argue that it is approaching a constant for sufficiently large Z. We demonstrate below that the latter is definitely not the case.

Note that the difference between HF energies^69,77 and Kohn-Sham exact-exchange (EXX) energies (as computed in our work) is smaller than other errors, so we ignore this difference here.

We now test the locality conjecture for correlation and use it to extract trends in Fig. 2 and Table I. The conjecture implies that

Thus we begin with the local density approximation, which uses the correlation energy of a uniform gas,

where $ϵ C unif ( n )$ is the correlation energy per electron. This is a non-trivial function of the density that is now well-known from many-body theory and quantum Monte Carlo simulations of the uniform gas.³⁹ There are several almost identical parametrizations^40,41 in common use. Our interest is the large-ζ limit, which is dominated by high densities. In a landmark of electronic structure theory, Gell-Mann and Brueckner⁷⁸ applied the random phase approximation (RPA) to the uniform gas to find

where r_S = (3/(4πn))^1/3 is the Wigner-Seitz radius of density n, c₀ = 0.031 091, and c₁^RPA = 0.070 82. In fact, Eq. (45) yields the exact high-density limit if c₁ = 0.046 64, the correction from RPA being due to the second-order exchange.⁷⁹ But, as is clear in Fig. 2, LDA¹ greatly overestimates the magnitude of the correlation energy of atoms (factor of 2 or more). For atoms with large Z, insert n^TF(r) into Eq. (44) to find

where A_C = 2c₀/3 = 0.020 73. $B C LDA$ is found by considering the contributions linear in Z generated by replacing Eq. (37) in the high-density limit of $E C LDA [ n ]$⁠. As a result,

where

and is known numerically to be about −3.331 462.³⁷ This yields

as reported in Ref. 67.

The locality conjecture, applied to correlation alone, implies that the leading term of Eq. (46) is exact. In fact, this has been proven for atoms relatively recently,²⁹ which shows definitively that the curves in Fig. 2 do not saturate. It also suggests that we define

as the asymptotically corrected correlation energy per electron. This is plotted in Fig. 3. Now the curves do appear to saturate, although the oscillations across open shells make accurate extrapolation of the large-Z limit very difficult. For any such curve, we define

and try to estimate this value as accurately as possible. This will be the subject of Sec. III D.

But to get a quick idea, using techniques that have worked for exchange,^21,37 we largely eliminate the effect of shell structure by taking only noble gas atoms and plot as a function of 1/n_HOMO, the principal quantum number of the highest occupied shell, as in Fig. 4. Clearly, the RPA and LDA values of B_C are quite different from the reference QC value. (The PBE correlation functional⁸ is included, as its asymptotic behavior is almost perfectly parallel to the benchmark data, and we use this in Sec. III D to improve our extrapolation.) Repeating the procedure for alkali earths shows that almost identical results occur, as shown in Fig. 5. Fitting noble gas data to the form

we find $B C nob =0.0371$ and $C C nob =−0.0536$⁠. Repeating the fit for the alkali earths, we find $B C AE =0.0378$ and $C C AE =−0.0644$⁠, reflecting a consistency in the value of B_C extracted from either series.

Finally we can convert our numerical extrapolation versus n_HOMO into a smooth function $e AC nob ( Z )$ in order to explore trends across the entire periodic table. To do so, first, notice that the data oscillate slightly around each fit, with every other point being above or below. This reflects the double-step in the periodic table between the appearance of new rows. To account for this, while still analyzing atoms down a given column, we define

which smoothly interpolates n_HOMO between even and odd noble gas atoms and is asymptotically correct as Z → ∞.⁸⁰ We can now write

an estimate of e_AC for all Z. To repeat the process for the alkali earths, we use $n ̃ nob ( Z − 2 ) +1$ instead of $n ̃ nob ( Z )$ in Eq. (54), and the AE coefficients $B C AE$ and $C C AE$⁠. The two asymptotic energies converge rapidly onto each other for large-Z, as shown in Fig. 8.

Having identified that B_C appears to approach a finite value, in Sec. III D, we make a variety of constructions in order to estimate its value as accurately and reliably as possible from the available data. We also use some of the same techniques to fill in large-Z non-spherical atomic correlation energies in Tables I and II.

As an obvious test of our numerical extrapolation methods, we consider self-consistent LDA calculations for neutral atoms. Because such calculations are of irrelevant computational cost, we enlarge the data set up to Z = 460. This is the 12th row of the non-relativistic periodic table, assuming Madelung’s rule continues indefinitely. Note that any calculations beyond Z = 100 are simply a mathematical device to reduce the range of Z over which we extrapolate.

As we go to larger Z, we find even/odd oscillations becoming less problematic; they can be largely eliminated by averaging over the noble gas and the alkali earth atom in each row. These are both spherical (and so appear in our original QC data set) and not spin polarized. In Fig. 6, we plot the averaged results of self-consistent LDA calculations (squares connected by dotted line), the linear extrapolation of these to the asymptotic limit (solid line), and mark the analytical form of this limit from Eq. (47). The LDA energy evaluated on the TF density for each value of Z (dashed line) is also shown.

One can see that the naïve linear trend of e_AC with 1/n_HOMO must eventually fail for 1/n_HOMO < 0.1. This occurs even when the TF density is used, due to the logarithmic terms in the high-density expansion of $ϵ C unif$⁠. The difference between the self-consistent and TF curves also appears to be non-linear in this region. To understand this, we expand the uniform gas energy in the high-density limit to the next two orders beyond Eq. (45),

where c₂ = 0.0066 and c₃ = 0.0104 are known from perturbation theory.⁴¹ Inserting the TF density yields

where

and

In these equations, r_s = a(3/f(x))^1/3 is the local Wigner-Seitz radius of the TF density for Z = 1, and

The important feature is that this expansion is ill-behaved. The culprit is the nature of the average over TF density [Eq. (59)], as x²f(x)[r_s(x)]ⁿ decays only as x²ⁿ⁻⁴ as x → ∞. As a consequence, these next two terms, of order n = 1, have much larger coefficients than the leading, order n = 0 terms of Eq. (45), and the following two terms diverge on the TF density. This poor behavior is inherited from the uniform gas and shows how unsuitable all but the leading terms of the LDA are for approximating the correlation energy of atoms. Furthermore, it causes the structure we see in the figure as Z → ∞. Thus, the LDA curve is especially difficult to extrapolate accurately. Nonetheless, performing a naive application of the linear extrapolation for n_HOMO ≥ 3 yields −3.0 mhartree, an underestimate of 1.5 mhartree from the analytic value. Thus we expect no larger error from our extrapolation process.

Finally, we perform the extrapolation that produces our best estimate for B_C. To achieve the maximum accuracy, we compare the accurate QC data with PBE results.⁸ The primary reason is that the PBE functional was constructed to yield an accurate finite result in this limit, by correcting the LDA energy. Our use here is simply that we can extrapolate the difference between QC and PBE, while extracting the PBE result analytically. The parallel trends in n_HOMO for the PBE and QC in Figs. 4 and 5 imply that their difference tends to a constant. Thus the resulting fit ought to be almost exclusively a measure of B_C, with minimum influence from the details of the fit.

The PBE beyond-local correlation energy for spin-unpolarized densities is⁸ 

where t = |∇n|/(2k_sn) is the dimensionless gradient for correlation, k_S = 2(3n/π)^1/6 is the TF screening length, and

where

and β = 0.066 725. This peculiar combination of functions results from requiring it to reduce to LDA for uniform densities, yield a finite value in the high-density limit of finite systems, and, when combined with PBE exchange, recover the LDA linear response for a uniform electron gas.⁸ It is the first two conditions which ensure that it yields a finite value for B_C that differs from LDA.

$B C PBE$ is found by considering the contributions linear in Z generated by replacing Eq. (37) in the high-density large-Z limit, where $A ̃ →1/r S$⁠, y → 0, and f_C → 1. As a result,

where

is a well-defined integral over the TF density.³⁷ For u = β/c₀ = 2.146 119, we find I₃ = 1.411 164, yielding

in agreement with Ref. 67.

In Fig. 7, we plot the energy difference in mhartree between QC and PBE per electron. Because the difference is so small, on this scale, the differences between noble gases and alkali earths are very visible, especially for odd rows. But their mean is seen to converge nicely. A fit yields

Combining the analytic result of Eq. (66) with this one, we conclude that the asymptotic value of B_C is 37.2 mhartree. Table III compares the different estimates, showing their consistency, and suggesting the extrapolation error is no more than 1.5 mhartree.

There is considerable current interest in using RPA and beyond-RPA methods for quantum chemistry.⁸¹ While many good features are known, a major drawback to pure RPA calculations for thermochemistry is the inaccuracy of RPA atomic energies and its impact on atomization energies.

Unfortunately, the quantum chemical data become sparse for large Z, being limited to spherical atoms. However, our RPA calculations are less limited, and their error relative to the quantum chemical data set is relatively simple. We fit the differences between the RPA correlation energies and the QC ones for noble-gas atoms using

where $n ̃ nob ( Z )$ is given by Eq. (53). We denote by acRPA the RPA energies with the correction of Eq. (68) added to them. The results are presented in Tables I and II.

In Fig. 8, we plot both the accurate correlation energies and the acRPA values, and the residual error. Clearly, the fit is extremely accurate and errors become negligible for large enough Z. Finally, we note that the constant contribution to the fit, 19.9 mhartree, is close to (but different from) the error the RPA makes for the correlation energy per electron in the high density limit of the uniform gas, 24.2 mhartree, as discussed in Section III B. We can understand one contribution to this difference by considering the PBE for RPA correlation energies developed by Kurth and Perdew.⁸² The expression of the corresponding enhancement factors is reported in Eq. (27) of Ref. 83, from which we readily obtain

the difference with PBE being again due to RPA correlation energy of the uniform gas. Since the PBE construction uses the value of the gradient expansion coefficient β from the high-density limit in the RPA approximation, the gradient contributions to $B C PBE$ and $B C RPA$ are the same. Assuming this procedure is highly accurate for approximating RPA, the remaining difference is only 3 mhartree—well within the numerical error bars for RPA correlation energies, estimated to be up to 4.6 mhartree (3% of 84 mhartree) for Z = 86.⁸⁴ 

Note that this correction formula concerns the relation between an approximation to the many-electron problem and the exact result and is entirely independent of any density functional analysis. Its value is that we can create an enhanced data-set of the QC data, plus RPA-corrected estimates for the larger open-shell atoms. We use acRPA for benchmark values wherever QC values are unavailable, i.e., the non-spherical atoms of Table II. Our simple atomic correction, Eq. (68), is remarkably accurate, as shown in Tables I and II, with few errors greater than 1 mhartree. A scheme which reproduces this correction, but also works well for molecules at equilibrium bond lengths, would overcome the limitations of pure RPA. We suggest testing existing schemes for their behavior in the large-Z limit, characterizing them by the value of B_C.

To motivate this final section, we return to the didactic example from the Introduction. In Fig. 9, we show (blue dotted line) the error in the asymptotic curve B₁(N) derived in that section and plotted in Fig. 1. Subtracting the mean value of b₂(N) across a shell, [ − 12(3N)^1/3]⁻¹, yields the error on $B 1 ∗$ shown as a blue solid line. This error is not a smooth function of N but has cusps at shell boundaries, is periodic in ν, the fraction of the shell that is filled, and only slowly decaying with shell number. Careful analysis³² yields

and addition of this term to form B₂(N) yields the error shown as red dashes. The errors between shells are reduced by an order of magnitude relative to B₁(N), and the error for N = 1 is a mere 3 × 10⁻⁵, i.e., 30 μH. Thus it behooves us to look at the structure of such terms.

Of course, exchange and correlation of real atoms is much more complicated, and we do not have access to the exact asymptotic expansion for these quantities. But as our example shows, we may numerically illustrate the nature of the underlying higher-order asymptotic corrections simply by showing the deviation of our data from smooth behavior. To do so, we first obtain smooth curves from our analysis of QC data from single columns of the periodic table, equivalent to the leading order terms of our didactic example. Second, we plot the difference between the smooth asymptotic curve and the QC data augmented by acRPA estimates to obtain a sense of the higher-order terms that determine shell structure along rows of the periodic table. We have already constructed such curves for correlation: $e AC nob ( Z )$ fitting values for the noble gases using Eq. (54) and an analogous curve, $e AC AE ( Z )$⁠, for the alkali earths. As shown in Fig. 8, the two curves have the same behavior at large Z, but the AE curve is more reliable for lower Z, so we adopt the latter as standard.

To create smooth functions for exchange, we fit exchange energy results with continuous functions of Z and apply these functions for every value of Z. We standardize on alkali earths to be consistent with correlation. As we know LDA becomes relatively exact for exchange for large Z, we define

and fit EXX energies for alkali earths with

where x = 1 − Z^−1/3. This gives B_X = − 0.260 in agreement with a previous estimate of −0.24.²¹ The value of −0.08 at Z = 1 is somewhat off the actual value of −0.11, but gives a close fit to e_AX for all Z > 2.

In Fig. 10, we plot the deviations of the exchange and correlation energies per electron from the smooth fits. This is where acRPA is important: To see the periodic pattern in the deviations in the correlation energy by including the large Z values. We have also included vertical lines at the noble gas atoms, to show the periodicity.

There are many interesting features in this plot.

• There is a clear (anti-)correlation between the correlation deviation and the exchange deviation, showing tremendous cancellation of “errors” in the separate components of the deviations. Thus, adding Eqs. (72) and $e AC AE ( Z )$ yields a much more accurate estimate of XC than either separately. Shell structure is much greater in either exchange or correlation alone than in XC.

• The greatest deviations occur about the middle of the odd rows, and deviations are almost negligible throughout the 4th row.

• Filled angular momentum subshells tend to be marked by cusp-like behavior, especially for correlation. These are more apparent if the irregularities in filling for d and f subshells are ignored.

• We observe the pattern that whenever the most loosely bound electron also has the highest angular momentum, the deviation grows as that angular momentum subshell is filled and drops once the next subshell is filled. For example, for Z = 21–29, the 3d-shell is being filled and the deviation grows and then drops rapidly as the 4p’s are filled. Note that this only happens when a given subshell is filled for the first time; there is little effect when the 4d’s are being filled. A similar pattern appears for the 4f’s and (to a lesser extent) even for the 2p’s.

The majority of the numerical work in this manuscript has been devoted to the extraction of an accurate numerical estimate for B_C, which characterizes the leading correction to LDA correlation for atoms, and is strong evidence for the locality conjecture. In this section, we assume that conjecture is true and discuss its consequences for DFT approximations. We use the mathematical prolog and its relevance to the kinetic energy of non-interacting atoms to illustrate this, by discussing T_S[n], E_X[n], and E_C[n].

In the math prolog, we showed that, for the kinetic energy of non-interacting atoms,

For real atoms,

where A_S = 0.768 745 is given exactly by TF theory, while