Improving results by improving densities: Density-corrected density functional theory

Density functional calculations have become ubiquitous in modern chemistry and materials science since the award of the 1998 Nobel prize in chemistry.[1] There are now many computer codes available for performing such calculations.[2, 3, 4, 5, 6, 7] It is a straightforward matter to choose a basis set and an approximate functional, and calculate an interesting property, such as a reaction barrier, bond length, or dipole moment. But it requires judgment and experience to choose wisely.[8] Ensuring the quantity is converged with respect to basis is relatively simple. Given hundreds of possible DFT approximations available in a code, the choice can be difficult.[9]

There are myriad approaches to constructing exchange-correlation (XC) approximations, varying from appeals to general principles of quantum mechanics to fits to large databases.[10, 11, 12, 13] Modern approximations include generalized gradient approximations (GGA), hybrids, range-separated functionals, the random phase approximation and variants thereof, dispersion corrections of at least three distinct flavors, double-hybrids, and many, many more.[11, 14, 15, 16] All over the world, theorists of many different backgrounds work at improving (or at least, expanding) on our current choices, either with improved accuracy, lower computational cost, or greater reliability.[17]

In each of the countless DFT calculations performed worldwide each year, the Kohn-Sham (KS) equations[18] are iterated to a self-consistent (SC) electronic density and orbitals, and the total energy of the system is reconstructed with these final quantities. By definition, this process finds the unique[19] density that minimizes the approximate energy. All components of that energy are exactly determined, apart from the notorious XC energy. It is that piece which is approximated in DFT and whose derivative appears in the KS equations as the XC potential.

Thus, whatever choice of XC is made, it is actually used twice in the calculation. Once in finding the density and a second time in finding the energy, so that neither is quite correct. As the foundation of DFT is to consider the energy as a functional of the density[19], we may write the error in any self-consistent KS calculation as:

where $n({\bf r})$ is the exact density and $E[n]$ is the exact functional, while tildes denote approximate quantities. In most practical calculations, modern XC approximations yield excellent approximate densities[20], so that the energy error would barely change if the approximation were evaluated on the exact density.

It is certainly extremely convenient to use the self-consistent solution density. It is easily computed from the KS equations. By being self-consistent, many important properties, especially those depending on derivatives of the energy, are much simpler and many additional terms need not be calculated. This is so convenient that essentially all modern codes use the self-consistent density in almost all circumstances. However, this was not always so. In the earliest days, the Hartree-Fock (HF) density was often used instead.[24, 25, 26, 27, 28] Mostly, it was used as a matter of convenience, so as to avoid needing to do a self-consistent calculation, in the belief that it mattered little. Later, the HF density was used as a matter of principle, to compare functionals against each other without having to worry about changes in the density.[26, 27, 28] It was even presciently noted that, in some cases, it really did seem to matter, and in those cases, it was often better to use HF densities.[29, 30, 31]

This article shows that, in fact, it really does matter, both theoretically and very practically. Until about 10 years ago, no careful, systematic analysis had been performed on this question. In fact, every single KS-DFT calculation ever run can be analyzed, to separate its functional error (energy error made on exact density) from its density-driven error (the remainder). Surprisingly large classes of calculation, such as typical reaction barriers, contain significant density-driven errors with standard functionals, such as B3LYP. One of the major reasons for this is the over-delocalization of charges and spins due to semilocal XC approximation.[32, 33] These errors are typically substantially reduced by using the HF density instead of the self-consistent density. Even highly accurate (and expensive) DFT approximations such as double-hybrids can be improved by separating out these two error sources in their design.

Figure 1 is a panoply of calculations where the density really matters. In every case, when self-consistent densities are replaced by HF densities, the energy errors drop by a substantial margin. In panel (a), we show an energy diagram for a textbook S_N2 reaction. Starting from either reactants or products, negatively charged complexes are formed barrierlessly, while the interconversion between the two involves a barrier of $\sim$35 kcal/mol in the backward and $\sim$5 kcal/mol in the forward directions. Standard DFT provides reasonable reaction energies, but fails badly for barrier heights in both directions for the complex interconversion. DFT underestimates the backward barrier height by about 10 kcal/mol, implying a reaction many orders of magnitude faster than reality. The barrier height is indeed smaller in the forward direction, but standard DFT yields no barrier at all. A long time ago it has been demonstrated that the use of HF densities fixes failures of DFT for barrier heights, and so does here.[31, 34] For the backward barrier, density-corrected DFT (DC-DFT) reduces the error of DFT by about 6 kcal/mol, whereas the DC-DFT forward barrier height matches the reference.

Panel (b) demonstrates the power of using DC-DFT to fix the failures of DFT for difficult torsional barrier heights, whose accurate predictions play a crucial role in describing a range of chemical processes (e.g., selectivity, protein folding, molecular electronics, etc.).[21] Most torsional barriers are very accurate with standard DFT (errors below 1 kJ/mol), but barriers of a single bond participating in $\pi$-conjugation are particularly problematic for DFT.[12, 35] For the oxalyl chloride shown, the standard DFT energy diagram is qualitatively wrong, incorrectly predicting that the perpendicular conformation is more stable ($\phi\sim 90^{\circ}$) than the trans conformation ($\phi=180^{\circ}$), where $\phi$ is the torsional angle depicted in panel (b). DFT also finds that there is no barrier upon conversion from trans to perpendicular. Using DC-DFT with HF densities, these barriers become far more accurate as shown.

For some weak interactions, such as halogen bonds, DC-DFT greatly improves over its self-consistent counterpart. The binding energy curve for one halogen bonded complex is shown in panel (c).[22] Standard DFT overbinds the complex by about 50 $\%$ at equilibrium, whereas the DC-DFT binding curve is almost indistinguishable from the reference. In contrast to DC-DFT, dispersion corrections (such as the commonly used D3[36]) cannot fix bad DFT densities, and their addition has almost no effect on the DFT binding curve in panel (c). Despite these large improvements in energetics when HF densities are used in place of self-consistent densities, their electrostatic potentials are almost identical, as shown. From the DC-DFT perspective, there is no need to stare at density or electrostatic potential plots to decide which density is better. DC-DFT measures the accuracy of densities directly in terms of their impact on the energy, the quantity that really matters. Even the tiny differences visible in the electrostatic potentials can be measured.

Finally, panel (d) reminds us of one of the earliest successes of DC-DFT – the description of odd-electron radical complexes, which play important roles in atmospheric and environmental chemistry, cell biology, etc.[23] In panel (d), we compare potential energy surface for the HO$\cdot$Cl^- complex by varying the $R$ distance and $\theta$ angle, as shown. Self-consistent DFT fails badly in simulating the potential energy surface: (i) it finds that the equilibrium structure is bent instead of linear ($\theta\sim 30^{\circ}$ instead of $0^{\circ}$), (ii) gives contour of the wrong shape leading to wrong forces, and (iii) gives too blue (too negative) potential energy surface. DC-DFT again saves its self-consistent counterpart by not only yielding the correct linear structure of HO$\cdot$Cl^- as the most stable, but also producing a far more accurate potential energy surface and equilibrium structure. In this way, we see that DC-DFT not only improves DFT energetics, but also gives more accurate geometries and force fields. For applications of the principles of DC-DFT to geometry optimization of any electronic structure method, see Ref. [37, 38], which contains many surprising results about functional performance for geometry.

For the cognoscenti, in Figure 1, all DFT calculations are with PBE, except in (a) which uses B3LYP, and all accurate reference calculations are DLPNO-CCSD(T)-F12, except in (d), which is simple CCSD(T). The rest of this article is about why the basic ideas of DFT do not imply always choosing the self-consistent density. This is followed by a discussion of practical DC-DFT, with many examples illustrating crucial aspects of density-sensitive systems and calculations. We next explore some of the finer points of theory, ending with a surprise: Although Diels-Alder reactions are not density sensitive, functionals (double-hybrids) designed to take advantage of DC-DFT perform better. We end with many challenges and potential of DC-DFT.

First we ask, why is this a question at all? Surely the self-consistent density is ‘best’ because it minimizes the (approximate) energy functional? It does, but because the functional is approximate, its minimum might well be below the true ground-state energy. Moreover, all useful properties are actually energy differences, and the difference between two minima does not obey a variational property. One of the most well-documented failings of most density functional approximations is that they are too smooth, especially as particle numbers pass through integer values.[39, 40] You might object that, in reality, all molecules have integer numbers of electrons. But as a bond is stretched, the exact functional develops sharp cusps that force integer numbers of electrons onto each fragment, while typical DFT approximations are smooth.[41] Figure 2 is a cartoon showing how function with a cusp can be well-approximated by a smooth one everywhere, but whose derivative is very wrong in the vicinity of that cusp. Semilocal XC approximations yield curves that are smooth everywhere, which causes the overdelocalization of charges when bonds are stretched.[42] The HF energy functional depends explicitly on only occupied orbitals, making it often even sharper than the exact functional. Relative to HF, correlation includes infinite sums over orbitals, which typically dampen the cusp as the particle number changes.

With some thought, these statements would appear paradoxical. If the functional is working well for the system you are calculating, how could its density be wrong? Well, this happens because its derivative, the XC potential, is sufficiently inaccurate as to produce a sufficiently flawed density as to mess up your energy evaluation. Return to Figure 2 to see a good approximation to a function whose derivative is lousy. Doesn’t a better functional automatically imply a better XC potential? No, it does not. Almost all modern XC approximations have very poor-looking XC potentials, often shifted by very large amounts relative to the exact XC potential.[47] Yet they still usually yield highly accurate densities in the regions where it matters. GGA approximations to XC often have worse looking potentials than their LDA counterparts, but nonetheless have much better energetics.[48, 49]

In the original work[20], the term ’abnormal’ was used to designate those KS-DFT calculations whose results were contaminated by density-driven errors, and this is a characteristic of the approximate functional, the property of interest, and the given system. By contaminated, one means that the error in the energy being calculated changes substantially if the exact density is used instead. There a small KS HOMO-LUMO gap in self-consistent DFT was identified as a signal of abnormality. But the use of the gap as the abnormality indicator is not ideal, as some calculations (e.g., those involving stretching of homonuclear bonds) have small HOMO-LUMO gaps without density-driven errors. More appropriate indicators of abnormality have been built and are detailed below.

In practice, much of the above is just so much theorizing, as, if we need to do a DFT calculation, we surely cannot afford to calculate the exact (or highly accurate) density. Fortunately, we show below that, in the cases where there is a significant density-driven error with a standard DFT calculation, very often using the HF density significantly reduces the density-driven error. This is presumably because HF, although yielding woefully inadequate energetics, suffers from the reverse of the errors of most density functional approximations. Essentially, DFT approximations almost always include some variety of semi-local density functional (i.e., depending on the density, its gradient, and/or its laplacian or kinetic energy density). These approximations tend to delocalize the density relative to the exact one, whereas the HF density is typically overlocalized. This is not to say that HF densities are somehow ’better’ than approximate self-consistent DFT densities. As discussed above, there is no well-defined meaning to being better. All this means is exactly what is stated: In cases where the self-consistent approximate DFT density is unusually poor, the HF density is often more accurate in the very precise and limited sense of yielding more accurate energetics.

But consider panel (b), where the results for forward barrier heights are plotted. Most functionals give about the same answer, except when the HF density is used. Now the differences between gray and purple are huge. Moreover, the answer often changes sign (i.e., goes from no barrier in SC-DFT to having a small one in HF-DFT). Also, the spread in the different self-consistent answers (shown by open circles) is now far greater than the variation in the purple bars. This is a pattern we often see: If a problem is density sensitive, often a standard bag of functionals which usually agree with one another show a wide disparity of results, when evaluated self-consistently. But on the HF density, their spread is smaller than usual.

Several years ago, a study was performed on pseudo octahedral Fe(II) complexes with various ligands.[65] All were wide spin gap cases, of order 1 eV (about 20 kcal/mol) but even so, different functionals yielded wildly different gaps. The spread in their results dropped by about a factor of 2 when HF densities were used instead. The average results differed significantly from those of CCSD or CCSD(T), but agreed (within error bars) with some very expensive, state-of-the-art QMC calculations. Since then, many authors have tackled these systems with many variations on many methods, so the jury is still out on whether or not DFT on HF densities yields accurate spin gaps here.[59, 66, 67]

Figure 5 illustrates some results for the Fe[(NCH)₆]²⁺ complex. Panel (a) shows energy differences between high- and low-spin states. A metaGGA called TPSS[68], when applied self-consistently, incorrectly yields the low-spin state as lower than the high-spin state, contradicting the QMC result. This is true of many semi-local functionals. Inclusion of a moderate fraction of exact exchange may bring the high-spin state slightly lower, but not enough (see blue curve in (c)). On the other hand, almost all functionals have the correct ordering when evaluated on HF densities, and most yield quite accurate spin gaps (red curve in (c)). Just as in the rainbow plot of Figure 4, there is a characteristic reduction in the spread of predictions when the HF density is used in a density-sensitive system. Finally, panel (b) shows the localized orbital locator (LOL)[69] for both calculations. One can see small differences in the bonding regions using LOL, because it is specifically designed to make such differences visible, but it is impossible to tell by visual inspection of densities or their differences which one is better and why.

In principle, DC-DFT is a much more general concept than those that appear in the literature: self-interaction, strong-correlation, delocalization, straight-line behavior of the energy as a functional of non-integer particle number, etc. It is based on a two-line decomposition of the error in any DFT calculation. Thus it can be applied to every functional approximation ever suggested and every DFT calculation ever performed, including the first ever Thomas-Fermi atomic calculations.[74, 75] Our focus here has been on KS calculations based on semi-local approximations, where the HF density typically works to cure significant density-driven errors, but DC-DFT can be applied much more broadly. DC-DFT analysis will usually provide different insights to these traditional analysis tools and may be less useful. But more importantly, DC-DFT can unite aspects of these other characterization tools.

With the tools provided by DC-DFT, it is straightforward to address this question from a pragmatic viewpoint. The first and foremost point is that, despite its name, the primary purpose of (ground-state) DFT is to produce ground-state energies for different molecular configurations, not densities. Few users ever output or examine the density closely, precisely because it is not what matters to their results. Thus the success of DFT in predicting those energies does not depend on how accurately approximations reproduce the density. Nevertheless, when the density or a property computed from it (e.g., electrostatic potential, partial atomic charges, etc.) is of interest to a user, it is usually better to use HF densities than SC ones provided that a calculation is density sensitive. This is illustrated later in Figure 9, where we compare HF and DFT atomic charges as we stretch NaCl.

Of course, the exact functional reproduces both the exact energy and the exact density but, as we have seen, a functional which yields usefully accurate energies need not yield accurate densities. This leads directly to a second important point. No matter how one might choose to measure density errors (and there are infinitely many choices, including infinitely many reasonable ones), there must be some sense of scale. If density differences are miniscule, why should anybody care, as they will have essentially no impact on predicted energies? Thus DC-DFT is the perfect tool for answering this question, as it measures the accuracy of densities directly in terms of their impact on the energy.

Figure 8 shows density errors for two very different two electron systems. The first is the He atom which is density-insensitive. Errors in these densities have little effect on energies. Furthermore, given the error profiles, which approximation has the ‘best’ density? The ranking depends entirely on one’s choice of measure. On the right, we have the same errors for H^-. This is a case of extreme density sensitivity, and was used as the prototype[20] for understanding density sensitivity. These are among the ‘largest’ density errors found in self-consistent DFT calculations, yet they are comparable in magnitude to those in He. But the observant will notice that the approximate functionals (PBE and M06 here) have errors that do not integrate to zero! How can this be? In fact, a correct self-consistent calculation[20] has about 0.3 electrons escape entirely from the system (and a HOMO of exactly zero).

There are two lessons here. First, in almost all interesting cases and all those discussed in this article, the density errors are small and subtle. We have never been able to understand energetics from studying these small differences directly. The relationship between densities and approximate XC is just far too complicated. On the other hand, we never need to do this, as our measures are all based on calculable energies, which speak for themselves.

The second lesson concerns electron affinities, and is more subtle. A popular method for calculating anions with self-consistent semi-local functionals is to use a basis set similar to that used for the neutral, and find electron affinities by subtraction. This method works surprising well, often yielding errors smaller than those of ionization potentials. But a sure sign that the anionic calculation is unconverged is the existence of a positive HOMO.[81, 82, 83] The basis-set is artificially binding the last electron. The fully converged calculation is like that of H^- shown, with a fraction of the last electron lost to the void. The beauty of the DC-DFT treatment is that it produces a well-bound density for the anion without a positive HOMO, and yields accuracies comparable to the artificial methods in common practice. In fact, it was studies of this issue[84, 85, 86], mostly couched in the language of self-interaction, that ultimately led to the more general concepts of DC-DFT[20].

In fact, this is not quite the same as evaluating on the HF density, as the HF kinetic energy is not quite the same as the KS kinetic energy.[96] However, this difference has been found to be much smaller than the improvement typically provided by using HF densities in cases where density-driven errors are large.[45] (See also Figure 11 below). In other words, HF densities do far more good than harm for density-sensitive calculations. Moreover, to the extent practical with finite basis sets, the differences with using the exact density have been found to also be very small, i.e., use of HF-DFT yields almost all the benefits that the exact density would confer.

The functional error of LDA is enormous on this scale, but the density-driven contribution is very small, illustrating the normalcy of this system and the high typical accuracy of even LDA densities. But note the accuracy of the total energy of PBE, i.e., at $\alpha=0$. That this is accidental can be seen both by the increased deviation as $\alpha$ grows, but more relevant is that the total error is small because the functional and density-driven errors are both much larger in magnitude, but mostly cancel. In fact, this cancellation is greatest at $\alpha=0$ and is less effective as more of the exact answer is added. Technically, this makes the H atom density sensitive for PBE, but this is entirely due to the accidental accuracy at self-consistency.

Finally, we see that SC-SCAN has a larger energy error than SC-PBE for the H atom, but in fact this is all density-driven. The SCAN functional error is zero for the H atom, by construction, but its error is non-zero when performed self-consistently. It can be perfectly possible for an approximate functional to be designed to be self-interaction error free for exponential densities, and yet produce a finite density-driven error for the H atom, because it’s XC potential will be incorrect. Ironically, SCAN is less accurate for the H atom than PBE is, despite SCAN using the H atom as an appropriate norm.

Figure 1(d) showed the PES of self-consistent and DC PBE calculations of OH$\cdot$Cl^-. Figure 11 shows what is happening as OH$\cdot$Cl^- dissociates, breaking down the PBE errors into their density-driven and functional-driven components, i.e., it is like Figure 10, but looking here at differences. The shaded regions indicate uncertainties due to the limitations of KS inversion methods with atomic basis sets, which we use to reverse-engineer the ’exact’ density and KS orbitals from a correlated wavefunction.[45] The density-driven component is sufficiently large as to be off-scale beyond about 2.5 Å. Also plotted are the errors of HF-PBE and the small error due to using the HF kinetic energy instead of the KS kinetic energy.

But, you might be surprised. The parameters in such functionals are chosen by minimizing errors on large curated databases, such as the GMTKN55 collection of 55 databases.[12] This process matches the calculated energetic errors relative the exact energies, not the functional errors. Moreover, the densities used in the calculation are typically not quite self-consistent, as it is more involved[99, 100] to optimize a functional involving MP2 (which depends on orbital energies). Thus the finding of the best parameters has (very small) inefficiencies. Because these functionals are so accurate, it only requires a very little density-driven error to make them suboptimal.

We recently used this insight to demonstrate such issues, creating our own 1-parameter double-hybrid, BL1p, but optimized to minimize functional errors rather than total energy errors.[46] For standard semilocal density functionals, DC-DFT tells us to use HF densities in density-sensitive cases. But because of the inclusion of approximate ab initio correlation in double-hybrids, the fraction of exact exchange is typically much higher, and it is fine to always use HF density as long as training is done with that density. The crucial step is to train on the functional error alone, i.e., subtracting the density-driven contribution from energy errors. Thus BL1p fixes the failures of standard double-hybrids for typical density sensitive calculations (e.g., dissociation of NaCl)[46, 52], but also provides improvements for density insensitive cases as will be illustrated below.

In Figure 12, we show results for the DARC database of Diels-Alder reactions. These reactions are long known as cases where standard functionals like B3LYP fail badly, and even double-hybrids. The figure shows B2PLYP, perhaps the most popular current double-hybrid, along with the more recent DSD-PBEP86.[101] On the other hand, BL1p uses the exact same ingredients, contains only one empirically determined parameter and uses the HF density. Furthermore, BL1p is trained on atomization energies of ony 6 molecules, while DSD-PBEP86 is trained on many more datapoints. The improvement of BL1p over B2PLYP is remarkable, further reinforcing the need to account for DC-DFT even with the latest, greatest approximations. We have recently shown that DFT calculations for DARC reactions are density insensitive making their errors almost entirely functional-driven.[46] Nevertheless, our BL1p still gives improvement for the DARC dataset because it is designed by the minimization of functional errors, while density-driven errors are taken care already by its construction.

As we have seen, the concept of density-driven errors is becoming widespread in the chemical literature and to a lesser extent, in the materials world.[102, 103, 104] Moreover, increasing numbers of authors are finding that the selective use of HF densities does indeed significantly reduce density-driven errors. In this section, we list some of the more obvious limitations of the current theory and also where it might be expanded.

In its current form, DC-DFT has nothing to say about how semi-local DFT can be improved for these systems, as they do not appear to be density sensitive. The errors made by the semi-local approximations on the stretched bonds are not much different if one uses exact densities or approximate restricted densities. And evaluation of the approximate functionals on the exact densities still produces the large errors. Thus these are functional errors by our current classification scheme.

But did we not say that a success of DC-DFT was to improve the dissociation limit of many molecules? Yes we did, but these are heteronuclear molecules whose stretched limit is not symmetric, and whose HF density is much more accurate in that limit, because of charge localization. We suspect that some generalization of DC-DFT ought to be able to include both stretched H${}_{2}^{+}$ and H₂ but we have not yet found it.

Of course, a much more satisfactory approach is to construct functionals that yield both good energetics and good densities when performed self-consistently.[110] Perhaps the foremost approach is that of Yang and co-workers in this direction.[33] Their approach is designed to reduce delocalization errors (both functional- and density-driven components) by explicitly imposing the well-known linearity condition with respect to particle number.[111] Many failures in DFT are attributed to deviations from this condition.[112] Cohen and co-workers have recently constructed DeepMind 21, a functional where machine learning has been used to address deviations from linearity condition.[113] This and other machine-learning approaches[114, 115] are also very promising when it comes to building functionals that give both good energetics and self-consistent densities.

We conclude by simply noting that DC-DFT is based on a simple one-line decomposition of DFT errors, based on the variational principle. In the past, many aspects of this decomposition had been noticed and mused over in understanding DFT results, but DC-DFT is a formal analysis that puts all these disparate pieces (and disparate sources of error) together. The concepts of DC-DFT are appearing more and more frequently in the chemistry and materials literature, and calculations using DC-DFT are being reported. As long as researchers continue to use KS-DFT as a standard tool for scientific discovery, DC-DFT will play an ever-expanding role in analyzing the inevitable errors.

Acknowledgment

ES and SS are grateful for support from the National Research Foundation of Korea (NRF-2020R1A2C2007468 and NRF-2020R1A4A1017737). KB acknowledges funding from NSF (CHEM 1856165). SV acknowledges funding from the Marie Sklodowska-Curie grant 101033630 (EU’s Horizon 2020 programme). ES and SS thank Prof. Soo Hyuk Choi for useful comments on the illustrations.