Approximate energy functionals for one-body reduced density matrix functional theory from many-body perturbation theory
Klaas J. H. Giesbertz, Anna-Maija Uimonen, Robert van Leeuwen
EEPJ manuscript No. (will be inserted by the editor)
Approximate energy functionals for one-body reduced densitymatrix functional theory from many-body perturbation theory.
Klaas J. H. Giesbertz , Anna–Maija Uimonen , and Robert van Leeuwen Theoretical Chemistry, Faculty of Exact Sciences, VU University, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands,e-mail: [email protected] Department of Physics, Nanoscience Center, University of Jyv¨askyl¨a, P.O. Box 35, 40014 Jyv¨askyl¨a, Survontie 9, Jyv¨askyl¨a,Finland, e-mail: [email protected]
July 4, 2018
Abstract
We develop a systematic approach to construct energy functionals of the one-particle reduceddensity matrix (1RDM) for equilibrium systems at finite temperature. The starting point of our formu-lation is the grand potential Ω [ G ] regarded as variational functional of the Green’s function G based ondiagrammatic many-body perturbation theory and for which we consider either the Klein or Luttinger–Ward form. By restricting the input Green’s function to be one-to-one related to a set on one-particlereduced density matrices (1RDM) this functional becomes a functional of the 1RDM. To establish theone-to-one mapping we use that, at any finite temperature and for a given 1RDM γ in a finite basis, thereexists a non-interacting system with a spatially non-local potential v [ γ ] which reproduces the given 1RDM.The corresponding set of non-interacting Green’s functions defines the variational domain of the functional Ω . In the zero temperature limit we obtain an energy functional E [ γ ] which by minimisation yields anapproximate ground state 1RDM and energy. As an application of the formalism we use the Klein andLuttinger–Ward functionals in the GW-approximation compute the binding curve of a model hydrogenmolecule using an extended Hubbard Hamiltonian. We compare further to the case in which we evaluatethe functionals on a Hartree–Fock and a Kohn–Sham Green’s function. We find that the Luttinger–Wardversion of the functionals performs the best and is able to reproduce energies close to the GW energywhich corresponds to the stationary point. One-body reduced density matrix (1RDM) functional the-ory provides an interesting alternative to density functiontheory (DFT), which holds the promise to alleviate manyof the current practical failures of DFT. Especially theinherent better capability of 1RDM functional theory todeal with strongly correlated systems, e.g. the breaking ofchemical bonds [1–4], Mott insulator transitions [5, 6] andthe fractional quantum Hall effect [7], are an encourage-ment to invest even further in the development of 1RDMfunctionals.Currently, two main strategies are followed in the de-velopment of new 1RDM functionals. The first strategyuses the fact that the exact 1RDM functional for two-electron systems is available which is used as a paradigmto device general N -electron functionals [1, 2, 8–10]. Theother main approach is to reconstruct the two-body re-duced density matrix (2RDM) from the 1RDM guided by N -representability conditions [3, 11–17]. A more extensiveoverview of different approaches can be found in Ref. [18]. Correspondence to : K.J.H. Giesbertz
A disadvantage of the current approaches is that theydo not provide a systematic route towards improved 1RDMfunctionals. In many-body perturbation theory (MBPT),however, such a route is available by including an increas-ing amount of terms in the perturbation expansion [19,20]. These terms in the perturbation expansion are conve-niently represented by Feynman diagrams. It would there-fore advantages to connect the MBPT framework to 1RDMfunctional theory, since this would allow one to use theMBPT functionals for the construction of increasingly moreaccurate 1RDM functionals in a systematic manner. Thisapproach has already been heavily pursued in DFT withsome considerable success [21–30], so the additional flex-ibility of 1RDM functional theory to a further improve-ment of the results, especially when strong correlationsplay an important role. Due to the peculiar propertiesof 1RDM functional theory, however, the connection toMBPT is not so trivial as for DFT. In DFT one canuse the Kohn–Sham (KS) system [31], from which a non-interacting Green’s function can be extracted to be in-serted into the MBPT functional. The construction ofa similar KS system in 1RDM functional theory — anon-interacting system with the same 1RDM as the in-teracting system — leads to a completely degenerate or- a r X i v : . [ c ond - m a t . o t h e r] J u l K.J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. bital spectrum, since the orbitals are fractionally occu-pied [32–34]. For Coulomb systems we know at least thatinfinitely many orbitals are fractionally occupied [35] andthere is strong evidence that they all are [36–38]. The cor-responding Green’s function of the non-interacting systemin 1RDM functional theory has therefore all its poles lo-cated at the same position. Such an unphysical Green’sfunction is bound to lead to nonsensical results when in-serted directly into an MBPT functional.The degeneracy of the orbital energies in 1RDM func-tional theory seems to pose a serious problem if one wantsto use MBPT for the construction of functionals. How-ever, the massive degeneracy of the orbital energies iscaused by the fact that we evaluate the system at zerotemperature and can be lifted by elevating the tempera-ture to any finite value,
T >
0. The degeneracy problemcan therefore circumvented by performing the calculationat a finite temperature. Within a finite basis set, 1RDMfunctional theory allows for a completely rigorous foun-dation in which all required functionals are nicely differ-entiable and no v -representability issues arise [39]. Thezero temperature limit, T →
0, can be taken at the endof the calculation. An additional advantage is that alsothe finite-temperature MBPT formalism is more reliable,since finite-temperature MBPT does not rely on adiabat-ically turning on the interaction, which leads to problemsfor zero-temperature MBPT in the case of level cross-ings [19, 20]. We will elaborate later in this article moreon this extraction process of the non-interacting Green’sfunction from the 1RDM.The extracted Green’s function can now be deployed inthe MBPT framework which can be done in several ways.The most straightforward option is to use the extractedGreen’s function as a zeroth order Green’s function inthe MBPT formalism. Unfortunately, this leads to compli-cated expressions: the effective potential needed to forcethe non-interacting system to have the same 1RDM as theinteracting system needs to be subtracted from all occur-rences of the self-energy [40–44]. This means that also theenergy expression depends on the effective potential it-self, so becomes a severe complication. More importantly,following this approach one only obtains a reformulationof MBPT in terms of a different (less convenient) zerothorder Green’s function. The results will therefore be thesame as for an MBPT calculation with the same diagrams,so the only achievement would be a more complicated for-malism, if one aims for fully self-consistent solutions.A more viable approach is to use variational MBPTfunctionals Ω [ G ] for the total energy [22, 23, 45–47]. Thetotal energy from these functionals are relatively stablewith respect to variations in the Green’s function, duetheir variational property δΩ/δ G = 0 when G is ob-tained self-consistently form the Dyson equation. This hasthe advantage that one does not need to solve the fullDyson equation, but that it is sufficient to use a rea-sonable approximate Green’s function to obtain a reason-ably accurate value for the total energy. Such an approxi-mate Green’s function can be the KS Green’s function or the non-interacting Green’s function extracted from the1RDM.We have different variational functionals at our dis-posal and popular forms are the Luttinger–Ward [48] andKlein [49] functionals. Especially the Klein functional isquite popular, as it yields the simplest expression for theenergy [22, 23]. The Klein functional can be reformulatedin terms of the exchange-correlation kernel to yield theadiabatic connection fluctuation-dissipation expression [26,27, 50, 51]. Though the Luttinger–Ward functional resultsin a more complicated expression for the energy than theKlein functional, it has the advantage of superior varia-tional properties. We will take both functionals in consid-eration in this article, but limit ourselves to one of thesimplest perturbative expansions: the GW. The perfor-mance of both the Klein and Luttinger–Ward functionalsin the GW approximation will be tested on an extensionsof the two-site Hubbard model which is capable to give agood representation of the ground state of the hydrogenmolecule. To achieve this, the extended two-site Hubbardmodel also takes the intersite interactions into accountand the interaction between the electrons and the nuclei.The paper is divided as follows. In Sec. 2 we presentthe general theory of the variational functionals that weuse. In Sec. 3 we show how to construct a non-interactingGreen’s function G s [ γ ] that yields a prescribed 1RDM γ .In Sec. 4 we present the details for the molecular modelthat we use to test our functionals. In Sec. 5 we describethe various input Green’s functions that we use as inputfor the variational functionals. In Sec. 6 we perform varia-tional calculations on the molecular model and finally, inSec. 7 we present our conclusions. The required theoretical framework for the use of varia-tional MBPT functionals has basically already been pre-sented in the framework of DFT in Ref. [22]. The only dif-ference is that in the framework, we now use non-local po-tentials to keep the complete 1RDM fixed at all interactionstrengths. The procedure is very analogous to the deriva-tion of the Luttinger–Ward functional [20, Sec. 11.4]. Forcompleteness, we will provide here a short overview. Con-sider theˆ H λ = (cid:88) ij (cid:0) k ij + v λij (cid:1) ˆ c † i ˆ c j + λ (cid:88) ijkl w ijkl ˆ c † i ˆ c † j ˆ c k ˆ c l , (1)where k ij and w ijkl are are one-and two-body matrix el-ements. The term v λ is a potential to keep the 1RDMidentical at all interaction strengths λ to the fully inter-acting 1RDM, γ λ = γ = γ . The potential at λ = 1is therefore the real external potential v = v ext . In theframework of KS-DFT the potential v λ would be local,i.e. diagonal in the coordinate representation or site basis,and only the density would be kept fixed at all interac-tion strengths, n λ = n . At λ = 0, the Hamiltonian be-comes non-interacting ˆ H s := ˆ H and the its correspondingGreen’s function G s has a particularly simple form (see .J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. 3 Sec. 3). The corresponding potential to keep the 1RDMor density of the non-interacting system identical to theone of the non-interacting system is denote by v s := v .At general λ , the Green’s function is related to thenon-interacting Green’s function via the Dyson equation G λ = G s + G s ˜ Σ λ G λ , (2)where ˜ Σ λ := Σ λ + v λ − v s . To relate the grand potentialsof the interacting and non-interacting system, we will usethe fundamental theorem of calculus. So we first work outout the derivative of the grand potential with respect tothe interaction strengthd Ω λ d λ = (cid:28) d ˆ H λ d λ (cid:29) = 1 β Tr (cid:26)(cid:18) d v λ d λ + Σ λ λ (cid:19) G λ (cid:27) , (3)where Tr indicates a summation over all indices as wellas a summation over the Matsubara frequencies. We nowuse the Φ -functional which can be constructed by sum-ming over all irreducible self-energy diagrams closed withan additional Green’s function [48]. Each of these closeddiagrams is multiplied by the pre-factor 1 / (2 n ), where n denotes the number of interaction lines Φ λ [ G ] = (cid:88) n,k n Tr (cid:8) Σ λ ( n ) k G (cid:9) . (4)The self-energy derived from the Φ functional as Σ [ G ] := δΦδ G , (5)leads to a conserving approximation for Σ [52, 53]. Itsderivative with respect to the coupling strength is readilyworked out asd Φ λ d λ = 12 λ Tr (cid:8) Σ λ G λ (cid:9) + Tr (cid:26) Σ λ d G λ d λ (cid:27) , (6)which allows us to rewrite the derivative of the grand po-tential as β d Ω λ d λ = dd λ (cid:18) Φ λ − Tr (cid:8) Σ λ G λ + ln (cid:0) − G s ˜ Σ λ (cid:1)(cid:9)(cid:19) . (7)The last term is readily verified by working out the deriva-tive of Tr (cid:8) ln (cid:0) − G s ˜ Σ λ (cid:1)(cid:9) [20]. Integrating from the non-interacting system to the fully interacting one, we findthat their grand potentials are related as Ω LW [ G ] = Ω s + tr (cid:8) ( v ext − v s ) γ (cid:9) (8)+ 1 β (cid:16) Φ [ G ] − Tr (cid:8) ˜ ΣG + ln (cid:0) − G s ˜ Σ (cid:1)(cid:9)(cid:17) , where ˜ Σ := ˜ Σ = Σ + v ext − v s (9)and tr only sum over matrix indices. By regarding thegrand potential as a functional of the one-body Green’s function we have retrieved the Luttinger–Ward (LW) func-tional Ω LW [ G ], for an arbitrary non-interacting referenceGreen’s function G s . Note that this expression assumesthat we are keeping the complete 1RDM constant for vary-ing interaction strength. If we were only to keep the den-sity fixed, γ should be replaced by n and the potential v ext and v s would be local. Note that apart from the al-lowance of an arbitrary spatially non-local potential in-stead of the true external potential to be able to workwith more general G s , the result in (8) is identical to thatobtained originally by Luttinger and Ward [48].When integrating over the interaction strength, we haveassumed that the chemical potential is also kept fixed atits interacting value. This condition fixes the constantin v λ , i.e. tr (cid:8) v λ (cid:9) , as it should be chosen such that theparticle number remains constant at varying interactionstrengths. This is readily guaranteed by assuring that theLuttinger–Ward functional yields correct number of par-ticles, N = − ∂Ω LW /∂µ . This is particularly convenient asthe term µN can now easily be eliminated to obtain thetotal energy.The Luttinger–Ward functional is variational in thesense that when G solves the Dyson equation (2) with Σ = δΦ/δ G , perturbations in Ω LW vanish to first order δΩ LW δ G = 0 . (10)This does not guarantee that the Luttinger–Ward func-tional achieves its minimum at the solution of the Dysonequation, but only that it is stationary.One can construct infinitely many different expressionsfor the grand potential with a different functional depen-dence of the Green’s function which both yield the correctvalue of the grand potential and are stationary at the so-lution of the Dyson equation [20]. A popular alternativevariational functional is the one due to Klein [49] Ω K [ G ] := Ω s + tr (cid:8) ( v ext − v s ) γ (cid:9) + 1 β (cid:104) Φ [ G ] + Tr (cid:8) ln( GG − s ) + − GG − s (cid:9)(cid:105) . (11)The Klein functional is simpler to evaluate, as there is noexplicit reference to the self-energy. The price we pay forthis simplification is that the Klein functional is less sta-ble to perturbations than the Luttinger–Ward functional.Though the first order variations vanish at the solutionof the Dyson equation, the higher order terms typicallyhave larger amplitudes. This higher sensitivity to the inputGreen’s function of the Klein functional has been demon-strated for atoms [22] and diatomic molecules [23]. Thesimplicity of the Klein functional is especially apparentif we insert the non-interacting Green’s function G s intothe Klein functional. All the terms within Tr {·} disappearand only the Φ -functional remains as a non-trivial part.Finally, we address some practical applications of thevariational functionals. First of all let us consider that thestationary point of the functional is attained for a Green’sfunction G but that we insert a different input Green’sfunction ˜ G . Due to the stationarity property (10) (and K.J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. similarly for the Klein functional) we have, with ∆ G =˜ G − G that Ω [ ˜ G ] = Ω [ G ]+ 12 (cid:90) δ Ωδ G δ G [ G ] ∆ G ∆ G + O (cid:0) ( ∆ G ) (cid:1) , (12)where the second order derivative in the integrand on theright hand side is a tensor contracted with terms ∆ G and we further imply a double integration of imaginarytimes. We therefore see that if we make an error ∆ G inthe Green’s function then the error in Ω is, due to itsvariational property, only quadratic in this error. This isa very useful property as it ensure that we may obtaingood values of the grand potential from rather simple in-put Green’s functions, provided they are close enough insome sense to the stationary point. The size of the erroralso depends on the value of second derivative in the inte-gral, which is different for the LW and Klein forms of theenergy functional and in practice the variational propertyof the LW form is found to be superior [22, 23]. Anotherway to use the energy functionals is by looking for thestationary point on a restricted domain. For example byinserting Green’s functions from a general non-interactingsystem with a local potential we obtain a density func-tional [22, 23] since the potential is in one-to-one corre-spondence with a density by means of the Hohenberg–Kohn (or Mermin) theorem. In this work we extend thisto non-local potentials since at finite temperature andin a finite basis it can be proven that any ensemble N-representable 1RDM is v -representable by a non-local po-tential [39]. Using this property, we can consider the vari-ational domain of Ω [ G ] to be the set of Green’s functionsof non-interacting systems at a finite temperature witha non-local potential. By varying over the non-local po-tentials we vary over a domain of 1RDMs and we are ef-fectively using a 1RDM functional the stationary pointof which yields an approximate 1RDM. This approachwill be used in the present work. Note that by restrict-ing the variational domain we will not recover the exactGreen’s function anymore at the stationary point of thevariational equations. This would even be the case hadwe used the exact Φ -functional. We further remark thata different way to construct density matrix functionals onthe basis of Green’s functions is given in Ref. [44] and theproposed approximate scheme in Section IV C of Ref. [44]coincides with our use of the Luttinger–Ward functional.The method that we present here is close in spirit to theone presented in Ref. [54]. G s [ γ ] In this Section we outline how we can construct a non-interacting system at finite temperature that produces agiven 1RDM. The Green’s function G s [ γ ] of this systemproduces functionals of the 1RDM when inserted in thevariational Luttinger–Ward and Klein functionals. These1RDM functional will be studied in greater detail in latersections. The Matsubara Green’s function of a non-interactingsystem consisting of bosonic/fermionic particles can beshown to be of the form [19, 20] G s,kl ( τ ) = − i δ kl (cid:2) θ ( τ ) ¯ f k ± θ ( − τ ) f k (cid:3) e − τ(cid:15) M k , (13)where f k := (e β(cid:15) M k ∓ − is the Bose/Fermi distributionfunction, ¯ f k := 1 ± f k and β := 1 /T is the inverse temper-ature. The Heaviside step function is defined as θ ( τ ) = (cid:40) τ >
00 if τ < . (14)The Matsubara energies are obtained by subtracting thechemical potential from the eigenvalues of the one-bodyHamiltonian, (cid:15) M k := (cid:15) k − µ , where the chemical poten-tial can be adjusted such that the system has the desirednumber of particles N = (cid:88) k f k = (cid:88) k β ( (cid:15) k − µ ) ∓ . (15)The corresponding non-interacting 1RDM is readily ex-tracted from the Green’s function as γ s,kl = ± i G s,kl (0 − ) = f k δ kl . (16)Since the 1RDM is diagonal in the eigenbasis of the one-body Hamiltonian, the Bose/Fermi distribution functionsare the natural occupation numbers.Now let us limit the discussion to electrons, i.e. fermions.It is clear from the form of the Fermi function that it sat-isfies the strict inequalities 0 < f k < T >
0. In the limit T →
0, however, the Fermi functioncollapses to a Heaviside function and the occupation num-bers become integers 0 or 1 depending on the sign of (cid:15) M k .For (cid:15) M k = 0 the value of the occupation number is unde-termined at T = 0 and can be anything between zero andone, 0 ≤ f k ≤
1, in principle. However, by considering howthe (cid:15) M k → T → (cid:15) M k → H M := ˆ H − µ ˆ N , which yields this1RDM. Since the eigenstates of the one-body Hamiltonianand the 1RDM coincide, it is clear that the eigenfunctionsof the one-body Hamiltonian should be the natural or-bitals (eigenfunctions of the 1RDM). The correspondingeigenvalues are readily obtained by inverting the Fermidistribution (cid:15) M k = 1 β ln (cid:18) − f k f k (cid:19) = 1 β ln (cid:18) ¯ f k f k (cid:19) . (17)From this expression it is clear that if the Matsubara en-ergies in the β → ∞ limit decay as (cid:15) M k ∼ /β that theoccupation number will have a definite value 0 < f k < .J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. 5 β -2 β -3/2 β -1 β -1/2 β β β -1/2 β -1 β -3/2 β -2 f k ( T ) β Figure 1.
Behaviour of the occupation numbers calculatedfor different temperature dependence of the Matsubara energy, (cid:15) M k = β − n ln(3). The β curve (black) corresponds to the moreaccustomed temperature independent behaviour. f k = 0.9f k = 0.8f k = 0.1 ε k ( T ) − − β Figure 2.
Behaviour of the Matsubara energies as a functionof the temperature for fractional occupation numbers. Theseoccupation numbers are 0 . values depending on its sign. If the Matsubara energy de-cays faster than 1 /β , the corresponding occupation num-ber converges to f k = / . The behaviour of the occupationnumber in the T → T → To keep the calculations as simpel as possible, we limitourselves to a two-orbital model for the hydrogen molecule.Let us start from the full interacting Hamiltonian for H ˆ H = (cid:88) ij,σ h ij ˆ c † iσ ˆ c jσ + (cid:88) ijkl,σσ (cid:48) w ijkl ˆ c † iσ ˆ c † jσ (cid:48) ˆ c kσ (cid:48) ˆ c lσ , (18)where the one-electron matrix elements, h ij := (cid:104) ϕ i | ˆ h | ϕ j (cid:105) = (cid:90) d r ϕ ∗ i ( r )ˆ hϕ j ( r ) , (19) contain the kinetic energy and external potential. Thetwo-electron matrix elements, w , describe the Coulombinteraction between the electrons and are given as w ijkl := (cid:90) d r (cid:90) d r (cid:48) ϕ ∗ i ( r ) ϕ ∗ j ( r (cid:48) ) ϕ k ( r (cid:48) ) ϕ l ( r ) | r − r (cid:48) | = [ il | jk ] . (20)In this equation we also used the quantum chemical no-tation [55] which regards the two-electron integral as aweighted overlap between charge distribution ϕ ∗ i ( r ) ϕ l ( r )and ϕ ∗ j ( r (cid:48) ) ϕ k ( r (cid:48) ). When we work with orthonormal or-bitals, we can make the tight-binding approximation tothe two-electron integrals[ il | kl ] ≈ δ il δ kl [ ii | jj ] . (21)The underlying idea is that the charge distributions ϕ ∗ i ( r ) ϕ l ( r )integrates to zero for i (cid:54) = l and therefore leads to a muchsmaller value of the two-electron integral than the diag-onal term i = l . If we further only use one orbital perhydrogen atom, the Hamiltonian simplifies in the tight-binding approximation toˆ H = α ˆ N + t (cid:88) σ (cid:0) ˆ c † σ ˆ c σ + ˆ c † σ ˆ c σ (cid:1) + w ˆ n ˆ n + U (cid:0) ˆ n ↑ ˆ n ↓ + ˆ n ↑ ˆ n ↓ (cid:1) + ¯ U (cid:88) i,σ ˆ c † iσ ˆ c † iσ ˆ c iσ ˆ c iσ , (22)where ˆ n iσ := ˆ c † iσ ˆ c iσ , ˆ n i := ˆ n i ↑ + ˆ n i ↓ and ˆ N = ˆ n + ˆ n .This Hamiltonian can be regarded as an extension of theHubbard model for two sites. The normal Hubbard Hamil-tonian only includes the hopping term between the twosites, whose strength is governed by t , and the on-siteinteraction, whose strength is set by U . We additionallyinclude a term depending on the number of particles withstrength α . This term does not have an influence the eigen-states, since they are eigenfunctions of the total numberof particles, but does change the corresponding eigenval-ues, however, so is important to recover the correct elec-tronic energy. Apart from the on-site interaction, we alsoinclude the interaction between the electrons if they re-side on different sites and its strength is determined bythe parameter w .The last term in the extended Hubbard Hamiltonian isa self-interaction term and does not contribute in princi-ple. However, when we construct energy functionals fromthe perturbation expansion, the expansion for the Φ func-tional cannot not always directly be associated to a properexpansion of an anti-symmetric 2-body Green’s function.In that case the energy functional is not guaranteed tobe self-interaction free and ¯ U typically does make a con-tribution. The GW (RPA) approximation considered inthis article is an example of such an expansion which isnot self-interaction free, hence the final result depends on¯ U . In this work, we will consider two different values for¯ U . The value ¯ U = U corresponds to a spin-independentinteraction, such as the non-relativistic Coulomb interac-tion used in molecular Hamiltonians. The other sensiblevalue to use is ¯ U = 0, which explicitly eliminates the self-interaction at the level of the Hamiltonian. This is theusual choice in the Hubbard model. K.J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. ζ op t [ a . u .] H-H [a.u.]1 2 3 4 5 6 7 8 9 10
Figure 3.
The exponent of the orbital, ζ , as a function of theinternuclear distance, R H–H . To determine the values of the extended Hubbard pa-rameters, we choose normalised 1s orbitals located at eachhydrogen atom to construct our basis χ ( r ) = N ζ e − ζ | r − R A | , χ ( r ) = N ζ e − ζ | r − R B | , (23)where the normalization factor is given by N ζ = (cid:112) ζ /π .The exponent ζ will be variationally optimized to ob-tain the lowest energy. In the dissociation limit R := | R A − R B | → ∞ we have that ζ = 1. A localised orthonor-mal basis is readily construct by L¨owdin orthogonaliza-tion [56], ϕ = χS − / , where S ij = (cid:104) χ i | χ j (cid:105) is the overlapmatrix of the non-orthogonal basis. Using this L¨owdin or-thogonalised basis, the one-electron matrix elements in theextended Hubbard Hamiltonian can be determined as α = (cid:104) ϕ | ˆ h | ϕ (cid:105) = (cid:104) χ | ˆ h | χ (cid:105) − s (cid:104) χ | ˆ h | χ (cid:105) − s , (24a) t = (cid:104) ϕ | ˆ h | ϕ (cid:105) = (cid:104) χ | ˆ h | χ (cid:105) − s (cid:104) χ | ˆ h | χ (cid:105) − s , (24b)where s := (cid:104) χ | χ (cid:105) denotes the overlap. For more detailsone these expressions and explicit forms in terms of the or-bital exponent ζ and the internuclear distance R , consultAppendix A.The two-electron matrix elements can be expressed as U = (11 |
11) + s − s ) (cid:2) (11 | − (11 | (cid:3) , (25a) w = (11 |
22) + s − s ) (cid:2) (11 | − (11 | (cid:3) , (25b)where we need the two-electron matrix elements in thenon-orthogonal 1s-basis( ij | kl ) := (cid:90) d r (cid:90) d r (cid:48) χ ∗ i ( r ) χ j ( r ) χ ∗ k ( r (cid:48) ) χ l ( r (cid:48) ) | r − r (cid:48) | . (26)These two-electron matrix elements for the 1s-basis areworked out explicitly in terms of ζ and R in Appendix A.Due to the limited dimension of the Hilbert space, theextended two-site Hubbard model is easy to solve exactly. In the two-particle sector we obtain the following tripletsolutionsˆ c † ↓ ˆ c † ↓ |(cid:105) , √ (cid:0) ˆ c † ↑ ˆ c † ↓ + ˆ c † ↓ ˆ c † ↑ (cid:1)(cid:12)(cid:12)(cid:11) , ˆ c † ↑ ˆ c † ↑ |(cid:105) , (27)with the eigenvalue E = 2 α + w . The singlet states can besubdivided according to their parity. There is only oneungerade singlet state, √ (cid:0) ˆ c † ↑ ˆ c † ↓ − ˆ c † ↑ ˆ c † ↓ (cid:1)(cid:12)(cid:12)(cid:11) , with theeigenvalue E = 2 α + U . There are two gerade singlet stateswhich can be expressed as √ (cid:2) cos( α ± ) (cid:0) ˆ c † ↑ ˆ c † ↓ − ˆ c † ↓ ˆ c † ↑ (cid:1) − sin( α ± ) (cid:0) ˆ c † ↑ ˆ c † ↓ + ˆ c † ↑ ˆ c † ↓ (cid:1)(cid:3)(cid:12)(cid:12)(cid:11) , (28)where the angles α ± can be determined from the equationtan( α ± ) = w − U t ± (cid:115)(cid:18) w − U t (cid:19) + 1 . (29)The corresponding electronic energies of the gerade singletstates are E ± = 2 α + w + U ± (cid:115)(cid:18) w − U (cid:19) + 4 t . (30)The energy E − is the ground state energy which is mini-mized by optimizing the value for the orbital exponent ζ .The optimal value, ζ opt , is shown in Fig. 3. The value of theexponent goes to 1 when stretching the bond. This is ex-pected, since in the dissociation limit the system consistsof two separated hydrogen atoms. The asymptotic valueis approached from below, since the 1s orbital needs tobecome more diffuse to facilitate binding. Binding wouldbe more efficiently achieved by allowing the 1s orbital topolarize towards the other hydrogen atom, for example bymixing in the 2 p z orbital. By allowing for polarization, thereduction in the exponent would be less.The values for the total energy from our extended two-site Hubbard model are compared to the accurate val-ues obtained for the non-relativistic hydrogen molecule inthe Born–Oppenheimer approximation by Ko(cid:32)los and Wol-niewicz [57–60]. Since the free parameter ζ is optimized forthe ground state, X Σ + g , it is reproduced quite accurately.The energy of the triplet b Σ + u state is of less qualityaround the equilibrium bond length R e = 1 . B Σ + u state is reasonable and correctly becomes degener-ate with the EF Σ + g state in the dissociation limit. Thedouble-well structure of the EF Σ + g is completely absent,which is no surprise, since the higher lying GK Σ + g re-quired for the necessary avoided crossing is not present inour simple two-orbital model. In the following sections we will assess the performance ofthe Klein and the Luttinger–Ward (LW) functionals using .J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. 7 X Σ +g EF Σ +g GK Σ +g b Σ +u B Σ +u E t o t [ a . u .] − − − − − − − − H-H [a.u.]1 2 3 4 5 6 7 8 9 10
Figure 4.
Comparison of the extended Hubbard energies (cir-cles) with the exact potential energy curves (triangles) [57–60].The energy of the (first) two/three Σ + g states are shown inblack, the energy of the (lowest) Σ + u state in blue and the(lowest) Σ + u state in red. the Green’s functions as input: the Kohn–Sham (KS), theHartree–Fock (HF) and 1RDM Green’s function. Thesenon-interacting Green’s functions are based on a one-bodyHamiltonian. Since we only use two basis functions in themodel, the symmetry of the system dictates that the eigen-functions of the one-body Hamiltonian need to be φ g/u ( r ) = ϕ ( r ) ± ϕ ( r ) √ χ ( r ) ± χ ( r ) (cid:112) ± s ) . (31)So the Matsubara non-interacting Green’s function (13)are also diagonal in this symmetry adapted basis. As werestrict the optimisation to symmetry adapted solutions,we do not allow for broken symmetry states to handlestrong correlation effects. The full burden is placed on thefunctional.We will limit the discussion to neutral H , so N =2( f g + f u ) = 2, which immediately implies that the chem-ical potential should be chosen such that (cid:15) M g + (cid:15) M u = 0,so the chemical potential in the non-interacting system isrequired to be µ s = (cid:15) g + (cid:15) u . (32) One of the simplest perturbation expansions involves onlythe Hartree–Fock (HF) diagrams − Φ HF = 12 + 12 = − β (cid:104)(cid:0) U + 2 w (cid:1) ˜ n − w ˜ f (cid:105) , (33)where ˜ n := f g + f u = N/ f := f g − f u denotesthe difference in occupation of the gerade and ungeradeorbital. The corresponding HF self-energy is particularlysimple, since it is local in time Σ HF kl ( z , z ) = z k l z + z k l z = v HF kl δ ( z , z ) , (34)where the HF potential for our simple model system isdiagonal in the symmetry-adapted basis and the non van-ishing elements are (see supplement) v HF gg = v HF¯ g ¯ g = (cid:0) U + 2 w (cid:1) ˜ n − w ˜ f , (35a) v HF uu = v HF¯ u ¯ u = (cid:0) U + 2 w (cid:1) ˜ n + w ˜ f . (35b)The lowest energy is obtained by fully occupying the σ g or-bitals, ˜ f = 1. As the HF potential differs by w between thegerade and ungrade orbitals, so the HF gap is increased by w with respect to the gap of the non-interacting system.Hence, we have (cid:15) HF = (cid:15) HF u − (cid:15) HF g = − t + w. (36a)As the HF potential is completely fixed by (35), the chem-ical potential in the HF approximation becomes µ HF (cid:0) ˜ n = 1 (cid:1) = 12 tr (cid:8) h + v HF (cid:9) = α + U + 2 w . (37)Later, we will demonstrate that the correlation part ofthe self-energy does not introduce any additional contribu-tions to the constant in the potential, so µ = µ HF . Hence,the effective potential in (8) becomes v HF s = v ext + v HF . (38a) The KS system has a particularly simple realisation inthis two-orbital model. As the KS potential is local, itscan only have diagonal elements in the site basis. As theKS should not break the symmetry of the system, it needsto be equal on both sites, so it can only be a constant shift.As the only constant contributions come from the externalpotential and the Hartree part, we need to set v KS s = µ = v ext + U + 2 w (38b)in order ot have µ KS = µ . Since t < R H–H ,occupying the σ g orbital will always lead to the lowest KSenergy. The gap in the KS system is now completely fixedby hopping matrix elements (cid:15) KS = (cid:15) KS u − (cid:15) KS g = − t. (36b) K.J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. β = ∞ β = 8 β = 2 β = 1 f = f g - f u − − ε [a.u.] − − Figure 5.
The difference between the (natural) occupationnumbers as a function of the gap for several temperatures.
The 1RDM Green’s function is defined in a similar man-ner as the KS Green’s function, except that we do notconstrain the effective potential to be local anymore. Al-lowing for non-local potential, gives more flexibility, butas the potential should retain the symmetry of the sys-tem, the non-local 1RDM potential remains diagonal inthe symmetry adapted basis of our model system. As thetrace of the potential is fixed by the particle number, thenon-local potential can only adjust the gap (cid:15) in the fol-lowing manner v s ( (cid:15) ) = v KS s + (cid:15) + 2 t (cid:18) − (cid:19) . (38c)So in the 1RDM setting, we have one variable to optimise:the gap (cid:15) . At finite inverse temperature, the relation be-tween the gap and the difference in occupation numbersis one-to-one and explicitly given by˜ f := f g − f u = tanh( β(cid:15)/ . (39)In the zero-temperature limit ( β → ∞ ), this function col-lapses into a step function as illustrated in Fig. 5. Thisdemonstrates explicitly that the occupation numbers canonly be fractional at zero temperature, if the gap closes [32–34, 61]. In particular, we have˜ f = (cid:15) > − ,
1] if (cid:15) = 0 − (cid:15) < . (40)As we loose the one-to-one relation between the gap andthe occupation number difference in the zero-temperaturelimit, we will need to optimise over the separate pieces ofthe β = ∞ curve: { (cid:15) ≤ , ˜ f = − } , { (cid:15) = 0 , − ≤ ˜ f ≤ } and { (cid:15) ≥ , ˜ f = 1 } . For simplicity we have limited the discussion on possibleinput Green’s functions to the simple two-orbital model for H , but all three possibilities can readily be used inmore general settings. The simplest case is the use of theHF Green’s function, since one only needs to do a (re-stricted) Hartree–Fock calculation and insert the resultingGreen’s function into either the Klein (11) or Luttinger–Ward functional (8) with some appropriate approximationfor Φ [ G ]. As we will use exclusively the GW approxima-tion for Φ in this work, we will denote these calculationsby K-GW@HF and LW-GW@HF respectively.The insertion of Kohn–Sham type Green’s function re-sults in some variational freedom, as the Kohn–Sham typeGreen’s function encompasses all non-interacting Green’sfunctions which can be generated via a local potential, i.e.potentials diagonal in spatial representation or site basis.So in a general setting, we can use the freedom of the lo-cal potential optimise the energy expression which resultsfrom inserting the KS type Green’s function into the Kleinor LW functional. Within the GW approximation, we willdenote such calculations as K-GW@KS and LW-GW@KSrespectively. As discussed before in Sec. 5.2, the symmetryconstraints in the limited two-orbital model (homogeneousdensity) effectively leaves no variational freedom in the KSpotential. This is clearly an artefact of our limited setting.The space of non-interacting Green’s functions to besearched over can be enlarged by allowing the potential tobe non-local, i.e. we only require the potential to be her-mitian ( v † s = v s ). This class of potentials is exactly thetype of potentials used in 1RDM functional theory. Theinsertion of these more general non-interacting Green’sfunctions in the energy functionals within the GW ap-proximation and subsequent optimisation, will thereforebe denoted as K-GW@1RDM and LW-GW@1RDM re-spectively. In this section we consider the evaluation of the LW andKlein functionals for various input Green’s functions.To benchmark the results we have also solved the Dysonequation in the GW approximation self-consistently to ob-tain the true stationary Green’s function of these function-als. The code has been implemented for spin-independentinteractions, so we only have these results for the spin-independent interaction, ¯ U = U . For bond distances largerthan 6 . . The Hartree and exchange part of the Φ functional werealready evaluated in the previous section (Hartree–Fock),so we only need to evaluate the correlation part. The cor-relation part of the Φ -functional in the GW approximation .J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. 9 can be written as − Φ GWc = 14 + 16 + · · · = ∞ (cid:88) n =2 n (cid:0) ¯ wP (cid:1) n = −
12 Tr (cid:8) ln (cid:0) − ¯ wP (cid:1) + ¯ wP (cid:9) , (41)where the trace Tr is both over matrix entries as well asthe Matsubara time/frequencies. The bar over the inter-action indicates that its indices are in chemical ordening,¯ w ijkl := [ il | jk ] = w iklj , and the polarization bubble isdefined as P abb (cid:48) a (cid:48) ( ω ) := i b b (cid:48) a a (cid:48) ω = 1 β (cid:88) m G a (cid:48) a ( ω m ) G bb (cid:48) ( ω + ω m ) . (42)in frequency space. Using that tr (cid:8) ln( M ) (cid:9) = ln (cid:12)(cid:12) M (cid:12)(cid:12) , i.e.the logarithm of the determinant, the expression for thecorrelation part of the Φ -functional in the GW approxi-mation can be further simplified to Φ GWc = 12 (cid:88) p e ηω p (cid:0) ln (cid:12)(cid:12) − ¯ wP ( ω p ) (cid:12)(cid:12) + tr { ¯ wP ( ω p ) } (cid:1) , (43)where the trace, tr, is now only over matrix entries. For oursimple non-interacting Green’s function, the determinantcan be worked out as (see supplement) (cid:12)(cid:12) − ¯ wP s ( ω ) (cid:12)(cid:12) = (cid:32) − (cid:15) ( u + v ) ˜ fω − (cid:15) (cid:33) × (cid:32) − (cid:15) ( u − v ) ˜ fω − (cid:15) (cid:33) , (44)where (cid:15) := (cid:15) u − (cid:15) g denotes the gap and we introduced theshort-hand notations u := (cid:0) ¯ U − w (cid:1) and v := (cid:0) U − w (cid:1) .Evaluating the remaining sum over the bosonic frequencies(see supplement) yields Φ GWc ,s = ln (cid:18) sinh (cid:16) βζ + (cid:17)(cid:19) + ln (cid:18) sinh (cid:16) βζ − (cid:17)(cid:19) (45) − (cid:18) sinh (cid:16) β | (cid:15) | (cid:17)(cid:19) − β u (cid:0) ˜ f + ˜ n (2 − ˜ n ) (cid:1) , It is tempting to remove the ¯ wP term from the trace, sinceit seems we would include exchange immediately. However, thetime limits in this integral do not correspond to the ones forexchange, so this does not work unfortunately. E K - G W t o t [ a . u .] − − − ε [a.u.] − − g -f u − ε [a.u.]0 1 2 ε = U = . ε = - U = - . Figure 6.
The total GW energy as a function of the gap andoccupation number difference at R H–H = 1 . U = U , yields a continuous de-pendence (red), whereas excluding the self-interaction term,¯ U = 0, results in a discontinuous energy dependence (blue)due to negative numbers in the square root in (47). where ζ ± := (cid:113) (cid:15) + 2 (cid:15) ( u ± v ) ˜ f . (46)Taking the T → n = 1 E K-GWc = lim β →∞ β Φ GWc = ζ + ζ − − | (cid:15) | − u (1 + ˜ f )2 . (47)Since the K-GW correlation energy (47) depends explicitlyon the self-interaction term ¯ U , we should investigate dif-ferent values. As mentioned in the introduction, we will in-vestigate a spin-independent interaction ( ¯ U = U ) and theexplicitly self-interaction free model ( ¯ U = 0). The total K-GW energy is plotted as a function of the gap/occupationnumbers for these two cases in Fig. 6 for R H–H = 1 . β = ∞ -path as depicted inFig. 5. In the case of a spin-independent interaction, theK-GW energy can be calculated for all values of the gap asis shown by the continuous line. If the self-interaction termis set to zero, ¯ U = 0, small values of the gap yield an imag-inary part to the energy, since we have u − v = − U/ < ζ − of the K-GW correlationenergy (47). Only when | (cid:15) | ≥ U or (cid:15) = 0, we obtain areal value for the K-GW energy as is apparent from thediscontinuous curve in Fig. 6.The relative positioning of the minima of the total GWenergy as depicted in Fig. 6 turns out to be the same at allbond distances. The minimum is always located at a fullyoccupied gerade orbital ˜ f = 1 for both choices of ¯ U . In thecase of a spin-independent interaction, ¯ U = U , the optimalgap is always located at (cid:15) = 0. Excluding self-interactionterm, ¯ U = 0, always yields the global minimum at thesmallest strictly positive gap, (cid:15) = U .Different methods to use the K-GW energy expres-sion (47) in the spin-independent interaction case are com-pared to the exact ground state in Fig. 7. The 1RDM ver- E K - G W t o t [ a . u .] − − − − − − − − H-H [a.u.]1 2 3 4 5 6 7 8 9 10 − − − − − − Figure 7.
The ground state energy calculated with the Klein-GW functional compared to the exact Hubbard results (blackcircles), HF (red crosses) and the self-consistent GW energy(orange stars). The self-interaction term is included: ¯ U = U .The following three different versions are shown: the 1RDMversion (blue triangles), the DFT version (green squares) andthe HF version (purple diamonds). sion of the K-GW functional (K-GW@1RDM) is consis-tently too low, but becomes exact in the dissociation limit.Restricting the potential to be local (K-GW@KS) raisesthe total energy and it becomes quite accurate around theequilibrium distance, though still slightly too low. In thedissociation regime, the upshift is too high and createsthe infamous RPA bump. Since the KS gap also vanishesin the dissociation limit, the correct dissociation limit isretained [23, 62]. Alternatively, one could use the HF gapas input (K-GW@HF). The HF gap is larger than the KSgap due to the additional non-local exchange in the po-tential (36a). This leads to a higher total energy, sincethe K-GW correlation energy increases monotonically in (cid:15) for positive gaps as can be seen from (47). The K-GW@HF energy gives an improvement over the KS gap re-sult around the equilibrium distance, in the sense that thetotal energy agrees even better with the exact result (seeFig. 7). Upon dissociation however, the artificial bump israised to even higher values, leading to a worse perfor-mance than the KS gap. Though not so clear from Fig. 7,the dissociation limit remains correct, since also the HFgap goes to zero when R H–H → ∞ . Albeit, the closurerate is much slower. Due to the additional exchange con-tribution to the gap, the HF gap only closes as 1 /R H–H ,whereas the KS gap closes exponentially.Since the GW approximation is not self-interaction free(depends on ¯ U ), one would expect that the results wouldimprove if the self-interaction term would be explicitly R H-H = 1.275R
H-H = 2.269 ε [ a . u .] E K - G W t o t [ a . u .] − − − − − − − − H-H [a.u.]1 2 3 4 5 6 7 8 9 10
Figure 8.
The upper panel is the same as in Fig. 7 but nowwithout the self-interaction term: ¯ U = 0. In the lower panel thevalue of the gap is shown which went into the GW correlationenergy expression (47). excluded from the Hamiltonian. The results for the totalenergy for ¯ U = 0 are shown in the upper panel of Fig. 8. Infact, elimination of the self-interaction does actually notlead to any improvement at all. Using the K-GW func-tional as a 1RDM functional leads to even lower energiesaround the equilibrium distance and upon dissociation theenergy becomes higher than the exact results. Even thedissociation limit is not correct anymore, since the gapneeds to remain finite, (cid:15) = U , to prevent the energy fromhaving a complex part, as mentioned before in connectionwith Fig. 6. The results are even more disastrous when theKS or HF gaps are used. Since both gaps decrease whenstretching the bond, there is always a point from wherethey become smaller than U and the K-GW correlationenergy becomes complex. This is illustrated in the lowerpanel in Fig. 8, where the different gaps of each calculationare plotted. Since the 1RDM gap is always the minimumgap, the failure of the K-GW functional with the KS/HFgap exactly occurs when it crosses the 1RDM gap (cid:15) = U .At these points, the K-GW energy from the KS/HF gapbecomes exactly equal to the K-GW@1RDM functional. .J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. 11 In an attempt to improve the results, one can use theLuttinger–Ward (LW) functional instead of the Klein func-tional, as the LW functional has superior variational prop-erties [22, 23]. The superior variational properties come atthe cost of a more complicated functional which requiresthe evaluation of the following additional terms Ω LW − Ω K = − β Tr (cid:8) ˜ ΣG + ln (cid:0) − G s ˜ Σ (cid:1)(cid:9) , (48)where the modified self-energy ˜ Σ was defined in (9).Since the HF part has already been evaluated before,cf. (34) and (35), we only need the correlation part of theself-energy in the GW approximation evaluated at G s .The evaluation is somewhat simplified by the fact thatthe self-energy is also diagonal in the symmetry-adaptedbasis. Its diagonal elements can be evaluated to be Σ GWc ,s ; gg [ G s ]( ω ) = R + ( ω − (cid:15) M u ) , (49a) Σ GWc ,s ; uu [ G s ]( ω ) = R − ( ω − (cid:15) M g ) , (49b)where we introduced R ± ( ω ) := (cid:15) ˜ f (cid:20) ( u + v ) ζ + ωω − ζ coth (cid:16) βζ + (cid:17) + ( u − v ) ζ − ωω − ζ − coth (cid:16) βζ − (cid:17) (50)+ ( u + v )( ω − (cid:15) ) − (cid:15) ˜ f ( u − v )( ω − ζ )( ω − ζ − ) × (cid:0) − ˜ n ± ˜ f (cid:1)(cid:21) . For ˜ n = f g + f u = 1, R ± satisfies R + ( − ω ) = − R − ( ω ).As in that case also (cid:15) Mg = − (cid:15) Mu , the components of theself-energy are simply related as Σ GWc ,s ; gg [ G s ]( − ω ) = − Σ GWc ,s ; uu [ G s ]( ω ) if ˜ n = 1. (51)This relation is convenient, as it implies that the numberof particles is not affected by the correlation part for ˜ n = 1(see Appendix B), proving the validity of the constant in(trace of) the potentials used to generate the input Green’sfunctions (38).The linear term of the LW correction in (48) can nowbe evaluated to beTr (cid:8) Σ GWc ,s G s (cid:9) = β(cid:15) ˜ f ( u + v ) ζ + coth (cid:16) βζ + (cid:17) + β(cid:15) ˜ f ( u − v ) ζ − coth (cid:16) βζ − (cid:17) − β (cid:0) ˜ n (2 − ˜ n ) + ˜ f (cid:1) u. (52)Unfortunately, the logarithmic term in (48) cannot beevaluated analytically. The integrand decays as 1 /ω , sonumerical evaluation is not a problem. In the gapless limit, (cid:15) →
0, however, an analytical solution is in reach. In the limit of vanishing gap, the self-energy simplifies consider-ably to lim (cid:15) → Σ GWc ,s ; gg/uu ( ω + (cid:15) M u/g ) = uβω . (53)In the β → ∞ limit its contribution will therefore van-ish, so only the contribution of the correlation potentialin the logarithmic term remains. Evaluation of the loga-rithmic term with only the correlation potential is ratherstraightforward and giveslim β →∞ β Tr (cid:8) ln (cid:0) + v c G s (cid:1)(cid:9) = 2 | v c gg | . (54)where v c := v s − v ext − v HF . (55)The values for the total energies evaluated with the LWfunctional in the GW approximation are plotted as a func-tion of the gap and the occupation number difference inFig. 9 for a bond distance of R H–H = 1 . U = U , the totalenergy can be calculated for any gap and the minimumis now located at a positive gap. In the case of the self-interaction free interaction, ¯ U = 0, the energy is only realwhen the gap is either sufficiently large or zero. This is ex-actly the same situation as described before for the Kleinfunctional. For short bond distances the minimum is lo-cated at ˜ f = 1 and (cid:15) = 0, but at larger bond distances theminimum shifts to a positive gap. This relocation of theminimum is demonstrated in Fig. 10 for R H–H = 5 . R H–H = 4 .
25 Bohr. The relocation ofthe minimum only occurs for self-interaction free case andnot for the spin-independent interaction case.In the top panel of Fig. 11 we show the results for thetotal energies from the Luttinger–Ward functional for thespin-free Hamiltonian, ¯ U = U . Compared to the Kleinfunctional we see a huge improvement of the total energywhen the LW functional is used as a 1RDM functional,but unfortunately, the dissociation limit is not correct any-more. As a matter of fact, the optimised 1RDM results arealmost identical to the results when the HF Green’s func-tion is used as input. Only in the dissociation limit theoptimisation of the non-local potential yields a somewhatlower result. This trend is corroborated by comparing theHF gap and the 1RDM gap visualised in the lower panelof Fig. 11. At equilibrium distance the gaps are nearlyidentical and they start to deviate somewhat when theinteratomic distance is increased. Somewhat surprisingly,increasing the gap slightly from its HF value yields a lowerenergy.We have also included the results when the KS Green’sfunction is used as input for the LW functional. The KSGreen’s function now gives the worst result of all Greensfunctions, which is somewhat surprising, as it performedso well in the Klein functional. The results for the Kleinfunctional are mainly due to error cancelations betweenthe inflexibility of the KS Green’s function and the badvariational properties of the Klein functional. A similarinferior behaviour of the KS Green’s function as input forthe LW functional has been obtained for the H moleculein a large Slater basis [23]. ε = U = . ε = - U = - . E L W - G W t o t [ a . u .] − − − ε [a.u.] − − g -f u − ε [a.u.]0 1 2 Figure 9.
The total LW-GW energy as a function of the gapand occupation number difference at R H–H = 1 . U = U , yields a continuousdependence (red), whereas excluding the self-interaction term,¯ U = 0, results in a discontinuous energy dependence (blue)due to negative numbers in the square root in (47). E L W - G W t o t [ a . u .] − − − − − − − ε [a.u.] − − g -f u − ε [a.u.]0 1 2 ε = U = . ε = - U = - . Figure 10.
The total LW-GW energy as a function of the gapand occupation number difference at R H–H = 5 Bohr. Again,the continuous (red) line corresponds to the spin-dependentinteraction, ¯ U = U , and the the discontinuous one (blue) theself-interaction free interaction ¯ U = 0. Now let us consider the results in the self-interactionfree case, ¯ U = 0. The results for the total energies aredepicted in Fig. 12. Most notable are similar catastro-phes for the KS and HF input Green’s functions as forthe Klein functional, cf. Fig. 8. When the bond is suffi-ciently stretched, both the KS and HF Green’s functionyield from some distance onwards a complex energy, dueto the square root in ζ − . These catastrophes therefore oc-cur at exactly the same bond distances as for the Kleinfunctional. The only difference is that for the HF Green’sfunction the catastrophe for the LW functional occurs in adifferent manner than for the Klein functional. In the caseof the KS functional the behaviour of the two functionalsis nearly identical.When using the full flexibility of a non-local potentialto generate a trial Green’s function, the LW functionaldoes not perform better than the Klein functional in the¯ U = 0 case. Around the equilibrium distance the total en- ε [ a . u .] E L W - G W t o t [ a . u .] − − − − − − − − − H-H [a.u.]1 2 3 4 5 6 7 8 9 10 − − − − − − Figure 11.
Upper panel: the ground state energy calcu-lated with the Luttinger–Ward-GW functional compared tothe exact Hubbard results (black circles), HF (red crosses)and the self-selfconsistent GW energy (orange stars). The self-interaction term is included: ¯ U = U . The following three dif-ferent versions are shown: the 1RDM version (blue triangles),the DFT version (green squares) and the HF version (purplediamonds). In the lower panel the value of the gap is shownwhich went into the LW-GW correlation energy expression. ergy is somewhat lower compared to the Klein functional,so the LW functional yields worse results. Approachingthe dissociation limit, R H–H = 10 . R H–H = 4 .
25 Bohr. This is the location where the globalminimum jumps from the (cid:15) = 0 solution to (cid:15) > U solution(see Figs. 9 and 10). When one of these solutions becomeshigher in energy, we have continued plotting the solutionwith a dotted line in Fig. 12. If now the (cid:15) = 0 solutionwould be disregarded altogether, we see that the finitegap solution actually yields a quite reasonable result. Wecould not think of a proper physical argument to discardthe gapless solution, so we do not promote this procedure.All these issues remain a topic for future research. .J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. 13 R H - H = . R H - H = . R H - H = . E L W - G W t o t [ a . u .] − − − − − − − − H-H [a.u.]1 2 3 4 5 6 7 8 9 10
Figure 12.
The same as the upper panel in Fig. 11, thoughnow with the self-interaction free interaction, ¯ U = 0. The nexthigher energy solutions when using the LW as a 1RDM func-tional are indicated by the (blue) dotted lines. In this work we considered two different energy functionalsof the 1RDM. These were obtained by restricting the do-main of the LW and Klein functionals to non-interactingGreen’s function for a Hamiltonian with a spatially non-local potential. These functionals were evaluated withinthe GW approximation for two different cases of self-interactionfor an extended Hubbard Hamiltonian representing a modelsystem for the hydrogen molecule. We compared the bind-ing curves of this system with exact results as well as withenergy functionals with HF and KS input Green’s func-tions. We found that the LW functional outperforms theKlein functional and gives results very close to the opti-mal self-consistent GW Green’s function which shows thatthere is room for future improvement of the method. Es-pecially the GW approximation does not correctly treatthe self-interaction correctly and a probable improvementcan be attained by including T-matrix and ladder diagramexpansions to introduce the necessary exchange-type di-agrams. The system that we considered is very rigid dueto its low dimensionality and systems in larger basis setsare likely to increase the variational freedom to attainmore accurate results. We further studied the system inthe zero-temperature limit which led to discontinuities inthe energy landscape. It is therefore worthwhile to furtherexplore systems at finite temperature to smoothen thesediscontinuities.
KJHG thanks B.C.E. Giesbertz for useful suggestions and grate-fully acknowledges a VENI grant by the Netherlands Founda-tion for Research NWO (722.012.013).
A Derivation of matrix elements of thetwo-orbital model for H Defining ρ := ζR , the one-electron integrals in the non-orthogonal 1s-basis become s := (cid:104) χ | χ (cid:105) = (cid:18) ρ + ρ (cid:19) e − ρ , (cid:104) χ | ˆ h | χ (cid:105) = ζ − ζ + ζρ (cid:0) e − ρ (1 + ρ ) − (cid:1) , (A.1) (cid:104) χ | ˆ h | χ (cid:105) = e − ρ (cid:20) ζ (cid:18) ρ − ρ (cid:19) − ζ (1 + ρ ) (cid:21) . The unique two-electron integrals in the non-orthogonal1s-basis are(11 |
11) = 5 ζ , (A.2a)(11 |
22) = ζ (cid:20) ρ − e − ρ (cid:18) ρ + 118 + 3 ρ ρ (cid:19)(cid:21) , (A.2b)(11 |
12) = ζ ρ (cid:2) e − ρ (5 + 2 ρ + 16 ρ ) − e − ρ (5 + 2 ρ ) (cid:3) , (A.2c)(12 |
12) = ζ e − ρ ρ (cid:104) ρ − ρ − ρ − ρ + 16 (cid:0) γ + ln( ρ ) (cid:1)(cid:0) ρ + 15 ρ + 6 ρ + ρ (cid:1) − (cid:0) − ρ + ρ (cid:1) e ρ Ei( − ρ )+ 16(3 − ρ + ρ ) e ρ Ei( − ρ ) (cid:105) , (A.2d)where γ ≈ . x ) := − (cid:90) ∞− x d y e − y y . (A.3)The formula (A.2d) has been derived using the procedurepresented in [63] and agrees to all digits with the numericalresults presented in [64]. In the calculations we use theMulliken approximation to the two-electron integrals [65]( il | jk ) ≈ S il S jk (cid:2) ( ii | jj )+( ll | jj )+( ii | kk )+( ll | kk ) (cid:3) , (A.4)where S ij := (cid:104) χ i | χ j (cid:105) is the overlap. That means that thefollowing integrals are approximated as(11 | ≈ s (cid:2) (11 |
11) + (11 | (cid:3) , (A.5a)(12 | ≈ s (cid:2) (11 |
11) + (11 | (cid:3) . (A.5b) (11|12)(12|12) ( e x a c t - M u lli k en ) [ m H ] − − − H-H [a.u.]0 1 2 3 4 5 6 7 8 9 10
Figure 13.
Difference between the exact integrals and theintegrals in the Mulliken approximation. Note that the error isin milliHartree.
The Mulliken approximation has been verified numericallyas a function of the bond length combined with the opti-mized exponent. The error of the Mulliken approximationis shown in Fig. 13 and is of the order of mH in around theequilibrium distance and decreases quickly for the (12 | |
11] = 12(1 − s ) (cid:2) (2 − s )(11 |
11) + s (11 | − s (11 |
12) + 2 s (12 | (cid:3) , (A.6a)[11 |
22] = 12(1 − s ) (cid:2) (2 − s )(11 |
22) + s (11 | − s (11 |
12) + 2 s (12 | (cid:3) , (A.6b)[11 |
12] = 12(1 − s ) (cid:2) s )(11 | − s (12 | − s (cid:0) (11 |
11) + (11 | (cid:1)(cid:3) , (A.6c)[12 |
12] = 12(1 − s ) (cid:2) | − s (11 | s (cid:0) (11 |
11) + (11 | (cid:1)(cid:3) . (A.6d)Now using the Mulliken approximation, the expressionsfor the two-electron integrals in the L¨owdin orthogonalisedbasis simply to[11 |
11] = (2 − s )(11 | − s (11 | − s ) , (A.7a)[11 |
22] = (2 − s )(11 | − s (11 | − s ) , (A.7b)[11 |
12] = 0 , (A.7c)[12 |
12] = 0 . (A.7d)Hence, we find that the Mulliken approximation to thetwo-electron integrals for the H molecule in a minimalbasis, automatically leads to the extended Hubbard model used in our paper, with U = [11 |
11] and w = [11 | B Particle number consistency
To have a consistent number of particles in the variousapproximations to the grand potential of the interactingsystem, we need to guarantee that we do not change theparticle number when going from Ω s to Ω , so we require N s = N . In other words, we need that ∂ ( Ω − Ω s ) ∂µ = 0 . (B.1)When inserting G s in the Klein functional, this amountsto the requirementtr (cid:8) ( v ext − v s ) ∂ γ ∂µ (cid:9) + 1 β Tr (cid:26) δΦδ G s ∂ G s ∂µ (cid:27) = 1 β Tr (cid:26) ˜ Σ ∂ G s ∂µ (cid:27) = 0 . (B.2)Let us first consider the HF approximation. The HF self-energy is local in time, so G s can be replaced by γ . Thederivative of the 1RDM with respect to the chemical po-tential is readily worked out as ∂ γ /∂µ = βf g f u , so theparticle conservation condition (B.2) puts the followingcondition on the trace of the effective potential of the non-interacting reference systemtr (cid:8) v s (cid:9) = tr (cid:8) Σ HF + v ext (cid:9) = 4 (cid:18) α + U + 2 w (cid:19) , (B.3)where we used the matrix elements of the HF self-energy (35)and have set ˜ n = 1. This exactly agrees with the chemicalpotential in the HF approximation (37), as this leads to (cid:15) M g + (cid:15) M u = 0.Including the correlation part of the GW functionaldoes not lead to a change in the trace of v s . This is mosteasily established by considering the derivative of Φ GW c in (45) directly, as it yields ∂Φ GWc ,s ∂µ (cid:12)(cid:12)(cid:12)(cid:12) ˜ n =1 = − βu − n ) ∂ ˜ n∂µ (cid:12)(cid:12)(cid:12)(cid:12) ˜ n =1 = 0 . (B.4)We find therefore, that the trace condition on v s remainsthe same as in (B.3) when using the GW approximationin the Klein functional.In the Luttinger–Ward functional we have some ad-ditional terms (48). For the derivative of the linear termwith respect to the chemical potential we find directlyfrom (52) ∂∂µ Tr (cid:8) ˜ ΣG s (cid:9)(cid:12)(cid:12)(cid:12)(cid:12) ˜ n =1 = − β (2 − n ) u ∂ ˜ n∂µ (cid:12)(cid:12)(cid:12)(cid:12) ˜ n =1 = 0 . (B.5)For the logarithmic term we do not have an explicit ex-pression, but we can work it out as follows .J.H. Giesbertz, A. Uimonen, R. van Leeuwen: Approximate energy functionals for 1RDM functional theory MBPT. 15 − ∂∂µ Tr (cid:8) ln (cid:0) − G s ˜ Σ (cid:1)(cid:9) = Tr (cid:26)(cid:0) − G s ˜ Σ (cid:1) − (cid:18) ∂ G s ∂µ ˜ Σ + G s ∂ Σ ∂µ (cid:19)(cid:27) (B.6)The derivatives of the non-interacting Green’s functionwith respect to the chemical potential is readily workedout in the frequency domain as ∂G s,kl ( ω ) ∂µ = − δ kl (cid:0) ω − (cid:15) M k (cid:1) − , (B.7)which has the following useful property ∂G s,gg ( ω ) ∂µ = ∂G s,uu ( − ω ) ∂µ . (B.8)Likewise, the derivative of the GW self-energy with G s inserted, can be calculated directly from its explicit ex-pression in (49). For the gerade component we get ∂∂µ Σ GWc ,s ; gg/uu [ G s ]( ω ) = ∂∂ω Σ GWc ,s ; gg/uu [ G s ]( ω ) − ( u + v )( ω u/g − (cid:15) ) − (cid:15) ˜ f ( u − v )( ω u/g − ζ )( ω u/g − ζ − ) ∂ ˜ n∂µ , (B.9)where we used ω k = ω − (cid:15) M k as an abbreviation. Due tothe relation between the GW self-energy in (51) for ˜ n = 1,the derivative with respect to the chemical potential obeysthe following relation ∂∂µ Σ GWc ,s ; gg [ G s ]( ω ) = ∂∂µ Σ GWc ,s ; uu [ G s ]( − ω ) if ˜ n = 1.(B.10)Note that (51) also implies that if the trace of the singleparticle potential v s satisfies (B.3), that we also have˜ Σ GWc ,s ; gg [ G s ]( ω ) = ˜ Σ GWc ,s ; uu [ G s ]( − ω ) if ˜ n = 1. (B.11)Let us introduce the following abbreviation for the geradecomponent of the trace in (B.6) as F ( ω ) := ∂G s,gg ( ω ) ∂µ ˜ Σ gg ( ω ) + G s,gg ( ω ) ∂ ˜ Σ gg ( ω ) ∂µ − G s,gg ( ω ) ˜ Σ gg ( ω ) . (B.12)Using G s,gg ( ω ) = − G s,uu ( − ω ) and the relations in (51),in (B.8) and in (B.10), we find that the ungerade compo-nent is given by − F ( − ω ). The derivative of the logarith-mic can therefore be written as − ∂∂µ Tr (cid:8) ln (cid:0) − G s ˜ Σ (cid:1)(cid:9) = (cid:88) p (cid:0) F ( ω p ) − F ( − ω p ) (cid:1) = 0 , (B.13)where we used that the direction of summation over theMatsubara frequencies ω p is immaterial. Note that wecould drop the convergence factor e ηω p as the integrantalready converges without it. So also for the Luttinger–Ward functional (B.3) is the appropriate trace of the sin-gle particle potential to keep N = 2. References
1. O.V. Gritsenko, K. Pernal, E.J. Baerends, J. Chem.Phys. , 204102 (2005)2. D.R. Rohr, K. Pernal, O.V. Gritsenko, E.J. Baerends,J. Chem. Phys. , 164105 (2008)3. M. Piris, X. Lopez, F. Ruip´erez, J.M. Matxain, J.M.Ugalde, J. Chem. Phys. , 164102 (2011)4. F. Ruiperez, M. Piris, J.M. Ugalde, J.M. Matxain,Phys. Chem. Chem. Phys. , 2055 (2013)5. S. Sharma, J.K. Dewhurst, N.N. Lathiotakis, E.K.U.Gross, Phys. Rev. B , 201103(R) (2008)6. S. Sharma, J.K. Dewhurst, S. Shallcross, E.K.U.Gross, Phys. Rev. Lett. , 116403 (2013)7. E. T¨ol¨o, A. Harju, Phys. Rev. B , 075321 (2010)8. M. Buijse, Ph.D. thesis, Vrije Universiteit, De Boele-laan 1105, Amsterdam, The Netherlands (1991), http://theochem.chem.rug.nl/publications/PDF/ft217.pdf
9. M. Buijse, E.J. Baerends, Mol. Phys. , 401 (2002)10. (cid:32)L.M. Mentel, R. van Meer, O.V. Gritsenko, E.J.Baerends, J. Chem. Phys. , 214105 (2014)11. D.A. Mazziotti, Chem. Phys. Lett. , 323 (2001)12. M. Piris, P. Otto, Int. J. Quantum Chem. , 317(2003)13. M. Piris, Int. J. Quantum Chem. , 1093 (2006)14. M. Piris, J.M. Matxain, X. Lopez, J.M. Ugalde, J.Chem. Phys. , 021102 (2009)15. M. Piris, J.M. Matxain, X. Lopez, J. Chem. Phys. , 234109 (2013)16. M. Piris, J. Chem. Phys. , 044107 (2014)17. M. Piris, Phys. Rev. Lett. , 063002 (2017)18. K. Pernal, K.J.H. Giesbertz, in Density-FunctionalMethods for Excited States , edited by N. Ferr´e, M. Fi-latov, M. Huix-Rotllant (Springer, Berlin Heidelberg,2015), Vol. 368 of
Topics in Current Chemistry ,chap. 4, pp. 125–183, ISBN 978-3-319-22080-219. A.L. Fetter, J.D. Walecka,
Quantum Theory of Many-Particle Systems (Dover Publiations, Inc., 2003)20. G. Stefanucci, R. van Leeuwen,
Nonequilibrium Many-Body Theory of Quantum Systems: A Modern In-troduction (Cambridge Univeristy Press, New York,2013), ISBN 978-0-521-76617-321. F. Furche, Phys. Rev. B , 195120 (2001)22. N.E. Dahlen, U. von Barth, Phys. Rev. A , 195102(2004)23. N.E. Dahlen, R. van Leeuwen, U. von Barth, Phys.Rev. A , 012511 (2006)24. F. Furche, J. Chem. Phys. , 114105 (2008)25. M. Hellgren, U. von Barth, J. Chem. Phys. ,044101 (2010)26. A. Heßelmann, A. G¨orling, Mol. Phys. , 359 (2010)27. H. Eshuis, J.E. Bates, F. Furche, Theor. Chem. Acc. , 1084 (2012)28. J.E. Bates, F. Furche, J. Chem. Phys. , 171103(2013)29. P. Bleiziffer, A. Heßelmann, A. G¨orling, J. Chem.Phys. , 084113 (2013)30. M. Hellgren, F. Caruso, D.R. Rohr, X. Ren, A. Rubio,M. Scheffler, P. Rinke, Phys. Rev. B , 165110 (2015)31. W. Kohn, L.J. Sham, Phys. Rev. , A1133 (1965)
32. T.L. Gilbert, Phys. Rev. B , 2111 (1975)33. K. Pernal, Phys. Rev. Lett. , 233002 (2005)34. K.J.H. Giesbertz, E.J. Baerends, J. Chem. Phys. ,194108 (2010)35. G. Friesecke, Proc. R. Soc. London A , 47 (2003)36. K.J.H. Giesbertz, R. van Leeuwen, J. Chem. Phys. , 104109 (2013)37. K.J.H. Giesbertz, R. van Leeuwen, J. Chem. Phys. , 104110 (2013)38. K.J.H. Giesbertz, R. van Leeuwen, J. Chem. Phys. , 184108 (2014)39. K.J.H. Giesbertz, M. Ruggenthaler (2017),
40. A. G¨orling, M. Levy, Phys. Rev. B , 13105 (1993)41. A. G¨orling, M. Levy, Phys. Rev. A , 196 (1994)42. T. Baldsiefen, Ph.D. thesis, Institut f¨ur TheoretischePhysik Freie Universit¨at Berlin (2012)43. T. Baldsiefen, E.K.U. Gross (2012),
44. P.E. Bl¨ochl, T. Pruschke, M. Potthoff, Phys. Rev. B , 205139 (2013)45. C.O. Almbladh, U. von Barth, R. van Leeuwen, Int.J. Mod. Phys. B , 535 (1999)46. F. Aryasetiawan, T. Miyake, K. Terakura, Phys. Rev.Lett. , 166401 (2002)47. T. Miyake, F. Aryasetiawan, T. Kotani, M. van Schil-fgaarde, M. Usuda, K. Terakura, Phys. Rev. B ,245103 (2002)48. J.M. Luttinger, J.C. Ward, Phys. Rev. , 1417(1960)49. A. Klein, Phys. Rev. , 950 (1961)50. Y.M. Niquet, M. Fuchs, X. Gonze, Phys. Rev. A ,032507 (2003)51. F. Caruso, D.R. Rohr, M. Hellgren, X. Ren, P. Rinke,A. Rubio, M. Scheffler, Phys. Rev. Lett. , 146403(2013)52. G. Baym, L.P. Kadanoff, Phys. Rev. , 287 (1961)53. G. Baym, Phys. Rev. , 1391 (1962)54. U. von Barth, N.E. Dahlen, R. van Leeuwen, G. Ste-fanucci, Phys. Rev. B , 235109 (2005)55. A. Szabo, N.S. Ostlund, Modern Quantum Chemistry:Introduction to Advanced Electronic Structure Theory (Dover Publiations, Inc., 31 East 2nd Street, Mineola,N.Y. 11501, 1989)56. P. L¨owdin, J. Chem. Phys. , 365 (1950)57. W. Ko(cid:32)los, L. Wolniewicz, J. Chem. Phys. , 2429(1965)58. W. Ko(cid:32)los, L. Wolniewicz, J. Chem. Phys. , 509(1966)59. W. Ko(cid:32)los, L. Wolniewicz, J. Chem. Phys. , 3228(1969)60. L. Wolniewicz, K. Dressler, J. Chem. Phys. , 3292(1985)61. R. Requist, O. Pankratov, Phys. Rev. B , 235121(2008)62. M. Fuchs, Y.M. Niquet, X. Gonze, K. Burke, J. Chem.Phys. , 094116 (2005)63. N. Rosen, Phys. Rev. , 2099 (1931)64. J.O. Hirschfelder, J.W. Linnett, J. Chem. Phys. ,130 (1950) 65. R.S. Mulliken, J. Chim. Phys.46