OOn connecting density functional approximations to theory
Andreas SavinLaboratoire de Chimie ThéoriqueCNRS and Sorbonne University4 place Jussieu, F-75252 Paris, [email protected]:02 Tuesday 10 th November, 2020 submitted for publication in:
Density Functional Theory
Eds. Eric Cancés, Lin Lin, Jianfeng Liu, and Gero FrieseckeSpringer series on Molecular Modeling and Simulation, Vol. 1
Abstract
Usually, density functional models are considered approximations to density func-tional theory, However, there is no systematic connection between the two, and this canmake us doubt about a linkage. This attitude can be further enforced by the vaguenessof the argumentation for using spin densities. Questioning the foundations of densityfunctional models leads to a search for alternative explanations. Seeing them as usingmodels for pair densities is one of them. Another is considering density functional ap-proximations as a way to extrapolate results obtained in a model system to those of acorresponding physical one.
Dedicated to Jean-Paul Malrieuon his 80th birthday
On approximations in DFT
Density functional theory is here. It has changed the way the computation of electronicsystems is seen by the scientific community. It has a sound theoretical foundation. How-1 a r X i v : . [ phy s i c s . c h e m - ph ] N ov ver, following exact theory is more complicated than solving the Schrödinger equation.Furthermore, it does not tell us how to produce systematically approximations. Usual ap-proximations are convenient and (to a large degree) successful, but how to improve them? Excuses
This is not a review. References are erratic and biased. Own publications dominate, notbecause they are more important, but because they are only given to complement argumen-tation.
Summary
After giving the notations, and reminding the Hohenberg-Kohn theorem, some practicalsolutions are recalled, such as the decomposition of the universal functional, in order tocomply with different physical requirements. It is argued that this does not necessarily solvethe problem. Refinements, such as using the spin density as a supplementary variable, arediscussed. It is argued that the need for these refinements may hide a different foundationfor the approximations. In order to introduce a “systematic” way to approach the physicalHamiltonian, model Hamiltonian are defined that via an adjustable parameter approach thephysical Hamiltonian. Finally, examples show that simple mathematical recipes provide aquality similar to that of density functional approximations.
We start with a Schrödinger equation: H Ψ = E Ψ (1)The wave function Ψ depends on the coordinates of the electrons r , r , . . . , r N and theirspins. We will be mainly concerned with ground state eigenvalues, E = E . We considerHamiltonians of the form H = T + V + W (2) T is the operator for the kinetic energy, T = − N (cid:88) i =1 ∇ i (3) V is a local one-particle potential V = N (cid:88) i =1 v ( r i ) (4) = (cid:90) R v ( r ) ˆ ρ ( r ) d r (5)2he density operator, ˆ ρ ( r ) , can be written using Dirac’s δ function, ˆ ρ ( r ) = N (cid:88) i =1 δ ( r − r i ) (6)its expectation value is the density ρ ( r ) = (cid:104) Ψ | ˆ ρ ( r ) | Ψ (cid:105) (7)Please notice that it integrates to N , N = (cid:90) R ρ ( r ) d r (8) W is a two-particle local potential, W = N (cid:88) i Most DFAs start with the so-called local density approximation (LDA). Within this model,a general functional G [ ρ ] is replaced by the ansatz ˜ G [ ρ ] = (cid:90) R g ( ρ ( r )) d r (23)The function g has to be defined in some way. Traditionally, it is fixed in the uniformelectron gas, a system with an infinite number of particles, and where ρ does not depend on5he position (cf. chapter []) . Typically g is either obtained analytically as a function of ρ , or computed for a series of values of ρ , and fitted to them satisfying asymptotic conditions.LDA has the important advantage of being (to a certain extent) size-consistent, i.e., satisfying E A...B = E A + E B (24)where E A...B is the system composed of two parts, A and B , at infinite separation, while E A and E B are the energies of these parts computed independently. For the violations ofsize-consistency by LDA, see, e.g., [21, 27]. Another, major, advantage is its computationalsimplicity (just a numerical integration to obtain ˜ G ), and its linear scaling with system size.Both result from the local character of g : if ρ can be decomposed into contributions fromtwo spatial parts, ρ ( r ) = { ρ A ( r ) for r ∈ Ω A ( r ) ρ B ( r ) for r ∈ Ω B ( r ) (25)so can be g ; ˜ G becomes the sum of the two contributions.LDA can be extended by making g depend on other local quantities such as derivatives ofthe density, giving generalized gradient approximations (GGAs),etc. (cf. chapter []) . Applying the LDA, eq. (23), to F [ ρ ] , eq. (22), does not provide the accuracy needed in mostelectronic structure calculations. The strategy chosen is to define some density functional F d [ ρ ] , and approximate only the remaining part, ¯ F d [ ρ ] = F [ ρ ] − F d [ ρ ] .In the following, some choices for the partitioning of F will be presented. In the classical limit, the electrostatic interaction is given by the nuclear repulsion, V nn = (cid:88) A,B ( >A ) Z A Z B | R A − R B | the interaction between the electron cloud and the positive charges of the nuclei, (cid:90) v ne ( r ) ρ ( r ) d r and the repulsion inside the electron clouds, the Hartree energy, E H [ ρ ] = 12 (cid:90) R (cid:90) R ρ ( r ) ρ ( r (cid:48) ) | r − r (cid:48) | d r d r (cid:48) (26)There is a balance between these contributions. For example, between distant neutral atomsthese compensate (there is no | R A − R B | − term in the limit | R A − R B | − → ∞ ). Thisbalance is destroyed if E H is approximated, e.g., by using LDA, eq. (23). An excess or deficit6f repulsion produces an unphysical repulsion, or attraction of neutral atoms. Furthermore,even if this balance is enforced by parametrization for a given system, it is not kept foranother, even closely related system (see, e.g., [29]). The solution to this problem wasalready proposed in the original Hohenberg-Kohn paper [12]: E H is treated exactly, andonly the remaining part approximated.Finding good models for F [ ρ ] − E H [ ρ ] is still an active field of research; there are alreadyapproximations that work well for classes of systems, but one does not have yet a universallyapplicable model. The Pauli principle is hidden in the wave function used for defining F [ ρ ] , eq. (22). A way toimpose it is to use a model system, with F [ ρ, w (cid:54) = v ee ] , where the Pauli principle is imposed,and use approximations for the remaining part. E = min Ψ (cid:0) (cid:104) Ψ | T + V ne + W | Ψ (cid:105) + ¯ E Hxc [ ρ Ψ , w ] (cid:1) (27)where the subscript Ψ indicates that ρ is obtained from this wave function, and ¯ E Hxc [ ρ, w ] = F [ ρ, v ee ] − F [ ρ, w ] (28)This expression is derived using equations (18),(19),(21). In general, one takes into accountthe remark made above about E H [ ρ ] , and defines ¯ E H [ ρ, w ] = 12 (cid:90) R (cid:90) R ρ ( r ) ρ ( r (cid:48) ) (cid:18) | r − r (cid:48) | − w ( r , r (cid:48) ) (cid:19) d r d r (cid:48) (29)The remaining part, ¯ E xc [ ρ, w ] = ¯ E Hxc [ ρ, w ] − ¯ E H [ ρ, w ] (30)is called exchange-correlation energy.With eq. (27) one is back to an unconstrained variation of a wave function that is chosen tobe anti-symmetric, thus satisfying the Pauli principle.The problem is made simpler by a proper choice of w . For the Kohn-Sham model, onechooses the simplest one, namely w = 0 .The Kohn-Sham model is usually presented as a modified Schrödinger equation that is ob-tained by the variation of Ψ in eq. (27), H ( w )Ψ( w ) = E ( w )Ψ( w ) (31)where H ( w ) = T + V ne + V Hxc [ ρ, w ] + W (32) V Hxc = N (cid:88) i =1 v Hxc ( r i ) (33) v Hxc ( r , w ) = δ ¯ E Hxc [ ρ, w ] δρ ( r ) (34)7lease notice that this step (with the extra problem of the existence of the functional deriva-tive) is not needed to obtain E . Furthermore, E ( w ) = E [ v ne + v Hxc , w, N ] , so that E = E ( w ) + ¯ E Hxc [ ρ , w ] + (cid:90) R ρ ( v ne ( r ) − v Hxc ( r , w )) d r (35)where ρ is a minimizing density. One can also use the minimizing model wave function, Ψ( w ) , and choose to approximate the correlation density functional ¯ E c [ ρ, w ] = F [ ρ, w ] − (cid:104) Ψ( w ) | T + W | Ψ( w ) (cid:105) (36) F Separating F into a part defined, F d , and a remainder to be approximated, ¯ F d , does notnecessarily guarantee that an approximation, like that given in eq. (23), works better.Separating the Hartree part, E H , analogously to what was done in eq. (30), removes aproblem, but introduces a new one. Take the limiting case of one-electron systems. Thereis no contribution of the interaction between electrons: E Hxc = 0 . Thus, calculating exactlythe Hartree part means that the remaining part has to cancel E H exactly. But obtainingapproximations for − E H is as difficult as obtaining them for E H , and this was considerednot to be reachable with approximations of LDA-type. This is known as the self-interaction problem.Another (not unrelated) problem is due to degeneracy. For example, this appears when weconsider two parts of the system far apart, and this even in the simplest molecules like H ,or H +2 when they are stretched (the internuclear distance goes to infinity). Then, somethingrelated to the Einstein-Podolsky-Rosen effect shows up: an infinitesimal perturbation canproduce a drastic change in the wave function, the density, etc., but not in the energy.Unfortunately, this gets in conflict with the general philosophy of constructing DFAs thatare aimed to produce significant changes in the energies for small changes in the density.One could imagine to detect degeneracy. However, the model systems, in special the mean-field models (such as Kohn-Sham) do not necessarily have ground states presenting the samedegeneracy as the physical system: while one can present some degeneracy, the other maynot. While the physical wave functions have the symmetry of the Hamiltonian, the modelwave function often breaks symmetry to reduce the energy. (Well-known is the breaking ofspin symmetry showing up when bonds are stretched.) The opposite can occur, too: Kohn-Sham system can produce degeneracy, while the latter does not show up when Coulombinteraction is present (see, e.g., Fig. 11 in [32]).Even more difficult is the case of near-degeneracy, i.e., when a small change in the parameterscharacterizing H can produce degeneracy. In this case, detecting degeneracy is not a trivial For the standard Kohn-Sham model, E ( w = 0) is a sum of orbital energies. By construction, the minimizing model Ψ( w ) gives an exact ground state density. Someother properties can be reproduced, too. Trivially, all the expectation values of local one-particle operators, as they need only the density to compute them. At first surprisingly, theexact ionization potential can be also obtained. However, this can be easily understood, asit can be related to the asymptotic decay of the density (see, e.g., [17], [6]).Often, quantities that are not proven to be reproduced exactly by the model system arenevertheless expected to be good approximations. However, there is the danger of over-stretching this analogy. For example, it is fashionable to judge DFAs by their ability toreproduce fundamental gaps (differences between the ionization potentials and the electronaffinities) from differences between orbital energies (of the lowest unoccupied and highestoccupied ones). This is wrong, however [22, 34]. Let us consider, for example, a system withzero electron affinity. For a neutral system the Kohn-Sham potential, v ne + v Hxc , eq. (32),decays at large distances as − /r (see, e.g.,[17]), we know that it supports Rydberg series.Thus, its gap (ionization potential) is necessarily larger than its first excitation energy. Infact, accurate Kohn-Sham orbital energy differences give good approximations to excitationenergies. Let us take the He atom as an example [33]. An extremely accurate Kohn-Shampotential can be obtained from an extremely accurate density. The Kohn-Sham one-particleHamiltonian lowest eigenvalue corresponds to the doubly occupied state (1s). However,higher eigenvalues exist. The next eigenvalue (2s) is ≈ . hartree about the lowest one. Itcan be compared to the excitation energies to the triplet and singlet ( ≈ . , and . hartree,respectively). However, the fundamental gap of the He atom is of ≈ . hartree. (Thiscomparison should not to be confused with potentials produced by DFAs, as LDA for E xc that generates a potential that does not support excited states, and has a ionization potentialof ≈ . hartree.) Thus, in general, a DFA that produces an orbital energy difference thatreproduces the exact fundamental gap can be expected not to be a good approximation tothe exact Kohn-Sham system.The preceding discussion leads to a slippery ground. Could it be that the Kohn-Shamapproximations are used because they produce convenient mean-field models? Could it bethat (for specific purposes) they may be better than the exact Kohn-Sham system would be? The quality of approximations improves considerably when spin densities are introduced,i.e., when the functional ˜ G is made to depend not only on ρ , but on its components, thespin-up, ρ ↑ ( r ) , and the spin-down ρ ↓ ( r ) densities, ρ ( r ) = ρ ↑ ( r ) + ρ ↓ ( r ) ρ that on the spin polarization ζ ( r ) = ρ ↑ ( r ) − ρ ↓ ( r ) ρ ( r ) (37)A justification is brought by the fact that the exchange term acts only for electrons of thesame spin, and that correlation is not the same for a pair of electrons of different spins asfor that between two electrons of the same spin (that are kept apart by the Pauli principle).An example of the importance of making the the functional depend on ρ ↑ and ρ ↓ is shownin fig. 1. According to the Hohenberg-Kohn theorem, neither the energy, nor the value of F , for the hydrogen atom should depend on ζ . However, for LDA where a dependence on ζ is introduced by adjusting the exchange-correlation of the spin-polarized uniform electrongas, there is a clear dependence on ζ , the best value being obtained when ζ = ± , i.e.,formaximal spin polarization. - - ζ =( ρ ↑ - ρ ↑ )/ ρ F [ ρ ↑ , ρ ↓ ] Figure 1: Dependence of the local density approximation of the functional F on the spinpolarization ζ , eq. (37) for the exact density of the hydrogen atom; the exact value of F is0.5 hartree (dotted line).In spite of contributing in a decisive way to the success of DFAs (most achievements inthermo-chemistry would be inexistent without using spin-densities), there is a problem: thetheoretical foundation of this approach has never been established. This affirmation shouldbe supported here by a few of several arguments. One hears that the spin-density shows upin a weak magnetic field, and wrongly assumes that1. a weak magnetic field should not affect the result,2. a linear magnetic field should be sufficient, because the field is weak,3. it is sufficient to take into account the interaction between the magnetic field and thespins (i.e., only a term B z S z ), 10. the magnetic fields used for spin-polarized systems are weak.The first point is wrong, because lifting degeneracy by a magnetic field can produce a differentground state. For example, putting the stretched H molecule in a weak uniform magneticfield produces a triplet ground state, while in absence of the magnetic field, it is a singlet. Thesecond point is wrong, because it ignores a general problem:“a small perturbation parameterdoes not mean a small perturbation” [24]. For the specific case we consider, we noticethat even a weak linear magnetic field stabilizes states with high angular momentum belowthe ground state in the absence of the magnetic field. The variational principle cannot beapplied, and the Hohenberg-Kohn theorem cannot be proven [28]. The third point is wrong,as we know from the elementary treatment of the Zeeman effect: the orbital momentum isas important as the spin, but if we introduce a dependence on it, we have a dependence onthe external potential, and this is not allowed for a universal density functional. Finally, theforth point is wrong, because in order to produce a spin-polarized electron gas (for densitiesof chemical interest, ρ ≈ / π , i.e., r s = 1 ) a strong electronic excitation is needed, and thiscan be produced only by a huge magnetic field (see Fig. 2). ζ B z (a.u.) r s =1 Figure 2: Strength of the magnetic field, B , needed to stabilize the the uniform electron gaswith polarization ζ with respect to the unpolarized electron gas with density ρ = 3 / π , i.e., r s = 1 . The strongest magnetic field ever produced on earth is indicated by a horizontaldotted line.There is, however, a different viewpoint: the spin-density stays for another quantity thatcan be related to the spin-density. It has been noticed long ago for unrestricted Hartree-Fock calculations [38] that spin-densities can be connected to on-top pair density, P ( r , r ) ,cf. eq. (12). Starting from ρ ( r ) = ρ ↑ ( r ) + ρ ↓ ( r ) (38) P ( r , r ) = ( ρ ↑ ( r ) + ρ ↓ ( r )) − (cid:0) ρ ↑ ( r ) + ρ ↓ ( r ) (cid:1) (single determinant) (39)11n alternative interpretation of the spin-density in DFT is obtained [2, 23]: | ρ ↑ − ρ ↓ | = (cid:112) ρ ( r ) − P ( r , r ) (single determinant) (40)Could it be that the theory behind DFAs is not DFT?Let us mention that a relationship can be found also between spin-densities and first-orderreduced density matrices ([36], [35]). The adiabatic connection was invoked in order to understand what a density approximationshould do [11, 15, 9, 39]. The basic idea is that one constructs a model Hamiltonian de-pending continuously on a parameter, H ( µ ) . The corresponding Schrödinger equation hasan eigenvalue E ( µ ) and an eigenfunction Ψ( µ ) . We require that for a certain value of thisparameter the model Hamiltonian becomes the physical one. Let us now choose µ = ∞ for it.Furthermore, we assume that the Hellmann-Feynman theorem (or first-order perturbationtheory) can be applied to this model: ddµ E ( µ ) = (cid:104) Ψ( µ ) | ∂ µ H ( µ ) | Ψ( µ ) (cid:105) (41)Suppose that the model system, say at µ , is accessible (for example, it is a Kohn-Shamcalculation). We want to know how to correct the model energy, E ( µ ) , to obtain E = E ( µ = ∞ ) . For the missing part, we use the notation, ¯ E ( µ ) : E = E ( µ ) + ¯ E ( µ ) (42)By integrating eq. (41) we get what is also called integrated Hellmann-Feynman formula [4] ¯ E ( µ ) = E − E ( µ ) = (cid:90) ∞ µ (cid:104) Ψ( µ ) | ∂ µ H ( µ ) | Ψ( µ ) (cid:105) dµ (43)If we consider (as above) that the model only changes V and W , we also write E − E ( µ ) = (cid:90) ∞ µ (cid:104) Ψ( µ ) | ∂ µ ( V ( µ ) + W ( µ )) | Ψ( µ ) (cid:105) dµ (44)In density functional theory, one furthermore assumes that one can choose V ( µ ) such thatthe density does not change with µ . Using eqs. (5),(7), and the convention used here that v → v ne when µ → ∞ , we can write E − (cid:104) Ψ( µ ) | T + W ( µ ) + V ne | Ψ( µ ) (cid:105) = (cid:90) ∞ µ (cid:104) Ψ( µ ) | ∂ µ W ( µ ) | Ψ( µ ) (cid:105) dµ (45)Using a relationship analogous to eq. (10), and eq. (12), (cid:104) Ψ( µ ) | ∂ µ W ( µ ) | Ψ( µ ) (cid:105) = 12 (cid:90) R (cid:90) R P µ ( r , r , µ ) ∂ µ w ( | r − r | , µ ) d r d r (46)12lease notice that as Ψ depends on µ , so does P . A comparison with eq. (27) (where thedependence on w is replaced by that on µ ) gives the correction to E ( µ ) : ¯ E Hxc ( µ ) = (cid:90) R d r (cid:90) ∞ µ dµ (cid:90) R d r P µ ( r , r , µ ) ∂ µ w ( | r − r | , µ ) (cid:124) (cid:123)(cid:122) (cid:125) e ( r ) (47)The integrand e ( r ) shows a superficial similarity with the function g appearing in LDA,eq. (23). However, unlike LDA, the connection with ρ is not evident.One can eliminate a known term from ¯ E Hxc , and correct correspondingly the r.h.s. Forexample, if we would like to have ¯ E xc , eq. (30), we eliminate the contribution of ¯ E H , bytaking the derivative w.r.t. µ in eq. (29), i.e., by subtracting ρ ( r ) ρ ( r (cid:48) ) from P on the r.h.s.of eq. (47). Starting from the eq. (47), one may ask whether one should not construct functionals of thepair density, P ( r , r (cid:48) ) , instead of one that depends on ρ ( r ) . One can first notice that the pairdensity, P ( r , r (cid:48) ) , yields, by integration over r (cid:48) , the density ρ ( r ) , up to a factor N − . Thealready mentioned relationship between spin-densities and the on-top pair density, eq. (40),presents itself as a further argument. However, the conditions to be imposed on P such thatit is a fermionic one are difficult, while those to be imposed on ρ are simple ( ρ should benon-negative, and integrate to N ).In fact, LDA can be seen as replacing, in each point of space r , P ( r , r ) in eq. (47) bythat obtained in the uniform electron gas with density ρ ( r ) (see, e.g., [9]). This idea canbe extended beyond LDA: many successful functionals were constructed starting from thisperspective (among them those developed by A.D. Becke, or J.P. Perdew and co-workers,see, e.g., [1]).Some people consider the random phase approximation (RPA) as a density functional model.It can also be seen as constructing a simplified form of P to be used in eq. (47) (see, e.g.,[5]).Recently, new approximations using the pair density showed up (see, e.g., [37]). Even if by miracle we had the exact Kohn-Sham determinant (and potential), we still wouldmiss information about the exact system (with Coulomb interaction). For example, we stillwould not have the exact energy. Unfortunately, the task of obtaining simple functionalscapable of dealing with cases when a single Slater determinant is not a good approximationis not solved. 13ometimes ensembles of Kohn-Sham states are discussed. A formula expressing the correla-tion energy in terms of weighted Kohn-Sham orbital energies exists [25]. However, we do notknow a simple expression for obtaining the weights, and it does not seem that they follow aBoltzmann distribution [30].A long experience in quantum chemistry shows that a single Slater determinant is often abad starting point for obtaining many properties such as the energy. There, it seems naturalto consider multi-reference methods, i.e., wave functions where more than one determinantdeserve a preferential treatment. The selection of determinants is an art, unless selectivemethods are used, such as CIPSI (configuration interaction by perturbation with multicon-figurational zeroth-order wave function selected by iterative process) [13]. In the following, aspecial way of generating a multi-determinant wave function will be discussed, namely usingsome (ideally) weak interaction operator W . Degenerate (and near-degenerate) states aredetected by such operators, and this automatically introduces more than one Slater deter-minant if needed. Using more complicated wave functions is a price to pay for getting formsthat make existing DFAs closer to a theoretically justifiable form. w Eq. (47) suggests that it may be more easy to obtain approximations for E hxc when w (cid:54) = 0 ,i.e., µ > . Indeed, if ∂ µ w is short-ranged, we can use some approximation of P ( r , r (cid:48) ) thatis valid only when r (cid:48) is close to r , and use an expansion around r . In the limit of zero-range ( δ -function) we obtain the on-top density P ( r , r ) that for a single Slater determinantproduces a connection to the spin-density (eq. (40)), i.e., a form that resembles LDA withspin-dependence. Furthermore, expanding P in r (cid:48) around r produces semi-local terms suchas density derivatives [7]. Finally, we can expect a better transferability between systemswhen electrons are close, justifying the transferability from other systems like the uniformelectron gas, in other words, expecting “universality”.Also, it seems advantageous to avoid using w that posses a singularity (like the Coulombinteraction), because this induces a strong dependence on the basis set used, a very slowconvergence to the exact results (cf. the difficulty of converging (cid:104) Ψ | δ ( r ) | Ψ (cid:105) with a finitebasis set [3]).A simple and computationally convenient form for w satisfying the requirements above isgiven by w ( r , µ ) = erf( µ | r | ) | r | (48)Its derivative is short-ranged, ∂ µ w ( r , µ ) = 2 √ π e − µ | r | (49)and, when µ is very large ∂ µ w ( r , µ ) → πµ δ ( r ) , for µ → ∞ (50)The interaction in eq. (48) that also has the properties:14 w → v ee when µ → ∞ ,• w = 0 when µ = 0 i.e., by changing µ it is possible to switch between the Kohn-Sham and the physical system.This allows considering this method systematically improvable, in the sense that increasing µ brings the model closer to the physical Hamiltonian.However, we do not know how far we have to get away from w = 0 to get reliable approxi-mations. This can be explored numerically. w > Below are results obtained with w given by Eq. (48) and the dependence on µ is analyzed.First, to construct a density functional approximations to ¯ E xc , Eq. (30), uniform electrongas calculations are used [26, 19]. Now, the LDA, Eq. (23) can be applied to ˜ G = ¯ E xc forany value of µ .The numerical results given below are for the 2-electron harmonium, a system with theHamiltonian H = T + (cid:88) i =1 ω r i + erf( µ | r − r | ) / | r − r | (51)The variables can be separated, and the solutions can be found for real values of µ and ω by solving numerically a one-dimensional differential equation (see, e.g., [14]). For ω = 1 / ,that is chosen below, analytical solutions are known for the non-interacting ( µ = 0 ) and thefully interacting ( µ = ∞ ) system.Fig. 3 shows the errors made for the harmonium as a function of the choice of the parameter µ . At µ = 0 , the error is that given by the usual LDA. It decreases steadily, and around µ = 0 . . . . a change of behavior occurs, quickly reaching chemical accuracy (1 kcal/mol ≈ w (cid:54) = 0 requires having more than a Slater determinant, the time required forcomputing the wave function rapidly increases with µ . However, as the convergence withthe basis set is faster for w having no singularity, the computational effort for obtainingthe wave function is smaller. Fig. 4 shows the error that can be achieved in a given time.Calculations were done first for spherically symmetric basis functions to saturation (s-limit).Next a new value was obtained for the p-limit ( l = 1 ), next for the d-limit, ( l = 2 ), etc. Forsuch a small system, there is no gain in computing the integrals. However, one can see thatone reaches much faster a high accuracy when µ = 1 than with µ = ∞ .For this system, choosing a value of µ between 0.5 and 1 seems to provide a good compromisebetween the supplementary effort needed to have w (cid:54) = 0 , and having a good density functionalapproximation. 15 .5 1.0 1.5 2.0 2.5 3.0 μ ( a.u. ) Δ E ( a.u. ) Figure 3: Errors made by the local density approximation for the ground state energy ofharmonium (Eq. (51)) as a function of the range separation parameter of the model, µ (Eq. (48)) Instead of using universal models for P in eq. (47), one can construct corrections for a given model Hamiltonian energies determined by some v and w , ¯ E = E [ v ne , v ee , N ] − E [ v, w, N ] ,Eq. (42). The role of the approximation is to correct for the difference between the energyof the exact and that of the model system. We explore whether standard techniques fromnumerical analysis could compete withe density functional approximations in estimatingthese corrections.Please notice that as v → v ne and w → v ee , the correction vanishes: ¯ E → . One can alsotry to improve the result by using a set of model Hamiltonians for which obtaining the modelenergy is simpler than finding E [ v ne , v ee , N ] .In a density functional context it is tempting to use v as given by some density functionalapproximation, or even to use the potential that yields the exact density ρ (to show theprinciple of the procedure). Below the simplest expression for the external potential ischosen, v = v ne . Of course, this brings the model system very far from the physical systemwhen the interaction w is weak : the errors of the model at w = 0 are a very important partof the total energy. For example, for the harmonium studied above, at µ = 0 the error is of0.5 hartree, as shown in Fig. 3.First, we analyze how the energy of the model system, E ( µ ) , approaches that of the Coulombsystem, i.e., how E ( µ ) approaches E ( µ = ∞ ) . From the large µ behavior of the interaction w of eq. (48) and eq. (50), we derive E = E ( µ ) + a − k µ − k + a − k − µ − k − + . . . (52)16 ■ ■●● ● ● ( s ) Δ E ( au ) Harmonium ω = / ■ μ = ● μ = ∞ Figure 4: Harmonium energy errors obtained by saturating the basis set with l = 0 , , , ,in a calculation with µ = 1 , blue, and for the Coulomb interaction, red.The coefficient a in the equation above is proportional to (cid:90) R P ( r , r , µ = ∞ ) d r The coefficient b is proportional to a , and given by the cusp condition, as Ψ( µ ) has toapproach Ψ = Ψ( µ = ∞ ) when µ gets large [8], k is equal to l + 2 , l being the power ofthe expansion of P ( r , r + u ) in | u | around zero. In particular, for a pair of singlet coupledelectrons (of anti-parallel spin), we have k = 2 and P ( r , r , µ ) = P ( r , r , µ = ∞ )(1 + 2 √ π µ − + . . . ) (53)yielding a − a − = 4 √ √ π We consider a Taylor series for large µ . First, we make a change of variable to x ( µ ) suchthat x monotonously approaches 0 as µ → ∞ , E ( x = 0) = E ( x ) − xE (cid:48) ( x ) + 12 x E (cid:48)(cid:48) ( x ) + . . . (54)Using the chain rule, we go back to the µ variable, E ( µ = ∞ ) = E ( µ ) − x ( µ ) E (cid:48) ( µ ) x (cid:48) ( µ ) + 12 x ( µ ) [ E (cid:48)(cid:48) ( µ ) − E (cid:48) ( µ ) x (cid:48)(cid:48) ( µ ) /x (cid:48) ( µ )] ( x (cid:48) ( µ )) − + . . . (55)17btaining the first derivative of E with respect to the energy is not expensive, because itdoes not require computing a new wave function. Of course, the cost increases with higherderivatives.A choice for change of variable that makes the expansion correct at large µ , is x = µ − , cf.Eq. (52). As we know the next term in this expansion, we can also choose x = µ − + κµ − .Not surprisingly, the latter choice is more reliable than the first (cf. Fig. 5). However, itis a surprise that the approximation works very well up to µ ≈ , while the expansion wasderived in the limit µ → ∞ .Comparing these results with Fig. 3, one notices that the range of models for which thedensity functional approximation works well is comparable to that for which the Taylorseries works well. - - - - - μ ( a.u. ) Δ E ( a . u . ) Taylor series - - × to ( ) , x = μ - to ( ) , x = μ - to ( ) , x = μ - + κ μ - to ( ) , x = μ - + κ μ - E ( μ ) Figure 5: Errors of different approximations for the ground state energy of harmonium: E ( µ ) (thin, black), and the Taylor series around µ to order 1 (blue), and to order 2 (red); the dashedcurves correspond to a transformation to x ( µ ) = µ − , the others to x ( µ ) = µ − + κmu − .The horizontal dotted lines indicate chemical accuracy ( ± kcal/mol). The inset shows azoom on the same curves. Most applications of density functional theory rely on the simplicity of using a single Slaterdeterminant. This chapter does not intent to discourage the traditional search of densityfunctionals. They are suuccesful in practice, and there still is romm for improvement. How-ever, using simple mathematical techniques as discussed in the preceding section allowsobtaining a good quality, and this is encouraging. There are many paths that could befollowed. One, of course, is to improve the mathematical techniques. Another is to changethe interaction w to a form for which the extrapolations considered here would work better.Finally, using density functional models in the connnection withy the extrapolation approachpresented here, although this is could be envisaged.18 Acknowledgment The comments of Eric Cancès on the first draft on the manuscript are gratefully acknowl-edged. References [1] A. D. Becke. Int. J. Quantum Chem. , 23:1915, 1983.[2] A. D. Becke, A. Savin, and H. Stoll. Theoret. Chim. Acta , 91:147, 1995.[3] E. R. Davidson. J. Chem. Phys. , 39:875, 1963.[4] Saul T. Epstein, Andrew C. Hurley, Robert E. Wyatt, and Robert G. Parr. J. Chem.Phys. , 47:1275, 1967.[5] Henk Eshuis, Jefferson Bates, and Filipp Furche. Electron correlation methods basedon the random phase approximation. Theor. Chem. Acc. , 131:1084, 2012.[6] S. Fournais, M. Hoffmann-Ostenof, Th. Hoffmann-Ostenhof, and Th. OstergaardSorensen. Analytic Structure of Many-Body Coulombic Wave Functions. Commun.Math. Phys. , 289:291, 2009.[7] P. M. W. Gill, R. D. Adamson, and J. A. Pople. Mol. Phys. , 88:1005, 1996.[8] P. Gori-Giorgi and A. Savin. Phys. Rev. A , 73:032506, 2006.[9] O. Gunnarsson and B. I. Lundqvist. Exchange and correlation in atoms, molecules, andsolids by the spin-density-functional formalism. Phys. Rev. B , 13:4274, 1976.[10] C. Gutlé and A. Savin. Phys. Rev. A , 75:032519, 2007.[11] J. Harris and R. O. Jones. The surface energy of a bounded electron gas-solid. J. Phys.F , 4:1170–1186, 1974.[12] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev. , 136:B 864, 1964.[13] B. Huron, J.-P. Malrieu, and P. Rancurel. J. Chem. Phys. , 58:5745, 1973.[14] J. Karwowski and L. Cyrnek. Ann. Phys. (Leipzig) , 13:181, 2004.[15] D. C. Langreth and J. P. Perdew. Solid State Commun. , 17:1425, 1975.[16] M. Levy. Proc. Natl. Acad. Sci. U.S.A. , 76:6062, 1979.[17] M. Levy, J. P. Perdew, and V. Sahni. Phys. Rev. A , 30:2745, 1984.[18] E. H. Lieb. Int. J. Quantum Chem. , 24:24, 1983.1919] Simone Paziani, Saverio Moroni, Paola Gori-Giorgi, and Giovanni B. Bachelet. Phys.Rev. B , 73:155111, 2006.[20] J. Percus. Int. J. Quantum Chem. , 13:89, 1978.[21] J. P. Perdew. In R. M. Dreizler and J. da Providencia, editors, Density FunctionalMethods in Physics , page 265. Plenum, New York, 1985.[22] J. P. Perdew and M. Levy. Phys. Rev. Lett. , 51:1884, 1983.[23] J. P. Perdew, A. Savin, and K. Burke. Escaping the symmetry dilemma through apair-density interpretation of spin-density functional theory. Phys. Rev. A , 51:4531,1995.[24] F. Rellich. Perturbation Theory of Eigenvalue Problems . Gordon and Breach, NewYork, 1969.[25] A. Savin. Phys. Rev. , 52:4531, 1995.[26] A. Savin. On degeneracy, near degeneracy and density functional theory. In J. M.Seminario, editor, Recent Developments of Modern Density Functional Theory , pages327–357. Elsevier, Amsterdam, 1996.[27] A. Savin. Chem. Phys. , 356:91, 2009.[28] A. Savin. Mol. Phys. , 115:13, 2017.[29] A. Savin and F. Colonna. J. Mol. Struct. (Theochem) , 501-502:39, 2000.[30] A. Savin and F. Colonna. J. Mol. Struct. (Theochem) , 527:121, 2000.[31] A. Savin, F. Colonna, and M. Allavena. J. Chem. Phys. , 115:6827, 2001.[32] A. Savin, F. Colonna, and R. Pollet. Adiabatic connection approach to density func-tional theory of electronic systems. Int. J. Quantum Chem. , 93:166, 2003.[33] A. Savin, C. J. Umrigar, and X. Gonze. Chem. Phys. Lett. , 288:391, 1998.[34] L. J. Sham and M. Schlüter. Phys. Rev. Letters , 51:1888, 1983.[35] V.N. Staroverov and E.R. Davidson. Chem. Phys. Lett. , 330:161, 2000.[36] K. Takatsuka, T. Fueno, and K. Yamaguchi. Theor. Chim. Acta , 48:175, 1978.[37] L. Wilbraham, P. Verma, D.G. Truhlar, L.Gagliardi, and I. Ciofini. J. Phys. Chem.Lett. , 8:2026, 2017.[38] K. Yamaguchi and T. Fueno. Chem. Phys. , 19:35, 1977.[39] W. Yang.