The duhka of DFT: a noble path to better functionals via a point electron approximation for the exchange–correlation hole^†,‡

Keywords: density functional theory, electron correlation, exchange correlation energy, exchange correlation functional, exchange correlation hole, Jacob’s Ladder, pair correlation function.

When it comes to predicting chemical properties, density functional theory (DFT) in its Kohn–Sham (KS) form is today unarguably the workhorse of quantum chemistry. This is due to its attractive accuracy-to-computational cost ratio. In short, DFT works well, and it will not break your bank of computers. DFT has this remarkable status because contrary to the methods based on wavefunctions that depend on every electronic coordinate, DFT depends on the electron density ρ – a function of only three coordinates, regardless of the number of electrons. Today, conventional DFT methods are routinely applicable to systems with hundreds of atoms and even to systems with millions of atoms using linear-scaling and orbital-free DFT techniques. (For a comprehensive review of large scale DFT calculations, see ref.¹). We also note that DFT is the method of choice for generating massive amounts of chemical data needed for training machine learning models.^[2]

However, everyone knows DFT’s dirty little secret. There are, in fact, two of them. The first is that the rarely-employed Hohenberg–Kohn (HK) version of DFT – the only version that depends explicitly on the electron density – is in practice a dismal failure when it comes to prediction of chemical properties. (This failure of HK-DFT was in fact the motivating reason for the more successful KS-DFT theory.)^[3] The second, more often mentioned, is that DFT involves an unknown exchange–correlation (XC) functional which must be modelled somehow.^[4–14] This incompleteness of DFT poses a fundamental limitation on the applicability and intrinsic accuracy of the theory and has led to a proliferation in the number of developed DFT methods over the past three decades, resulting in hundreds of different models for the XC functional.^[8] From the point of view of the proverbial used-car salesperson then, DFT is an iconic example: it does not depend explicitly on its namesake, the electron density, and its contract is written without any guarantee of repair.^[15]

Into this density functional zoo, at the turn of the millennium, John Perdew shrewdly proposed his ‘Jacob’s Ladder of DFT’. In a loose analogy to John Pople’s two-dimensional diagram of systematically improved wavefunction methods,^[16] Jacob’s Ladder provides a route where, with suitable effort, one might toil upwards towards the nirvana of chemically accurate functionals.^[17]

The first rung of this ladder includes functionals that depend only on the value of the electron density through an energy density at a point in space, which itself depends only on the electron density at that same point, the so-called local density approximation (LDA). On the next rung is the generalised gradient approximation (GGA), which employs in addition the (reduced) gradient of the electron density at that point. By now you have the idea: the third rung, the meta-GGA (mGGA) functionals, employ in addition the second derivatives of the density at that point in the form of the Laplacian or the local kinetic energy density (the latter requiring not only the electron density but the KS orbitals themselves). This third rung is the last which can be regarded as ‘pure’ DFT, that is, DFT functionals that depend only on the electronic density and its gradients.

The fourth rung on the ladder involves admixing some component of the exact Hartree–Fock (HF) exchange. These functionals are therefore referred to as hybrid functionals. The fifth rung additionally involves admixing second-order Møller–Plesset perturbation theory (MP2) correlation energy. These functionals are referred to as double-hybrid DFT (DHDFT) functionals. The justification for these steps within the framework of DFT is the adiabatic connection formula of Harris and Jones.^[18] Notwithstanding, it is hard to say that a combination of DFT with post-Hartree–Fock wavefunction theory really constitutes a DFT method at all; as Becke puts it, ‘I think KS-DFT is about occupied orbitals only’.^[19]

Formally the computational cost of DHDFT methods is higher by 1–2 orders of magnitude relative to semi-local and hybrid DFT methods. The increase in computational cost of DHDFT relative to conventional DFT methods is further exacerbated by the need to use larger basis sets due to the slower basis set convergence of the MP2-like correlation term. Most chemical properties require at least a quadruple-zeta (ζ) quality basis set for obtaining results sufficiently close to the basis set limit.^[10,20] Thus, a major disadvantage of functionals from the fifth rung of Jacob’s Ladder is that they are applicable to substantially smaller systems relative to functionals from the lower rungs.

Another useful way to classify DFT functionals is as being nonempirical or empirical. Nonempirical methods attempt to model the XC functional by satisfying known exact conditions.^[21,22] Examples of popular nonempirical XC functionals are the GGA PBE,^[23] meta-GGAs, TPSS^[24] and SCAN,^[25] hybrid-GGAs PBE0^[26] and PBE0-1/3,^[27] and hybrid-meta-GGA TPSSh.^[28] Empirical methods, on the other hand, introduce parameters, which are determined by fitting to relevant thermodynamic, kinetic, and spectroscopic data obtained from experiment or more often from high-level ab initio wavefunction theory. Examples of popular highly empirical functionals with 10–40 adjustable parameters are the GGA HCTH407,^[29] meta-GGAs τ-HCTH^[30] and M06-L,^[31] hybrid GGA B97-1, and hybrid-meta-GGAs BMK^[32] and M06-2X.^[33] As expected, empirical DFT methods tend to outperform nonempirical methods for chemical systems and properties that are covered in the training sets.^[34–37] Nevertheless, this approach may lead to overfitting^[38–40] and it has been found that heavily parameterised functionals may exhibit erratic basis set dependencies and grid sensitivity.^[34,41–44] Thus, there is now a trend to develop empirical DFT functionals with a reduced number of parameters.^[34,45–47] Interestingly, the ‘hyperparametric disorder’ of DFT was foreseen by Gill in an Essay published in the present journal two decades ago ‘Obituary: Density Functional Theory (1927–1993)’.^[48] Still, we would caution, as von Barth does, regarding parameterisation, ‘There are certainly, theoreticians who would be prone to scorn the described [parameterisation] procedure but we find such attitudes unwise’.^[5]

Considering the empirical nature of many generally applicable DFT methods along with the limitations of Jacob’s Ladder outlined above, such as the jump in computational cost of methods on the fifth rung relative to the lower rungs and lack of guidelines for systematic improvements within each of the rungs, it seems to us that DFT is in a state of existential dukha.^[49] Is there a solution out of this fraught situation, a middle way?

The concept of the exchange–correlation hole plays an important role in DFT. It is essentially the depression in the probability of finding an electron (with coordinates ‘r₂’) around the position of a reference electron (labelled ‘r₁’).

where by definition

and where ρ₂ is the second-order reduced density matrix.^[3] From the above, we see that ρ(r₁)ρ_XC(r₁, r₂) is the difference between the exact electron pair distribution and the classical distribution comprising a product of the electron densities of the two electrons. In the early days, Gunnarson and Lundqvist explain the success of LDA, the first rung, by casting it as an approximation to the exchange–correlation hole with the correct sum rule.^[50] Importantly, the DFT exchange correlation energy E_XC can be obtained very easily from a knowledge of the XC hole

Thus, the functional of the exchange–correlation energy in terms of the exchange–correlation hole is known – a starting point not shared by the usual DFT functional of the electron density (where r₁₂ is the interelectronic distance). Gunnarson et al. later developed more refined models based on this hole approximation idea, which were much better than the LDA.^[51] Colle and Salvetti developed an approximate formula for the correlation energy for a special kind of wavefunction with the correct electron–electron cusp.^[52] This was later simplified by Lee, Yang, and Parr by avoiding the dependence on the on-top pair density to form the widely used LYP exchange correlation function of the density and its gradient.^[53] The hole-modelling idea has also been used by Becke and Roussel to make a very simple two-parameter exchange energy functional,^[54] by Becke and Johnson to develop elegant van-der-Waals density functionals,^[55] and by Baerends for orbital-dependent DFT functionals.^[56]

Obtaining the exchange correlation hole is challenging because it requires knowledge of the pair electron density, which, although obtainable from the wavefunction or quantum Monte Carlo calculations, is nevertheless computationally demanding, if not intractable for larger systems. There seem to be very few pictures of these in the literature, for example for H₂ and Si(s).^[57–59] From the existing images one learns that the exchange correlation hole is not spherical, and it extends onto nearby nuclei. This explains why a Taylor-series idea like the Jacob’s ladder scheme must necessarily be limited in its scope. We present here an idea to obtain a proxy for this quantity which is easier to calculate.

Suppose we have a program that can calculate a wavefunction for a given set of nuclei placed at some positions in space. The wavefunction may be as accurate as one desires HF, MP2, or CCSD (i.e. coupled cluster with single and double excitations). Consider now removing one electron from the system and placing it at a fixed position in space r_fixed, just like the other nuclei in the system. If the program makes use of atom-centered basis sets, then place a basis set onto this electron. Next, solve the Schrödinger equation for the wavefunction of the remaining electrons. Then clearly, the electron density ρ^N–1(r₁, r_fixed) from this wavefunction is a very good proxy for ρ₂(r₁, r_fixed) from which the exchange–correlation hole can be derived (Eqns 1, 2) (where N is the number of electrons in the original system).

Of course, this is not the exchange-correlation hole for the original system because one electron has been removed, but if the basis set on the electron is large, then one at least expects that the electron–electron cusp condition is well modelled – at least as well as the electron–nuclear cusp condition is on the remaining nuclei. One may collect as many of these pseudo exchange–correlation holes as one likes, simply by placing the fixed electron at distinct positions in space. From this, one can even calculate a pseudo exchange correlation energy via Eqn 3, integrating only over the coordinate of the first (unfixed) electron. Furthermore, placing the electron at the positions of a DFT integration grid would make it possible to integrate this quantity over all space, and calculate a pseudo electron–electron exchange correlation energy.

It is important to keep in mind that the pseudo exchange-correlation energy of the system with n–1 electrons calculated via a fixed electron as described above, will approach that of the real system with n electrons in the limit of a large number of electrons. Regardless of this limiting situation, we think that the pseudo exchange–correlation holes generated by this procedure form a simple, effective, and practical way to collect information pertaining to the real exchange correlation hole.

One idea immediately comes to mind. One can collect information concerning the exchange-correlation energy, the pointwise pseudo exchange-correlation energy from the procedure above, information about the electron density of the other electrons at the point of the reference electron (this is the ‘on-top’ pair density) as well as other quantities, such as the gradient of the electron density, the second derivative, the kinetic energy density, the asphericity of the exchange-correlation hole (parameterised by, say a spherical harmonic expansion), some notion of how the nuclei are arranged around the reference point, and the values of the electron density at the nuclear positions. A practical approach for processing these copious amounts of data would be through machine learning models. One may use this data to either improve DFT functionals using the standard Jacob’s ladder scheme (and then it will truly be at its limit) or one could model the actual exchange-correlation hole following the DFT pioneers.

All data used for this article is available through the referenced literature.

Amir Karton is an Associate Editor of the Australian Journal of Chemistry but was blinded from the peer-review process for this paper.

A. K. gratefully acknowledges an Australian Research Council (ARC) Future Fellowship (FT170100373). The authors thank Professor Jan Martin from the Weizmann Institute of Science for helpful discussions.