Towards non-local density functionals by explicit modeling of the exchange-correlation hole in inhomogeneous systems

# Towards non-local density functionals by explicit modeling of the exchange-correlation hole in inhomogeneous systems

K.J.H. Giesbertz    R. van Leeuwen Department of Physics, Nanoscience Center, University of Jyväskylä, P.O. Box 35, 40014 Jyväskylä, Survontie 9, Jyväskylä, Finland    U. von Barth Department of Mathematical Physics, Institute of Physics, Lund University, Sölvegatan 14A, Lund S-22362, Sweden
July 27, 2019
###### Abstract

We put forward new approach for the development of a non-local density functional by a direct modeling of the shape of exchange-correlation (xc) hole in inhomogeneous systems. The functional is aimed at giving an accurate xc-energy and an accurate corresponding xc-potential even in difficult near-degeneracy situations such as molecular bond breaking. In particular we demand that: (1) the xc hole properly contains electron, (2) the xc-potential has the asymptotic behavior outside finite systems and (3) the xc-potential has the correct step structure related to the derivative discontinuities of the xc-energy functional. None of the currently existing functionals satisfies all these requirements. These demands are achieved by screening the exchange hole in such a way that the pair-correlation function is symmetric and satisfies the sum-rule. These two features immediately imply (1) and (2) while the explicit dependence of the exchange hole on the Kohn–Sham orbitals implies (3). Preliminary calculations show an improved physical description of the dissociating hydrogen molecule. Though the total energy is still far from perfect, the binding curve from our non-local density functional provides a significant improvement over the local density approximation.

## I Introduction

The local density approximation (LDA) is the simplest functional in density functional theory (DFT) and has been around since the advent of DFT Hohenberg and Kohn (1964); Kohn and Sham (1965). Although the LDA has only a local dependence on the density, it has been tremendously successful in describing the ground-state properties of solids, surfaces and large molecules. Its shortcomings have been obvious from the beginning and an enormous effort has gone into the search for better approximations. This task has proved to be exceedingly difficult. One should, however, not be surprised or disappointed. Density functional theory provides a way to reduce the full, interacting many-body problem to a simple non-interacting one. Therefore, an accurate density functional for the total energy would provide a surprisingly simple way to solve the complicated many-body problem — at least as far as static properties are concerned. Nevertheless, in the past decades considerable progress has been made. With the advent of the generalized gradient approximations (GGAs) Becke (1988); Perdew (1991), bond lengths and atomization energies were greatly improved as compared to those of the LDA. Unfortunately, the exchange-correlation (xc) potentials of the GGAs have several undesirable features. In particular, they decay too fast outside finite systems Lee and Zhou (1991); Engel et al. (1992), unlike a correct decay. Consequently, neither the LDA nor the GGAs produce proper Rydberg levels van Leeuwen and Baerends (1994) and ionization potentials are too low.

An important class of extensions to the GGAs came from the observation that exchange energies are much larger than correlation energies. This suggests that exchange should be treated exactly, while using a GGA only for the correlation energy. Full inclusion of “exact exchange” does however not work well in practice, since there is a large cancellation of errors between the LDA or the GGA versions of the exchange and correlation energies. On the other hand, using only a portion of exact exchange combined with a GGA does yield quite accurate bond energies in molecules Becke (1993). An important advantage of “exact exchange” is that it provides some necessary improvements of the xc potential. For instance, the inclusion of “exact exchange” gives a stepped structure in the xc potential van Leeuwen et al. (1995, 1996) thus improving the description of the atomic shells as well as enabling a neutral dissociation of heteronuclear molecules Perdew et al. (1982); Gritsenko and Baerends (1996); Baerends and Gritsenko (1997). It also gives an xc-potential with the proper asymptotic () behavior which gives rise to Rydberg levels and a highest occupied KS eigenvalue in better agreement with the negative of the ionization potential. Usually, however, only a fraction of full exchange is incorporated in the so called hybrid functionals meaning that the desirable features mentioned above are only partially obtained. It seems that only functionals with a massive amount of fitting parameters are able to handle 100% “exact exchange” Gill (2001); Zhao and Truhlar (2008), though they will always suffer from weak singularities in the response functions of metals. We mention here that inclusion of “exact exchange” is not mandatory in order to have good properties of the xc potential. The proper step structure as well as the correct asymptotics away from finite systems can also be obtained by modeling the potential directly van Leeuwen and Baerends (1994); Gritsenko et al. (1995). Such model potentials can indeed provide good response properties Gritsenko et al. (1999); Schipper et al. (2000) but they cannot easily be written as the functional derivative of some accurate functional for the xc energy. As a result, their implementations have been limited. Further improvements to the energy functional are also sought by adding the kinetic-energy density to the functional arguments Ghosh and Parr (1986); Tao et al. (2003) and including parts from many-body perturbation theory Görling and Levy (1994); van Leeuwen (1996); Gatti et al. (2007); Romaniello et al. (2009); Hellgren and von Barth (2010).

In this paper we would like to consider an alternative approach by looking directly at electron correlations in real space. In a correlated system we can consider the conditional probability of finding an electron at some point in space, when the position of another reference electron is given. The difference between this function and the unconditional probability (which is simply the density) is defined to be the xc hole and the knowledge of this function 111To be precise we need to know xc hole for different strengths of the electron interaction. is sufficient to calculate the xc-energy.

The effect of exchange and correlation is to dig a hole in the density around each electron, so as to remove one electron in that region. We can say that the task of a good xc functional is to provide an accurate description of the xc hole. The LDA and GGA assume that this hole has a spherical shape and an extent given by the Wigner–Seitz radius (, being the electron density) at the reference electron. Unfortunately for the LDA, xc holes are definitely not spherical Buijse (1991); Buijse and Baerends (2002). Although hybrid functionals lift this restriction, their blunt use of “exact exchange” actually worsens the description of the xc hole for stretched bonds compared to the LDA and GGA functionals.

A more sophisticated class of functionals which aims to build a non-spherical model of the xc hole is the so-called weighted density approximations (WDA) Alonso and Girifalco (1978); Gunnarsson et al. (1979); Sadd and Teter (1996); Rushton et al. (2002). These functionals avoid the spherical xc hole by digging a hole out of the real density rather than in the density at the position of the reference electron. A nice feature of the WDA is that the asymptotic behavior of its xc potential has a Coulombic asymptotic decay instead of an exponential behavior as in the LDA and the GGA. An important symmetry of the pair-density () is, however, broken. This causes an additional factor in the asymptotic decay of the xc potential, so that it decays too fast (as Ossicini and Bertoni (1985).

Here, we will advocate an approach in a similar spirit as the weighted density approximation von Barth (2004). However, we will take care not to destroy the symmetry of the pair-density and therefore, the xc potential will automatically have the correct asymptotic behavior. Furthermore, important information on the physics of the xc hole is provided by the dissociation of molecules. In particular, a proper localization of the xc hole around the reference electron for a dissociated molecule is a challenging task. The failure of current approximations to achieve this is reflected in their consistent inability to properly describe the breaking of chemical bonds. The most important aim of our new functional will therefore be a bold one: the new functional has to be able to describe molecular dissociation.

The paper is outlined as follows. First (Sec. II) we will give a more detailed discussion on the background to motivate the construction of the screened exchange (SX) functional. The actual construction of the SX functional is discussed in section III. In section IV we show preliminary results and finally conclude in section V.

## Ii Motivation

### ii.1 Symmetry and asymptotics

We will start from an exact expression Almbladh (1972); Langreth and Perdew (1975); Gunnarsson and Lundqvist (1976) for the exchange-correlation energy

 Exc=12∫dr∫dr′n(r)n(r′)(¯g(r,r′)−1)|r−r′|, (1)

where is the density. This expression gives the exact xc contribution to the interaction energy of the system, provided is the exact pair-correlation function of the system, defined as

 g(r,r′)≡Γ(r,r′)n(r)n(r′).

Here the diagonal of the two-body reduced density matrix is defined as

 Γ(r,r′) ≡∑σ,σ′Γ(rσ,r′σ′) ≡∑σ,σ′⟨Ψ|^ψ†(rσ)^ψ†(r′σ′)^ψ(r′σ′)^ψ(rσ)|Ψ⟩,

where and are the usual field operators. The xc energy, however, also contains a correlation contribution to the kinetic energy which is most conveniently included by integrating the interaction energy with respect to the strength of the Coulomb interaction — while keeping the density fixed at the fully interaction one (Almbladh (1972); Langreth and Perdew (1975); Gunnarsson and Lundqvist (1976); von Barth (2004), i.e.

 ¯g(r,r′)≡∫10dλgλ(r,r′).

The physical picture of representing the xc energy in this manner is that an electron does not interact with the full density, but depending on its position, , it sees an effective density with electrons. Therefore, the part

 ¯ρxc(r′|r)≡n(r′)(¯g(r,r′)−1) (2)

in the xc energy (1) exactly describes the density of minus one electron, a hole, which is reflected in the sum-rule of the xc hole

 (3)

The local density approximation (LDA) proceeds by using the xc hole from the homogeneous electron gas evaluated for the density at the position of the reference electron, so the pair-correlation function is approximated as . Furthermore, if the distance is large the pair-correlation function only differs slightly from one, so it is quite reasonable to replace by in the xc-energy (1). Combining both approximations, we obtain the expression for the xc energy of the LDA

 ELDAxc =12∫dr∫dr′n(r)n(r){¯gh(|r−r′|;n(r))−1}|r−r′| =12∫dr∫dr′n(r)n(r){¯gh(|r′|;n(r))−1}|r′| =∫drn(r)εxc(n(r)), (4)

where the function is just the exchange-correlation energy per electron of the homogeneous electron gas. This expression for the xc energy of the LDA shows explicitly that its hole is approximated by a spherical one. As mentioned above xc holes are not spherical Buijse (1991); Buijse and Baerends (2002). The original weighted density approximation (WDA) Alonso and Girifalco (1978); Gunnarsson et al. (1979) improves on the shape of the xc hole by not replacing by in the xc-energy (1). Since the sum-rule (3) for the xc hole is no longer trivially satisfied, an effective density, , is used as input into the pair-correlation function instead of the density at the reference position. The xc energy in the WDA is therefore

 (5)

where the effective density should be found by satisfying the sum-rule for the xc hole (3) at every point

 ∫dr′n(r′){¯gh(|r−r′|;¯n(r))−1}=−1.

Unfortunately, it was found that the pair-correlation functional of the homogeneous electron gas is not a good approximation to the pair-correlation function of inhomogeneous systems as molecules and surfaces. Using a more localized function for , the results were significantly improved Gunnarsson and Jones (1980); Sadd and Teter (1996); Rushton et al. (2002).

From the expression for the WDA xc energy (5) we immediately notice that the integral kernel is still asymmetric in the coordinates and and therefore breaks the important symmetry . This is not so important for the value of the xc energy. But for the corresponding xc potential, we obtain

The first term actually decays as as is obvious from the sum-rule. The long-range behavior of the other two terms is not so obvious, but in practice they turn out to decay exponentially Sadd and Teter (1996); Rushton et al. (2002). Therefore, the xc potential decays too slowly (as ) compared to the proper decay Gunnarsson et al. (1979). The incorrect long-range behavior of the potential is expected to have a significant effect on properties which depend strongly on a proper xc potential such as the ionization energy and Rydberg excitations Almbladh and von Barth (1985); van Leeuwen et al. (1995); Schipper et al. (2000). This failure can be attributed to the broken symmetry of the pair-correlation function. Consequently, one of the requirements of our new functional will be that it should satisfy the proper symmetry of the pair-correlation function, i.e. .

### ii.2 Step structure from exchange

As mentioned in the introduction, the steps in the KS potential are important for a proper description of the atomic shell structure van Leeuwen et al. (1995, 1996) and also the neutral dissociation of hetero-nuclear molecules Perdew et al. (1982); Gritsenko and Baerends (1996). It has been shown that the necessary steps for the atomic shell structure are already featured by the exchange energy van Leeuwen et al. (1995). However, the necessary step in the dissociation of hetero diatomic molecules is a less clear case, since the spin-symmetry has to be broken to provide the necessary localization Fuks et al. (2011); Makmal et al. (2011).

Closely related, the exchange energy also has the necessary features for the integer discontinuity van Leeuwen et al. (1995); Gritsenko and Baerends (1996), since the exchange term often changes radically when crossing an integer number of electrons due to the usual idempotency of the KS density matrix. The corresponding hole, the exchange hole, can be expressed in spin-integrated quantities as

 ¯ρx(r|rref)≡−12|γs(r,rref)|2n(rref), (6)

where the spin-integrated KS density matrix is defined as

 γs(r,r′) ≡∑σγ(rσ,r′σ) ≡∑σ⟨Ψs|^ψ†(r′σ′)^ψ(rσ)|Ψs⟩,

with as the KS wavefunction. The KS density matrix can alternatively be expressed directly in terms of the KS orbitals as

 γs(r,r′)=∑knkϕk(r)ϕ∗k(r′), (7)

where are the occupation numbers, being simply 0 or 2 in the non-degenerate case. The exchange hole has the convenient property that it already satisfies the hole xc hole sum-rule (3).

Since exchange already satisfies a number of important properties, it is often used as a starting point to model the full exchange-correlation effects. Traditionally, one adds a correction term, correlation, defined as the difference between the exact exchange-correlation quantities and the ones with one exchange taken into account. For example the correlation hole is simply defined as

 ¯ρc(r|rref)≡¯ρxc(r|rref)−¯ρx% (r|rref).

Although this approach has some appeal, it is inconvenient in practice, especially in bond-breaking situations. We show that explicitly in the next section.

### ii.3 Bond breaking

In Fig. 2 we have plotted the different holes, , and for the H molecule at equilibrium distance Bohr. Note that the quantities with correlation are at full coupling strength () and not the integrated ones. Ideally we would like to have shown the integrated ones, but obtaining them is a rather difficult task. Since the effect of the kinetic energy is not very large on the total energy van Leeuwen (1994), we expect that the holes do not differ too much as well. The reference electron is fixed at 0.3 Bohr to the left from the right nucleus. The holes were calculated from full configurations interaction (CI) results using the 1s, 2s, 3s, 2p and 3d hydrogen wavefunctions on each atom as a basis set 222The 3p orbitals caused linear dependency problems at Bohr, so they could not be included.. The exchange hole (x hole) can be worked out to be

 ¯ρx(r|rref)=−2|σg(r)|2|σg(rref)|2n(rref)=−|σg(r)|2,

thus the exchange hole is actually independent of the position of the reference electron. However, the real hole is deeper around the reference electron and therefore, depending on the position of the reference electron, the correlation hole (c hole) has to add and remove an equal amount of the hole to deepen it around the reference electron. Although the xc hole is more localized around the reference position, it still shows the two-peak structure of the x hole.

The localization effect becomes more prominent if we consider the hydrogen molecule at an elongated bond distance of Bohr (Fig. 2). Again, the x hole is independent of the position of the reference electron and is completely delocalized over the molecule. However, the real hole is completely localized around the reference electron, which is again located at 0.3 Bohr to the left from the right nucleus. Therefore, the c hole has to completely remove the hole from the left side of the molecule at put it back at the right side such that it integrates still to electron. The c hole can not be regarded as a “small” correction to the x hole anymore, since it is actually equal in magnitude. The lack of the “small” c hole in the Hartree–Fock (HF) description correcting the x hole is the main reason for the bad performance of HF for the dissociation of molecules.

It is instructive to consider the LDA holes, since they explain why DFT and its predecessor, X Slater (1951), are so successful. The LDA holes are actually the -integrated ones, so a direct comparison is strictly not correct. Luckily, however, in the dissociation limit the -integrated and the exact hole at are identical van Leeuwen et al. (1996); Perdew et al. (1995). The reason is that at long bond distances the interaction term for will favor wavefunction configurations in which the electrons are residing on different atoms, i.e. the Heitler–London wavefunction. The density corresponding to this ground state wavefunction will be exactly equal to the wavefunction at full coupling strength. It thus follows that for at infinite bond distance. At the wavefunction simply remains a KS determinant with a doubly occupied orbital. However, the region becomes unimportant Teale et al. (2009) in the -integration for the pair-correlation function and hence . Therefore at the longer bond distance of Bohr the exact xc hole should be close to its -integrated counterpart.

In Fig. 3 we show the LDA holes for the hydrogen molecule at Bohr. From our discussion in Sec. II.1 it is clear what the definition of the xc hole of the LDA should be. The following expression Gunnarsson et al. (1979); Dreizler and Gross (1990) is consistent with the LDA energy expression (II.1)

 ¯ρLDAxc(r|r% ref)=n(rref)(¯gh(n(rref),|r−rref|)−1). (8)

Gori–Giorgi and Perdew made an accurate model for the pair-correlation function of the homogeneous electron gas Gori-Giorgi and Perdew (2002), which we used to calculate . The x hole can be calculated by using the exchange part of the electron gas pair-correlation function in this expression and the c hole is simply defined as the difference between the other two. The most striking feature of the LDA holes in Fig. 3 is that the x hole is localized instead of delocalized over the two atoms, just as is the case for the exact x hole. Although the LDA x hole does not resemble the exact x hole at all, the full xc hole (the one of interest) is actually modeled quite well. Especially, if we consider the exact x hole which is the hole employed in HF theory, the LDA xc hole provides a large improvement. We also see that the LDA c hole only provides a small correction to the x hole: it removes the outward oscillations and localizes the LDA hole a bit further. Since the correction from the c hole is so small, we can understand why the old X method already outperformed HF so much. In particular, the deepening of the hole was empirically taken into account by scaling the exchange prefactor, , from its theoretical value, , to . These observations concerning the LDA holes also make it clear that it does not make sense at all to add “exact exchange” to LDA correlation. The same holds for GGA holes, since they are quite similar to LDA holes, only having slightly more wild oscillations and a discontinuity due to their cut-off to satisfy the sum-rule (3) as additional features Slamet and Sahni (1991); Burke et al. (1998).

## Iii The screened exchange Model

Considering these facts about the x hole in dissociating H, it seems to be unwise to add a correlation hole to an exact exchange hole. It will be hard to build a model for the correlation hole with the proper strongly varying features. But we also know that exchange effects often dominate and that correlation effects only provide a modification of the former. This observation suggests that we should not add correlation to exchange but rather modify the shape of the x hole by some correlation factor. From the holes for the hydrogen molecule (Figs 2 and 2) we observe that the main effect of correlation is to localize the x hole. This is not special for the H molecule, but is the main feature of correlation in any system, even in the electron gas where the correlation reduces the long range behavior of the x hole from to  Gori-Giorgi et al. (2000a, b). A further example was provided long ago by J.C. Slater when he pointed out that atomic term energies were often accurately described by term dependent Hartree–Fock theory (“exact exchange ”) by simply reducing Slater’s and integrals by  Slater (1960).

Following the discussion above it seems rather natural to use the following Ansatz for the xc energy of an inhomogeneous system

 Exc=−14∫dr∫dr′|γs(r,r′)|2|r−r′|F(r,r′). (9)

Here, is the spin-integrated non-interacting density matrix of the KS orbitals (7) and the “screening”-factor is intended to provide the necessary modification of exact exchange which will take care of the effects of correlations. The exact expression for is

and by multiplying the numerator and denominator by we immediately see that the screening function is symmetric in its arguments and , cf. (2) and (6).

Our task is thus to find a reasonable model for and we stress again the importance of keeping the symmetry of in order to have an ensuing xc potential with the correct tail outside finite systems. We also believe it to be essential to have a model which satisfies the sum-rule for the xc hole and in terms of , our model should thus obey

 ∫dr′|γs(r,r′)|2F(r,r′)=2n(r),

If we, like the founding fathers (KS), first turn to the electron gas, we realize that , and also for that matter, must be spherical functions of only . In the spirit of the older WDAs we could thus attempt such an Ansatz for our function. It is, however, important here to stress that in the original WDAs it is almost the entire xc hole which is modeled in this way, whereas in our case we just model a modification of the full, in general non-spherical x hole. As mentioned previously, we would also like to model the -function in the case of the dissociation of H. In the dissociation limit of the hydrogen molecule the -function actually takes the form

 FHL(r,r′)≈fa(r)fa(r′)+fb(r)fb(r′)

in terms of two one-point functions and located on the different hydrogen atoms. This follows from the fact that in the limit of large separation between the nuclei, the Heitler–London wavefunction becomes exact—but we will not present the details here. But it means that in this limit, the -function is very far from spherical.

When the two electrons are on different nuclei, the -factor should vanish, because we are then dealing with two non-interacting subsystems and there is neither exchange nor correlation. When both electrons are on the same atom, the -funcion should actually be 2 to make the xc hole equal to the negative of the local density. In this way the xc hole will precisely remove the full self-interaction on each atom—not just half the self-interaction as Hartree–Fock does—and we recover the correct limit of two separate and non-interacting atoms.

As suggested by the above observations, the following Ansatz for the screening function might stand a chance of carrying us from the limit of a homogeneous system to that of the complete breaking up of the H bond

 FSX(r,r′)≡A(r)A(r′)h(|r−r′|,¯n(r,r′)). (10)

Here, the spherical function has an effective screening radius given by . This form was also inspired by the success of a wavefunction for the helium atom by Hylleraas Hylleraas (1929) having precisely this form.

We are then left with the choice of satisfying the hole sum-rule either by adjusting the one-point function , the hole-depth function, or by varying the effective density . We stress that the latter must be a symmetric two-point function in order not to jeopardize the asymptotics of the potential. In terms of the Ansatz (10) the sum-rule reads

 A(r)∫dr′|γs(r,r′)|2A(r′)h(|r−r′|,¯n(r,r′))=2n(r). (11)

The effective density is expected to vary in a similar way as the density itself (from very small to very large values) and it is a two-point function. The hole-depth function on the other hand is a one-point function of limited variation—typically between zero and two depending on the normalization of the function . Consequently, it appears to be numerically more stable to use the -function for the purpose of satisfying the sum-rule. Indeed, we have encountered no difficulties in solving (11) in our applications. A further argument in favor of this choice is a lack of guidance from the electron gas when determining the -factor, should we have chosen to satisfy the sum-rule by varying the effective density . Most WDA models have used the latter procedure, but again, they have not considered an -factor.

In our case we are thus free to choose the effective density. Typical choices are or . in previous attempts to construct a symmetric version of an WDA-like model we have found that the first of these suggestions gave rise to numerical difficulties, whereas Gunnarsson et al. Gunnarsson et al. (1979) encountered difficulties with the second choice. A choice which seems to work in our previous attempts is

 ¯n(r,r′)≡√n(r)n(r′), (12)

which is the definition of the effective density which we will use here. We stress, however, that there is no compelling reason for this choice. As a matter of fact, this is one part of our new model which we might have to revise in future attempts to improve the accuracy of the model.

It is not clear what properties and shape the screening function should have. However, since its main task is to screen the x hole, we will use the very simple form inspired by the screened Coulomb (Yukawa) interaction

 hHEG(r,¯n)≡e−D(¯n)r12. (13)

The function is fitted to the homogeneous electron gas such that it yields the exact xc energies for all homogeneous densities. In this way also the kinetic energy contribution to the xc energy is included. More details on the fit can be found in Appendix B. This form for the screening function is definitely too simplistic. More knowledge is required to build more accurate SX models. This will be the subject of future research.

## Iv Illustrative results

One of the most severe tests for our new SX functional will be if it performs well for dissociating molecules. To keep matters simple, we have limited ourselves to the hydrogen molecule. The most natural grid for calculations on a diatomic molecule is based on a prolate spheroidal coordinate system. A planar elliptical grid with foci at the two nuclei and -axis joining those foci is rotated about that axis to generate the full grid. More details on the grid used are in Appendix C. As a further simplification we used the density from a full CI calculation in the same basis as used before (1s, 2s, 3s, 2p and 3d hydrogen wavefunctions). This CI expansion is not very good for obtaining an accurate total energy. However, we believe that the density will be accurate enough for the SX model.

The first step in evaluating the SX model is to solve for the hole-depth function, , by solving the integral equation (11) on the grid. Once the hole-depth function is obtained, we performed the double integral (9) with  (10) to obtain the xc energy according to our simple SX model. To calculate the total energy, we note that the non-interacting kinetic energy can be directly obtained from the density for the two-electron system, , and the Hartree and nuclear contributions are already known from the full CI calculation. In Fig. 4 we compare the total energies from the SX model with the ones from a full CI calculation as a function of the distance, , between the hydrogen atoms. As mentioned before the Slater basis is too poor to obtain a good total energy, so we performed a full CI calculation with the dalton package Dal (2005) in an aug-cc-pVQZ basis Dunning Jr. (1989) as a reference and additionally the corresponding HF result is shown as well. Unfortunately our SX model with the simplistic screening function is not performing much better than the HF method. It follows the HF curve rather closely and only around Bohr does the curve start to bend downwards to the correct total energy. The most important feature of the SX model is that it directly models the xc hole and therefore we can look at it what is going wrong. In Fig. 6 in the left panel, we show the hole at Bohr. If we compare with the exact holes in Fig. 2, we see immediately that our current SX model is not localizing the x hole sufficiently. Actually, the SX hole is almost identical to the x hole . In the left panel of Fig. 6 we show the hole the elongated distance Bohr. Comparing with the exact holes in Fig. 2 we see that the SX model actually does localize the x hole, but not sufficiently. The lack of localization of the hole explains exactly why the total energy in the SX model is consistently too high (Fig. 4). Due to the term in the expression for the xc energy (1), a too delocalized hole does not stabilize the energy enough, which results in a too high energy.

In retrospect we should not be surprised that parametrizing the screening function by the homogeneous electron gas did not work out. The main task of the screening function is to localize the hole which is most important in inhomogeneous systems like the hydrogen molecule. Therefore, its actual form and localization strength should be modeled by these systems and not the homogeneous electron gas. A detailed study and parametrization are beyond the aims of this article, but to reinforce our arguments for this statement, we also like to show some results for the following two heuristic screening functions

 h1(r12,¯n) =exp(−c1(r12¯rs)), (14a) h2(r12,¯n) =exp(−c2(r12¯rs)2). (14b)

The first one, , keeps the Slater like form, but replaces the parametrization by the electron gas by a simple division by to make the total dimensionless and a constant that we can choose to our liking. We found that gave a nice dissociation behavior for the energy. The second one is mainly included to emphasize that the required shape of the screening function is also unclear at the moment. Choosing the constant in the same manner as before, we found to be sufficient for our purposes. Note that both screening functions are not fitted to the electron gas anymore and are therefore not expected to give the correct xc energy density, , for the gas.

The results for the energy from these screening functions are shown also in Fig. 4. The situation is now rather different. The total energy is consistently underestimated. However, the shape of the curve is definitely an improvement. The total energy at elongated distances Bohr remains rather constant as we choose the constants to do so. From the figure it is not immediately clear if the shape is also an improvement over the simple LDA functional. However, shifting the curves such that their minima coincide with the full CI minimum (Fig. 7), we see that the energy from the SX models follow the full CI energy much closer. In particular the Gaussian, seems to reduce the overbinding most. Even at equilibrium distance the position and shape of the minimum seems to be somewhat improved.

Again, since we have built a direct model for the hole we can explain both features. Considering the holes from the heuristic and screening functions in Figs 6 and 6, we see that both screening function are too powerful: they screen the x hole too much. This explains why the total energy is consistently lower than the full CI result: the hole becomes too compact. Since the hole is even more compact for screening function than the screening function, the energy of the screening function is even lower than the energy of the screening function. However, if we consider the shapes of the holes, then they are definitely an improvement over the previous screening function. Especially at the elongated bond distance Bohr, we see that both heuristic screening functions nicely localize the hole completely at the correct side, which explains their improved energy in the dissociation limit over the erroneous of HF.

The LDA hole is also shown in Figs 6 and 6. Compared to the HF, cf.  in Fig. 2 and 2, the LDA hole is a big improvement since it localizes around the reference electron. This improvement is also visible in the total energy where the LDA does not exhibit the erroneous behavior as HF does (Fig. 4). However, the shape of the LDA hole is ridiculous. It is spherical by construction and it does not have the peaked features at the nuclei. Instead, the LDA hole attains its minimum at the position of the reference electron. The SX holes, especially with the screening function, are an improvement over both the HF hole and the LDA hole. It correctly localizes the hole around the reference electron and still retains the distinct peaked features of the hole. Hence, the binding curve is much improved over HF and LDA (see Fig. 7).

## V Conclusion

The aim of the present work is to construct a functional for the exchange-correlation (xc) energy of DFT which applies to such diverse systems as the electron gas and the dissociation of the hydrogen molecule (H). To this end we try to find a model for the xc hole of these and intermediate systems in real space. It has been known for long that the exchange (x) hole captures many important features of the full one. For instance, in atoms term energies are often well described by reducing the exchange effects by and in the electron gas correlations have the effect of reducing the range of the pure x hole from a decay to a decay Gori-Giorgi et al. (2000a, b), being the distance from the the center of the hole. Thus unlike previous models that have sought to model the xc hole in real space, our present model aims at modifying the full x hole. Our investigations have shown that in the case of the dissociation of H the xc hole is qualitatively different from the x hole. Therefore, it is an unwise strategy to add a correlation (c) hole to a full x hole. The former would have to replace the latter with a full xc hole with appropriate features. We show here that this can be achieved in a more natural manner by multiplying the x hole by an appropriate correlation factor. Judging from the electron gas, the correlation factor might be chosen as a function of the distance between an electron and a reference electron. (We here again remind the reader that the xc hole describes the depletion of negative charge around an electron known to be at the reference position .) Unfortunately, our investigations have shown that such a simple correlation factor will have difficulties in moving half an x hole from one atom to the other—which is the appropriate effect of correlation in H at large nuclear separation. Instead, we have chosen to include in our correlation factor an additional factor depending on only one coordinate—a modulation of the depth of the xc hole. Such a factor is by symmetry just a constant in the electron gas and irrelevant to the shape of the hole in that case. Consequently, we have no guidance from the gas in choosing the -factor. Instead, we have decided to determine this factor at every point in space from the sum-rule for the xc hole. This sum-rule is of course of utmost importance for obtaining an xc potential with the correct asymptotics outside finite systems () from our model. This important property also requires the full symmetry in and in the density multiplied xc hole, , something that we have emphasized throughout the paper.

The decision to use the -factor for satisfying the sum-rule leaves us with great freedom in choosing a screening factor depending only on . Thinking about the electron gas it is natural to let this screening function have a screening length determined by an effective density according to . In the present work we have made the somewhat arbitrary choice . We are, however, aware of that we might be forced to abandon this simple choice in future refinements of our model.

For the actual shape of the screening function we have simply made a couple of reasonable choices for illustrational purposes. They are rapidly decaying, analytic functions with a few parameters with density dependence. One of the screening functions had its parameters specifically chosen such to reproduce the “exact” xc energies of the electron gas in the homogeneous limit, whereas two other screening functions had more ad-hoc parameters to illustrate the effect of selecting different forms of the screening function.

We could, however, also have attempted to find a screening function with a shape that would have allowed us to obtain accurate results for one- and two-electron systems. We would then have had an approximation which is able to properly dissociate H, which would be exact for the electron gas and quite accurate for the one- and two-electron cases. Such a functional is likely to give good total energies in a large number of systems.

In the present work we have, however, refrained from going down this road. Instead, our aim here was to present the basic ideas and to leave further refinements to future investigations. In order to just illustrate our new approach, we thus chose to present results obtained from two simple, but rather arbitrary screening functions given by the equations (14a) and (14b). We have seen that the shorter ranged choice () gives better results, but they are still not very good. It is however, seen from both Fig. 4 and Fig. 7 that the errors are mainly associated with an inadequate description of a simple one-electron system. (At a bond distance of 10 Bohr we basically have two separate hydrogen atoms.)

The most important result of the present work is that we managed to design a model which is able to describe the proper behavior of the xc hole of a hydrogen molecule as it dissociates. The details are not overly accurate, but we nurture real hope that they may fall in place by a better choice of the screening function. For now, however, we have left the search for such a function to future work.

## Appendix A The xc potential

The xc potential is obtained by straightforward differentiation of with respect to the density

 vSXxc(r)=−14∫dr1∫dr21|r1−r2|×(δ|γs(r1,r2)|2δn(r)A(r1)A(r2)h(r12,¯n)+2|γs(r1,r2)|2δA(r1)δn(r)A(r2)h(r12,¯n)+|γs(r1,r2)|2A(r1)A(r2)δh(r12,¯n)δn(r)). (15)

The functional derivative of the hole-depth function can be obtained by differentiating its hole sum-rule (11) with respect to the density. Its derivative can be worked out as

Unfortunately this equation does not give a closed expression for the functional derivative of , but just like its original counterpart (11) has to be solved iteratively to self-consistency.

Since the hole correctly integrates to electron (11), the potential has the correct asymptotic behavior . However, due to all the functional derivatives this is not directly visible in the expression for the xc potential (15). The main complication arises from the functional derivative of the KS density matrix whose evaluation requires the application of the chain-rule multiple times. Fortunately, in the case of a two-electron system, the Kohn–Sham density matrix can be expressed directly in the density as allowing for direct differentiation. Working out the first part of the potential gives

 vSXxc(r)=−A(r)2∫dr′n(r′)A(r2)h(r12,¯n)|r−r′|+⋯, (16)

so we find the correct Coulombic behavior. We only have to check the charge. Using in the sum-rule for the hole-depth function (11), we find

 A(r)2∫dr′n(r′)A(r′)h(|r−r′|,¯n)=1, (17)

so indeed, the Coulombic part of the potential also has the correct charge of .

## Appendix B The electron gas screening function

One of the requirements of the new functional is that it should work form a broad class of systems. Not only for molecules, but also for extended systems and in particular the homogeneous electron gas. Since the homogeneous electron gas is well studied and much of its properties are known, it provides the ideal system to fit the screening function (13) such that the SX functional will be exact for the homogeneous electron gas, in particular, it should give the exact xc energy for the electron gas. Before we can start fitting the screening function, we have to solve for the hole-depth, . For the KS density matrix of the homogeneous electron gas one can straightforwardly calculate that

 γs(r12,kF)=k3Fπ2sin(kFr12)−kFr12cos(kFr12)(kFr12)3,

where the Fermi wavevector . Using this expression for the KS density matrix in the sum-rule for the hole-depth function (11) and the Ansatz for the screening function (13), we find

 1=6πA2F4(~D), (18)

where we defined and the integrals

 Fn(β)≡∫∞0dy(sin(y)−ycos(y))2y−ne−βy.

More details about these functions and explicit expressions are given in Appendix D. Now we have an explicit relation how to calculate the hole-depth function for the homogeneous electron gas (18), we will fix by requiring our functional to produce the exact xc energy for the electron gas. In terms of our model, the xc energy density becomes

 εxc=−3kFπA2F5(~D). (19)

Since the equations (18) and (19) are linear in , the hole-depth function can simply be eliminated by dividing the equations, which gives

 F5(~D)F4(~D)=32πεxc(rs)εx(rs), (20)

where we used the Seitz radius and that for the homogeneous electron gas. To find , we need an explicit expression for for which we used an accurate fit to the random-phase approximation (RPA) and Green’s function Monte Carlo data by Perdew and Wang Perdew and Wang (1992). Unfortunately, the expression for cannot be inverted analytically. However, one can show that is a monotonic increasing function over , so at least the solution to (20) is unique. Further, in the low density limit we have

 limrs→∞32πεxc(rs)εx(rs)=32π(1+4π33√49πα1β4)=0.92925,

where the parameters and are from the low-density fit of Perdew and Wang Perdew and Wang (1992). Inverting (20) numerically, we find that in the low density limit

 D(rs)≈D∞rs,

where . For the asymptotic behavior in the high density limit, , we find for the ratio of the xc and x energies

 32πεxc(rs)εx(rs)=32π−23√49πc0rsln(rs)+⋯,

where is a constant from the high-density fit by Perdew and Wang to the RPA Perdew and Wang (1992). Similarly for small arguments of we have

 F5(~D)F4(~D)=32π−92π23√49πD(rs)rsln(rs)+⋯.

Comparing these two high density limits, we find that for high density to lowest order

 D(rs)≈D0≡49(1−ln(2)).

To obtain a workable expression for , we simply solved (20) numerically for several and made a least square fit using the following the following Padé approximant

 D(rs)=a0+a1rs+b3D∞r2s1+b1rs+b2r2s+b3r3s,

for which we found the coefficients , , , and . Note that we did not enforce the proper limit for , since the low density region is more important and relaxing this constraint significantly increased the accuracy for Bohr, which are more relevant densities in non-relativistic molecules and solids.

## Appendix C The prolate spheroidal grid

In this appendix we very briefly introduce the prolate spheroidal grid and give details on the exact parameters used for the grid in the calculations. The prolate spheroidal grid is created from an elliptic grid by rotating it around the axis connecting the two foci, the -axis. The elliptic coordinates are defined as

 ξ ≡r1+r22ρ and η ≡r1−r22ρ,

where is the distance to nucleus , which are placed at from the origin on the -axis. One readily sees that the ranges of the elliptic coordinates are and . The curves describe ellipses and the curves describe hyperbolae. The intersection points of the ellipses and hyperbolae for different constants will be used as grid points. Also taking the revolution around the -axis into account one can derive expressions for the cartesian coordinates in therms of the prolate spheroidal ones

 x =ρ√(ξ2−1)(1−η2)cos(ϕ), y =ρ√(ξ2−1)(1−η2)sin(ϕ), z =ρξη.

The main advantage in using ellipses and hyperbolae is that they are orthogonal to each other. Therefore, the metric is simply diagonal, , with diagonal elements

 gξξ =⟨∂r∂ξ∣∣∣∂r∂ξ⟩=ρ2ξ2−η2ξ2−1, gηη =⟨∂r∂η∣∣∣∂r∂η⟩=ρ2ξ2−η21−η2, gϕϕ =⟨∂r∂ϕ∣∣∣∂r∂ϕ⟩=ρ2(ξ2−1)(1−η2).

The scaling factors are defined as , from which we immediately find the volume element

 Ω=λξληλϕ=ρ3(ξ2−η2)

 ∇

We could also write down an explicit expression for the Laplacian, but we do not need it. The disadvantage of discretizing the Laplacian directly is that it is not symmetric anymore. Instead we use that for functions vanishing sufficiently fast at the boundary Becke (1982)

 −∫drf(r)∇2g(r)=∫dr∇f(r)⋅∇g(r),

so formally we can write

 −Ω(r)∇2=\raisebox6.45pt\makebox[−0.5pt][l]$↼$∇√Ω(r)⋅√Ω(r)\raisebox6.45pt\makebox[−0.3pt][l]$⇀$∇,

which is inherently symmetric and its discretization can directly be constructed from the discretized gradient.

The -grid points are generated by starting from an equidistant grid, , and using a simplified version of the transformation used by Becke Becke (1982) to generate grid points and weights for the -grid

 ξ(u)=1(1−u2)pξ.

The parameter can be used to affect the distribution of the grid points between the inner and outer region. For a given maximum value of , , is easily determined as

 pξ=−ln(ξmax)ln(Δu(2−Δu)),

where is the distance between the points in the -grid. An important aspect is the which has the benefit that the cusps at the nuclei are transformed into smooth Gaussians in -grid. The disadvantage is that the weight becomes zero at the line between the nuclei (), so the solution cannot be calculated directly at these points and has to be extrapolated.

For the -grid we started from an equidistant grid, , and use the following transformation to generate the -grid points

 η(v) =−cos(v+pηsin(2v)),

where the parameter can be used to modify the point distribution between the intermediate region and the nuclei. Becke found to be optimal Becke (1982). Also this transformation has the property that the cusps at the nuclei are transformed into Gaussians in the -grid. As one might already expect, the disadvantage of this feature is that the weights at the boundaries vanish (), so the results will have to be extrapolated also to these points. The points correspond along the bond axis outside the molecule.

We found that 80, 81 and 40 points in the , and -directions respectively gave numerically sufficient converged results. For the self-consistent solution of the hole-depth equation from its sum-rule (11) it was important that the density did not become too small. Hence the grid should not have points far in the asymptotic region. Using was sufficient to prevent points with numerically vanishing density, though still including a sufficient part of the relevant space.

In the equidistant grids it is straightforward to define numerical derivatives. For example one can use cubic B-splines Becke (1982) or simple central finite differences Davstad (1990) which we used in our calculations. In practical calculations, 4th order derivatives already gave sufficient accuracy.

## Appendix D The functions Fn

In this appendix we show how the integrals defining the functions can be evaluated. The functions satisfy

 dFn+1(β)dβ=−Fn(β). (21)

Using the boundary condition , we find

 Fn+1(β)=∫∞βduFn(u), (22)

which can be used to obtain successively higher order functions. The integral most suitable for direct evaluation is , albeit the evaluation is still rather tedious. The final result is

 F0(β)=165β2+4β3(β2+4)3. (23)

By applying successively the integral formula (22) we can obtain the higher order integrals required for our model. Their evaluation is straightforward, but takes the necessary amount of time. The final results are

 F1(β) =ln(1+4/β2)4+12β2−4(β2+4)2−3/2β2+4, F2(β) =14arctan(β/2)+12ββ2+4, F3(β) =(β2+2)ln(1+4/β2)8−12, F4(β) =arctan(2/β)3−(β3+6β)ln(1+4/β2)24+β6, F5(β) =14+β2(β2+12)ln(1+4/β2)96−β224−βarctan(2/β)3, F6(β) =(5β2+4)arctan(2/β)30−β3(β2+20)ln(1+4/β2)480+β3120−11β60.

The integrals for do not converge anymore, so is actually the last function in this series.

###### Acknowledgements.
K.J.H.G. would like to thank N.E. Dahlen for writing the code to calculate the Slater integrals needed in the full CI program. K.J.H.G. and R.v.L. acknowledge the Academy of Finland for research funding under Grant No. 127739.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters