Maria Hellgren

Nicola Colonna

Stefano De Gironcoli

Beyond the random phase approximation with a local exchange vertex

With the aim of constructing an electronic structure approach that systematically goes beyond the GW and random phase approximation (RPA) we introduce a vertex correction based on the exact-exchange (EXX) potential of time-dependent density functional theory. The EXX vertex function is constrained to be local but is expected to capture similar physics as the Hartree-Fock vertex. With the EXX vertex we then unify different beyond-RPA approaches such as the various re-summations of RPA with exchange (RPAx) and the second order screened exchange (SOSEX) approximation. The theoretical analysis is supported by numerical studies on the hydrogen dimer and the electron gas, and we discuss the role of including the vertex correction in both the screened interaction and the self energy. Finally, we give details on our implementation within the planewave pseudo potential framework and demonstrate the excellent performance of the different RPAx methods in describing the energetics of hydrogen and van der Waals bonds.

I. INTRODUCTION

Kohn-Sham (KS) density functional theory (DFT) 1,2 has been extremely successful in predicting structural and vibrational properties of a wide range of systems. However, for sparse systems, or systems in reduced dimensions, the standard local or semi-local approximations are unreliable and sometimes fail qualitatively. Examples are the various silica polymorphs which require a proper description of the van der Waals forces, 3,4 or the layered materials which often exhibit effects of strong electron correlation. 5 Moreover, in reduced dimensions, excitonic effects become increasingly important and may be strong enough to influence ground-state properties. 6 Although highly accurate methods such as coupled cluster or quantum Monte Carlo (QMC) could overcome these limitations high computational cost limits them to the study of very selected systems. A cheaper solution to at least some of the problems is to include a fraction of Hartree-Fock (HF) exchange via the so-called hybrid functionals. 7,8 These functionals are well-known to improve the electronic structure of, e.g., transition metal compounds. 9 On the other hand, the hybrid functionals miss the van der Waals forces and rely on a parametrisation which often turns out to be system dependent.

Another route is to use many-body perturbation theory (MBPT) based on Green's functions which is a complete theoretical framework for studying both ground and excited state properties. A popular approximation to the many-body self-energy Σ is the GW approximation 10 (GWA) which simultaneously captures the exchange interaction and the van der Waals forces. 11,12 The GWA is widely used for band-structure calculations (i.e. quasiparticle energies) and total energy calculations are becoming feasible, at least within the random phase approximation (RPA), which can be seen as the variation-ally best GW total energy within a restricted space of KS Green's functions. [13][14][15] Approximations beyond the GW self-energy can be constructed using diagrammatic techniques by, for example, including the second order screened exchange (SOSEX) diagram. [16][17][18][19][20] Since, in general, it is not obvious which set of diagrams to include different approximation schemes have been developed. In the so-called Φ/Ψ-derivable schemes 21 diagrams are chosen such that only conserving approximations are generated, i.e., approximations that conserve energy, momentum and particle number. 22,23 Another approach is Hedin's iterative procedure, in which approximations are generated from a vertex function Γ = 1 + δΣ/δV where Σ = iGW Γ and V is the single particle potential. 10 By, e.g., starting from GW (i.e. Γ = 1) the self-energy can be updated iteratively. It is also possible to start from a simpler approximation to the self-energy, such as the HF approximation, generating a non-local vertex function similar to the one used in the Bethe-Salpeter equation.

In this work we study an approximate local HF vertex derived from the exact-exchange (EXX) potential of time-dependent density functional theory (TDDFT). The EXX potential is obtained by minimizing the HF energy in a constrained space of local potentials. [START_REF] Engel | A Primer in Density Functional Theory[END_REF] Local vertex corrections derived from LDA and GGA's have been studied in the context of quasi-particle energies, but their effect on the GW gaps was quite small. [START_REF] Del Sole | [END_REF][26][27] The HF or the local EXX vertex should more accurately capture the electron-hole interaction, expected to be important for correcting the deficiencies found within the selfconsistent GWA. 28 The main focus of this paper will be on the total energy and we will show how one can consistently generate different approximations beyond the RPA total energy using the EXX vertex. In this way, we establish theoretical connections between the different variants of RPA with exchange (RPAx), usually defined in terms of the density response function, [29][30][31][32][33][34][35][36][37][38][39] as well as between the RPAx and the SOSEX total energy approximation. The analysis also provides a many-body perspective on approximations usually defined within DFT.

The paper is organized as follows. In Sec. II we derive the basic equations for the adiabatic connection total energy formula within MBPT and DFT. We then define the vertex function, which establishes a connection between the two frameworks. We also review the comparison between GW and RPA. In Sec. III we define the local HF or EXX vertex and use it to derive previously defined approximations based on RPAx and a slightly new variant of SOSEX. In Sec. IV we present numerical results on H 2 and the electron gas to support the theoretical discussion in Sec. III. We also present results on the A24 test-set 40 to demonstrate the performance of the different RPAx methods in describing hydrogen and van der Waals bonds. In the end, we present numerical details from our implementation within the plane-wave and pseudo potential framework. Finally, in Sec. V, we present our conclusions.

II. CORRELATION ENERGY

We start by deriving the exact adiabatic connection formula for the correlation energy in terms of the singleparticle Green's function G. Similar derivations can be found in Refs. 41 and 42.

The Hamiltonian of the interacting electronic system is given by

Ĥ = T + Vext + Ŵ ( 1
)
where T is the kinetic energy operator, Vext is the external nuclear potential and Ŵ is Coulomb interaction between the electrons. The adiabatic connection path is defined in terms of a parameter λ ∈ [0, 1] that linearly scales the Coulomb interaction Ŵ from zero to full interaction strength. A single-particle potential V λ is added to the Hamiltonian such that the density is kept fixed to its interacting value (λ = 1) at every value of λ. Introducing the Hartree (H) and exchange-correlation (xc) potential VHxc , as defined within KS DFT, we thus write the scaled Hamiltonian on the adiabatic connection path as follows

Ĥλ = T + Vext + VHxc -V λ + λ Ŵ (2)
where

V λ=1 = VHxc , V λ=0 = 0. ( 3
)
With the Hamiltonian defined in Eq. ( 2) we derive the corresponding single-particle Green's function and write it in terms of a Dyson equation

G λ = G s + G s [Σ λ [G λ ] -V λ ]G λ (4)
where Σ λ is the irreducible self-energy and G s is the Green's function of the λ = 0 system (i.e. the KS Green's function). Using standard tricks [START_REF] Fetter | Quantum Theory of Many-Particle Systems[END_REF] we can write the total energy at full interaction strength as an integral over the coupling constant λ. Defining the single-particle energy

E s = T s + nv ext + E H (5)
the total energy is written as

E = E s - i 2 dω 2π dλ λ Tr{Σ λ [G λ ]G λ } (6)
where the trace is defined as Tr {AB} = drdr A(r, r )B(r , r) and the convergence factor (e iωη ) has been suppressed. This is an exact formula from which one can generate approximate total energies from any approximate Σ.

An alternative way to write the correlation energy is in terms of the density correlation function (or the reducible polarization propagator) χ

E = E s + i 2 dω 2π dλ Tr{v[χ λ -δn]} (7)
where δn = δ(r, r )n(r). In this form, approximations are generated diagrammatically via the Bethe-Salpeter equation, or via the linear response TDDFT Dyson equation.

In order to connect the two expressions, Eq. ( 6) and Eq. ( 7), we define the irreducible polarization propagator P in terms of the vertex function Γ

P = -iGGΓ, Γ = - δG -1 δV = 1 + δΣ δV (8)
where V = v ext + v H . With these quantities it is possible to rewrite the scaled density correlation function using the relation

χ λ v = 1 λ P λ W λ (9)
in which W = v + vP W is the screened Coulomb interaction. Since P is given by Eq. ( 8) and Σ = iGW Γ we can rewrite Eq. ( 7) in the form of Eq. ( 6). Within TDDFT the density correlation function is determined from the time-dependent KS system which is defined to reproduce the exact interacting density with a time-dependent local single-particle potential. 43 The TDDFT Dyson equation reads

χ λ = χ s + χ s [λv + f λ xc ]χ λ ( 10
)
where χ s = -iG s G s is the KS polarization propagator and f λ xc is the functional derivative of the timedependent xc potential with respect to the density. The irreducible polarization propagator P is then determined from P = χ s + χ s f xc P . But, this expression can easily be rearranged to P = -iG s G s Γ xc , that is, in terms of a local vertex function defined as

Γ xc = 1 + δv xc δV . (11)
In this way we can rewrite Eqs. ( 7) as

E = E s - i 2 dω 2π dλ λ Tr{Σ GW Γxc λ G s } (12)
where Σ GW Γxc = iG s W Γ xc , which can be compared to Eq. ( 6). Notice that if the exact xc kernel is used this is an exact expression for the correlation energy. This, however, does not mean that Σ GW Γxc is the exact selfenergy which should be defined in terms of the fully nonlocal vertex function. The extent to which a local vertex can approximate the exact vertex was studied in Ref. 27 and 44.

A. RPA and GW

The GWA and RPA have already been thoroughly compared 45,46 but for the sake of completeness we review this comparison here. Within the GWA we can calculate the total energy via Eq. ( 6) by simply setting

Σ GW λ = iG λ W λ , W λ = λv + λvP 0 W λ . ( 13
)
The screened interaction is calculated in the timedependent Hartree approximation (or RPA) for which P 0 = -iGG.

The RPA for the total energy is given by

E RPA = E s + i 2 dω 2π dλTr{v[χ RPA λ -δn]} (14)
where

χ RPA λ = χ s + λχ s vχ RPA λ . ( 15
)
Following the steps in the previous section [Eqs. ( 8)-( 12)] we see that this response function corresponds to Γ xc = 1. We can thus write

E RPA = E s - i 2 dω 2π dλ λ Tr{Σ GW λ [G s ]G s } (16)
and we see that this exactly corresponds to Eq. ( 6) in the GWA but with G λ replaced by G s . We notice that the diagrammatic structure within the GWA is preserved in the RPA. Furthermore, since only the explicit λ-dependence is retained, the λ-integral is easily carried out leading to the following expression

E RPA = E s - i 2
dω 2π Tr{ln[1 + ivG s G s ]}. ( 17
)
Since the GW self-energy is Φ-derivable, meaning that the self-energy can be generated from a functional Φ[G], i.e., Σ = δΦ[G]/δG, also the λ integral in Eq. ( 6) can be performed analytically. This leads to a total energy of the form This is a functional of G (also known as the Klein functional) which is stationary when G obeys the Dyson equation. If we replace G with G s it is easy to show that this expression is equivalent to Eq. ( 17). We thus have a second way to compare the GWA and RPA. 45,46 The RPA is the GW Klein energy constrained to the space of KS Green's functions.

E GW = - i 2 dω 2π Tr{ln[1 + ivGG]} + E H + i dω 2π Tr [GG -1 s -1 + ln(-G -

III. LOCAL HARTREE-FOCK VERTEX

We will now improve the GWA and RPA by constructing an approximate vertex function (Eq. ( 8)) from the HF self-energy (Σ HF = iGv)

Γ HF = 1 + δΣ HF δV . (19)
We note that since the self-energy depends on two arguments the vertex depends on three arguments (three space and time variables). The diagrammatic expansion of the HF vertex is given in Fig. 1. From the diagrammatic structure it is easy to see that the electronhole interaction is accounted for. If we replace the bare Coulomb interaction with a statically screened Coulomb interaction this vertex generates what is usually referred to as the Bethe-Salpeter equation for the polarization propagator.

An approximate local HF vertex can be defined in terms of the local time-dependent exact-exchange (EXX) potential v x

Γ x = 1 + δv x δV = 1 1 -Γ 1 x (20)
where

Γ 1 x = δv x δV s , V s = v ext + v Hx . (21)
Since the local approximation v x ∼ Σ HF depends only on one argument the corresponding local vertex function Γ x has only two arguments as opposed to three. The timedependent EXX potential is determined by the linearized Sham-Schlüter equation 47 and has been studied in several previous work. It can, e.g., be shown that v x exactly reproduce the HF density to first order in the Coulomb interaction. The equation for Γ 1 x is illustrated diagrammatically in Fig. 2. We note that the local EXX vertex reproduce the response function of the full HF vertex to first order in the Coulomb interaction. 48 Using the chain rule Eq. ( 21) also defines the standard EXX kernel

f x = δv x δn = Γ 1 x χ -1 s . (22)
With the EXX kernel we can generate the TDDFT density correlation function which is usually referred to as the RPAx response function

χ RPAx = χ s + χ s [v + f x ]χ RPAx . ( 23
)
This response function has been applied to atoms and molecules as well as to silicon for the calculation of excitons. 30,32,48,49 Although some excitation energies are sensitive to the local approximation 30 integrated quantities such as polarizabilities and van der Waals coefficients give results similar to the nonlocal time-dependent HF (TDHF) approach. 32 In the following we will use the EXX vertex to generate a set of approximations to the total energy.

A. GWΓx and GWΓ 1

HF

If we use the RPAx response function defined in Eq. ( 23) and insert it into the total energy of Eq. ( 7) we find what we have previously called the RPAx for the total energy

E RPAx = E s + i 2 dω 2π dλTr{v[χ RPAx λ -δn]}. ( 24
)
With the definition of the EXX vertex in Eq. ( 20) we can follow the same steps that lead to Eq. ( 12) and rewrite Eq. ( 24) as

E RPAx = E s - i 2 dω 2π dλ λ Tr{Σ GWΓx λ G s } (25)
with Σ GW Γx = iG s W Γ x . Since the vertex is local this self-energy is approximate but is expected to mimic the self-energy with the full HF vertex. If the latter is used, the approximation to the correlation energy corresponds to the RPAx defined in Ref. results of Refs. 29 and 50 to those of Ref. 51 we see similar trends. Only in the case of the Be dimer there is a qualitative difference and it is found that the use of a nonlocal potential for generating G s strongly improves the results. We notice that neither the HF nor the EXX vertex produces Φ-derivable self-energies due to the lack of symmetry.

If we expand the EXX vertex to first order only (Γ 1 x ) all terms in the expansion of the correlation energy have a diagrammatic representation. We can therefore write the total energy of the Γ 1

x -approximation [RPAx (1)] in terms of the full HF vertex to first order (Γ 1 HF ) plus a self-energy correction

E RPAx(1) = E s - i 2 dω 2π dλ λ Tr{Σ GWΓ 1 HF λ [G s ]G s }+E Σs . (26)
The diagrammatic representation of the first correlation term is given in Fig. 3(a). The second, self-energy correction, is of the form

E Σs = - i 2 dω 2π dλ Tr{G s W λ G 1 + G 1 W λ G s } (27)
with

G 1 = G s [Σ HF -v x ]G s . (28)
One of the terms in Eq. ( 27) is illustrated in Fig. 3(c).

Alternatively we can write the correlation energy of Eq. ( 26) as in Eq. (7), in terms of an approximate density correlation function obtained by expanding the irreducible polarizability to first order χ RPAx(1)

λ = P 1 λ + λP 1 λ vχ RPAx(1) λ ⇒ χ RPAx(1) λ = [1 -λP 1 λ v] -1 P 1 λ ( 29
)
where P 1 λ = χ s + λχ s Γ 1 x = χ s + λχ s f x χ s . This approximation was studied previously for the electron gas in Ref. 35 and for the total energy of molecules in Ref. 51. By keeping the vertex to first order only cured certain pathologies due to an unscreened Coulomb interaction. We note here that in Refs. 35 and 51 we employed an additional alternative re-summation for χ based on the first order exchange kernel, called t'RPAx, which is similar to the one proposed in Ref. 34. This approximation has been shown to produce results of similar quality as RPAx( 1), but since it does not have a counterpart in MBPT we will not study this approximation here.

B. SOSEX

In the RPAx(1) approximation defined above the vertex is treated to first order in both screened interaction and in the self-energy. If we set the vertex to one in the screened interaction (i.e. keeping it at the RPA level) we generate the so-called SOSEX approximation. For a diagrammatic representation see Fig. 3(b). The total energy is given by

E AC-SOSEX = E s - i 2 dω 2π dλ λ Tr{Σ SOSEX λ [G s ]G s } + ẼΣs (30)
where ẼΣ is a self-energy correction with a RPA screened interaction. This expression can also be exactly rewritten in terms of an approximate density correlation function given by

χ SOSEX λ = P 1 λ + λχ s vχ SOSEX λ ⇒ χ SOSEX λ = [1 -λχ s v] -1 P 1 λ . (31)
The SOSEX has been studied previously, both for total energies, and quasi-particle excitations of finite systems. [16][17][18][19] Total energies are obtained either from the adiabatic connection expression (Eq. ( 6)) or from the perturbative Galitskii-Migdal expression. It has been shown that so-called single-excitations improve the results 18 and our self-energy correction ẼΣs can be interpreted as an approximate single-excitation correction. In our formalism this correction naturally arises from the definition of the local EXX vertex. The so generated set of approximations RPAx, RPAx(1) and AC-SOSEX can thus be seen as a consistent way to do many-body approximations based on a local DFT framework. Self-energy corrections to the same order as the vertex naturally follows and compensate to some extent for a lack of selfconsistency within the MBPT scheme. Once the EXX vertex (or the EXX kernel) is calculated all approximations follow straightforwardly (see also Appendix A). In the next section we will compare these different approximations on H 2 and the electron gas and assess their performance in describing hydrogen and van der Waals bonds.

IV. NUMERICAL RESULTS

In this section we present numerical results on H 2 and the electron gas, in order to compare the different approximations defined in the previous section, and to compare to SOSEX results obtained by other groups. We also present calculations on the A24 test-set 40 containing molecules with hydrogen and van der Waals bonds. All calculations (except for those on the electron gas) have been done with a recently released RPA/RPAximplementation within the Quantum Espresso package. 52 See Sec. IVC, Appendix A and Refs. 35, 53, and 54 for details regarding the implementation.

A. H2 and the electron gas

In Fig. 4 we present H 2 dissociation curves within RPA, RPAx, RPAx(1) and AC-SOSEX, and compare the results to accurate results from Ref. 55. We also compare our AC-SOSEX to previous SOSEX results in the literature. 56 RPA, RPAx and RPAx(1) all include the vertex correction in both the screened interaction and the self energy (Γ = 1, Γ x and Γ 1 HF respectively). At the same time they are all able to capture the dissociation region, or static correlation, reasonably well. We note that the remaining "bump" seen in all approaches has been shown to disappear when employing a statically screened vertex. 57 Our AC-SOSEX results are essentially identical to previous results. SOSEX includes the vertex correction (Γ 1 HF ) in the self-energy only, while it treats the screened interaction at the RPA level. This inconsistency could explain why it fails to describe dissociation into open-shell atoms.

In the dissociation limit of H 2 the one-electron Hatoms are represented as a mixture between a spin-up and a spin-down state. Such an ensemble state requires a correlated method to be described. If, instead, the Hatom is isolated it can be treated as spin-polarized and HF gives the exact non-interacting result. The RPA is based on the Hartree response and, therefore, introduces self-interaction in the correlation functional. The RPAx naturally cures this error since it is based on the full HF response. 58 It can straightforwardly be proven that also AC-SOSEX is one-electron self-interaction free (see e.g. Ref. 59). RPAx(1) will, however, contain terms to 3rd order and beyond that are non-vanishing, but expected to be small.

In Fig. 5 we present the correlation energy per particle in the homogeneous electron gas (HEG). We compare our AC-SOSEX energies to previous SOSEX results by Freeman 16 and to accurate QMC results. 60,61 Again, we find that our AC-SOSEX agrees very well with previous SOSEX results. Both coincides with QMC results at a given density (AC-SOSEX at r s = 6.3 and SO- 11). A similar behavior has also been reported by the authors of Ref. 62 who observed imaginary frequencies eigenmodes for the polarization propagator computed by solving the Bethe Salpeter equation. Also in their case disregarding selected diagrams removes the instability in the low density regime.

1 2 3 4 R(H-H) [Å] -2.
In the inset of Fig. 5 we have plotted the static response functions at r s = 5, obtained from Eq. ( 23), Eq. ( 29) and Eq. (31). Again, RPAx(1) systematically improves upon RPA and has no pathological behavior as compared to RPAx. SOSEX strongly overestimates the static response up to its maximum suggesting that the good energies are subject to effects of error cancellation.

B. Dispersion interactions: The A24 test-set

In order to test the performance of the beyond-RPA methods in describing dispersion forces we calculated the binding energies of molecular dimers in the A24 test set. 40 The results are presented in Fig. 6 and in Table 1. The RPA tends to underestimate the binding energies with a Mean Absolute Error (MAE) of 0.44 kcal/mol. SOSEX has been shown to improve the description of dispersion interactions as compared to RPA and we find similar results in this work, with a MAE of 0.19 kcal/mol. Overall, the SOSEX energies are overestimated, in particular for the dimers with mixed bonds where the errors only slightly reduce with respect to RPA. RPAx and RPAx(1) give similar average performances (MAE of 0.14 kcal/mol and 0.13 kcal/mol, respectively) but RPAx (1) gives a more reliable improvement without exceptions.

Comparing to other beyond-RPA results in the literature we find that our approximations perform better than those of TDHF based approximations, unless a rangeseparation parameter is introduced. 38,50,63 As expected, our RPAx results are comparable to the RPAx results of Ref. 64.

C. Technical aspects

All molecular calculations of the ACFDT expression for the total energy were performed in a plane-wave and pseudo-potential formalism as implemented in a separate module of the Quantum ESPRESSO distribution. 52,65 We devote this section to review the basic aspects of our implementation of the exchange and correlation energy at all the levels of the theory discussed in Sec. II and Sec. III (for a more detailed description we refer to Refs. 35, 53, and 54). We complement the discussion with convergence tests on the relevant parameters of our implementation for a subset composed of three complexes, namely the water-ammonia, the methane-HF and the borane-methane, each one representative for a different type of bond in the A24 set.

All the ACFDT calculations have been performed in a post-DFT fashion. The KS single-particle wavefunc- 1.

tions and energies needed as input for the ACFDT total energy calculation have been calculated at the PBE 66 level using an energy cut-off of 80 Ry; Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials [67][68][69] have been used to model the electron-ion interaction. The molecules have been placed in an orthorhombic cell with 12 A of vacuum in each direction, sufficient to suppress the spurious interactions between periodic replica. As illustrated on the first and second row of Fig. 7, with this setup all the components of the binding energy (BE) converge within 0.05 kcal/mol. A slightly larger deviation is seen for the RPAx correlation BE with respects to energy cut-off (identical in all beyond-RPA methods). We note, however, that this is not due to a worse convergence with respect to the energy cut-off, but to the numerical precision associated with the EXX optimized-effective potential problem, needed for the evaluation of the RPAx correlation energy. 35 The error bars on the final beyond-RPA binding energies are estimated to be < 3% for the A24 test-set, corresponding to 0.006 kcal/mol per electron. Half of this value is estimated for the RPA binding energies. The exact xc-energy [second term on the rhs of Eq. ( 7)] can be further separated into the KS EXX energy

E x = - 1 2 drdr | occ i φ i (r)φ * i (r )| 2 |r -r | ( 32
)
and the ACFDT correlation energy determined by the difference between the interacting and non-interacting KS response functions

E c = - 1 0 dλ ∞ 0 du 2π Tr {v [χ λ (iu) -χ s (iu)]} . ( 33
)
The integrable divergence appearing in a plane-wave implementation of the EXX energy would lead to a slow convergence with respect to the size of the supercell. In this work this issue has been dealt with using the method proposed by Gygi and Baldereschi 70 plus the extrapolation scheme of Nguyen and de Gironcoli, 53 although other strategies are also possible. 71,72 For all cases analyzed, this correction scheme allows us to converge the EXX contribution to the BE within 0.05 kcal/mol with 10 A of vacuum separating the molecular replicas. The strategy we use to solve the Dyson equation for χ λ and efficiently compute the trace over the spatial coordinates in Eq. ( 33) is a generalization of the projective dielectric eigenpotential method 53,[73][74][75] and requires the solution of a well defined generalized eigenvalue problem (GEP) involving the non-interacting response function χ s and its first order correction in the limit of vanishing coupling-constant H Hx = χ s (v + f x )χ s . The set of eigenvectors of the GEP defines an optimal basis set on which i) the response matrices have a compact representation and ii) the Dyson equation has a straightforward solution, as detailed in Appendix A. Most importantly, only the lowest-lying eigenvectors and eigenvalues of the GEP at hand are relevant for the correlation energy calculation, and an efficient iterative diagonalization techniques can be used to compute the first N eig lowest-lying eigen- . Convergence of the binding energy, at all levels of theory presented in the text, with respect to kinetic energy cut-off (first row), vacuum (second row), number of eigenvalues Neig (third row) and number of frequencies Niu (last row) for three representative molecules in the A24 test-set. In each plot and for each approximation, the zero is set by the most converged calculation.

values/eigenvectors of the problem. Moreover, the matrix elements of χ s and of H x = χ s f x χ s , needed for the solution of the GEPs, can be efficiently computed 35,53 resorting to the linear response machinery of densityfunctional perturbation theory 76 (DFPT) thus avoiding explicit summations over empty states.

Once the GEP is solved, the trace over spatial coordi-nates can be written as a sum over a simple function E c of the GEP eigenvalue e α (and eventually of the diagonal matrix elements of χ s and H x ; see Appendix A for the actual expression of E c ), and a compact expression for the correlation energy is obtained:

E c = - 1 0 dλ ∞ 0 du 2π Neig α=1 E c [e α (iu, λ), λ]. (34)
Here the frequency integral has been moved from the real to the imaginary axis of the complex plane 77 to take advantage of the smooth behavior of the response functions for an efficient numerical quadrature (see below). The integration over the coupling constant is analytic for all the approximations but the RPAx(1) for which the eigenvalues e

RPAx(1) α inherit a non trivial dependence on the coupling constant from the λ-dependent RPAx(1) GEP (see Appendix A 3). In this case a numerical integration can be efficiently performed using a Gauss-Legendre quadrature. Usually very few points are needed to converge the integral within 0.01 kcal/mol (see Fig. 8). For all the other approximations the GEPs do not depend on λ and only the explicit dependence of E c on the coupling constant has to be considered. In these cases a straightforward analytic integration is always possible.

The number of eigenvalues N eig set the dimension of the optimal basis set used to represent the response functions, and, therefore, it is a convergence parameter. The third row of Fig. 7 shows the dependence of the correlation binding energy on N eig for the selected subset of molecules. In all the cases and for all the approximations studied, a relatively small number of eigenvalues, never larger than 20 times the number of electrons in the system at hand, is sufficient to converge the ACFDT correlation energy contribution to the BE within 0.05 kcal/mol.

Finally, the integral over the imaginary frequency iu can be efficiently performed using a Gauss-Legendre quadrature. A standard mesh of point z i ∈ [0, π/2] is mapped on an imaginary-frequency grid between [0, +∞) using the transformation u i = (ε HL -cz i ) tan(z i ) where ε HL is the DFT HOMO-LUMO gap of the systems and the parameter c is set by the upper integration limit (u max = 200 Ry). This transformation takes into account the typical dependence of the response functions on the imaginary frequency and the fact that at small frequency the scale of the excitation energies is given by the HOMO-LUMO gap. With this strategy a grid with N iu = 10 points is usually sufficient to converge the ACFDT energy within 0.05 kcal/mol as illustrated on the last row of Fig. 7.

V. CONCLUSIONS

In this work we have defined a local EXX vertex based on the EXX kernel of TDDFT. This allowed us to unify different beyond-RPA approaches such as the various resummations of RPAx and the SOSEX approximation, as well as to give a many-body perspective on approximations normally defined within DFT.

We have tested the theory on the H 2 molecule and the electron gas and we find that our AC-SOSEX agrees well with previous results in the literature. Although AC-SOSEX gives excellent total energies for the electron gas, it cannot be well represented as an approximation to the density response function. This could explain the poor performance in the dissociation region of H 2 . On the contrary, approximations that incorporate the vertex correction in both the screened interaction and in the self-energy (RPA, RPAx and RPAx(1)) dissociate H 2 correctly.

In terms of an overall performance we find that the RPAx(1) gives the most reliable results. It is as accurate as RPAx but shows no pathologies for any of the systems studied so far. High quality results within a reasonable computational cost are obtained on the A24 test-set, paving the way for accurate ab-initio treatment of sparse systems.

(iu) -χ s (iu)]} = i e i (iu) 1 -λe i (iu)
-e i (iu) .

(A2) The integration over the coupling constant can be performed analytically and the final result for the RPA correlation energy becomes:

E RPA c = 1 2π ∞ 0 du i [ln[1 -e i (iu)] + e i (iu)] . (A3)
The spectrum of the GEP in Eq. (A1) is bounded from above by zero, and only a small number of the lowest lying eigenvalues are significantly different from zero. 74,78 This implies that only a small fraction of the {e i } contributes significantly to the correlation energy, and an efficient iterative diagonalization scheme can be used to evaluate those. 53,[START_REF] Nguyen | Efficient calculation of RPA correlation energies in the Adiabatic-Connection Fluctuation-Dissipation theory[END_REF] The basic operation involved in the iterative solution of the GEP in Eq. (A1) is the calculation of the noninteracting response to a trial potential, and this is done resorting to the linear-response techniques of density functional perturbation theory 76 (DFPT), generalized to imaginary frequencies.

The RPAx correlation energy

A strategy analogous to the one adopted for the RPA problem is used to compute the RPAx correlation energy. The GEP to solve in this case is the following

-χ s [v + f x ] χ s |z RPAx α = e RPAx α [-χ s ]|z RPAx α . (A4)
The expression for H x = χ s f x χ s in terms of KS orbitals and eigenvalues has been obtained by Görling, [START_REF] Görling | [END_REF]81 and Hellgren and von Barth. 48 In this case the diagonalization was carried out in the basis set of the eigenvectors of the RPA GEP, Eq. (A1). This was done in order to i) avoid possible instabilities that may occur in the inversion of the non-interacting response function (i.e. the overlap matrix in the GEP) and ii) to speed up the calculations. On this basis set χ s and H H = χ s vχ s are diagonal and their matrix elements are readily available in terms of the eigenvalues of the RPA GEP (Eq. (A1)), i.e. χ ij s = w i |χ s |w j = δ ij e i and H ij H = w i |χ s vχ s |w j = δ ij e 2 i . The only additional operation required to solve the RPAx problem is the evaluation of the matrix elements H ij x = w i |χ s f x χ s |w j . These can be efficiently computed using DFPT as detailed in Ref. 35 At variance with previous cases, for the RPAx(1) problem the integration over the coupling constant has to be performed numerically because of the non-trivial dependence of the eigenvalues and eigenvectors on coupling constant inherited from the λ-dependent GEP Eq. (A7). However the integrand is a smooth function and a Gauss-Legendre quadrature with less than five points ensures converged results within 0.01 kcal/mol (see Fig. 8). We stress here that the most computationally intensive part of the whole procedure is the evaluation of the response functions matrix elements. Once χ s (iu) and H Hx (iu) have been represented, only linear algebra operations with matrices of dimension N eig × N eig are needed to solve the GEP for each λ on the grids, meaning that the additional cost associated to the numerical coupling-constant integration is only a tiny fraction of the total computational cost.

The AC-SOSEX correlation energy

For the AC-SOSEX approximation the RPA eigenvectors {w i } from Eq. (A1) can be used to explicitly solve the SOSEX dyson equation [Eq. (31)] and compute the action of the SOSEX response function

Figure 1 .

 1 Figure 1. Diagrammatic expansion of the exchange vertex in the HF approximation. Wiggly lines correspond to the bare Coulomb interaction v.

Figure 2 .

 2 Figure 2. Equation to obtain the first order local exchange vertex Γ 1x .

Figure 3 .

 3 Figure 3. Diagrammatic representation of the GW Γ 1 HF (a) and SOSEX (b) self-energies. The self-energy correction in Eq. (27) is illustrated in (c).

Figure 5 .

 5 Figure 5. The energy per particle of the homogeneous electron gas within different beyond-RPA methods compared to QMC results of Ref. 61 and to SOSEX results of Ref. 16. Symbol and line coding as in Fig. 4. The inset shows the corresponding static response functions at rs = 5.

Figure 7

 7 Figure 7. Convergence of the binding energy, at all levels of theory presented in the text, with respect to kinetic energy cut-off (first row), vacuum (second row), number of eigenvalues Neig (third row) and number of frequencies Niu (last row) for three representative molecules in the A24 test-set. In each plot and for each approximation, the zero is set by the most converged calculation.

  and 82. Once the solution of the GEP in Eq. (A4) is available the action of vχ RPAx λ on the eigenvectors {z RPAx α } can be explicitly expressed making use of the RPAx Dyson equation Eq. (23), i.e. vχ RPAx λ |z RPAx α = vχ s |z RPAx α /(1-λe RPAx α ), and the trace over spatial coordinates in the RPAx correlation energy can be written as Tr {v[χ RPAx λ (iu) -χ s (iu)]} solved in the basis set of the RPA eigenvectors [Eq. (A1)]. Also in this case the only additional information needed is the representation of the H x operator on this basis (this may be already available if an RPAx calculation has been previously performed). The solution of the RPAx(1) GEP together with the RPAx(1) Dyson equation [Eq.(29)] allows to write explicitly vχ

Figure 8 . 2 i×

 82 Figure 8. Convergence of the RPAx(1) binding energy with respect to the number of points in the λ grid for three different complexes.

  Figure 4. Dissociation curves for the hydrogen molecule within different beyond-RPA methods. Approximations derived in this work are compared to 'exact' results presented in Ref. 55 and to SOSEX results obtained by Ren et al. in Ref. 56 .

		H 2
		-1.8
	Total Energy (Ry)	-2.2 -2.0	Accurate AC-SOSEX RPAx(1) RPAx
			RPA
		4	SOSEX (Lit.)

SEX at r s = 4.8). Overall, SOSEX produces energies in very good agreement with QMC even up to low densities. RPAx(1) represents a systematic improvement upon RPA, performing very well in the metallic range (r s < 5) although slightly worse than RPAx. For a more detailed discussion of RPAx and RPAx(1), see Ref.

35

. But, we mention here that expanding the EXX vertex to first order only, i.e. going from RPAx to RPAx(1), removes an instability of the response function at low electronic densities

35

(r s >

Table I .

 I Binding energies (in kcal/mol) for the A24 test-set as obtained with the RPA, RPAx, RPAx(1) and AC-SOSEX methods Compared with CCSDT(Q) reference values40

			Symm	RPA	RPAx	RPAx(1)	AC-SOSEX	CCSDT(Q)
		Hydrogen-bonded systems
	01	water-ammonia	Cs	-5.71	-6.20	-6.14	-6.51	-6.492
	02	water dimer	Cs	-4.27	-4.89	-4.82	-5.16	-4.994
	03	HCN dimer	Cs	-4.20	-4.72	-4.50	-4.94	-4.738
	04	HF dimer	Cs	-3.86	-4.47	-4.42	-4.78	-4.564
	05	ammonia dimer	C 2h	-2.65	-2.92	-3.00	-3.12	-3.141
		ME		0.64	0.14	0.21	-0.12	-
		MAE		0.64	0.14	0.21	0.13	-
		MA%E		13.8%	3.2%	4.3%	2.7%	-
		Mixed-type systems
	06	methane-HF	C3v	-1.33	-1.70	-1.71	-1.91	-1.660
	07	ammonia-methane	C3v	-0.54	-0.86	-0.86	-0.92	-0.771
	08	methane-water	Cs	-0.54	-0.76	-0.76	-0.82	-0.665
	09	formaldehyde dimer	Cs	-2.78	-3.82	-4.33	-5.23	-4.479
	10	ethene-water	Cs	-2.11	-2.57	-2.59	-2.92	-2.564
	11	ethene-formaldehyde	Cs	-1.32	-1.78	-1.73	-1.94	-1.623
	12	ethyne dimer	C2v	-1.28	-1.66	-1.59	-1.79	-1.529
	13	ethene-ammonia	Cs	-1.05	-1.43	-1.44	-1.58	-1.382
	14	ethene-dimer	C2v	-0.71	-1.13	-1.15	-1.24	-1.106
	15	methane-ethene	Cs	-0.34	-0.64	-0.61	-0.66	-0.509
		ME		0.43	-0.01	-0.05	-0.27	-
		MAE		0.43	0.14	0.08	0.27	-
		MA%E		25.2%	9.3%	7.2%	18%	-
		Dispersion-dominated bonds
	16	borane-methane	Cs	-0.57	-1.23	-1.39	-1.58	-1.513
	17	methane-ethane	Cs	-0.58	-1.02	-1.00	-1.03	-0.836
	18	methane-ethane	C3	-0.44	-0.78	-0.75	-0.78	-0.614
	19	methane-dimer	C 3d	-0.38	-0.53	-0.52	-0.52	-0.539
	20	methane-Ar	C3v	-0.32	-0.42	-0.41	-0.42	-0.405
	21	ethene-Ar	C2v	-0.27	-0.35	-0.35	-0.35	-0.365
	22	ethene-ethyne	C2v	+1.00	+1.04	+0.96	+1.15	+0.794
	23	ethene dimer	D 2h	+1.11	+1.18	+1.09	+1.35	+0.909
	24	ethyne dimer	D 2h	+1.30	+1.30	+1.20	+1.38	+1.084
		ME		0.26	0.08	0.03	0.08	-
		MAE		0.26	0.15	0.10	0.17	-
		MA%E		29.5%	17.4%	12.3%	20.1%	-
	RPA correlation energy can be written as
	Tr {v[χ RPA λ

ACKNOWLEDGEMENT

We acknowledge IDRIS, CINES, TGCC and the computing center of EPFL for high performance computing resources.

Appendix A: The auxiliary basis sets

We provide in this section additional information on how the auxiliary basis set, used to solve the Dyson equation and calculate the trace over spatial coordinates, is computed for each approximation described in the text.

The RPA correlation energy

The evaluation of the RPA correlation energy is based on an eigenvalue decomposition of the noninteracting response function and has been carefully detailed in Ref. 53; we summarize here the main features of the RPA implementation because the solution of the RPA problem is a prerequisite for all the exchange-corrected approximations.

For each point on the imaginary-frequency grid the generalized eigenvalue problem (GEP)

is solved (the dependence on the imaginary frequency is implicitly assumed). Once the solution of the GEP is available the RPA Dyson equation, Eq. ( 16), can be readily solved since {w i } are also eigenvectors of vχ RPA λ with eigenvalues e RPA i = e i /(1 -λe i ). The trace over spatial coordinates appearing in the expression for the


Maria Hellgren

Jacopo Baima

Raffaello Bianco

Matteo Calandra

Francesco Mauri

Ludger Wirtz

Critical role of the exchange interaction for the electronic structure and charge-density-wave formation in TiSe 2

We show that the inclusion of screened exchange via hybrid functionals provides a unified description of the electronic and vibrational properties of TiSe2. In contrast to local approximations in density functional theory, the explicit inclusion of exact, non-local exchange captures the effects of the electron-electron interaction needed to both separate the Ti-d states from the Se-p states and stabilize the charge-density-wave (CDW) (or low-T) phase via the formation of a p -d hybridized state. We further show that this leads to an enhanced electron-phonon coupling that can drive the transition even if a small gap opens in the high-T phase. Finally, we demonstrate that the hybrid functionals can generate a CDW phase where the electronic bands, the geometry, and the phonon frequencies are in agreement with experiments.

The charge density wave (CDW) instability is a common phenomenon in layered semi-metallic transition metal dichalcogenides (TMDs) [1] and has attracted considerable interest over the years, from both the experimental and the theoretical side. The CDW phase is often found to compete with superconductivity and thus plays a similar role as the anti-ferromagnetic phase in strongly correlated heavy fermion systems or in high-T c cuprates [2][3][4][5]. This intriguing similarity has stimulated the search for a better understanding of the physical mechanism behind the CDW instability in TMDs [6].

The CDW instability in TiSe 2 is one of the most studied and debated. On the experimental side, neutron diffraction [7] and X-ray scattering [8] have established the existence of a commensurate 2×2×2 structural transition at 200 K. This is confirmed by angle resolved photo emission spectroscopy (ARPES) as well as by transport measurements [7], where an abrupt increase in resistivity is found at the same temperature. However, upon further cooling the resistivity reaches an anomalous maximum [7], after which a weak metallic behaviour is observed. By contrast, ARPES finds an insulating low-T phase, with a gap of approximately 0.15 eV [9][10][11][12][13][14]. In the high-T phase ARPES has not been able to conclude whether the system is semi-metallic or semi-conducting due to the very small indirect (possibly negative) gap. Theoretically, this fact makes TiSe 2 an ideal candidate to exhibit an excitonic insulator phase [9,[15][16][17] for which the CDW transition is driven by a purely electronic instability. Some recent experiments [18] partly support this scenario. On the other hand, excitonic correlations alone are insufficient as demonstrated in Ref. 19.

Additional experimental evidence for the CDW instability has been provided by vibrational spectra as a func-tion of temperature. A complete softening of an optical phonon at the L-point has been observed in inelastic Xray scattering experiments [20]. In Raman and infrared (IR) spectroscopy the transition is detected by the appearance of a large number of new modes [21,22], some of which can be related to the CDW transition due to their strong temperature dependence.

On the theoretical side, density functional theory (DFT) within the standard local density approximation (LDA) or the semi-local PBE approximation predict a structural instability at the correct wave vector [20,23,24] emphasising the role of the lattice distortion. However, orbital occupations obtained from the electronic band structure are in disagreement with ARPES measurements and CDW phonon frequencies strongly deviate from experiments. This suggests that electronic correlations cannot be neglected and may be crucial to describe the CDW instability. It has been shown that the inclusion of nonlocal exchange already gives a much better description of the electronic bands in both the highand low-T phase [13,14,[25][26][27]. In Ref. 24 it was further shown that also the DFT+U approach improves the electronic band structure, bringing it into good agreement with ARPES spectra. This seeming improvement, however, was accompanied by a complete loss of the CDW instability.

In this work we provide a unified description of both the electronic band structure and the lattice dynamics of TiSe 2 , and give a physical explanation for the CDW instability based on first-principle calculations. By an inclusion of exact-exchange within the hybrid functionals we capture the strong Coulomb repulsion due to the localised Ti-d-states, thus correcting the orbital occupations. We also capture the long-range exchange interac- We show also the results obtained by varying the α and µ parameters such that the band structure around conduction band minimum and valence band maximum is the same as the one of HSE06 and PBE+U.

tion which enhances the electron-phonon coupling (EPC) and stabilize the CDW phase through the formation of a Ti-d-Se-p hybridized state.

TiSe 2 is a layered material which in the high-T phase (>200 K) crystallises in the space group P 3m1. The CDW phase is characterised by a 2×2×2 superstructure with a different space group P 3c1. The PBE phonon dispersion of the high-T phase was analysed in Refs. 23 and 24. The distortion pattern d 3L (see Eq. 9 in Ref. 24) associated with the imaginary phonon frequencies at the three equivalent L points is of A u symmetry and its projection on a single TiSe 2 layer is shown in the inset of Fig. 1(b). It is completely identified by its symmetry, the magnitude of the displacement δTi and the ratio δTi/δSe. The layers are held together by van der Waals forces and in the CDW phase adjacent layers are rotated by 60 degrees. In all calculations we use the experimentally determined a and c lattice parameters. This choice indirectly accounts for the van der Waals interaction due to the fixing of the interlayer distance and does not sig-nificantly change the in-plane structure as compared to the theoretical one [24,28]. Although the spin-orbit coupling (SOC) splits the Se-p energy levels [27] it only gives a small contribution to the energetics [28]. Due to the increased computational cost we have therefore omitted SOC in calculations presented here. Local functionals in DFT such as PBE give rise to an excess of electron occupation on the Ti-d-states compared to ARPES, mainly due to the large p zd band overlap but also due to the overestimated p x,yd band overlap (see black curve in Fig. 1 (a)). To cure this problem it is necessary to include the Coulomb repulsion or 'U' due to the localised Ti-d states. A simple corrective procedure to incorporate this effect is the DFT+U approach [29]. Indeed, in Ref. 24 it was shown that the pd band overlap reduced significantly (see orange curve Fig. 1 (a)). However, it was also found that the energy gain in the distortion montonically decreased by increasing 'U' [28]. The hybrid functionals contain a fraction of exact-exchange and therefore naturally incorporate the repulsive 'U' interaction [30]. In addition, they exhibit the long-range exchange interaction, which generates the attractive electron-hole interaction at the linear response level [31,32]. The description of screening is rather simplified using a pre-factor α and a range-separated Coulomb potential defined by µ. These parameters can be fitted in various ways but often the HSE06 parameters (α = 0.25, µ = 0.2 Å-1 ) give good results [33]. In this work we have used the freedom of varying (α, µ) to reveal the critical role played by exchange, and to shed light on the mechanism that drives the CDW transition in TiSe 2 .

In Fig. 1 (a) we plot the results for the band structure in the high-T phase with PBE, PBE+U and HSE06 [28, [35][36][37][38][39]. In Fig. 1 (b) we plot the corresponding energy gain as a function of δTi by displacing the atoms according to the d 3L CDW pattern (using the experimentally determined ratio δTi/δSe = 3 [7]). We found HSE06 to give an enhanced negative curvature and energy gain of approximately a factor of 2 compared to PBE. This difference is mainly attributed to an overscreening at the PBE level, i.e., due to the overestimation of the metallicity of the system. However, this is not the complete picture. We see that PBE+U and HSE06 have almost identical orbital occupations and gaps -both in rather good agreement with experiment [13, 14, 28]but in PBE+U the instability is absent. This suggest that the long-range (or nonlocal) Coulomb exchange interaction is important for the formation of the CDW phase. That this is, indeed, the case is demonstrated explicitly in Figure 1 (b). We have calculated the energy gain using a set of (α, µ) parameters all tuned to give similar band-structures to HSE06 and PBE+U. By reducing the range, i.e., by increasing µ, thus approaching the PBE+U limit, the instability quickly vanishes (green curves). Instead, with an infinite range (HSE(17,0)) the energy gain increases substantially as does the amplitude of the distortion. After relaxing the atomic positions in the CDW phase, we find δTi = 0.055, 0.061, 0.082 Å with PBE, HSE06 and HSE(17,0) respectively. In experiment, δTi = 0.085 ± 0.014 Å [7], thus supporting the long range approximation. Based on the above analysis we will continue to discuss the mechanism for the CDW transition, but first we analyze the properties of the CDW phase.

The electronic band structure in comparison with ARPES suggests that the CDW phase is accurately captured with hybrid functionals [14,28]. However, as seen above, a good description of the band structure alone is not sufficient to judge the quality of the XC approximation. The vibrational frequencies are proportional to the second derivative of the energy and thus provide an additional stringent test for the validity of different approximations. We have calculated the harmonic Raman and IR active phonon frequencies in the CDW phase within both HSE06 and the infinite range HSE(17,0) using a finite difference approach as implemented in the CRYSTAL code [40].

In Fig. 2 we compare the theoretical frequencies with low-T experimental Raman [34] and IR [21] frequencies (at 11 K and 20 K, respectively). We analyze the spectra of the CDW phase by back-folding of the phonons at A, L and M of the high-T phase Brillouin zone (see HSE06 phonon dispersion [41] in Fig. 2 (a)) onto the Γ point of the low-T phase Brillouin zone [28]. In the following, the high-symmetry points of the 2×2×2 supercell are marked with a bar. For all three approximations, we find an overall good agreement with the experimental results. However, there are three notable exceptions of particular interest for this work: these are the modes related to the instabilities of the high-T structure (marked by vertical dashed lines). The three-fold degenerate imaginary phonon frequency at L splits into one A 1g and a two-fold degenerate E g mode, which are both Raman active. These have been previously identi-fied as being related to the CDW transition from their experimental temperature dependence [21,22]. In fact, the A 1g (amplitude) mode corresponds to the oscillations of the CDW order parameter. The imaginary phonon frequency at M similarly gives rise to an E u IR active phonon and to one A 1u inactive mode. The frequencies of these modes are systematically underestimated by PBE (20% and 30% for the A 1g and E g modes, respectively). HSE06 brings these frequencies in much better agreement with experiment (now underestimated by 4% and 7%, respectively). HSE(17,0), which gave the best geometry and overall phonon frequencies, also gives excellent CDW mode frequencies (overestimated by 5% and 2%, respectively). Since the frequencies of these modes depend very sensitively on the CDW potential energy surface (compare Fig. 1 (b)), the reproduction of these modes by the hybrid functionals is an important confirmation that they reliably reproduce the electronic structure of TiSe 2 .

With a reliable description of the vibrational and electronic properties, we can now analyze in detail the physical mechanism that is responsible for the CDW instability in TiSe 2 . Figure 3 (a) shows the HSE06 band structure in both the low-T and high-T phase (the latter one is back-folded into the Brillouin zone that corresponds to the 2×2×2 supercell of the low-T phase). At the CDW distortion, the three-fold degenerate Ti-d band splits into a band with dominant d z 2 character derived from the Ti atom in the supercell that does not move with the distortion, and a two-fold degenerate band with Ti-d and Se-p hybridization in a small region around Γ due to the interaction between the p and d states. Similarly, the Se-p bands move to lower energy and hybridize with Tid states around Γ. Such pd hybridization, discussed also in Refs. 11, 42, and 43, is observed by studying the site-projected orbitals onto spherical harmonics before and after the transition. Panel (b) presents a zoom of the region around Γ, showing the effect of a small d 3L distortion on the electronic structure. We see that the distortion only changes the pd gap while leaving the rest of the electronic structure unaffected. We can therefore assume that it is mainly the pd states around Γ which are involved in the transition and extract the EPC, or, more precisely, the deformation potential D for the different functionals. The deformation potential is defined as D = Γd|δV /δd 3L | Γp , where the variation of the total self-consistent potential V is taken with respect to the magnitude d 3L of the d 3L distortion. The states | Γd and | Γp with energy difference E 0 are obtained from the Ti-d and Se-p states at L and Γ in the undistorted unit cell by back-folding them onto Γ in the supercell. Without loss of generality, all quantities can be taken as real in the supercell. The p and d states are superimposed and illustrated in Fig. 3 (c) together with the pd hybridized state in the CDW phase. We also note from first order perturbation theory that D = δ( Γd)/δd 3L | Γp × E 0 and hence (at fixed E 0 ) D measures the overlap between the variation of the d state with the distortion and the undistorted p state. We then obtain D using a finite difference approach [44]. At small but finite d 3L we construct the 2 × 2 submatrix of the Hamiltonian H:

H dd = -E 0 /2 + O(d 2 3L
),

H pp = E 0 /2 + O(d 2 3L
) and

H pd = H dp = Dd 3L + O(d 2 3L ) and diagonalize it to obtain D 2 = 1 4 (E 2 -E 2 0 )/d 2 3L
, where E is the pd gap upon the distortion and we have discarded the higher order terms in d 3L . The results obtained from this procedure can be found in Fig. 3 (d). Although some coupling already exists between p and d states at the PBE level, the deformation potential (and so the EPC) is strongly enhanced via the nonlocal exchange potential in V provided by the inclusion of longrange exchange. Indeed, a comparison between PBE and PBE+U shows that U has no effect on the EPC. For the approximations with the same band structure we see the same trend in Figs. 1 (b) and 3 (d) and hence we conclude that it is the nonlocal exchange interaction that determines the strength of the instability. The inclusion of SOC splits and slightly changes the dispersion of the Se-p bands close to the Fermi level but has no qualitative and only a small quantitatve effect on the EPC [28]. Our result goes in line with some previous findings that the EPC can be strongly enhanced by nonlocal exchange, notably in weakly doped 2D materials [44,[46][47][48] and some high-T c superconductors [49]. Although the band structure in HSE06, which we used as a reference, has a small negative pd gap (-0.1 eV) in the high-T phase, the instability does not crucially depend on the existence of a Fermi surface. In fact, due to the enhanced EPC, the CDW phase can exist even if the gap in the high-T phase is positive. This is easily demonstrated by increasing the amount of exact-exchange (α). With µ = 0.2 Å-1 fixed, the instability persists up to α = 0.35 when the gap is as large as 0.2 eV. At larger values of the gap (or α) the hybridization of Ti-d -and Se-p x,y bands is suppressed and the instability disappears. If we instead decrease α < 0.25 the energy gain reaches a maximum around α = 0.15. At smaller values of α the CDW phase becomes metallic and the energy gain drops. This unusual behaviour as a function of α (or the gap) shares some features with the excitonic insulator transition (EIT) first proposed in the works of Kohn et al. [15,50]. The EIT arises in a two band model (similar to the p and d bands in TiSe 2 ) when tuned until the value of the pd gap is smaller than the exciton binding energy. Since the symmetry of the lattice is kept fixed the EIT is a purely electronic effect. In our calculations it was not possible to generate a symmetry lowering of the electronic density alone. The mechanism that we found relies instead on a simultaneous lattice distortion and strong EPC, enhanced by the electron-hole interaction. If the system is still close to a true excitonic instability or influenced by excitonic effects can only be determined by studying the excitation spectrum, but this we leave for future work.

M.H. and L.W. gratefully acknowledge financial support from the FNR (projects RPA-phonon and IN-TER/ANR/13/20/NANOTMD). The authors acknowledge support from the European Union Horizon 2020 research and innovation programme under Grant agree-ment No. 696656-GrapheneCore1. and by Agence Nationale de la Recherche under the reference n. ANR-13-IS10-0003-01 and by Prace. Calculations were performed at IDRIS, CINES, CEA and BSC TGCC. The authors are grateful to K. Rossnagel, for having supplied to them the image data of ARPES experiment.

Figure 1 .

 1 Figure 1. The upper panel (a) shows the band structure in the high-T phase in PBE (black solid line), HSE06 (green dashdotted line and circles), HSE(17,0) (blue solid line and circles) and PBE+U (U=3.5 eV) (orange solid line). The lower panel (b) shows the corresponding energy gain as a function of δTi.We show also the results obtained by varying the α and µ parameters such that the band structure around conduction band minimum and valence band maximum is the same as the one of HSE06 and PBE+U.

Figure 2 .

 2 Figure 2. (a) Phonon dispersion of the high-T phase (HSE06). (b-c) Experimental[21,34] versus calculated phonon frequencies of Raman and IR modes in the low-T phase. Modes related to the CDW distortion are marked by vertical dashed lines. The labels refer to the largest components of the phonon eigenvector with respect to the eigenvectors of the high-T structure.

2 /D 2 PBEFigure 3 .

 223 Figure 3. (a) HSE06 band-structure of undistorted (high-T) and CDW (low-T) phases of TiSe2. Panel (b) shows a zoom around the Fermi level and the EPC between the p-d states is calculated in panel (d)). The value for D 2 PBE = 8 eV 2 / Å2 (calculated with respect to δTi displacement). Panel (c) demonstrates the mixing (hybridization) of these orbitals [45] in the CDW phase.


Sebastien François

Rumen Andonov
email: randonov@irisa.fr

Hristo Djidjev
email: djidjev@lanl.gov

Metodi Traikov
email: metodi.gt@gmail.com

Nicola Yanev

Mixed Integer Linear Programming Approach for a Distance-Constrained Elementary Path Problem

Keywords: Graphs, Optimization, Longest Path Problem, Mixed Integer Programming

Given a directed graph G = (V, E, l) with weights l e ≥ 0 associated with arcs e ∈ E and a set of vertex pairs with distances between them (called distance constraints), the problem is to find an elementary path in G that satisfies a maximum number of distance constraints. We describe two MIP formulations for this problem and discuss their advantages.

Introduction

We are interested in the following problem. Consider a directed graph G = (V, E, l) where l e ≥ 0 is the weight associated with each arc e ∈ E. It is also given a set DC of vertex pairs (called distance constraints) such that with any couple (u, v) ∈ DC a distance interval [d(u, v), d(u, v)] is associated (i.e. (d(u, v) ≤ d(u, v) are lower (upper) bounds for the distance between u and v along the solution path). Let P be a path in G. We say that P satisfies a given distance constraint dc(u, v) = [d(u, v), d(u, v)] if both u and v are on P and the subpath of P between u and v has length in dc (u, v). Recall that the length of a path p is defined as the sum of the weights of the arcs in p. Then the problem we are considering is to find an elementary1 path in G that satisfies a maximum number of distance constraints. We call it Distance-Constrained Elementary Path (DCEP) problem. The DCEP is N P-hard due to a simple reduction from the Hamiltonian path problem. To see this, consider a set DC containing all ordered pairs of vertices and such that (

d(u, v), d(u, v)) = (-∞, ∞) for any (u, v) ∈ DC.
As far as we know, the DCEP problem has not been previously discussed in the combinatorial optimization community. The problem is motivated by applications in genome assembly in bioinformatics, and a variation of the problem has been originally described2 in two of our previous publications [START_REF] François | Global optimization methods for genome scaffolding[END_REF][START_REF] François | Global optimization approach for circular and chloroplast genome assemblies[END_REF]. The Genome assembly problem is a challenging computational task aiming at reconstructing the full genome of an organism from short DNA sequences (reads) [START_REF] Daniel | The greedy path-merging algorithm for contig scaffolding[END_REF]. Note that in spite of the efforts and the progress done by the bioinformatics community, no satisfactory solution is available today. The method proposed in [START_REF] François | Global optimization methods for genome scaffolding[END_REF][START_REF] François | Global optimization approach for circular and chloroplast genome assemblies[END_REF] is based on integer programming model for solving genome assembly as a problem of finding a long simple path in a specific graph, which satisfies additional constraints encoding the insert-size (distance) information. This methodology significantly differs from the heuristics described in the literature like BESST [START_REF] Sahlin | BESST -efficient scaffolding of large fragmented assemblies[END_REF] and SPAdes [START_REF] Bankevich | Spades: a new genome assembly algorithm and its applications to single-cell sequencing[END_REF]. The computations performed on several chloroplast genomes demonstrate that such an approach outperforms these widely-used assembly solvers by the accuracy of the results.

While the paper [START_REF] François | Global optimization approach for circular and chloroplast genome assemblies[END_REF] is application oriented, here we revisit the mixed integer linear programming formulation proposed there from a combinatorial optimization viewpoint. Furthermore, we show how to adapt the well known Miller, Tucker and Zemlin (MTZ) formulation for solving the longest path problem [START_REF] Miller | Integer programming formulation of traveling salesman problems[END_REF] in order to solve the DCEP problem. Note that the challenges in DCEP problem are somehow similar, but harder in practice, to the ones in solving longest/shortest (with real weights) elementary path problems [START_REF] Bui | Exact methods for solving the elementary shortest and longest path problems[END_REF][START_REF] Taccari | Integer programming formulations for the elementary shortest path problem[END_REF]).

Standard integer programming formulation related to elementary path description in a directed graph

One of the specificities of problem is that the beginning and the end of the path are unknown. The most natural approach to overcome this problem is to introduce two dummy vertices s and t (source and target) and to connect them with the rest of the vertices. The extended graph has 2|V | more edges. Hence we consider a directed graph Ĝ = { V , Ê, l} with set of nodes V and set of arcs Ê where

• V = V ∪ {s, t} and s and t are such that s / ∈ V and t / ∈ V .

• Ê = E ∪ {(s, v)|v ∈ V } ∪ {(u, t)|u ∈ V } and l(s, v) = l(u, t) = 0, (s, v) ∈ Ê, (u, t) ∈ Ê.
In all formulations below we assume that |δ -(s

)| = |δ + (t)| = ∅, where by δ -(v)/δ + (v)
we denote the set of ingoing/outgoing edges for a vertex v. The standard integer programming formulation for constructing an elementary path from s to t is to find values of binary variables x e for all e ∈ Ê (whose meaning is that x e = 1 if e belongs to the solution path, and x e = 0, otherwise) such that ∀u ∈ V,

e∈δ + (u) x e ≤ 1 (1)
∀u ∈ V :

e∈δ + (u)
x e -e∈δ -(u)

x e =    1 if u = s -1 if u = t 0 else (2)
Constraint (1) ensure that the outgoing degree of each node is at most one, and constraints (2) are flow conservation constraints.

In the sequel we describe two different sets of constraints and variables that can be added to formulation (1)-( 2) for sub-tour elimination in case of cycles.

Sequential formulation (MTZ) and its adaptation to the DCEP

To derive the well known MTZ formulation ( [START_REF] Miller | Integer programming formulation of traveling salesman problems[END_REF]), it is enough to introduce, for each vertex v, an auxiliary variable y v ≥ 0 that allows us to number/label the vertices along the path in an increasing order. The constraint (3) allows then to prevent sub-tours.

∀(u, v) ∈ E, (y v -y u ) ≥ x u,v -(1 -x u,v )|V |. (3)
Set W = e∈E l e and y s = W . To adapt (3) to the DCEP we replace (3) by the constraint

∀(u, v) ∈ Ê, y v ≤ y u -x u,v l u,v + (1 -x u,v )W, (4)
where W plays the role of a big constant. Since l e ≥ 0, the labels y v are decreasing along the path. This avoids cycles. Another advantage of ( 4) is that, for any couple (u, v) from the path, the gap y u -y v measures the distance between vertices u and v.

In order to manage the distance constraints, we apply a technique similar to the one in [START_REF] François | Global optimization methods for genome scaffolding[END_REF][START_REF] François | Global optimization approach for circular and chloroplast genome assemblies[END_REF]. We introduce a new variable g e ∈ {0, 1}, e ∈ DC, and we set to 1 the value of g (u,v) if and only if both vertices u and v belong to the selected path and its length between them is in the given interval [d (u,v) , d (u,v) ]. We also add the constraint

∀v ∈ V, y v ≤ W ( e∈δ -(v)
x e ), [START_REF] François | Global optimization methods for genome scaffolding[END_REF] which sets to zero the values of y v outside the selected path. The other constraints related to the distances DC are ∀e ∈ DC : g (u,v) ≤ y u and g (u,v) ≤ y v [START_REF] Daniel | The greedy path-merging algorithm for contig scaffolding[END_REF] as well as

∀(u, v) ∈ DC : d(u, v)g (u,v) ) + W (1 -g (u,v) ) ≥ y u -y v ≥ d(u, v)g (u,v) -W (1 -g (u,v)
). ( 7)

We search for a path that satisfies a maximum number of distance constraints max W (

e∈DC g e ) + v∈V y v . ( 8
)
The last term in ( 8) forces the labels of the vertices on the path to take their minimal values. We denote by MTZDC the obtained in this manner formulation.

GAT formulation ([4])

To any vertex v ∈ V we associate a variable i v , 0 ≤ i v ≤ 1, encoding whether v is on the solution path. The two possibles states for a vertex v (to be an intermediate vertex in the path or not) are enforced by the following constraints

i v = e∈δ + (v) x e = e∈δ -(v)
x e . (

One can then show [START_REF] François | Global optimization methods for genome scaffolding[END_REF] that the real variables i v , ∀v ∈ V take binary values. We introduce a continuous variable f e ∈ R + to express the quantity of the flow circulating along the edge e ∈ E. Without this variable, the solution found may contain some loops and hence may not be a simple path. We put a requirement that no flow can use an edge e when x e = 0, which can be encoded as

∀e ∈ E : 0 ≤ f e ≤ W x e , ( 10
)
where W is as defined above (W = e∈E l e ). We use the flows f e in the following constraint ∀v ∈ V :

e∈δ -(v) f e - e∈δ + (v) f e = e∈δ + (v) l e x e , (11)
while for the source vertex we require

e∈δ + (s)
f e = W . The flow then decreases along the path and this feature forbids cycles. Furthermore, a variable g e ∈ {0, 1} is associated with any distance constraint. The value of g (u,v) is set to 1 only if both vertices u and v belong to the selected path and its length between them is in the given interval [d (u,v) , d (u,v) ]. The constraints are as follows : ∀e ∈ DC : g e ∈ {0, 1} and g (u,v) ≤ i u and g (u,v) ≤ i v (12) as well as :

∀(u, v) ∈ DC : d(u, v)g (u,v) )+M (1-g (u,v) ) ≥ e∈δ -(u) f e - e∈δ -(v) f e ≥ d(u, v)g (u,v) -M (1-g (u,v) ) (13
) We search for a path that satisfies a maximum number of distance constraints. max e∈DC g e .

(14)

Conclusion and perspectives

In Table 1 we summarize the presented formulations. We observe that, on theory, they are very similar-both formulations have almost the same number of variables and the same number of constraints. However, they are dual-like in the sense that in MTZDC the distances are computed in the vertex variables y v , while in GAT they are in the flow variables f e which are associated in the arcs. The sub-tour elimination constraint in MTZDC (4) requires |E| inequalities, while the same constraints in GAT (11) requires |V | equations. We have implemented both formulations with AMPL language. On small instances they behave similarly. We are currently getting statistics on huge instances towards more precise performance analysis.

MTZDC

An elementary (also called simple) is a path that visits each vertex at most once.

However, the name Distance-Constrained Elementary Path is firstly used here.