Optimal power series expansions of the Kohn–Sham potential

A fundamental weakness of density functional theory (DFT) is the difficulty in making systematic improvements to approximations for the exchange and correlation functionals. In this paper, we follow a wave-function-based approach [N.I. Gidopoulos, Phys. Rev. A 83, 040502 (2011)] to develop perturbative expansions of the Kohn–Sham (KS) potential. Our method is not impeded by the problem of variational collapse of the second-order correlation energy functional. Arguing physically that a small magnitude of the correlation energy implies weak perturbation and hence fast convergence of the perturbative expansion for the interacting state and for the KS potential, we discuss several choices for the zeroth-order Hamiltonian in such expansions. Our first two choices yield KS potentials containing only Hartree and exchange terms: the exchange-only optimized effective potential (xOEP), also known as the exact-exchange potential (EXX), and the Local Fock exchange (LFX) potential. Finally, we choose the zeroth order Hamiltonian that corresponds to minimum magnitude of the second order correlation energy, aiming to obtain at first order the most accurate approximation for the KS potential with Hartree, exchange and correlation character.


I. INTRODUCTION
Electronic structure calculations are becoming indispensable in many areas of modern science, with applications spanning fields from drug discovery [1] to superconductivity [2].This change has been largely driven by the continuing development of density functional theory (DFT) over 50 years, which enjoys extraordinary and growing popularity [3]; and by the robustness of modern computational codes combined with the increasing speed of modern computers.
The importance of DFT in the theory of electronic structure was reflected in the 1998 divided Nobel Prize in Chemistry between W. Kohn [4] and J. Pople [5] for their developments of DFT and of computational methods in quantum chemistry respectively.The shared prize also reflected the importance of the synthesis of the two theories, since on their own both are limited, either by the difficulty to improve systematically on the approximations (DFT), or by poor scaling (wave function theory -WFT).In order to overcome these limitations, and to satisfy the growing demands for more accurate electronic structure calculations, on larger and more complicated systems, it is important to gain new insights.Such insights can be obtained from the integration of DFT with WFT [6][7][8][9][10][11][12][13].
In the scientific community, a dichotomy is perceived between DFT and WFT.The different emphasis of the two theories, density vs wave function, appears to hinder their smooth integration.The constrained search formulation by Levy [14] and by Lieb [15] and the adiabatic connection path construction [16][17][18] are seminal works in this area.Put together, they enabled Görling and Levy to formulate DFT perturbation theory (PT) [19,20], and Bartlett and co-workers to develop ab initio DFT [8].In these approaches, the correlation energy is approximated from second-order PT (or higher) and the KS potential is then determined using the optimized effective potential (OEP) method [21,22].
However, as the correlation energy from second-order PT is unbound from below, any minimization of the subsequent total energy functional is variationally unstable, tending to yield unphysically low total energies [10,23].Of course, there are physical situations (such as molecular dissociation) where the second order correlation energy term from a single reference Slater determinant will necessarily diverge to negative infinity; but the point is that using a correlation energy functional from second order PT, the tendency to diverge is inherent for all systems.In practice this divergence indeed turns out to be far more common than in quantum chemical methods employing perturbation theory, such as second order Møller-Plesset PT (MP2) [23,24].Whilst ways to alleviate the variational collapse have been put forward, for example by using Fock exchange energies instead of the true KS orbital energies [9,11], the absence of a rigorous solution to this issue has hindered progress in the cross-fertilization between WFT and DFT.Nevertheless, many-body perturbation theory (MBPT) has been employed successfully to yield an accurate correlation energy functional for DFT in the random phase approximation (RPA) [25][26][27][28][29][30][31][32][33], by combining the adiabatic connection path construction with the fluctuation dissipation theorem [34], but also using the Sham-Schlüter equation [35,36].
A few years ago, Gidopoulos proposed a natural way to integrate DFT and WFT, by constructing a pure WFT method whose solution happens to be the Kohn-Sham (KS) system of DFT [37].In that theory, it is no longer necessary to fix the electron density along the adiabatic connection path or elsewhere, making it a straightforward task to employ techniques from WFT in order to determine key quantities of DFT.The aim of the present paper is to demonstrate how the new formalism works, by constructing perturbative expansions of the KS potential that can be expected to converge optimally.
The paper is structured as follows.In section 2, we review the WFT method from [37], which is based on the minimization of the energy difference given by, where Ψ is the ground-state (g.s.) of the physical (interacting) system, and H v is an effective Hamiltonian, for some local potential v(r), which simulates the electron-electron repulsion.The g.s. of H v is Φ v and the g.s.energy is E v , The energy difference T Ψ [v] is strictly positive due to the Rayleigh-Ritz inequality; the positivity of the energy difference is preserved even when it is expanded with PT and an approximation up to second order is kept.Hence when T Ψ [v] is minimized there is no possibility of incurring the variational collapse of DFT with a correlation energy functional from second order PT.The relation between Eq. ( 1) and Lieb's functional [38] is explored in Ref. [39].
In section 3, we see how the the optimization over the total energy in the traditional OEP manner is equivalent to optimizing over the magnitude of the correlation energy.We then compare in section 4 our method with the traditional DFT perturbation theory (DFT PT) approach.In section 5, we discuss three different expansions for the KS potential: the first two yield at first order the exact exchange and local Fock exchange potentials respectively, and the final one already at first order includes correlation and has not been considered in the literature so far.Finally, we draw conclusions in section 6.

II. POWER SERIES EXPANSIONS OF THE KS POTENTIAL
In this section, we review the key results from the WFT approach developed by Gidopoulos in [37]; namely, how minimization of the energy difference in Eq. (1) yields the KS potential, and how to derive power series expansions of the KS potential from perturbation theory.
Ineq. (1) holds because the interacting state Ψ cannot be the exact g.s. of a non-interacting Hamiltonian H v ; however, we can view Ψ as an approximate g.s. of H v .Then, choosing v(r) to minimize T Ψ [v] amounts to selecting the Hamiltonian H v in the class (2) which optimally adopts Ψ as its approximate ground state.It transpires that the minimizing potential v s of T Ψ [v] is the KS potential, since setting the functional derivative of where ρ Ψ is the density of Ψ and ρ s is the density of v s .By the definition of the KS potential and the Hohenberg-Kohn theorem, the potential v s must be the KS potential (a detailed proof can be found in [37]).
With the variational principle (1), the problem of constructing a power series expansion of the KS potential is simplified, as it is no longer necessary to employ the adiabatic connection path formalism, where the local potential varies in an unknown manner along the path.Instead, we may substitute any power series expansion of Ψ in T Ψ [v], truncating the energy difference T Ψ [v] at a finite order.Optimization over v for a given expansion of T Ψ [v] then yields a corresponding expansion for the KS potential.
Of course, for a specific power series expansion of Ψ, it was always possible to truncate the expansion at any order and thus obtain its density; numerically inverting the density then leads to a (numerical) power series expansion of the KS potential.The difference with the present theory is that this procedure can be formally carried out for a whole class of Taylor series expansions of Ψ, characterized by the choice of zeroth-order Hamiltonian.It is then possible to consider the corresponding class of Taylor series expansions of the KS potential and search in that class for those expansions that converge faster than others.In other words, our method allows us to construct and then search a wide space of power-series expansions for the KS potential, to find those expansions which are expected to be the most accurate when truncated at some finite order.
In the following, we review from [37] the way to construct the lowest order in such expansions.In order to expand the energy difference, we use the interacting state Ψ u (α), g.s. of the perturbative Hamiltonian H u (α): The zeroth-order Hamiltonian is H u ; it belongs to the class of effective Hamiltonians (2) but with a local potential v en (r) + u(r) instead of v en (r) + v(r).Similarly to v(r), the effective potential u(r) mimics the electronic repulsion in a mean-field way.The fully interacting Hamiltonian H is obtained for α = 1, H u (1) = H.Obviously, for α = 0, Ψ u (0) = Φ u .If we substitute Ψ u (0) in place of Ψ in T Ψ [v] and search for the potential that minimizes T Ψu(0) [v], the minimizing potential will be v = u obviously.Hence, for small α, we expect that the potential which minimizes T Ψu(α) [v] will be close to u. Setting v(r) = u(r) + αv ′ (r), (7) the leading term in the energy difference T Ψu(α) [u + αv ′ ] turns out to be of second order: where Φ u,n , E u,n are the n-th eigenstate and energy eigenvalue of the effective Hamiltonian H u .
The second-order energy difference T u [w] is a functional of both the potentials u and w, but for now we take u to be fixed and focus on its dependence on w.
In the following, we seek to minimize T u [u + v ′ ] over v ′ : this is equivalent to minimizing T u [w] over w, because w = u + v ′ and u is fixed.In [37] the same symbol v was used for the potential appearing as the argument of the functional T Ψ in (1) and for the argument of T u in (9).Here, we use different symbols v and w to avoid confusion.
The functional derivative of T u [w] with respect to w, at fixed u, is given by [53] δT J u (r) is the direct Coulomb (or Hartree) local potential operator and K u is the Coulomb exchange non-local operator.φ u,i and φ u,a are respectively occupied and unoccupied orbitals in the Slater determinant Φ u , with ǫ u,i and ǫ u,a their corresponding eigenvalues.The functional derivative in Eq. ( 10) represents a charge density with zero net charge, Optimization over w in Eq. ( 9), by setting the functional derivative (10) equal to zero (at fixed u), yields the first order KS potential.We denote by w 0 [u] the minimizing potential of T u [w] for fixed u, From ( 8), the first-order term v ′ [u] in the KS expansion can be obtained from The desired expansion of the KS potential to first order is (2,7) The exact KS potential does not depend on u, but when the expansion is truncated at a finite order, the KS potential up to that order will depend on u.Hence, we write v s [u] to denote the KS potential up to first-order, and v s to denote the exact KS potential.We also denote by Φ s [u] the g.s. of v s [u], i.e. the KS determinant of the first-order KS potential v s [u].
In the Taylor expansion of the KS potential (14), the zeroth-order term, v en (r) + u(r), is the same as the potential in H u .The first-order term in the expansion of v s is v ′ [u].We may construct as many expansions for the KS potential as there are choices for u, and more besides using an altogether different expansion for Ψ, such as Møller-Plesset.
It is interesting to note that, by setting w = u in the functional derivative (10), we retrieve the equation for the exchange-only OEP (xOEP), also known in the literature as (exchange-only) exact exchange potential (EXX).This particular choice of u will be discussed in more detail in section 5; for now, we see how it also arises from an alternative perspective.
The density ρ Ψu(α) (r) of the weakly interacting state Ψ u (α) is given by where ρ u (r) is the density of the zeroth-order state Φ u .The density ρ Ψu(α) (r) of the weakly interacting system differs from the zeroth-order density ρ u (r) by a charge density equal (up to first order) to the functional derivative (10), where the latter is evaluated at w = u.Therefore, the search for the zeroth order potential u for which the g.s.density does not change to first order yields the exchange-only OEP (xOEP), as observed by Bartlett and coworkers [24].Furthermore, the density ρ Ψu(α) is related to the density ρ u+αv ′ (r) as follows: Hence, the density ρ Ψu(α) of the weakly interacting state differs from the density ρ u+αv ′ (r) of the non-interacting state by a charge density equal (up to first order) to the functional derivative (10), where the latter is evaluated at w = u + v ′ .Therefore, these densities are equal if the potential w is equal to the minimizing potential w 0 [u] (13); this minimizing potential defines the KS potential v s [u] (14).In other words, for any u, the density of the KS state is equal to the density of the weaklyinteracting state (to first order), where Although there can be several u that yield a converging expansion for Ψ and for the KS potential v s , we want to find those u whose expansions converge faster than others.We investigate this in section III.

A. Relation with the Sham-Schlüter method
Before proceeding to section III, we make contact with MBPT and the formalism of Green's functions.In MBPT, the requirement by Kohn and Sham that the density of the auxiliary noninteracting (KS) system be equal to the density of the interacting system leads to the Sham-Schlüter equation [35,36,40], in which G(r, r ′ ; ω) and G s (r, r ′ ; ω) are respectively the one-particle Green's functions for the interacting and the noninteracting (KS) systems.Eq. ( 18) determines the approximate exchange and correlation (xc) potential v xc in terms of an approximate xc self-energy Σ xc (r, r ′ , ω).Following Engel and Dreizler [41], who derive the OEP equation for the xc potential from the Sham-Schlüter equation, we point out the relation between Eq. ( 18) and Eqs. ( 16) and (17).Using (16) and requiring that the densities of the noninteracting and interacting systems be equal up to first order, i.e. requiring the validity of ( 17), yields the OEP equation, which determines the first-order KS potential v ′ [u] ( 13).This equation for v ′ [u] is equivalent to the Sham-Schlüter equation ( 18) with v ′ [u] in place of v xc and the modified self-energy, Σ − u, in place of Σ xc .Of course, in our theory, we do not impose the validity of (17), since the equality of the two densities comes out naturally from the optimisation of the second-order energy difference T u [w] (9).

III. REFERENCE DETERMINANTS WITH MINIMUM CORRELATION ENERGY
Historically, the xOEP is found by a minimization of the total energy Φ v |H|Φ v , where the Slater determinant Φ v depends on the effective potential v(r) (see Eq. 2).Since the exact energy Ψ|H|Ψ does not depend on v, the minimization of the energy is equivalent to the minimization over v of the magnitude of the correlation energy from the reference Slater determinant Φ v , we have explicitly shown the dependence of the correlation energy on the interacting Hamiltonian H of the system and on v. Hence, another interpretation of the xOEP follows: Corollary.xOEP is that effective potential v(r) with weakest correlation energy from its ground state Φ v .
The implication is that if we want to treat the interacting Hamiltonian perturbatively to all orders, then the effective Hamiltonian with the xOEP potential is the best zeroth-order Hamiltonian, as the remaining correlation energy to be treated perturbatively is smallest.
Often, we are interested in the lowest orders of perturbative expansions either because we want to study the limit of weak interactions or because we can only access the lowest orders numerically.Hence, we consider the partially interacting system described by the perturbative Hamiltonian H u (α) in (6) where the zeroth-order potential u(r) is meant to be determined later on in an optimal way.We make the following statement for the weakly interacting system described by the Hamiltonian H u (α), in the limit α → 0 and for any u: In this statement, the KS potential v s [u] is given to first order and the lowest (dominant) order in the correlation energy is second.
Proof.The correlation energy for the partially interacting system using as reference the g.s.Φ v (see Eq. 3) of an effective local potential v(r) in the class of Hamiltonians ( 2) is: For fixed u, the potential that minimizes the magnitude of the correlation energy E c Hu(α) [v] over v is the same as the potential that minimizes the expectation value Φ v |H u (α)|Φ v over v, since E u (α) does not depend on v.This optimal effective potential is different in general from the xOEP, due to the dependence of the former on u and on α.
Let us expand the correlation energy (21) in powers of α and obtain the dominant term.Obviously, when α = 0, the potential v that minimizes the energy Φ v |H u (0)|Φ v (or minimizes the magnitude of E c Hu(0) [v]) is v = u.Hence, for small α, we substitute Eq. 7 in (21) and we expand the correlation energy to second order in α to obtain where T u [w] is given by (9).Up to second order in α, the correlation energy ( 23) is thus equal to minus the energy difference ( 8): The KS potential v s [u] in (14) is that potential which minimizes the energy difference, and hence the statement follows.
It follows that when we minimize T u [w] over w to obtain the first order KS potential v s [u], the resulting potential not only has the same density as Ψ u (α) to first order, but it also has the following unique properties among other effective local potentials: • it best adopts Ψ u (α) (to first order) as its own approximate ground state and • its KS ground state Φ s [u] has the lowest magnitude of correlation energy to second order.
Let us denote by E c u [w] the negative of the energy difference T u [w], is a second order correlation energy expression.It is useful to use this notation to represent the total energy of the weakly interacting systems described by H u (α) using three different references: the zeroth-order state Φ u , the perturbative state Φ u+αv ′ , and the KS determinant Φ s [u].Keeping up to second order, we have in the limit α → 0: where Φ s [u] is the ground state of the first order KS potential v s [u] in (14).In general, for a given u, the optimal potential v ′ [u] (13) does not vanish and therefore the first-order KS potential v s [u] is different to the zeroth-order potential u.
In the following, we shall determine u optimally by selecting the one that makes E c u u + v ′ [u] small.

IV. COMPARISON OF DFT PT AND THE PRESENT WFT A. Traditional DFT PT method
In traditional DFT PT the KS potential is obtained from a perturbative expansion of the total energy functional, thus they are of the same order.The first order term in the expansion of the total energy is the exact exchange energy functional, which yields through functional differentiation the exchange potential, the first order term in the expansion of KS potential.Similarly, the correlation energy functional (truncated at second order) yields the correlation potential, which is the second order term in the expansion of the KS potential.The familiar scheme is summarized below: Because the exact exchange energy cannot be written explicitly in terms of the density, its functional derivative (the exact exchange potential) cannot be obtained directly from the density but only after solving an integral equation (Fredholm equation of the first kind), known as the equation for the optimized effective potential method.Although we are solving an OEP equation, the exchange potential is still the functional derivative of the exchange energy functional w.r.t. the density.

B. Present WFT method
In the present WFT method, which happens to have the KS potential as its solution, the xc potential is not the functional derivative of the xc energy w.r.t. the density (since the various quantities are not density functionals) and cannot be obtained directly.Instead, minimization of the magnitude of the second-order correlation energy functional E c u [v], Eq. ( 25), yields the minimizing potential w 0 [u], which emulates the Hartree exchange and correlation potential (Hxc) for the KS system with density ρ s [u] (17).The sum v en + w 0 [u] gives the KS potential up to first order, Eq. ( 14), for α = 1.
The xc-potential term in v s [u] is obtained by subtracting the Hartree potential from the optimal potential w 0 [u].The scheme is summarized below: We emphasize again the conceptual shift between the two theories: in DFT PT, the KS potential is obtained by minimizing the total energy of the system, while in the present WFT method the KS potential is obtained by minimizing the energy difference T Ψ [v].To dominant order, the latter optimization amounts to minimizing the magnitude of the correlation energy from the KS determinant.

V. OPTIMAL CHOICES FOR u
In the following we shall explore some choices for approximations to the interacting state Ψ.Based on the expansion Ψ u (α) discussed so far, this amounts to making a suitable choice for the potential u.However, we are free to pick any Ψ which might be expected to yield an accurate approximation to the exact KS potential; in addition, we shall also consider a Møller-Plesset expansion for Ψ.In any case, since we shall only consider perturbative expansions for Ψ, we wish to find expansions for the KS potential v s which are expected to give accurate results when the expansion is truncated at the lowest (meaningful) order: first order for Ψ and v s , and second order for the correlation energy.
In traditional DFT PT, the first order KS potential is restricted to Hartree and exact exchange, in fact in DFT PT the Hartree and exchange potential is defined as the first order term in the expansion of the KS potential.We shall discuss two choices for H u for which the first order KS potential indeed corresponds to Hartree and exchange only.Finally, we shall introduce a third choice for H u , which is expected to yield a first order KS potential with accurate Hartree and exchange and correlation character.

A. Exchange optimized effective potential
We anticipate that a good choice for u is such that the magnitude of the second-order correlation energy (25,26) is small, but we shall not discuss here how to find the global minimum of T u [u].An energetically better choice will be investigated in section V C.However, we present below an alternative argument which allows us to pick a u for which T u [u] is small.That choice of u yields xOEP.
For all zero-order potentials u, it holds that: The inequality holds because the search for the minimum over the potential w includes the value w = u.Inequality (29) states that for any potential u, the magnitude of its correlation energy is always larger or at most equal to the minimum of T u [w].It follows that the potential u Hx , for which equality holds in (29), will have correlation energy with small magnitude (but not the smallest possible).Equality in (30) holds when the first-order term in the expansion of the KS potential vanishes, The potential u Hx is then determined by setting w = u Hx in Eq. ( 10) and finding the potential u Hx which makes this functional derivative vanish, i,a φ * uHx,a (r)φ uHx ,i (r) + c.c. = 0 (32) This is the well-known equation for the xOEP.Hence the KS potential is We note two differences between our method and DFT PT, which also yields the xOEP: (a) In DFT PT the xOEP is the functional derivative of the exchange energy functional, which appears as the first-order term in the DFT PT expansion of the xc energy functional.
(b) In DFT PT, the total energy that gives rise to xOEP is truncated to first order and includes exchange energy and no correlation energy.There is no way to pair xOEP with a correlation energy functional without self-consistently altering the exchange potential away from its exchange only character.
In the current WFT, the first order KS potential is always paired, naturally, with a second order correlation energy, even when the first-order potential is xOEP.Specifically, the correlation energy corresponding to xOEP is given by Finally, we remark that the xOEP can be obtained [36] from the Sham-Schlüter equation ( 18), when we keep only the exchange term in the self-energy and approximate the interacting Green's function G with G s (linear Sham-Schlüter equation).

B. Local Fock exchange potential
So far, we have approximated the interacting state Ψ with the partially interacting state Ψ u (α), and considered which local potentials u(r) will give accurate approximations to the KS potential.In the prior section, we saw how one particular choice of u yields the wellknown xOEP.However, we now consider an altogether different approximation to Ψ, the Møller-Plesset (MP) expansion Ψ MP .
We initially consider only the zeroth-order term in the MP expansion, the Hartree-Fock (HF) determinant Φ HF .Following the approach in [42], we search for the effective Hamiltonian H v , with local potential v(r) (Eq.2), which optimally adopts Φ HF as its ground-state.We therefore minimize the energy difference T HF [v], given by over v(r) to determine the optimal H v .The functional derivative of T HF [v] is equal to where ρ HF (r) is the HF density, i.e. the density of Φ HF .
The ground-state whose potential minimizes this energy difference thus has the same density as the HF determinant.Denoting this optimal potential as v MP0 , the local Fock-exchange (LFX) potential is defined as and the MP expansion of the KS potential is The local potential with the HF density has been considered previously in the literature as an accurate approximation to xOEP and EXX [43][44][45].Much like for Ψ u (α) and subsequent expansions of the KS potential v s [u], we can consider higher order terms in the MP expansion of Ψ MP , which give rise to MP expansions of the KS potential.However, from Brillouin's theorem [46], singly excited Slater determinants do not couple directly with their zeroth-order HF state, which means the density of Ψ MP does not change to first order.Therefore, the potential which optimizes the energy difference, where Ψ MP1 is the first-order MP state, is the same potential as that which minimizes the energy difference in Eq. (34).Including first-order corrections to the MP expansion thus leaves the density and hence the expansion of the KS potential unchanged up to first-order.This is entirely analogous to our derivation of the xOEP: the xOEP is that zero-order effective potential (u Hx ), such that when we switch on the Coulomb interaction, the g.s.density of the weakly interacting state does not change to first-order (15), and whose corresponding power series expansion for the KS potential also has vanishing firstorder correction (31).Both potentials (LFX and xOEP) have exchange character.In DFT PT, the xOEP is the exact exchange potential (EXX) as it is the functional derivative of the exchange energy functional with respect to the density.The LFX potential cannot be expressed exactly but only approximately [42] as the functional derivative of the exchange energy functional.
Similarly to the xOEP, the LFX potential as well can be obtained from the Sham-Schlüter equation ( 18) when we keep the exchange part of the self-energy and omit correlation [40].However, unlike xOEP, the linear-response approximation (the replacement of G by G s ) is not employed to determine the LFX potential, and hence, from the point of view of the Sham-Schlüter method, the xOEP is an approximation of the LFX potential.On the other hand, from the DFT point of view, the LFX potential is instead an approximation of the exact exchange potential since only the latter is the functional derivative of the exact exchange energy w.r.t. the density [42,45].
As discussed in [42], the LFX and xOEP potentials are mathematically distinct, but share many physical properties, and would thus be expected to yield similar results when exchange dominates over correlation.Indeed, this was demonstrated to be the case, and it was theorized that the difference in results between the two methods is likely to indicate the correlation strength for a given system.Although the two methods are very similar, one advantage of the LFX method is that the functional derivative ( 35) is easier to compute, as there is no need to calculate the KS orbital shifts [47,48].

C. First order exchange and correlation potential
We previously saw that making the magnitude of the correlation energy E c u [u] (26) small gave rise to the wellknown Hartree and exact exchange potential in the first order KS potential (33).Whilst it is interesting to reproduce this result with our method, we want to develop a new expression that will give accurate results for systems where correlation is important.
As mentioned, finding the absolute minimizing potential of E c u [u] ( 26) is mathematically complex.The argument which gave rise to the Hartree and exact exchange potential does not fully optimize E c u [u] , and thus the expansion of v s [u] is not expected to converge as fast as desired.Let us instead try to minimize the magnitude of the correlation energy In principle, this involves a coupled minimization over u and v ′ which is even more complicated than minimizing E c u [u] over u.However, in practice the two minimizations of E c u [u + v ′ ] can be approximately decoupled, simplifying significantly the minimization scheme.To proceed, we split T u [w] into two terms, with and The first term S u [w] is a sum is over singly excited determinants from Φ u , while the second term D[u] is a sum over doubly excited determinants.
The potential w appears only in In practice [49], we have found that for any reasonable u, the minimization of T u [w] over w reduces S u [w] to very small values, compared with D[u]: Therefore, the dominant term is D[u], and the minimum of the energy difference T u [w] over w is given by D[u] to a good approximation, We conclude that in order to pick the best u, so as to minimize the minimum T u w 0 [u] , it is sufficient to choose u to minimize the double-excitations term D[u].This optimal u 0 , with minimum D[u 0 ] will correspond to the best zeroth-order effective Hamiltonian H u0 for a perturbative expansion of the interacting Hamiltonian ( 6).This dominance of D[u] over S u [w 0 [u]] also reinforces that u = u Hx is not energetically the best choice of u, since D[u] is not optimised in any way by this choice and can be quite large.
To minimize the double-excitations term, we first need to derive the functional derivative of D[u], so we need to determine how D[u] changes due to a perturbation u → u + λδu.Given that the ground and excited state wavefunctions, Φ u and Φ u,n , as well as their respective energy levels, E u and E u,n , are affected by the perturbation, D[u + λδu] to first order is where the dependence on u is now assumed and Φ 0 labels the ground state.We use the notation To write D[u + λδu] explicitly to first order in λ, we multiply it by the denominator in Eq. ( 45) and then write both sides of the subsequent expression as a power series in λ.We then expand the squared term which yields the following expression for the r.h.s. of Eq. ( 44), We must now determine the perturbed states and energies.We begin with the perturbed state |Φ (1) δu,n ; from Rayleigh-Schrödinger perturbation theory, this is given by where δU .= i δu(r i ).Since |Φ n is a doubly excited state, we can write it in the form |Φ ab ij , where i, j denote occupied orbitals in the ground state and a, b denote unoccupied orbitals.The matrix element, |Φ δu,n , is evaluated using Slater-Cordon rules and is given by where k ∈ |Φ ab ij , and c ∈ |Φ ab ij .The possible combinations for the pair (k, c) are therefore (a, i); (a, j); (b, i); (b, j); (µ, i); (µ, j); (a, ν); (b, ν), (49) where µ = (i, j), |µ ∈ |Φ and ν = (a, b), |ν ∈ |Φ .Any other permissible combination of (k, c) represents a triple excitation which will vanish in the final expression.We now determine the state |Φ abc ijk based on these possible combinations.We write where ĉ † and ĉ are fermion creation and annihilation operators.Using the anticommutator properties of these operators, namely and the fact that we get the following possible combinations for the state |Φ abc ijk : We are now able to compute the matrix element Φ 0 Φ δu,n in Eq. ( 46) (relabelling µ as k and ν as c), where ∆ αβ = ǫ α − ǫ β .The matrix element Φ δu,0 Φ n is determined in a similar manner and is given by Finally, we compute the perturbed energy levels E (1) δu,n and E δu,0 and hence the difference We collate these terms to determine the r.h.s. of Eq. ( 46).Let us first consider what happens to the first four terms in each of Eqs. ( 54) and (55) in the context of Eq. ( 46).The contribution from the very first term in each expression is given by with ij||ab = ij|V ee |ab − ij|V ee |ba .The other three terms in Eqs. ( 54) and ( 55) which involve a single-orbital substitution can be evaluated in a similar manner, and by relabelling dummy indices it can be shown that each of terms is equal.The total contribution from these terms is therefore It is noted that several of the other terms in Eqs. ( 54) and ( 55) are duplicates of each other, which again can be seen by relabelling dummy indices.After expanding all the remaining terms in Eqs.(54), ( 55) and (57) in terms of KS orbitals, the functional derivative of the double excitations term is found to be equal to The above expression is equal to zero at the minimizing potential, u(r) = u 0 (r).This result is reminiscent of the derivative of the double-excitations part of the second-order correlation energy in traditional DFT PT.In Ref. [10], in which part of the KS potential is expanded in terms of a Gaussian basis set {g t (r)} with coefficients b σ t , the derivative of the doubly-excited correlation energy term with respect to b σ t can be expressed as with δD[u]/δu(r) given by Eq. (61).However, as previously stressed, in Ref. [10] and other works in DFT PT, the minimization is carried out over the total energy, which is unbound from below.We discuss at some length the issues with a total energy minimization using a second-order correlation energy functional in section V D.
We can further simplify Eq. ( 61) in a manner which is also beneficial if we want to employ the Unsöld approximation [50] (common energy denominator approximation) [50? ????, 51].We note that some terms contain a denominator of mixed sign, which yields less accurate results if we approximate the denominators with a constant.Consider the complex conjugate of the expression where in the last step we have just swapped the labels of the dummy indices b and c.This term plus its complex conjugate is therefore equal to where the denominator is now of fixed (positive) sign.We can perform a similar procedure for the term with denominator ∆ ki , which with its complex conjugate becomes Using Eqs. ( 64) and (65), we can rewrite Eq. (61) as If desired, it is now straightforward to use the Unsöld approximation [50] (set all denominators ∆ equal to a constant) and remove the summations over the unoccupied orbitals using the completeness relation.The minimizing potential u 0 is determined by setting the functional derivative (61, 66) to zero: Eq. ( 67) must be solved iteratively with an energy minimization algorithm such as steepest descent.At the n th iteration, the potential will be u (n) (r).Substituting the single-particle orbitals φ u,p of u (n) (r) into (66), we obtain δD[u]/δu(r) at u (n) .Using this functional derivative we correct the potential, u (n) → u (n+1) , so as to lower D[u].Finally, we iterate until the functional derivative (67) vanishes.
Once the optimal potential u 0 has been found, together with its single-particle orbitals φ u0,p and energies ǫ u0,p , we may proceed to determine the first-order KS potential by minimizing S u0 [w] over w, keeping u 0 fixed.
The minimizing potential w 13) is given by (for fixed u 0 ): Eq. ( 68) is a standard OEP equation for the potential v ′ [u 0 ] with the simplification that during the solution of the OEP equation the orbitals φ u0,p and their energies ǫ u0,p remain fixed and independent of v ′ [u 0 ].The first-order correction v ′ [u 0 ] does not vanish.Finally, the KS potential to first order is given by Hence, the total energy (75) will have a stationary point (but not a minimum) at a potential lying somewhere between u Hx and u 0 .That potential will be the Hxc potential of DFT PT.It is intriguing to investigate the relation of the latter potential with the Hxc potential of the present theory, u 0 + αv ′ [u 0 ] (69).
The minimization of the total energy (75) over u amounts to a balanced search to achieve two goals: to minimize the expectation value Φ u |H|Φ u (well behaved) and to maximize the second-order difference T u w 0 [u] .Although bound from below, T u w 0 [u] is not bound from above and the search will be biased towards the maximization of T u w 0 [u] .During the iterations the potential is expected to move away from the minimum of T u w 0 [u] .Hence, the second term on the r.h.s. of (83), which we had omitted based on ( 42), (43), can no longer be neglected as it is prone to diverge, similarly to the third term.
The f.d. of the total energy (76) has only two terms (84) because the f.d. of Φ u |H|Φ u cancels with part of the f.d. of S u [u] (80).Thus, fully self-consistently and without risk of variational collapse, the Hxc potential (solution of δE[ρ u ]/ δu(r) = 0) can be obtained by searching for the potential ũ (dependent on w) that minimizes the (positive) second-order quantity S u [w] + D[u] and then choosing w so that ũ = w.From (77) and (84), it is evident that an algorithm to minimize S u [w] + D[u] will effectively maximize rather than minimize the total energy (77).Even more strongly than the previous case, the minimization of the total energy (76) does little to lower the value of Φ u |H|Φ u (since the f.d. of this term cancels) while it leads to the divergence of S u [w] + D[u].

VI. SUMMARY AND DISCUSSION
The research reported in this paper builds on previous work at the interface between wave function theory (WFT) and Kohn-Sham (KS) density functional theory (DFT) [37].The link between WFT and KS-DFT, established in [37], is that among all non-interacting Hamiltonians H v with an effective potential v(r), the KS effective Hamiltonian adopts energetically optimally the interacting ground state as its approximate ground state.Specifically, the KS potential turns out to be optimal in that it minimizes an appropriate energy difference T Ψ [v] (1) over all effective potentials v(r).This energy difference depends on the interacting state Ψ and is strictly positive, T Ψ [v] > 0 (1).
There is a large number of partially interacting Hamiltonians H u (α) (6), with 0 ≤ α ≤ 1, that yield the interacting Hamiltonian of interest H for α = 1; they differ in the choice of effective potential u(r) appearing in the zero-order Hamiltonian H u (2).For any of these partially/weakly interacting systems of electrons, their ground state Ψ u (α) can be expanded in a power series in the small perturbation α V ee − i u(r i ) .When we replace Ψ in the energy difference T Ψ [v], with an expansion of any of the partially interacting ground states Ψ u (α), truncated at a finite order, we obtain a corresponding power series expansion of the energy difference.Minimizing order-by-order the expansion of the energy difference w.r.t. the effective potential v, we obtain a corresponding power series expansion of the KS potential in powers of α.
There are at least as many expansions of the KS potential in powers of α as there are choices for the zero-order potential u.For any of the weakly interacting ground states Ψ u (α), and for small α, the dominant term in the expansion of the energy difference T Ψ [v] is second order: T Ψu(α) [u + αv ′ ] = α 2 T u [w] + O(α 3 ), with w = u + v ′ and T u [w] > 0.
Minimization of the second order energy difference T u [w] over w(r) gives an expansion of the KS potential up to first order.The aim is to choose optimally the zero-order effective potential u in order to obtain fast converging expansions for Ψ u (α) and for the KS potential.
The link between WFT and KS-DFT is explored further in the present work: We consider the correlation energy E c H [v] of the interacting system, with noninteracting reference the ground state Φ v of the Hamiltonian with effective potential v.The potential that minimizes the magnitude of the correlation energy E c H [v] over all effective potentials v is xOEP.
When we expand the ground state energy of the partially interacting system in powers of α, we obtain a power series expansion of the correlation energy.We consider the correlation energy E c Hu(α) [u + αv ′ ] of the partially interacting system with reference to the ground state of the effective potential u(r)+αv ′ (r).In the weakly interacting limit, α → 0, the dominant term in the expansion of the correlation energy is second order and it is equal to minus the second-order energy difference: E c Hu(α) [u + αv ′ ] = −α 2 T u [w] + O(α 3 ), with w = u + v ′ , see Eq. 23.This is the first important result of the paper.
We recall that the optimization of the energy difference T u [w] over all effective potentials w (i.e. over all reference ground states of w) yields the KS potential up to firstorder.We conclude that, for any u, the ground state of the KS potential (up to first-order) is the optimum reference for the correlation energy, since the magnitude of the second order correlation energy is minimum for that reference.We extend this reasoning by seeking the effective potential u for which the correlation energy from the KS reference state E c u w 0 [u] = −T u w 0 [u] (already a quantity with minimum magnitude over w) also has small or minimum magnitude over the zero-order potential u.
Intuitively, small magnitude of correlation energy implies weak perturbation and hence fast convergence of the perturbative expansion for Ψ u (α) and for the KS potential.
We consider three choices for the zeroth-order potential u.In the first two, the density of the zero order state is equal to the density of the weakly interacting state, within first order.In both cases, the first order term in the expansion of the KS potential vanishes.These two choices yield the Hartree and exact exchange potential of DFT (xOEP) and the Hartree and LFX potential [42].
By minimizing the magnitude of the correlation energy over w and over u (our third choice) we hope to obtain the fastest converging power series expansion for Ψ u (α) and for the KS potential, with the latter having exchange and correlation character.Since our second order expressions are bound from below, their minimization is mathematically well posed.We claim then that we have derived for the first time well behaved equations determining in an ab initio manner the KS potential with Hartree, exchange and correlation character, in a power series expansion of the potential up to first order.

VII. ACKNOWLEDGEMENT
NIG acknowledges financial support by The Leverhulme Trust, through a Research Project Grant with number RPG-2016-005.