The derivative discontinuity of the exchange–correlation functional

Abstract

The derivative discontinuity is a key concept in electronic structure theory in general and density functional theory in particular. The electronic energy of a quantum system exhibits derivative discontinuities with respect to different degrees of freedom that are a consequence of the integer nature of electrons. The classical understanding refers to the derivative discontinuity of the total energy as a function of the total number of electrons (N), but it can also manifest at constant N. Examples are shown in models including several hydrogen systems with varying numbers of electrons or nuclear charge (Z), as well as the 1-dimensional Hubbard model (1DHM). Two sides of the problem are investigated: first, the failure of currently used approximate exchange–correlation functionals in DFT and, second, the importance of the derivative discontinuity in the exact electronic structure of molecules, as revealed by full configuration interaction (FCI). Currently, all approximate functionals, including hybrids, miss the derivative discontinuity, leading to basic errors that can be seen in many ways: from the complete failure to give the total energy of H₂ and H₂⁺, to the missing gap in Mott insulators such as stretched H₂ and the thermodynamic limit of the 1DHM, or a qualitatively incorrect density in the HZ molecule with two electrons and incorrect electron transfer processes. Description of the exact particle behaviour of electrons is emphasised, which is key to many important physical processes in real systems, especially those involving electron transfer, and offers a challenge for the development of new exchange–correlation functionals.

---

I. Introduction

The total energy of a system of electrons moving in an external potential, v_ext(r), in density functional theory (DFT)^1,2 is given by

E[ρ] = T_s[ρ] + ∫ρ(r)v_ext(r)dr + J[ρ] + E_xc[ρ] (1)

with explicit expressions for the non-interacting kinetic energy,and Coulomb energyAll the unknown complexity and many-body physics are in the remaining term, the exchange–correlation functional, E_xc[ρ]. The orbitals used to construct the density, , are solutions of the Kohn–Sham equation

(−½∇² + v_ext(r) + v_J(r) + v_xc(r))ϕ_i(r) = ε_iϕ_i(r) (4)

where the Coulomb potential is and the exchange–correlation potential is given by the functional derivative of the exchange–correlation energy, . This is exact Kohn–Sham DFT. With the exact functional, the solution of the Kohn–Sham equation (eqn (4)) to minimise the total energy (eqn (1)) yields the exact energy and density of the Schrödinger equation. However, the exact form of E_xc[ρ] is unknown and it is necessary to use density functional approximations (DFA). DFAs have both an approximate energy expression, E^DFA_xc, and approximate Kohn–Sham potential, v^DFA_xc(r), giving rise to approximate density, ρ^DFA(r), and eigenvalues, {ε^DFA_i}.

A. Density functional approximations

There are many different functional forms, starting with semilocal functionals that range from the local density approximation (LDA)^3–5 to the generalized gradient approximation (GGA)^6–10 and meta-GGA functionals.^11–14 There are also many functionals that mix in Hartree–Fock exact exchange in some manner, such as hybrid functionals with a varying degree of constant admixture of exact exchange, from B3LYP (20% HF) to PBE0 (25%) to M06-2X (54%) and M06-HF (100%).^15–21 In the last decade many functionals have emerged that examine the idea of range-separation, with functionals that include all the long-range part of HF exchange (LC-BLYP, LC-ωPBE) to CAMB3LYP and ωB97xd^22–25 to mixing in only the short range part of Hartree–Fock.^26,27 All these functionals use only the occupied orbitals and fit in a general sense to the first four rungs of Jacob's ladder of density functional approximations.²⁸ On the fifth rung there have been some ideas that use the unoccupied orbitals and eigenvalues in functionals such as B2PLYP,²⁹ which mix in some MP2-like terms or the random phase approximation (RPA), for example direct RPA (dRPA)^30,31 which uses the Coulomb only response in the adiabatic fluctuation dissipation theorem.

Fig. 1 presents a similar figure to Fig. 8 of ref. 32 to illustrate the performance of a large range of functionals for a set of thermochemical data (the heats of formation of the G3 set³³) and a set of reaction barrier heights.^34,35 For each of the functionals the calculated error of an energetic quantity for every individual molecules is represented by a single dot, in Fig. 1 each red dot is the error in the heat of formation of a molecule in the G3 set and each blue dot is one of the errors in a reaction barrier height. The results for dRPA for thermochemistry are only for the G2 set³⁶ and the barriers are the DHB24 set³⁷ and the results for these are taken from the paper and ESI of ref. 38 all other calculations are post-B3LYP.³² This allows one to see the performance of many different functionals in a global manner. In addition to the usual thermochemistry and barriers for each individual functional we also plot the errors for the two simplest molecules in the whole of chemistry infinitely stretched H₂⁺ and infinitely stretched H₂. Individually these molecules are known to be difficult for functionals to describe with stretched H₂⁺ (ref. 39) epitomising self-interaction error^40,41 and stretched H₂ the problem of static correlation.⁴² The errors for these two molecules, as one can see in Fig. 1, are very large but importantly connected. One can see that in changing functionals it is only possible to improve one error but with a corresponding failure on the other, the two seem connected. No functional is able to describe correctly these two simple molecules. It is this connection between different systems that epitomises the challenge of making one functional that can act discontinuously for different particle numbers, which is a markedly different challenge to the usual atomization energies of the G3 set or barrier heights.

Fig. 1 Calculated errors in thermochemical heats of formation (red dots), reaction barrier heights (blue dots), and also the errors of stretched H₂⁺ (black circle) and stretched H₂ (black dots) for a wide variety of functionals.

While there has been much improvement in the prediction of thermochemistry and reaction barriers over many years, using many different ideas in functionals, there is no functional that can reproduce the energy of these two simple systems. This can be viewed in two ways: one, is that the challenges of chemistry are not so related to the electronic structure of stretched H₂/H₂⁺, which has lead to the concept that DFT (and more specifically DFAs) works well as long as one does not stretch bonds. It is hoped that these errors do not cause a problem in the systems under study. However, the other view is that if current approximations are not able to correctly describe these two simple systems, then it should not be expected that for an unknown chemical they will give the correct answer. The key is not to just focus on the system, but on the behaviour of the electrons themselves. The errors in stretched H₂/H₂⁺ show a fundamental failure to correctly describe the electrons in those molecules and, as such, the description of similar electronic structure in many other systems will also fail. If these errors are not corrected, the inconsistencies of functionals will continue to dominate over the true behaviour of electrons.

B. Newer ideas in functionals

There are several groups working on new ideas in DFT, which is greatly needed to address the qualitative problems that can be seen in simple model systems. For example, Gori-Giorgi and coworkers are looking at the strictly-interacting λ → ∞ limit of the adiabatic connection for ideas based on strictly correlated electrons (SCE).^43–45 The concept of SCE can deal with problem of stretched H₂. Other notions beyond DFT, such as partition DFT, have been developed to attempt to tackle some of these problems.⁴⁶ Burke and coworkers have also looked at exact Kohn–Sham calculations in 1D using the exact functional by doing a DMRG calculation via a density perspective.⁴⁷ The particle–particle RPA (ppRPA) has been developed for electronic structure theory by van Aggelen, Yang and Yang,⁴⁸ showing a relation with ladder coupled-cluster.^49,50 Becke also has used real space ideas to address non-dynamical correlation and delocalization error using inverse hole models.^51,52 There are schemes to combat this issue beyond DFT, from combining multiconfigurational methods with DFT^53–55 to embedding methods in quantum chemistry^56–59 as well as ideas in density-matrix functional theory. Hopefully, the culmination of these theories will provide new functionals to be able to apply to a large spectrum of chemistry and physics without the drawbacks of many of the currently used functionals.

C. Challenge for DFAs

In this paper we highlight the difficult question of the derivative discontinuity of the exchange–correlation functional. The issue of describing the energies of stretched H₂⁺ and stretched H₂ with the same functional epitomises this challenge in a clear manner, but there are many other ways to the view the problem. This is a much larger issue than just that of stretching molecules, it is to correctly describe the energy of the electrons in all situations. An improved functional should be able to accurately describe the general behaviour of electrons, especially in interesting physical processes where competing effects act equally. This is the key problem to tackle, so that DFT calculations can describe the important chemical reactions and responses to electric and magnetic fields that are needed for the correct understanding of the behaviour of electrons in enzyme catalysis, Li-ion batteries, solar cells and many other technological applications.

II. Derivative discontinuity of the energy versus number of electrons

The famous paper from Perdew, Parr, Levy and Balduz in 1982⁶⁰ showed that the energy for a system with fractional electron number is given by a straight line connecting integer electron numbers

E(N + δ) = (1 − δ)E(N) + δE(N + 1) (5)

ρ_N+δ(r) = (1 − δ)ρ_N(r) + δρ_N+1(r). (6)

The energy and density are piecewise linear with straight lines connecting the integer points. This means that, at the integers, both the energy and density show (or can show) derivative discontinuities. In most situations there is a large discontinuity in the density on changing electron number. This is especially true in closed shell molecules where the density difference between last electron added (given by ρ_N − ρ_N−1) is very spatially different from where the next electron is added (ρ_N+1 − ρ_N). Or, in other words, when the frontier orbitals are spatially (and energetically) different.

We will focus on understanding manifestations of the derivative discontinuity in the total energy of integer systems. However, much of the understanding of the DD is often related to the exact Kohn–Sham potential, which was shown by Levy and Perdew⁶¹ to undergo a jump by a constant when passing through the integer. This constant, C, is the derivative discontinuity, as ε → 0

v^N−ε_xc(r) = v^N_xc(r)

v^N+ε_xc(r) = v^N_xc(r) + C.

This can be confusing to understand. For example, for a functional such as LDA, what does it matter if the potential is shifted by a constant? If that shifted potential is put in to the Kohn–Sham equations eqn (4), it will give rise to identical orbitals and density, however the eigenvalues are shifted by a constant. If those orbitals and density are put in to the energy expression eqn (1) an identical energy will be obtained. This means that there is a discontinuous change only in our eigenvalues, not in the total energy. This question is part of the challenge of understanding the importance of exact conditions in DFT, in this case the relation between the derivative discontinuity and the energy expression in relation to the potential. This whole discussion is fraught with problems, but in this paper we will cement our understanding by finding molecules (or model systems) where these questions become clarified. In this work we will elaborate on the implications of the derivative discontinuity for the energies of systems with integer number of electrons, where there is no need to invoke eigenvalues, fractional numbers of electrons or ensemble densities. Nevertheless, it should be noted that orbital energies are of course useful for many purposes,^62–66 and ensemble densities with fractional electron numbers have developed many constraints on the exact exchange–correlation functional that are essential to approximate it accurately as well as provoking many stimulating ideas.^67–71

A. Hydrogen atom and flat plane condition

Consider the energy of a hydrogen atom, going from H⁺ to H⁻, but with a possibly fractional number of electrons, N = n_α + n_β, 0 ≤ N ≤ 2 and also n_α ≤ 1 and n_β ≤ 1. This is a very simple system to help understand fractional numbers of electrons. The behaviour of the exact energy is given by the flat plane condition, which is a very stringent test of approximate functionals.⁷² Especially important is the understanding of the behaviour at N = 1, where two planes intersect giving an energy derivative discontinuity, and to consider adding and subtracting fractional numbers of electrons. More complicated surfaces were also investigated in the work of Gal and Geerlings.⁷³

DFAs really struggle to describe the flat plane and they completely fail to recover the discontinuous behaviour seen in the exact behaviour of the energy at N = 1. It is this failure that is connected to (or the root of) all the subsequent failures that we see. For the n_α = n_β line, approximate functionals massively fail, as they need to know that on going from E[0.4α,0.4β] → E[0.5α,0.5β] → E[0.6α,0.6β] they have passed through one electron. Compare this with the edge of the flat plane, E[0.8α,0β] → E[1.0α,0β] → E[1.0α,0.2β], where there is a change of the orbital being occupied (from α to β) on passing through the integer. It is clear from the energy expressions of eqn (1) that without a change of orbitals, the only term that can give the discontinuity is the exchange–correlation term. This is where the failure of current functionals and the challenge for future functionals lie. There is also a large discontinuity in the density with electron addition f₊ = ρ_H⁻ − ρ_H very spatially different from electron removal f₋ = ρ_H − ρ_H⁺.

1. Hydrogen atom with 1 basis function. To simplify the argument, it is useful to consider the calculation of the hydrogen atom just using a single basis function. The main reason to do this is that now the density is constrained by the basis function to be completely determined up to a factor such that it is no longer discontinuous on passing through N = 1, with ρ_H⁻ = 2ρ_H. A secondary point is that the discontinuity in the energy is in this case greatly enhanced if the basis function is chosen as (note the discontinuity could be reduced slightly if the basis function is chosen to be give the correct density at N = 2, i.e.). Additionally, the use of one basis function provides a direct connection to the Hubbard model where there is also one basis function per site. One of the features of the flat plane condition is that it considers fractional numbers of electrons, which can lead to conceptual confusion as well as technical challenges in extending methods to fractional numbers of electrons.^74,75 However, the flat plane condition was developed to explain the root cause in functionals of a general problem that affects real systems and hence can be equally seen in integer systems.

Let us consider a cube of 8 hydrogen atoms (each with one basis function). For this system FCI calculations can be easily carried out with different numbers of electrons and spins as in Fig. 2, where we have used a very large (1000 Å) distance along the edge of the cube just for simplicity. The density at each H atom is constrained by the basis set so an increase in energy happens when more than one electron per site is added; this is exactly the same as in the Hubbard model, where the on-site repulsion U causes an increase in energy past half-filling. Of course if the basis set allowed, these electrons would be unbound from the molecule (or as in a real H atom they would be much more diffuse to slightly lower in energy). These issues could be circumvented by changing the nucleus to be a He atom so that the attraction to the nucleus makes it much more favourable in energy. The real challenge for functionals it is to give the line of discontinuity crossing 8 electrons (one electron per H atom).

Fig. 2 The total energy for different number of electrons of an H₈ cube separated by 1000 Å (sketched in inset (a)), with one basis function per hydrogen, calculated by FCI with up to 8α and 8β electrons.

Let us just consider a single line in the flat plane with N_α = N_β, i.e. only closed shell systems. For this line, calculations with any functional can be easily carried out as the density and orbitals are known. Fig. 3 shows the performance of HF, dRPA, PBE, B3LYP and FCI. First, the FCI energy with 16 electrons (i.e. two electrons per site) is the same as HF, as the basis set does not allow for any electron correlation. As DFT functionals such as PBE and B3LYP treat correlation in a completely different manner, they have a slightly lower energy at 16 electrons. We could have considered just exchange functionals but have left in the correlation part so that one can see the relatively small effect of dynamic correlation functionals.

Fig. 3 Energy of closed shell H₈ with a minimal basis set with 0, 2, 4, 6, 8, 10, 12, 14 and 16 electrons. Methods such as RHF, dRPA, BLYP and B3LYP completely miss any discontinuous behaviour in the total energy and the exchange–correlation part, as shown in the inset.

Overall, all the approximate methods completely fail to reproduce any discontinuous behaviour of the total energy, and have a smooth behaviour in contrast to the correct answer of FCI. For less than one electron per site the energy decreases by −0.5E_h per electron, and for more than one electron per site the energy increases by 0.13E_h per electron. This H₈ molecule clearly illustrates the same physics as fractional electron numbers in one single H atom and also the same error of functionals, the missing derivative discontinuity. This example illustrates the important point that the missing derivative discontinuity in functionals can manifest itself as an error in the energy of integer systems.

B. The 1-dimensional Hubbard model

The Hubbard model⁷⁶ has a very simple Hamiltonian; for a set of sites, i, with creation operators () and annihilation operators (c_i) and the number operator , it contains only hopping terms to nearest neighbour sites and an on site repulsion when occupied by two electrons. It is specified by two parameters t and U

The physics depends only on the ratio U/t so we will fix t = 1 and vary U. It is a much studied system in strongly correlated condensed matter physics as it is a very simple model which describes interacting electrons in narrow energy bands, and which has been applied to problems as diverse as superconductivity, band magnetism, and the metal–insulator transition. The interplay between the delocalised hopping, t, and the localised repulsion, U, can lead to interesting balance in physical behaviours. Here we will examine the 1-dimensional Hubbard model (1D-HM) to highlight the connection with hydrogen atoms and derivative discontinuities, especially that seen at half filling. For the 1D-HM the exact answer is known in the thermodynamic limit using the Bethe-Ansatz, for example, the gap of the 1D-HM is given by⁷⁷

where J₁ is the first order Bessel function.

The Hubbard model is symmetric around half filling (except that above half-filling an electron-interaction term is included) i.e. for a system with 2N sites and doping fraction of M/N

E[N + M] = E[N − M] + MU. (7)

This leads to a clear picture of the derivative discontinuity that exists in the Hubbard model at half-filling, and raises the question of how to include this physics in a functional. This is exactly the same as the question of how to get a gap at the whole line in the flat plane condition of the hydrogen atom. The key of how to generally put this in a functional is distinct from how to predict the gap of the Hubbard model, where of course the knowledge of the system and the property in eqn (7) can be more specifically used to get the gap. For example, Capelle and coworkers^78–80 have several specific functionals (BA-LDA (LSOC, FVC)) and different parameterizations that are able to give the gap of the Hubbard model.

We consider finite Hubbard rings (chains with periodic boundary conditions) where exact results can be computed using FCI. The size of the chain can be increased and the large number of site limit corresponds to the thermodynamic limit. The two site Hubbard model has a very clear connection to infinitely stretched hydrogen molecule, with one electron this corresponds to H₂⁺ and with two electrons to H₂. Fig. 4 shows the total energy for two different Hubbard models, both with t = 1 and a large values of U = 10. The two site model is examined in Fig. 4a. The performance of methods is exactly as expected from calculations on stretched H₂⁺/H₂; HF and MP2 are good for one electron but fail for two electrons (HF is too high due to static correlation error and MP2 diverges to −∞ as U is increased). dRPA performs better for two electrons at the cost of a massive error for one electron.⁸¹ Here, dRPA would still give a gap at half-filling (N = 2). However, as shown for the 50 site model in Fig. 4b, the real problem arises in approaching the thermodynamic limit, where methods such as HF, dRPA or even MP2, smooth out and all semblance of a gap at half-filling disappears. However, the exact Bethe-Ansatz results give a very large gap. The complete failure of functionals to give a derivative discontinuity means that their results on the Hubbard model are completely physically incorrect as they miss one of the key behaviours of true Hubbard electrons.

Fig. 4 (a) Energy of the two site Hubbard model with different numbers of electrons calculated with FCI, HF, MP2 and dRPA. N = 1 is like a stretched H₂⁺ molecule and N = 2 is like a stretched H₂ and we can see how the errors of methods such as HF, dRPA and MP2 lead to an incorrect gap at N = 2. (b) Energy of a 50 site 1-d periodic Hubbard model with between 0 and 100 electrons calculated with HF, MP2 and dRPA compared with the exact Bethe-Ansatz result. In contrast to the exact Bethe-Ansatz result the DFA methods completely fail to give any sort of discontinuous behaviour at half-filling (50 electrons).

C. Fractional nuclei: HZ^{2e}

A picture of the discontinuous behaviour of the electrons is beautifully illustrated in the two electron example of an H₂ like molecule changing the nuclear charge of one of the protons to be fractional, giving HZ^{2e}.⁸² This encompasses a set of systems connected by a very simple change to the one-electron potential. This is a smooth and continuous change to the Hamiltonian, however, how does the electronic structure behave on these small changes, does it also change smoothly? As demonstrated in Fig. 5 the answer is, of course, that it depends. In some cases, as illustrated in short bond distances of HZ^{2e}, the electron also moves smoothly on this change. However, at long distances, the electron moves discontinuously, being either on the H or the Z, it is not shared between them.

Fig. 5 Occupation of Z atom, 〈n_Z〉, as the nuclear charge on the Z atom is varied at four different bond lengths comparing FCI and several different approximate DFT functionals.

The true behaviour of electrons, as given by FCI calculations, is simple to understand for HZ^{2e} at stretched bond lengths. When Z = 0 there are two electrons on the H atom (i.e. H⁻) and as Z increases an electron moves from the H to the Z when the energy of putting one electron on the Z atom gives a lower energy. This occurs when the energy of one electron on the Z atom, , is lower than the negative of the electron affinity of the H atom (EA(H) = 0.0277 a.u.). This happens at Z = 0.235. The atoms are too far apart to be bonded and no fractional electron transfer happens, as can be understood from the PPLB straight line interpolation of the true FCI energy of eqn (5). Something similar happens around Z = 1.67, where it becomes more energetically favourable to have two electrons on the Z atom and none on the H (at that point the electron affinity of the Z atom crosses 0.5). This is simple and clear, it is just the qualitative failure of functionals that is surprising.

Consider Z = 1, i.e. the H₂ molecule, all functionals get a qualitatively correct density due to the symmetry of the problem. However, they respond completely incorrectly to a small change in one of the atoms. DFAs smoothly move electrons when in fact they should not do so; no fractional charge is seen from FCI calculations. This turns the well known static-correlation error in the energy into a qualitative failure in the density. Approximate energy functionals are not able to describe the integer nature of electrons and they do not penalise correctly the splitting up of an electron in these stretched cases.

The HZ^{2e} system is very interesting and has advantages over very similar physics that can be seen in asymmetric Hubbard models or Anderson models in that functionals can be simply tested and their performance directly analyzed in real space. The conclusion is that they all fail for this problem. The reasons can be traced back to the failure of functionals for the closed-shell line of the flat plane (see Fig. 3). Along that line these functionals have no clear knowledge that they have passed through one electron and this translates into the problems seen in the density for these systems. This example is very illustrative but it is a slightly hard test as it requires self-consistent determination of the density. For methods such as dRPA, which is usually evaluated with PBE (sometimes HF) orbitals, it requires a self-consistent calculation⁸³ to highlight a problem in the density rather than energy, as carried out for Fig. 5.

III. The derivative discontinuity as a challenge for v_xc(r)?

Our main view is that the derivative discontinuity must be understood as a challenge for the energy, E_xc[ρ]. For example, in the HZ^{2e} problem, a functional that correctly gives the energy for the transfer of an electron would also give a qualitatively correct electron density and energy. For the set of HZ^{2e} systems, if the energy is wrong then the density will follow, giving rise to a qualitatively incorrect transfer. An equivalent, but in our view, confusing way to phrase this challenge is about the Kohn–Sham potential. To illustrate this point take the case Z = 1.5 for a stretched geometry as an example. Let us ask the question what is the exact Kohn–Sham potential that gives rise to integer number of electrons on each atom. This is clear in the HZ^{2e} system, as we have access to the exact Kohn–Sham potential, v_s(r) = v_ext(r) + v_J(r) + v_xc(r). This comes from a very simple rearrangement of eqn (4), substituting in the fact that the density is only made of one orbital, , giving rise toIn this case, the exact Kohn–Sham potential is directly available from an exact density in a FCI calculation.

Many plots have been seen in the literature, in both ground state and time-dependent analysis of the potential, for example.^84–89 The way to understand any features of the potential is that they are present as a consequence given a particular electron density, i.e. the exact Kohn–Sham potential (eqn (8)) is just a restatement of the density. This is shown in Fig. 6, where we just evaluate eqn (8) with an already minimised FCI density. If the potential shows any bumps it is because it comes from a density that gives rise to that structure, not the other way around. Furthermore, what gives rises to such a density is what is energetically favourable (for example, to have one electron each end), so it is the energy that is key.

Fig. 6 Kohn–Sham potential of HZ^{2e} for Z = 1.5 with a bond length of 10 Å. The features of the potential given by labels (1)–(5) are explained in the text and the inset shows the number of electrons on the Z atom for all HZ^{2e} systems.

Fig. 6 is produced by FCI calculation in a large basis set, but it is completely understandable and the same as that given by a density of the form ρ(r) = n₁e^−2r + n₂e^−3(r−10). The divergences at the nuclei (labels (1)) are because ∇²ρ goes as . For an exponent e^−2αr, the potential at large r goes to , so on the side of the H atom, e.g. label (2), it goes to a value of +½ and on the side of the Z = 1.5 atom it goes to a value of at label (3). The bump in the middle at label (4) is where ∇ρ goes to zero because of the overlapping densities, and the second term on the RHS of eqn (8) disappears, so a value of roughly the average of 1² and 1.5² is obtained. Finally, the change at label (5) is because at long range the e^−2r from the H atom dominates over the e^−3r behaviour from the Z atom, and from label (5) onwards the structure is a continuation of the line approaching the bump at label (4). Understanding all these features in the potential is of course relatively simple and known, and it should help to dispel any mysterious nature of them.

We want to stress that the bump in the middle (at label (4)) which has repeatedly been related to the derivative discontinuity, is just because ∇ρ goes to zero at some point in between the atoms where the density is very small, and this can even be thought of as a non-covalent interaction (NCI).⁹⁰ It is not related to any deeper physics and in particular is nothing to do with stopping electrons moving from one side to another. It should also be noted that the one thing that is not determined by this potential is the constants in front of the density (n₁ and n₂), this always cancels out. As such, it could be possible to have 0.8 electrons on the H atom and 1.2 electrons on the Z atom with an identical potential (note that the Z atom density would not respond to having more than 1 electron for the potential to remain unchanged). In general, these steps and bumps do not stop the electrons moving, however, having a correct energy functional does.

IV. Unrestricted calculations

Most of the problems we have highlighted to do with the integer nature of electrons and the derivative discontinuity are in fact well captured by unrestricted type methods. For example, the energies of infinitely stretched H₂⁺ and infinitely stretched H₂ are both given by UHF. As H₂ is stretched, there is at some bond distance a symmetry breaking (at the Coulson–Fischer point) beyond which the alpha and beta densities can be found each on one of the atoms. Therefore, the spin density is incorrect but the total density and energy are very good. Also for HZ^{2e} UHF recovers a discontinuity, and even for the Hubbard model UHF does very well, giving qualitatively correct gaps for all values of U. Of course, there are other known problems for UHF, such as stretching of F₂, for which RHF gives a reasonable minimum in terms of geometry but the minimum is actually above the dissociated atoms. RHF also dissociates incorrectly due to static correlation error and in this case the change to UHF does more than just correct the infinite limit, it actually means that UHF curve has no minimum⁹¹ (the Coulson–Fischer point is before the minimum in the RHF curve). For F₂ one may argue that a method such as UB3LYP gives a reasonable representation of the stretching. Another similar example is given by the stretching of O₂²⁺.⁹²

Another system that shows up a qualitative problem of the UHF method is that of stretching of He₂⁺. Here, the problem is in symmetry breaking, as illustrated in Fig. 7. Infinitely stretched He₂⁺ with symmetry gives two He^½+ atoms. UHF would give too high an energy due to its localization error (a concave behaviour of the energy for fractional systems). However, this too high energy is avoided by breaking the symmetry to give dissociation products of He and He⁺. This itself does not seem wrong, but the examination of the dissociation curve indicates that the symmetry breaking occurs at too short bond lengths, leading to an incorrect smooth transfer of electrons around 1.8 Å. Also the UHF curve falls off far too quickly compared to the FCI curve. This, of course, is very similar to the symmetry breaking in the spin densities seen in the UHF stretching of H₂, which is often argued to be acceptable⁹³ as in H₂ the total density is not qualitatively wrong, whereas in the case of He₂⁺ it is incorrect. For He₂⁺, DFT methods such as UB3LYP dissociate with half an electron on each atom but give a completely wrong and much too low energy due to the delocalisation error. Similar symmetry breaking by UHF can be seen in many systems, for example.⁹⁴

Fig. 7 Unrestricted Hartree–Fock stretching of He₂⁺ compared with FCI. Shown at the top is the difference in charge between the two He atoms. FCI has ΔCharge = 0 for all bond distances, however UHF incorrectly breaks the symmetry at around 1.8 Å, with an incorrect smooth transfer of half an electron from one atom to the other.

V. Fractional electron transfer coordinate

To understand in more detail the problem of the derivative discontinuity in HZ, let us consider the simplest case of HZ stretched to large distance (like 1000a₀). We now examine the transfer of electrons from one end of the molecule to another, in just one particular HZ system. In a single basis function per atom calculation this is very easy, as FCI gives just a trivial number of states. Let us first look at the system with one electron, where the FCI states are either one electron on the H, Ψ_H with energy E_H, or one electron on the Z, Ψ_Z with energy E_Z. Let us now consider a state that is a general coherent sum of these, , 0 ≤ α ≤ 1. As Ψ_H and Ψ_Z are FCI wavefunctions (one is the ground state the other is the first excited state), they are orthogonal and eigenfunctions of the many-body Hamiltonian, so it is trivially obtained that E[Ψ_α] = αE_H + (1 − α)E_Z and ρ_α(r) = αρ_H(r) + (1 − α)ρ_Z(r). This is very akin to fractional numbers of electrons as given by PPLB (eqn (5)), such that when α is varied the energy varies linearly and the electron moves smoothly from the H to the Z. In contrast to PPLB all possible values of electron transfer, α, correspond to an integer system with one electron and are represented by a wavefunction.

HZ with two electrons can be analyzed similarly. In this case, at the stretching limit, there are three singlet states with first order density matrices , and . The FCI wavefunctions that reduce down to give these density matrices are eigenfunctions of the Hamiltonian and orthogonal, so a linear combination of them will give a linear combination of both density matrices and energies. With more basis functions the idea is similar but the analysis more complex as there are many more possible states (one basis function excludes any excited states of the atoms) (Fig. 9).

Consider the case of HZ with Z = 1.2 from FCI and compare it with a functional such as PBE, as shown in Fig. 8. First, for one electron (HZ^{1e}), where Hartree–Fock and FCI are equivalent, there is a straight line interpolation between the energies of the two atoms. Of course, a minimization of the FCI leads to one electron wholly on one side or the other, in this case the Z atom, as it is much lower in energy. The behaviour of the energy with a functional such as PBE is qualitatively incorrect for one electron, as it does not have the correct linear straight line interpolation, but instead the energy varies smoothly (almost parabolically) with electron transfer. This leads to an incorrect minimum at around 0.26 electron on the H atom and 0.74 electron on the Z atom. The same result could be found with the compatible fractional calculations on each atom and piecing them back together, however, it is very good to see it in an integer electron system with a corresponding wavefunction. To summarise, in Fig. 8, PBE gives a good energy for the ground state Ψ = Ψ⁽⁰⁾ (and also for the first excited state Ψ = Ψ⁽¹⁾) but incorrectly gives a much lower energy for a wavefunction .

Fig. 8 (a) The HZ^1e system with Z = 1.2 for the wavefunctions , with −½ ≤ α ≤ ½. PBE incorrectly gives a minimum at α = 0.24, which corresponds to a qualitatively incorrect wavefunction with fractional charges on both atoms. (b) HZ^{2e} with Z = 1.2 along the electron transfer coordinate that moves both electrons from the H atom at α = −1 to both electrons on the Z atom at α = 1. For α > 0, , and α < 0, . PBE and RHF both transfer fractional numbers of electrons because of the incorrect shape of their energy surfaces compared with the FCI energy surface for electron transfer.

Fig. 9 The total energy of HZ with Z = 1.2 along the coordinate of transferring electrons for 0 ≤ N ≤ 2. The line in blue shows the energy at one electron and the line in black at two electrons, the red lines correspond to fractional number of electrons and are given by PPLB. The left hand side of the picture is when all the electron(s) are on the H atom and the right-hand side when all the electron(s) are on the Z atom.

For the two electron system (HZ^{2e}), FCI has a minimum with one electron on each atom (Ψ⁽⁰⁾) and two excited states, one with two electrons on the Z atom (Ψ⁽¹⁾) and the other with two electrons on the H atom (Ψ⁽²⁾). The straight line interpolations to be considered are between the ground state and first excited state and between the ground state and the second excited state. These are the lowest energy ones, a higher energy would be given by the interpolation between the first and second excited states. For any of these wavefunctions, the HF and PBE energies can be trivially evaluated from the density matrix, given that T_s[ρ] for any two electron system is the von-Weizsäcker expression . It is observed in Fig. 8 that both HF and PBE qualitatively fail to describe the electron transfer in this system. The energy on electron transfer is incorrectly smooth for both methods and with minima at the wrong values, HF at 0.33 electron on the H and 1.67 electron on the Z atom, and PBE with 0.4 electron on the H and 1.6 electron on the Z atom. In terms of the wavefunction, PBE is able to give a good energy for the excited states Ψ = Ψ⁽¹⁾ and Ψ⁽²⁾, but due to its static correlation error PBE is unable to give a good energy for Ψ = Ψ⁽⁰⁾. This means that PBE gives an incorrect minimum at . In general, density functional approximations favour an incorrect fractional charge transfer. Understanding all the possible states of this one system is perhaps a much simpler challenge than understanding many different systems (i.e. changing the molecule by changing Z) and yet captures the same physics.

VI. Perspectives for the future

The key picture of the derivative discontinuity of the total energy is shown by FCI calculations in Fig. 2 and 3, where the density increases smoothly but the derivative of the energy is discontinuous on passing through one electron per site. Therefore, there is an intrinsic discontinuity in the exchange–correlation term that is a consequence of the particle nature of electrons. Approximate functionals in the literature completely miss this behaviour and the failure in the total energies is clear as, for example, in the qualitative breakdown to describe the energies of stretched H₂ and stretched H₂⁺. This error can be transformed into an error in the density as shown in the two electron example, HZ^{2e} (Fig. 5). This leads to qualitative failures to describe charge transfer, with an artificial bias to fractionally transfer electrons. Remarkably, this is seen in systems with integer number of electrons characterised by a wavefunction, in contrast to the delocalization error typical of fractionally charged systems.

It is clear that the problems caused by missing the derivative discontinuity are not just about stretching molecules, but it is the same physics that occurs in transition metal complexes, chemical reactions, or especially in electron transfer processes. Our hope is that the use of simple chemical model systems gives a more complete understanding of the nature of all types of electronic behaviour that occur in the intricate nature of electronic structure.

The bumps of the exchange–correlation potential at the bond regions, where the density is very small, have been shown to have no effect on how electrons move and, therefore, do not capture the challenge of the derivative discontinuity. This illustrates the point that when the understanding appears quite paradoxical, it is just a clue that a deeper comprehension of the problem is needed.

We have highlighted the importance of the derivative discontinuity as a challenge for the energy functional at an integer number of electrons. We hope that understanding the examples given in this work can help highlight the avenue for development of new exchange–correlation functionals that contain the physics of the derivative discontinuity and represent the integer nature of electrons. This is needed so that DFT can play its role in helping tackle many important technological applications by correctly describing the movements of electrons in systems such as batteries and solar cells, chemical reactions in proteins, and transition metal compounds.

Acknowledgements

We gratefully acknowledge funding from the Royal Society (AJC) and Ramon y Cajal (PMS). PMS also acknowledges grant FIS2012-37549 from the Spanish Ministry of Science.

References

1. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864,  DOI:10.1103/PhysRev.136.B864.
2. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133,  DOI:10.1103/PhysRev.140.A1133.
3. P. A. M. Dirac, Proc. Cambridge Philos. Soc., 1930, 26, 376.
4. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200.
5. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 13244.
6. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098.
7. C. T. Lee, W. T. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
8. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
9. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8822.
10. N. C. Handy and A. J. Cohen, Mol. Phys., 2001, 99, 403.
11. A. D. Becke, J. Chem. Phys., 1988, 88, 1053.
12. A. D. Becke and M. R. Roussel, Phys. Rev. A: At., Mol., Opt. Phys., 1989, 39, 3761.
13. J. M. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
14. Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101.
15. A. D. Becke, J. Chem. Phys., 1993, 98, 1372.
16. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
17. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623.
18. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158.
19. A. D. Becke, J. Chem. Phys., 1997, 107, 8554.
20. F. A. Hamprecht, A. J. Cohen, D. J. Tozer and N. C. Handy, J. Chem. Phys., 1998, 109, 6264.
21. T. W. Keal and D. J. Tozer, J. Chem. Phys., 2005, 123, 121103.
22. H. Iikura, T. Tsuneda, T. Yanai and K. Hirao, J. Chem. Phys., 2001, 115, 3540.
23. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51.
24. O. A. Vydrov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 234109.
25. J. D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615.
26. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207.
27. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906.
28. J. P. Perdew and K. Schmidt, in Density Functional Theory and its Application to Materials, AIP Conference Proceedings, ed. V. VanDoren, C. VanAlsenoy and P. Geerlings, 2001, vol. 577, pp. 1–20.
29. S. Grimme, J. Chem. Phys., 2006, 124, 034108.
30. D. Bohm and D. Pines, Phys. Rev., 1952, 85, 332.
31. F. Furche, J. Chem. Phys., 2008, 129, 114105.
32. A. J. Cohen, P. Mori-Sánchez and W. T. Yang, Chem. Rev., 2012, 112, 289.
33. L. A. Curtiss, K. Raghavachari, P. C. Redfern and J. A. Pople, J. Chem. Phys., 2000, 112, 7374.
34. Y. Zhao, B. J. Lynch and D. G. Truhlar, J. Phys. Chem. A, 2004, 108, 2715.
35. Y. Zhao, N. Gonzalez-Garcia and D. G. Truhlar, J. Phys. Chem. A, 2005, 109, 2012.
36. L. A. Curtiss, K. Raghavachari, P. C. Redfern and J. A. Pople, J. Chem. Phys., 1997, 106, 1063.
37. J. J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2009, 5, 808.
38. Y. Yang, H. van Aggelen, S. N. Steinmann, D. Peng and W. Yang, J. Chem. Phys., 2013, 139, 174110.
39. R. Merkle, A. Savin and H. Preuss, J. Chem. Phys., 1992, 97, 9216.
40. J. P. Perdew and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048.
41. P. Mori-Sánchez, A. J. Cohen and W. T. Yang, Phys. Rev. Lett., 2008, 100, 146401.
42. A. J. Cohen, P. Mori-Sánchez and W. T. Yang, J. Chem. Phys., 2008, 129, 121104.
43. P. Gori-Giorgi, M. Seidl and G. Vignale, Phys. Rev. Lett., 2009, 103, 166402.
44. P. Gori-Giorgi and M. Seidl, Phys. Chem. Chem. Phys., 2010, 12, 14405.
45. A. Mirtschink, M. Seidl and P. Gori-Giorgi, Phys. Rev. Lett., 2013, 111, 126402.
46. J. Nafziger and A. Wasserman, arXiv preprint arXiv:1305.4966, 2013.
47. L. O. Wagner, E. Stoudenmire, K. Burke and S. R. White, Phys. Chem. Chem. Phys., 2012, 14, 8581.
48. H. van Aggelen, Y. Yang and W. Yang, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 88, 030501.
49. D. Peng, S. N. Steinmann, H. van Aggelen and W. Yang, J. Chem. Phys., 2013, 139, 104112.
50. G. E. Scuseria, T. M. Henderson and I. W. Bulik, J. Chem. Phys., 2013, 139, 104113.
51. A. D. Becke, J. Chem. Phys., 2005, 122, 064101.
52. A. D. Becke, J. Chem. Phys., 2013, 138, 074109.
53. A. Savin, Int. J. Quantum Chem., Symp., 1988, 22, 59.
54. J. Grafenstein and D. Cremer, Chem. Phys. Lett., 2000, 316, 569.
55. A. Stoyanova, A. M. Teale, J. Toulouse, T. Helgaker and E. Fromager, J. Chem. Phys., 2013, 139, 134113.
56. G. Knizia and G. K. L. Chan, Phys. Rev. Lett., 2012, 109, 186404.
57. G. Knizia and G. K. L. Chan, J. Chem. Theory Comput., 2013, 9, 1428.
58. C. Huang, M. Pavone and E. A. Carter, J. Chem. Phys., 2011, 134, 154110.
59. F. Libisch, C. Huang, P. Liao, M. Pavone and E. A. Carter, Phys. Rev. Lett., 2012, 109, 198303.
60. J. P. Perdew, R. G. Parr, M. Levy and J. L. Balduz, Jr., Phys. Rev. Lett., 1982, 49, 1691.
61. J. P. Perdew and M. Levy, Phys. Rev. Lett., 1983, 51, 1884.
62. E. J. Baerends, O. V. Gritsenko and R. van Meer, Phys. Chem. Chem. Phys., 2013, 15, 16408.
63. T. Stein, H. Eisenberg, L. Kronik and R. Baer, Phys. Rev. Lett., 2010, 105, 266802.
64. S. Refaely-Abramson, R. Baer and L. Kronik, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 075144.
65. L. Kronik, T. Stein, S. Refaely-Abramson and R. Baer, J. Chem. Theory Comput., 2012, 8, 1515.
66. T. Stein, J. Autschbach, N. Govind, L. Kronik and R. Baer, J. Phys. Chem. Lett., 2012, 3, 3740.
67. I. Dabo, A. Ferretti, N. Poilvert, Y. Li, N. Marzari and M. Cococcioni, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 115121.
68. E. Kraisler and L. Kronik, Phys. Rev. Lett., 2013, 110, 126403.
69. E. Kraisler and L. Kronik, J. Chem. Phys., 2014, 140, 18A540.
70. M. A. Mosquera and A. Wasserman, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 052506.
71. J.-D. Chai and P.-T. Chen, Phys. Rev. Lett., 2013, 110, 033002,  DOI:10.1103/PhysRevLett.110.033002.
72. P. Mori-Sánchez, A. J. Cohen and W. T. Yang, Phys. Rev. Lett., 2009, 102, 066403.
73. T. Gal and P. Geerlings, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 81, 032512.
74. W. Yang, P. Mori-Sánchez and A. J. Cohen, J. Chem. Phys., 2013, 139, 104114.
75. S. N. Steinmann and W. Yang, J. Chem. Phys., 2013, 139, 074107.
76. J. Hubbard, Proc. R. Soc. A, 1963, 276, 238.
77. E. H. Lieb and F. Y. Wu, Phys. Rev. Lett., 1968, 20, 1445.
78. N. A. Lima, M. F. Silva, L. N. Oliveira and K. Capelle, Phys. Rev. Lett., 2003, 90, 146402.
79. V. V. França, D. Vieira and K. Capelle, New J. Phys., 2012, 14, 073021.
80. K. Capelle and V. L. Campo Jr., Phys. Rep., 2013, 528, 91.
81. P. Mori-Sánchez, A. J. Cohen and W. T. Yang, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 85, 042507.
82. A. J. Cohen and P. Mori-Sánchez, J. Chem. Phys., 2014, 140, 044110.
83. P. Mori-Sánchez, W. T. Yang and A. J. Cohen, 2014, manuscript in preparation.
84. J. P. Perdew, in Density Functional Methods in Physics, ed. R. M. Dreizler, Providencia, 1985, pp. 265–308.
85. O. V. Gritsenko, R. van Leeuwen and E. J. Baerends, Phys. Rev. A: At., Mol., Opt. Phys., 1995, 52, 1870.
86. N. Helbig, I. V. Tokatly and A. Rubio, J. Chem. Phys., 2009, 131, 224105.
87. D. G. Tempel, T. J. Martinez and N. T. Maitra, J. Chem. Theory Comput., 2009, 5, 770.
88. M. Hellgren and E. K. U. Gross, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 85, 022514.
89. J. I. Fuks, P. Elliott, A. Rubio and N. T. Maitra, J. Phys. Chem. Lett., 2013, 4, 735.
90. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-Garcia, A. J. Cohen and W. T. Yang, J. Am. Chem. Soc., 2010, 132, 6498.
91. V. Klimo and J. Tiňo, Collect. Czech. Chem. Commun., 1981, 46, 1365.
92. R. H. Nobes, D. Moncrieff, M. W. Wong, L. Radom, P. M. W. Gill and J. A. Pople, Chem. Phys. Lett., 1991, 182, 216.
93. J. P. Perdew, A. Savin and K. Burke, Phys. Rev. A: At., Mol., Opt. Phys., 1995, 51, 4531.
94. T. Stein, C. A. Jiménez-Hoyos and G. E. Scuseria, J. Phys. Chem. A, 2014 DOI:10.1021/jp410713d.

---