J ' Aimerais Remercier

Éric Cancès Ainsi Que

Gabriel Stoltz

Antoine Levitt

David Gontier

Irek Grabowski

Paola Gori-Giorgi

Trond Saue

Stéphane Carniato

Elisa Stéphanie

Inès Etienne Samira

Aixiao, Olivier, Mohamed, Ozge, Eleonore, Emanuele, Sehr Daniel

Chantal

Chiara, Silvia, Jérémy, Loïc Roman Lydie

Odile Franck
email: odile.franck@etu.upmc.fr

Bastien Mussard
email: bastien.mussard@upmc.fr

Eleonora Luppi
email: eleonora.luppi@upmc.fr

Julien Toulouse
email: julien.toulouse@upmc.fr

Szymon Śmiga
email: szsmiga@fizyka.umk.pl

Adam Buksztel

Ireneusz Grabowski

Basis convergence of range-separated density-functional theory

J'aimerais remercier toutes les personnes qui sont venues assister à la soutenance.

Introduction

Early attempts of building a density-functional theory can be found in the work of Thomas [1] and Fermi [2] in 1927. In their model the energy is expressed with respect to the density with a simple expression for the kinetic energy based on the uniform electron gas and nuclearelectron and electron-electron interactions are described classically. This model was extended to the Thomas-Fermi-Dirac model by including an exchange energy formula for the uniform electron gas introduced by Dirac [3] in 1930. Another early model was given by Slater [4] in 1951 who proposed an approximation to the non-local exchange in the Hartree-Fock method that depends only on the local electron density. Density-functional theory as known today was first introduced in 1964 by Hohenberg and Kohn [5] as an alternative to solving the Schrödinger equation. Density-functional theory applied in the Kohn-Sham [6] scheme is an exact method if the exact exchange-correlation density functional is known, unfortunately it is not and a major topic of research is to define better and better approximated functionals. The first approximated functional that was proposed and that was used as a starting point to further development is based on an model system, an hypothetical uniform electron gas. In this approximation at each point of an inhomogeneous system the local exchange-correlation energy per particle is taken as the exchange-correlation energy per particle of the uniform electron gas of the same density. This approximation is called the local-density approximation (LDA), the most known parametrization was given by Vosko et al. [7]. Even if LDA is based on a simple approximation, realistic systems being different from an uniform electron gas, it performs surprisingly good and can even be comparable to or better in accuracy than Hartree-Fock, showing good accuracy for molecular properties such as equilibrium structure but failing to describe energetical quantities such as bonding energies with an overbinding tendency. An extension to improve the performance of LDA is to take information from the gradient of the density to take into account the nonhomogeneity of the true electronic density which leads to a new family of approximations: generalized-gradient approximations (GGAs). An example of such functional is given by B88 [8] for the exchange and LYP [9] for the correlation. A way to improve the performance of GGAs is to take into account the Laplacian of the density and/or the kinetic energy density, leading to a new family of approximations: meta-generalized-gradient approximations (meta-GGAs). One of the most used meta-GGA approximation is TPSS defined by Tao et al. [10].

LDA, GGAs and meta-GGAs are often referred to as semilocal approximations because they INTRODUCTION only depend on the density at a point or on the derivatives of the density at this point (or the derivatives of the orbitals for meta-GGAs). These semilocal approximations often give an accurate description of short-range dynamical correlation but fails to describe long-range or static correlation. These semilocal approximations present typically a self-interaction error which tends to favor delocalization of the electrons and induces too low total energies.

A way to improve the performance of approximated functionals, in particular by reducing the self-interaction error, is to combine density-functional theory with wave-function theory and create hybrid approximations. Combining both theories can be done in different ways. One of the simplest way of doing it is by doing a linear separation of the electron-electron Coulomb interaction into two parts

1 r = λ r W F T + (1 -λ) r DF T
, the first term being treated with using wave-function theory (WFT) while the second term being treated with density-functional theory (DFT) and λ is the parameter of this hybridation.

A first realization of this was done in 1993 by Becke with the half-and-half combination [11] (i.e. λ = 0.5) of Hartree-Fock exchange and a density-functional approximation. Most of the time this fraction of Hartree-Fock exchange was too important and some error compensation between exchange and correlation was lost. Later this hybrid scheme was extended using several empirical coefficients with a smaller coefficient of Hartree-Fock exchange [12]. Common hybrid approximations nowadays use a fraction λ 0.2 -0.25 of Hartree-Fock exchange [13].

An extension of such hybrid approximations is achieved by introducing a fraction of correlation energy calculated using second-order Møller-Plesset perturbation theory (MP2) and is known as double-hybrid approximations. It was originally introduced by Grimme [14]. Double-hybrid approximations allow us to use a more important fraction of Hartree-Fock exchange (λ 0.5 -0.7) than for hybrid approximations without loosing too much the benefit of the error compensation but the method fails to describe phenomena that cannot be treated with MP2, for example static correlation. The fraction of correlation energy calculated with wave-function theory can be treated using other approximations such as random-phase approximations [15]. To improve the description of static correlation, density-functional theory can also be combined with the multiconfiguration self-consistent-field (MCSCF) method [16].

Combining density-functional theory and wave-function theory can be done going beyond the linear combination with the range-separated approach introduced by Savin [17] in 1996 by decomposing the electron-electron Coulomb interaction into a long-range part and a short-range part using the error function , where µ is a parameter controlling the range of the separation. The long-range interaction is described using wave-function theory and the short-range interaction is described using densityfunctional theory. A version limited to the range separation of the exchange, named long-range correction scheme (LC), was proposed by Iikura et al. [18] by introducing long-range Hartree-Fock exchange while short-range exchange and correlation are treated using density-functional approximations. The range separation can also be done on the correlation, for example using MP2 for the long-range correlation with long-range Hartree-Fock exchange and a short-range exchange-correlation functional [19]. Such a decomposition can be performed using other methods to describe the long-range correlation such as random-phase approximations [20] or coupledcluster methods [21], which are well adapted to describe van der Waals dispersion interactions.

A multiconfigurational treatment of the long-range correlation can be used to improve the description of static correlation such as long-range MCSCF [22] or long-range density-matrixfunctional theory (DMFT) [23].

The density-functional theory formalism has been extended to describe excited states, with linear-response time-dependent density-functional theory. The key quantity in this approximation is the exchange-correlation kernel. The spatial and frequency dependence of this kernel needs to be approximated. The simplest approximation is the adiabatic semilocal approximation in which the kernel is local in time (i.e., independent of the frequency) and local in space. This approximation gives reasonably good results for low-lying valence electronic excitation energies of molecular systems but it fails to describe some phenomena such as multiple excitations, charge-transfer excitation energies and Rydberg excitation energies. To overcome these limitations, time-dependent density-functional theory has been extended to range separation. The decomposition of the exchange kernel into a long-range Hartree-Fock exchange kernel and a short-range exchange kernel described by a density-functional approximation has first been performed by Tawada et al. [24] and was able to correct some problems from the semilocal approximation such as the description of Rydberg excitation energies and charge-transfer excitation energies. The range-separated scheme can be extended by using a short-range correlation kernel calculated with a density-functional approximation and using long-range linear-response MCSCF [25] or long-range linear-response DMFT [26] approaches. The use of multiconfigurational methods to describe the long-range response improves the description of static correlation and it also allows one to calculate double excitations. The short-range correlation kernel can also be combined with a long-range correlation kernel calculated using the many-body Green function formalism [27] which gives us a frequency-dependent long-range correlation kernel.

INTRODUCTION

In this thesis we will investigate several aspects of hybrid methods combining wave-function theory and density-functional theory. The first chapter will give a brief overview of densityfunctional theory and these hybrid approaches.

In the second chapter the study is centered on the basis-set convergence of range-separated hybrid methods. The basis-set convergence has been studied a lot for wave-function theory and range-separated hybrid methods have been shown to converge faster, but the convergence rate had not been explored yet. In this chapter we first studied the convergence in a partial-wave expansion of the long-range wave function with respect to the maximal angular momentum. We then studied the convergence of the long-range second-order Møller-Plesset correlation energy with respect to the cardinal number of the Duning basis sets (cc-p(C)VXZ). The obtained results allowed us to propose a three-point extrapolation scheme for the complete basis set energy of range-separated hybrid density-functional theory.

The third chapter will be focused on double-hybrid density-functional methods combining density-functional theory with second-order Møller-Plesset perturbation theory (MP2). These methods give accurate results for thermochemical properties. Commonly the orbitals are evaluated without the MP2 term, which is added a posteriori. Recently Peverati and Head-Gordon [28] proposed an orbital-optimized double-hybrid method where the orbitals are self-consistently optimized in the presence of the MP2 correlation term. This orbital-optimized double-hybrid method has shown an improvement in the spin-unrestricted calculations for symmetry breaking and open-shell situations. In this study we will consider an alternative orbital-optimized doublehybrid method based on the optimized-effective-potential (OEP) method that could bring advantages for calculations of excitation energy and response properties and a better description of the LUMO orbital energy. We will compare the results for such an OEP double-hybrid method to the standard double-hybrid method for the calculation of atomic and molecular properties such as ionization potentials and electronic affinities.

In the fourth chapter we will consider range-separated linear-response time-dependent densityfunctional theory and we will study the short-range exchange and correlation kernels. We started by generalizing the exact-exchange kernel [29,30] to range-separated time-dependent density-functional theory. We then studied the behavior of the kernel with respect to the rangeseparation parameter (µ) and we compared the behavior of the short-range exchange kernel with the adiabatic LDA for He and H 2 . Finally we studied the frequency-dependent short-range correlation kernel for a model system: H 2 in a minimal basis set.

Finally in the last chapter we will give some general concluding remarks and outline.

Chapter 1

Review of density-functional theory and hybrid methods

In this chapter we will first briefly recall the many-body problem in Sec. 1.1 before introducing density-functional theory in Sec. 1.2 and finally an introduction of the range-separated hybrid approximations in Sec. 1. 3. Further details can be found in Refs. [1,2].

Schrödinger equation

The time-independent non-relativistic Schrödinger equation allowing us to describe atomic, molecular or solid-state systems is given by

Ĥ|Ψ = E|Ψ , (1.1)
with the energy E, the wave function Ψ and the Hamiltonian operator Ĥ. For a system with M nuclei and N electrons the Hamiltonian in position representation is

Ĥ = - N i=1 1 2 ∇ 2 i - M A=1 1 2M A ∇ 2 A - N i=1 M A=1 Z A r iA + N i=1 N j>i 1 r ij + M A=1 M B>A Z A Z B R AB . (1.
2

)
The equation is given in atomic units as will be all the equations in the following. In this equation the sum for A runs over all the nuclei up to M and the sum for i runs over all electrons up to N . The first two terms of the right-hand of the equation are the kinetic energy for the electrons and the nuclei, respectively. The Laplacian ∇ 2 q is the sum of the second-order partial derivatives; in Cartesian coordinates it is

∇ 2 q = ∂ 2 ∂x 2 q + ∂ 2 ∂y 2 q + ∂ 2 ∂z 2 q .
The three last terms are the nuclei-electron interaction, the electron-electron interaction and the nuclei-nuclei interaction. We can simplify this Hamiltonian by considering the Born-Oppenheimer

Density-functional theory

The aim of density-functional theory is to solve the many-body problem by expressing the energy as a functional of the one-electron density n(r) = N . . . |Ψ(x, x 2 , . . . , x N )| 2 ds dx 2 . . . dx N (1.4) with the integration over s is the sum over both s = +1/2 and s = -1/2. The density is normalized as n(r)dr = N .

Hohenberg-Kohn theorems

First Hohenberg-Kohn theorem

The external potential v(r) is (to within a constant) a unique functional of n(r); since, in turn v(r) fixes Ĥ we see that everything including the full many-particle ground-state energy is a unique functional of n(r).

The proof for this theorem [3] is given if we consider two external potentials v(r) and v (r) that differ from more than one constant and that both lead to the same density n(r) for an N -electron system. Each potential leads to different Hamiltonians Ĥ and Ĥ , respectively, and to different ground-state wave functions Ψ and Ψ , respectively. If we consider Ψ as a trial function for Hamiltonian Ĥ we have

E 0 < Ψ | Ĥ|Ψ = Ψ | Ĥ |Ψ + Ψ | Ĥ -Ĥ |Ψ = E 0 + n(r) v(r) -v (r) dr, (1.5)
with E 0 the ground-state energy of Ĥ and E 0 the ground-state energy of Ĥ . If we now consider Ψ as a trial function for Hamiltonian Ĥ E 0 < Ψ| Ĥ |Ψ = Ψ| Ĥ|Ψ + Ψ| Ĥ -Ĥ|Ψ = E 0n(r) v(r)v (r) dr. (1.6) Combining Eq. (1.5) and Eq. (1.6) we obtain E 0 + E 0 < E 0 + E 0 which is nonsense and shows that there can only be one potential up to a constant leading to the same ground-state density.

We now introduce the Hohenberg-Kohn functional F such that the energy becomes, for the specific external potential v ne (r),

E[n] = F [n]
+ n(r)v ne (r)dr.

(1.7)

The functional F includes the terms of energy that are universal (i.e., independent from the external potential),

F [n] = T [n] + E ee [n] = Ψ [n] | T + Ŵee |Ψ [n] , (1.8)
where

Second Hohenberg-Kohn theorem

We have shown that the energy is a functional of the density and that the ground-state energy is then obtained by using the ground-state density. The second Hohenberg-Kohn theorem allows us to use the variational principle for the Hohenberg-Kohn functional.

F [n], the functional that delivers the ground-state energy of the system, delivers the lowest energy if and only if the input density is the true ground-state density, n 0 . This is analogous to the variational principle applied to wave functions

E 0 ≤ E [ñ] = F [ñ] + ñ(r)v ne (r)dr.
(1.9)

The proof for this theorem is based on the variational principle. We consider a trial density ñ that implies the Hamiltonian H which defines the ground-state wave function Ψ

Ψ| Ĥ| Ψ = T [ñ] + V ee [ñ] + ñ(r)v ne (r)dr = E [ñ] ≥ E 0 [n 0 ] = Ψ 0 | Ĥ|Ψ 0 .
(1.10)

Finally the ground-state energy can simply be expressed as Considering the definition of the universal functional in Eq. (1.12) we obtain

E 0 =
E 0 = min n F [n] + n(r)v ne (r)dr . (1.16)

Kohn-Sham approach

We saw previously that the energy can be expressed as a functional of the density but we still have no expression for the Hohenberg-Kohn functional F . This functional should include the kinetic energy, the classical Coulomb interaction (Hartree) and the non-classical contributions (exchange and correlation).

Decomposition of the universal functional

We first consider the kinetic energy. We can exactly define the kinetic energy of a non-interacting system at a given density

T s [n] = min Φ→n Φ| T |Φ = Φ [n] | T |Φ [n] (1.17)
where Φ is a single determinant and Φ[n] is the wave function minimizing T and yielding n.

This kinetic energy T s is different from the exact kinetic energy T and we will need to take care of this difference in the functional. We decompose F [n] as

F [n] = T s [n] + E Hxc [n] (1.18)
with the Hartree-exchange-correlation energy E Hxc which can be decomposed in two contributions: the Hartree energy

E H [n] = 1 2 n(r)n(r ) |r -r | drdr , (1.19)
and the exchange-correlation energy E xc . The non-interacting kinetic energy and the Hartree energy are known exactly while the remaining terms are included in the exchange-correlation energy

E xc [n] = (T [n] -T s [n]) + (E ee [n] -E H [n]). (1.20)
The exchange-correlation energy thus contains the remaining of the kinetic energy and the nonclassical contribution to the electron-electron interaction energy E ee . The ground-state energy [5] for a given potential v ne is then

E = min n min Φ→n Φ| T |Φ + E Hxc [n Φ ] + v ne (r)n Φ (r)dr = min Φ Φ| T + Vne |Φ + E Hxc [n Φ ] . (1.21)
The minimizing single-determinant KS wave function giving the exact ground-state density is written as

Φ s = 1 √ N ! |χ 1 (x 1 )χ 2 (x 2 ) . . . χ N (x N )| . (1.22)
The function χ i (x) with x = (r, s) is a spin orbital composed of a product of spatial orbital ϕ i (r) and one of the two spin functions α(s) or β(s).The spatial orbitals fulfills

- 1 2 ∇ 2 + v s (r) ϕ i (r) = i ϕ i (r). (1.23)
This potential v s (r) is such that the density of the reference system is the density of the real system and is

v s (r) = v ne (r) + v H (r) + v xc (r), (1.24)
with the Hartree potential corresponding to the derivative of the Hartree energy with respect to the density v H (r) = n(r )/|rr |dr and the exchange-correlation potential

v xc (r) = δE xc [n] δn(r) . (1.25)
If the exchange-correlation energy was known exactly the method would be exact. Actually it is unknown and approximations are used for practical calculations. The exchange-correlation energy functional can be splitted in an exchange and a correlation energy functional. The exchange functional is known

E x [n] = Φ[n]| Ŵee |Φ[n] -E H [n] (1.26)
while the correlation energy functional is

E c [n] = F [n] -(T s [n] + E H [n] + E x [n]) = Ψ[n]| T + Ŵee |Ψ[n] -Φ[n]| T + Ŵee |Φ[n] (1.27)
with Ψ[n] the wave functions minimizing T + Vee and yielding n.

Some approximated functionals

LDA

The local-density approximation (LDA) is based on the simple model of the uniform electron gas. In this model the electrons are moving on a background of positive charge distribution such that the complete system is electrically neutral and defined by its number of electrons N and its volume V that are infinite and its density n = N/V that is finite. In this approximation, the exchange-correlation functional is expressed as

E LDA xc [n] = n(r) xc (n(r)) dr, (1.28)
whith xc (n) the exchange-correlation energy per particle of a uniform electron gas of density n that can be split into exchange ( x ) and correlation ( c ) contributions

xc (n) = x (n) + c (n) . (1.29)
The exchange energy per particle of a uniform electron gas was given by Dirac [6] and Slater [7] x (n) = -

3 4 3n π 1/3 , (1.30)
while the correlation energy per particle of a uniform gas is obtained by analysis and interpolation of highly accurate quantum Monte-Carlo simulations of the homogeneous electron gas [8]. The most known parametrization was proposed by Vosko et al. [9].

Another local approximation includes the spin resolution, where the energy becomes a functional of the spin densities n α and n β with n α + n β = n. Even if the functional does not have to depend on the separate spin densities (if the potential is spin-independent) it may improve

the approximation for open-shell systems. The local-spin-density approximation (LSD) is then simply defined as

E LSD xc [n α , n β ] = n(r) xc (n α (r), n β (r)) dr, (1.31)
where xc (n α , n β ) is the spin-resolved exchange-correlation energy per particle of the uniform electron gas.

The LDA shows a good performance, comparable or even better to the Hartree-Fock approximation for properties such as equilibrium structures or harmonic frequencies but problems remain to describe properties such as bond energies or atomization energies. For example the LDA overestimates systematically the atomization energies. Finally the LDA was mostly used in solid-state physics and less for computational chemistry.

GGAs and meta-GGAs

An extension of LDA can be found in the generalized-gradient approximation. In this approximation the functional at a point r depends not only on the density n(r), but also on the gradient of the density ∇n(r) in order to introduce inhomogeneity in the electron density of the model.

The first attempt was made with the gradient-expansion approximation where the exchangecorrelation energy functional is defined as a Taylor expansion with respect to the density and its derivatives. The first term of this expansion corresponds to the LDA approximation. This expansion does not actually improve the performance of LDA because the exchange-correlation hole defined by the expansion does not reproduce the constraints of the physical exchange-correlation hole.

The generalized-gradient approximation is then defined as follows

E GGA xc [n] = f (n(r), ∇n(r)) dr, (1.32)
where the integrand f is a function depending both on the density and the gradient of the density. The most important and most used GGA functionals are BLYP [10,11] and PBE [12].

The Laplacian of the density ∇ 2 n(r) and/or the kinetic energy density τ (r) can also be used in the definition of functionals to improve the performance of GGAs. This defines a new family of approximations: the meta-generalized-gradient approximations (meta-GGAs ou mGGAs). This family is defined as

E mGGA xc
[n] = f n(r), ∇n(r), ∇ 2 n(r), τ (r) dr, (1.33) and the most important meta-GGA functional is TPSS [13].

Hybrid approximations

To improve the performance of GGA and meta-GGA functionals an hybrid approximation was proposed, first by Becke [14] where the exchange was decomposed in two contributions, a fraction calculated using a density-functional approximation (DFA) and a fraction of Hartree-Fock (HF) exchange

E hybrid xc = a x E HF x + (1 -a x )E DFA x + E DFA c (1.34)
with the Hartree-Fock exchange energy

E HF x = - 1 2 occ. i,j χ * i (x 1 )χ j (x 1 )χ * j (x 2 )χ i (x 2 ) |r 1 -r 2 | dx 1 dx 2 . (1.35)
The ratio can be modified and it was further shown that the optimal fraction of Hartree-Fock exchange should be around a x 0.2 -0.3. This type of functionals brings improvement to the description of some properties (such as thermodynamic properties) at a reasonable computational cost but there is sometimes a loss with respect to GGA and meta-GGA functionals due to the error compensation. Different options of development can then be considered to improve the description of the correlation effects and particularly the non-local correlation effects. One way to do this is to hybridize the correlation.

Double-hybrid approximations

Another type of hybrid approximation includes a fraction of correlation calculated with secondorder Møller-Plesser perturbation theory (MP2), namely the double hybrid (DH) approximation

E DH xc = a x E HF x + (1 -a x )E DFA x [n] + a c E MP2 c + (1 -a c )E DFA c [n] (1.36)
where a x is the fraction of Hartree-Fock exchange and a c is the fraction of MP2 correlation given by

E MP2 c = occ. i<j unocc. a<b | χ i χ j | ŵee |χ a χ b -χ i χ j | ŵee |χ b χ a | 2 i + j -a -b , (1.37)
where χ i χ j | ŵee |χ a χ b are the two-electron integrals with ŵee the electron-electron interaction.

A rigorous formulation of these double-hybrid approximation was given by Sharkas et al. [15], in which the approximation has one parameter (a c = a 2 x ) and a density scaling. A rigorous formulation of the two-parameter double-hybrid approximation was given by Fromager [16].

Range-separated hybrid approximations

The range-separated hybrid approximations are obtained by decomposing the electronic interaction in a short-range and a long-range contribution [17] This separation scheme is motivated by the idea of using the best of density-functional theory and wave-function theory: first the good description at short-range given by density functional theory and second the good performance to describe static correlation effects with wave-function theory, while keeping the computational cost reasonably low. The decomposition of the electronic interaction is 1 r 12 = w lr,µ ee (r 12 ) + w sr,µ ee (r 12 ) (1.38) where w lr,µ ee is the long-range interaction and w sr,µ ee is the complement short-range interaction. The transition between the two interactions is made by the use of the error function, the longrange interaction being

w lr,µ ee (r 12 ) = erf(µr 12 ) r 12 , (1.39)
where µ is the range-separation parameter.

Considering these decomposition the universal functional F [n] becomes

F [n] = F lr,µ [n] + E sr,µ Hxc [n] (1.40)
with the long-range universal functional F lr,µ [n] and the complement short-range Hartreeexchange-correlation functional E sr,µ Hxc [n]. The long-range universal functional is given by

F lr,µ [n] = min Ψ→n Ψ| T + Ŵ lr,µ ee |Ψ (1.41)
where T is the kinetic operator, Ŵ lr,µ ee is the long-range interaction and Ψ is a multideterminant wave function. The ground-state energy for a given potential v ne is then

E 0 = min n F lr,µ [n] + E sr,µ Hxc [n] + n(r)v ne (r)dr = min Ψ Ψ| T + Ŵ lr ee + Vne |Ψ + E sr,µ Hxc [n Ψ ] = Ψ µ | T + Ŵ lr,µ ee |Ψ µ + E sr,µ Hxc [n Ψ µ ] + n Ψ µ (r)v ne (r)dr. (1.42)
The minimizing ground-state multi-determinantal wave function Ψ µ fulfills

T + Ŵ lr,µ ee + V sr,µ |Ψ µ = E µ |Ψ µ (1.43)
where the short-range potential V sr,µ = i v sr,µ (r i ) with

v sr,µ (r) = v ne (r) + δE sr,µ Hxc [n Ψ µ ] δn(r) . (1.44)
The v sr,µ is unique up to a constant as shown by the first Hohenberg-Kohn theorem. To express the second term we need to decompose the short-range Hartree-exchange-correlation functional

E sr,µ Hxc [n] = E sr,µ H [n] + E sr,µ xc [n] (1.45)
the first term is the complement short-range Hartree energy functional

E sr,µ H [n] = E H [n] - 1 2 n(r)n(r )w lr,µ ee (|r -r |)drdr , (1.46)
and E sr,µ xc [n] the unknown short-range exchange-correlation energy.

If E sr,µ xc [n] is known the method is exact but in practice approximations need to be introduced. The first level of approximation is the range-separated hybrid (RSH) [18]

using a N -electron normalized single-determinant wave-function instead of Ψ µ E µ RSH = min Φ Φ| T + Vne + Ŵ lr,µ ee |Φ + E sr,µ Hxc [n Φ ] (1.47)
where the minimizing Φ µ fulfills

T + Vne + V lr,µ Hx,HF [Φ µ ] + V sr,µ Hxc [n Φ µ ] |Φ µ = E µ 0 |Φ µ . (1.48)
The long range correlation energy is then added a posteriori

E = E µ RSH + E lr,µ c,MP2 (1.49)
where the long-range MP2 correlation energy is given by

E lr,µ c,MP2 = occ. i<j unocc. a<b χ µ i χ µ j | ŵlr,µ ee |χ µ a χ µ b -χ µ i χ µ j | ŵlr,µ ee |χ µ b χ µ a 2 µ i + µ j -µ a -µ b , (1.50)
with the set of RSH spin orbitals {χ µ k } and the RSH orbital energies µ k . Different methods can be used to evaluate this correlation energy: coupled-cluster theory [19,20], RPA approximations [21,22] Chapter 2

Basis convergence of range-separated density-functional theory

In this chapter we focused on the basis convergence of range-separated density-functional theory (DFT). This work has been published in The Journal of Chemical Physics [J. Chem.

Phys. 142, 074107 (2015)] and included in this chapter. The results presented in this article are the summary of two studies on the basis convergence in two different contexts: the partial wave expansion (where the basis is constructed by adding at each step all the orbitals corresponding to a given angular momentum ) and the principal expansion with the Dunning basis sets (where the cardinal number of the basis X that can be linked to the principal quantum number).

The starting point of the work on the partial-wave expansion was based on the study of the second-order energy (E 2 ) for a two-electron atom proposed by Schwartz [1]. In this study the author expressed E 2 with respect to , starting from perturbation theory and expanding the first-order wave-function Ψ 1 and E 2 in a basis of Legendre polynomials

E 2 ( ) = - 45 256λ 2 + 105 256λ 3 - 25965 65536λ 4 + O 1 λ 5 ,
where λ = ( + 1/2) 2 . We then considered a similar work of Kutzelnigg and Morgan [2] which proposed a similar study based on the form of the wave-function proposed by Kutzelnigg [3] to reproduce the correlation cusp condition defined by Kato [4] by imposing linearity with respect to the inter-electronic distance. The first step of our study was to reproduce those proofs. The next step was to extend this work to range separation. The long-range second-order energy is then given by

E lr,µ 2 = Ψ 0 | Ŵ lr,µ ee -E 1 |Ψ lr,µ 1
where the long-range second-order energy E lr,µ 2 , the long-range first-order wave-function Ψ lr,µ 1 and the long-range interaction Ŵ lr,µ ee need to be expanded in the basis of Legendre polynomials.

BIBLIOGRAPHY

Such work shows a lot of complexity and the method may not be the best to describe the convergence (the convergence of the first terms of a series may not be sufficient to make a statement on the convergence of the series). A way to overcome this limitation was to rather focus on the convergence of the wave function in the region of electron coalescence, as we know that in the Coulomb case there is a singularity that appears at coalescence which is the limiting factor of the convergence. We choose to consider a spherical model where two electrons are on a 1a 0 sphere and we compared the convergence of the Coulomb and long-range wave functions using the form proposed by Kutzelnigg for the Coulomb interaction and the wave function proposed by Gori-Giorgi and Savin [5] for the long-range interaction. We observed a change of convergence rate with range separation that converges exponentially.

The second part of this work was focused on the convergence with respect to the cardinal number X of the Dunning basis sets (cc-pVXZ). To connect with the part on partial-wave expansion we began with a study of the convergence of the wave function of the helium atom. This work was done in collaboration with Bastien Mussard. We then wanted to extend previous works on the convergence of correlated calculations [6] to range separation. The basis convergence of the second-order energy with respect to X becomes exponential with range separation. Finally a three-point extrapolation scheme was proposed for the complete basis set limit. Supplementary results that were not included in the article are presented in Appendix A.

Discussions with other researchers pointed out that the computational cost of the three calculations needed for the extrapolation was too high and an extension of this work could be to find a way to simplify this extrapolation scheme for only two points so that it could be used in practical calculations. Another point to discuss is to know whether similar results could be expected if the calculations were performed in a self-consistent way (a preliminary study on the path to a self-consistent RSH approach is presented in Chapter 3). We expect the results on the convergence to be the same in this situation. coupled-cluster theory [17][18][19][20][21], or random-phase approximations [22][23][24][25][26][27][28][29][30][31][32][33][34].

Bibliography

Among the advantages of using such many-body methods for the long-range part only of the electron-electron interaction is that they are much less sensitive to the one-electron atomic basis compared to the case of the standard Coulomb interaction. This has been repeatedly observed in calculations using Dunning correlation-consistent basis sets [35] for second-order perturbation theory [4,6,10,15,16], coupled-cluster theory [17] and randomphase approximations [22,23,25,27,31]. The physical reason for this reduced sensitivity to the basis is easy to understand. In the standard Coulomb-interaction case, the many-body wave-function method must describe the short-range part of the correlation hole around the electron-electron coalescence which requires a lot of one-electron basis functions with high angular momentum. In the range-separation case, the many-body method is relieved from describing the short-range part of the correlation hole, which is instead built in the densityfunctional approximation. The basis set is thus only used to describe a wave function with simply long-range electron-electron correlations (and the one-electron density) which does not require basis functions with very high angular momentum.

In the case of the Coulomb interaction, the rate of convergence of the many-body methods with respect to the size of the basis has been well studied. It has been theoretically shown that, for the ground-state of the helium atom, the partial-wave expansion of the energy calculated by second-order perturbation theory or by full configuration interaction (FCI) converges as L -3 where L in the maximal angular momentum of the expansion [36][37][38][39][40].

Furthermore, this result has been extended to arbitrary atoms in second-order perturbation theory [41,42]. This has motivated the proposal of a scheme for extrapolating the correlation energy to the complete-basis-set (CBS) limit based on a X -3 power-law dependence of the correlation energy on the cardinal number X of the Dunning hierarchical basis sets [43,44].

This extrapolation scheme is widely used, together with other more empirical extrapolation schemes [45][46][47][48][49][50][51]. In the case of range-separated DFT the rate of convergence of the manybody methods with respect to the size of the basis has never been carefully studied, even though the reduced sensitivity to the basis is one of the most appealing feature of this approach.

In this work, we provide a detailed study of the basis convergence of range-separated DFT.

First, we review the theory of range-separated DFT methods (Section II) and we study the convergence of the partial-wave expansion of the long-range wave function near the electronelectron coalescence. We show that the rate of convergence is exponential with respect to the maximal angular momentum L (Section III). Second, we study the convergence of the long-range second-order Møller-Plesset (MP2) correlation energy of four systems (He, Ne, N 2 , and H 2 O) with the cardinal number X of the Dunning basis sets, and find that the error in the correlation energy is best fitted by an exponential in X. This leads us to propose a three-point CBS extrapolation scheme for range-separated DFT based on an exponential formula (Section IV).

Hartree atomic units are used throughout this work.

II. RANGE-SEPARATED DENSITY-FUNCTIONAL THEORY

In range-separated DFT, the exact ground-state energy of an electronic system is expressed as a minimization over multideterminantal wave functions Ψ (see, e.g., Ref. 1)

E = min Ψ Ψ| T + Vne + Ŵ lr,µ ee |Ψ + E sr,µ Hxc [n Ψ ] , (1)
where T is the kinetic-energy operator, Vne is the nuclear-electron interaction operator, E sr,µ Hxc [n Ψ ] is the short-range Hartree-exchange-correlation density functional (evaluated at the density of Ψ), and Ŵ lr,µ ee = (1/2) w lr,µ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 is the long-range electronelectron interaction operator written in terms of the pair-density operator n2 (r 1 , r 2 ). In this work, we define the long-range interaction w lr,µ ee (r 12 ) with the error function

w lr,µ ee (r 12 ) = erf(µr 12 ) r 12 , (2)
where r 12 is the distance between two electrons and µ (in bohr -1 ) controls the range of the separation, with r c = 1/µ acting as a smooth cutoff radius. For µ = 0, the long-range interaction vanishes and range-separated DFT reduces to standard Kohn-Sham DFT. In the opposite limit µ → ∞, the long-range interaction becomes the Coulomb interaction and range-separated DFT reduces to standard wave-function theory. In practical applications, one often uses µ ≈ 0.5 bohr -1 [52,53].

The minimizing wave function Ψ lr,µ in Eq. ( 1) satisfies the Schrödinger-like equation

T + Ŵ lr,µ ee + Vne + V sr,µ Hxc [n Ψ lr,µ ] |Ψ lr,µ = E lr,µ |Ψ lr,µ , (3)
where V sr,µ Hxc is the short-range Hartree-exchange-correlation potential operator (obtained by taking the functional derivative of E sr,µ Hxc ), and E lr,µ is the eigenvalue associated with Ψ lr,µ . In practice, many-body perturbation theory can be used to solve Eq. ( 3). An appropriate reference for perturbation theory is the range-separated hybrid (RSH) approximation [4] which is obtained by limiting the search in Eq. (1) to single-determinant wave functions Φ

E µ RSH = min Φ Φ| T + Vne + Ŵ lr,µ ee |Φ + E sr,µ Hxc [n Φ ] . (4)
The corresponding minimizing wave function will be denoted by Φ µ . The exact ground-state energy is then expressed as

E = E µ RSH + E lr,µ c , (5)
where E lr,µ c is the long-range correlation energy which is to be approximated by perturbation theory. For example, in the long-range variant of MP2 perturbation theory, the long-range

correlation energy is [4] E lr,µ c,MP2 = Φ µ | Ŵ lr,µ ee |Ψ lr,µ 1 , (6)
where Ψ lr,µ 1 is the first-order correction to the wave function Ψ lr,µ (with intermediate normalization). In the basis of RSH spin orbitals {φ µ k }, E lr,µ c,MP2 takes a standard MP2 form Convergence of the truncated partial-wave expansion 1 + f L (r 12 ) for the Coulomb interaction (left) and 1 + f lr,µ L (r 12 ) for the long-range interaction using µ = 0.5 bohr -1 (right) for different values of the maximal angular momentum L. The functions are plotted with respect to the relative angle θ between the position vectors r 1 and r 2 of the two electrons, using r 12 = r 2 1 + r 2 2 -2r 1 r 2 cos θ. We have chosen r 1 = r 2 = 1 bohr, giving r 12 = √ 2 -2 cos θ. In the insert plot on the right, the curves for L = 2, 3, and 4 are superimposed.

E lr,µ c,MP2 = occ i<j vir a<b φ µ i φ µ j | ŵlr,µ ee |φ µ a φ µ b -φ µ i φ µ j | ŵlr,µ ee |φ µ b φ µ a 2 ε µ i + ε µ j -ε µ a -ε µ b , (7)
where φ µ i φ µ j | ŵlr,µ ee |φ µ a φ µ b are the long-range two-electron integrals and ε µ k are the RSH orbital energies. The long-range correlation energy can also be approximated beyond second-order perturbation theory by coupled-cluster [17] or random-phase [22,23,[28][29][30] approximations. Beyond perturbation theory approaches, Eq. (3) can be (approximately) solved using configuration interaction [1,54,55] or multiconfigurational self-consistent field [53,56,57] methods. Alternatively, it has also been proposed to use density-matrix functional approximations for the long-range part of the calculation [58,59].

Since the RSH scheme of Eq. ( 4) simply corresponds to a single-determinant hybrid DFT calculation with long-range Hartree-Fock (HF) exchange, it is clear that the energy E µ RSH has an exponential basis convergence, just as standard HF theory [60]. We will thus focus our study on the basis convergence of the long-range wave function Ψ lr,µ and the long-range MP2 correlation energy E lr,µ c,MP2 .

III. PARTIAL-WAVE EXPANSION OF THE WAVE FUNCTION NEAR ELECTRON-ELECTRON COALESCENCE

In this section, we study the convergence of the partial-wave expansion of the wave function at small interelectronic distances, i.e. near the electron-electron coalescence, which for the case of the Coulomb interaction determines the convergence of the correlation energy.

We first briefly review the well-known case of the Coulomb interaction and then consider the case of the long-range interaction.

A. Coulomb interaction

For systems with Coulomb electron-electron interaction w ee (r 12 ) = 1/r 12 , the electronelectron cusp condition [61] imposes the wave function to be linear with respect to r 12 when

r 12 → 0 [62] Ψ(r 12 ) Ψ(0) = 1 + 1 2 r 12 + O(r 2 12 ). (8)
Here and in the rest of this section, we consider only the dependence of the wave function on r 12 and we restrict ourselves to the most common case of the two electrons being in a natural-parity singlet state [41] for which Ψ(0) = 0. The function

f (r 12 ) = 1 2 r 12 (9)
thus gives the behavior of the wave function at small interelectronic distances. Writing

r 12 = ||r 2 -r 1 || = r 2 1 + r 2 2 -2r 1 r 2 cos
θ where θ is the relative angle between the position vectors r 1 and r 2 of the two electrons, the function f (r 12 ) can be written as a partial-wave expansion

f (r 12 ) = ∞ =0 f P (cos θ), (10)
where P are the Legendre polynomials and the coefficients f are

f = 1 2 1 2 + 3 r +2 < r +1 > - 1 2 -1 r < r -1 > , (11)
with r < = min(r 1 , r 2 ) and r > = max(r 1 , r 2 ). The coefficients f decrease slowly with when r 1 and r 2 are similar. In particular, for r 1 = r 2 , we have f ∼ -2 as → ∞ [START_REF] Stanton | Bartlett Integral packages included are VMOL[END_REF]. Therefore, the approximation of f (r 12 ) by a truncated partial-wave expansion, ≤ L, also converges slowly with L near r 12 = 0. This is illustrated in Figure 1 (left) which shows 1+f L (r 12 ) as a function of θ for r 1 = r 2 = 1 bohr for increasing values of the maximal angular momentum L. Comparing with the converged value corresponding to L → ∞ [Eq. (10)], it is clear that the convergence near the singularity at θ = 0 is indeed painstakingly slow.

f L (r 12 ) = L =0 f P (cos θ), (12)
This slow convergence of the wave function near the electron-electron coalescence leads to the slow L -4 power-law convergence of the partial-wave increments to the correlation energy [36-38, 41, 42] or, equivalently, to the L -3 power-law convergence of the truncation error in the correlation energy [39,40].

B. Long-range interaction

For systems with the long-range electron-electron interaction w lr,µ ee (r 12 ) = erf(µr 12 )/r 12 , the behavior of the wave function for small interelectronic distances r 12 was determined by Gori-Giorgi and Savin [64] Ψ lr,µ (r 12 ) Ψ lr,µ (0) = 1 + r 12 p 1 (µr 12 ) + O(r 4 12 ), (13) where the function p 1 (y) is given by

p 1 (y) = e -y 2 -2 2 √ πy + 1 2 + 1 4y 2 erf(y). ( 14
)
We thus need to study the function

f lr,µ (r 12 ) = r 12 p 1 (µr 12 ). ( 15
)
For a fixed value of µ, and for r 12 1/µ, it yields

f lr,µ (r 12 ) = µ 3 √ π r 2 12 + O(r 4 12 ), (16)
which exhibits no linear term in r 12 , i.e. no electron-electron cusp. On the other hand, for µ → ∞ and r 12 1/µ, we obtain

f lr,µ→∞ (r 12 ) = 1 2 r 12 + O(r 2 12 ), (17)
i.e. the Coulomb electron-electron cusp is recovered. The function f lr,µ (r 12 ) thus makes the transition between the cuspless long-range wave function and the Coulomb wave function.

As for the Coulomb case, we write f lr,µ (r 12 ) as a partial-wave expansion

f lr,µ (r 12 ) = ∞ =0 f lr,µ P (cos θ), (18)
and calculate with Mathematica [65] the coefficients f lr,µ for each

f lr,µ = 2 + 1 2 1 -1
f lr,µ (r 12 )P (x)dx, (19) with x = cos θ, r 12 = r 2 1 + r 2 2 -2r 1 r 2 x, and using the following explicit expression for P (x)

P (x) = 2 k=0 k +k-1 2 x k . ( 20
)
Since the partial-wave expansion of the first term in r 2 12 in Eq. ( 16) terminates at = 1, we expect a fast convergence with of f lr,µ , for µ small enough, and thus also a fast convergence of the truncated partial-wave expansion

f lr,µ L (r 12 ) = L =0 f lr,µ P (cos θ). (21)
Plots of this truncated partial-wave expansion for µ = 0.5 in Figure 1 (right) confirm this expectation. The Coulomb singularity at θ = 0 has disappeared and the approximation f lr,µ L (r 12 ) converges indeed very fast with L, being converged to better than 0.001 a.u. already at L = 2.

We study in detail the dependence of the coefficients of the partial-wave expansion f lr,µ on . We compare two possible convergence behaviors, a power-law form

f lr,µ = A -α , (22)
and an exponential-law form

f lr,µ = B exp(-β ), (23)
where A, B, α and β are (µ-dependent) parameters.

To determine which form best represents f lr,µ , in Figure 2 we plot ln |f lr,µ | for r 1 = r 2 = 1 as a function of ln (left) and as a function of (right), for several values of µ, as well as for the Coulomb case (µ → ∞) [66]. A straight line on the plot of ln |f lr,µ | vs. ln indicates a power-law dependence, whereas a straight line on the plot of ln |f lr,µ | vs. indicates an exponential-law dependence.

For the Coulomb case (black curve nearly superimposed with the curve for µ = 10), we observe that the plot of ln |f | vs. ln is linear, whereas the plot of ln |f | vs. is curved upward. This is expected for a power law A -α form. Moreover, we find α ≈ 2 as expected from Section III A. When going from large to small values of µ, we observe that the plot of ln |f lr,µ | vs. ln becomes more and more curved downward, and the plot of ln |f lr,µ | vs.

becomes more and more linear. We thus go from a power-law dependence to an exponentiallaw dependence. Already for µ ≤ 2, the exponential law is a better description than the power law.

When µ decreases, the absolute value of the slope of the plot of ln |f lr,µ | vs. increases, i.e. the convergence becomes increasingly fast. More precisely, we have found β ≈ 2.598 -

1.918 ln µ for µ ≤ 2.
The exponential convergence of the partial-wave expansion of the long-range wave function near the electron-electron coalescence implies a similar exponential convergence for the partial-wave expansion of the corresponding energy. The present study is thus consistent with the approximate exponential convergence of the partial-wave expansion of the energy of the helium atom in the presence of a long-range electron-electron interaction reported in Refs. 67 and 68. However, no quantitative comparison can be made between the latter work and the present work since the form of the long-range interaction is different.

IV. CONVERGENCE IN ONE-ELECTRON ATOMIC BASIS SETS

In this section, we study the convergence of the long-range wave function and correlation energy with respect to the size of the one-particle atomic basis. This problem is closely related to the convergence of the partial-wave expansion studied in the previous section.

Indeed, for a two-electron atom in a singlet S state, it is possible to use the sphericalharmonic addition theorem to obtain the partial-wave expansion in terms of the relative angle θ between two electrons by products of the spherical harmonic part Y ,m of the oneparticle atomic basis functions

P (cos θ) = 4π 2 + 1 m=- (-1) m Y ,m (θ 1 , φ 1 )Y ,-m (θ 2 , φ 2 ), (24)
where cos θ = cos θ 1 cos θ 2 +sin θ 1 sin θ 2 cos(φ 1 -φ 2 ) with spherical angles θ 1 ,φ 1 and θ 2 ,φ 2 . The partial-wave expansion can thus be obtained from a one-particle atomic basis, provided that the basis saturates the radial degree of freedom for each angular momentum . In practice, of course, for the basis sets that we use, this latter condition is not satisfied. Nevertheless, one can expect the convergence with the maximal angular momentum L of the basis to be similar to the convergence of the partial-wave expansion. (right). For the case of the Coulomb interaction, an essentially exact curve has been calculated

with a highly accurate 418-term Hylleraas-type wave function [73][74][75]. For comparison, we also show the results obtained with the single-determinant HF and RSH wave functions (with the cc-pV6Z basis) which just give horizontal lines since they do not depend on θ. In the insert plot on the right, the V5Z and V6Z curves are superimposed.

A. Convergence of the wave function

We start by analyzing the convergence of the FCI ground-state wave function of the He atom with respect to the cardinal number X of the cc-pVXZ basis sets. We perform a FCI calculation with the long-range Hamiltonian in Eq. ( 3) using a fixed RSH density, calculated from the orbitals obtained in Eq. ( 4

Ψ lr,µ (r 1 , r 2 ) = i≥1 c µ i ϕ µ i (r 1 )ϕ µ i (r 2 ), (25)
where the coefficients c µ i are related to the NO occupation numbers n µ i by the relation

n µ i = 2|c µ i | 2 .
As the signs of c µ i are undetermined we have chosen a positive leading coefficient c µ 1 = n µ 1 /2, and we assumed that all the other coefficients are negative c µ i =n µ i /2 for i ≥ 2 [78]. Even though it has been shown that, for the case of the Coulomb interaction, there are in fact positive coefficients in the expansion in addition to the leading one, for a weakly correlated system such as the He atom, these positive coefficients appear only in larger basis sets than the ones that we consider here and have negligible magnitude [79][80][81].

In Figure 3 we show the convergence of the FCI wave function Ψ lr,µ (r 1 , r 2 ) with the cardinal number X for µ → ∞ which corresponds to the Coulomb interaction (left) and for µ = 0.5 (right). The first electron is fixed at the Cartesian coordinates r 1 = (0.5, 0., 0.) (measured from the nucleus) and the position of the second electron is varied on a circle at the same distance of the nucleus, r 2 = (0.5 cos θ, 0.5 sin θ, 0.). For the Coulomb interaction, we compare with the essentially exact curve obtained with a highly accurate 418-term Hylleraastype wave function [73][74][75]. The curve of Ψ lr,µ (r 1 , r 2 ) as a function of θ reveals the angular correlation between the electrons [82]. Clearly, correlation is much weaker for the longrange interaction. Note that a single-determinant wave function Φ

(r 1 , r 2 ) = ϕ 1 (r 1 )ϕ 1 (r 2 ),
where ϕ 1 is a spherically symmetric 1s orbital, does not depend on θ, and the HF and RSH single-determinant wave functions indeed just give horizontal lines in Figure 3.

The fact that Figure 3 resembles Figure 1 confirms that the convergence with respect to X is similar to the convergence of the partial-wave expansion with respect to L, and thus corroborates the relevance of the study of Section III for practical calculations. As for the partial-wave expansion, the convergence with X of the Coulomb wave function near the electron-electron cusp is exceedingly slow. The long-range wave function does not have an electron-electron cusp and converges much faster with X, the differences between the curves obtained with the cc-pV5Z and cc-pV6Z basis being smaller than 0.03 mhartree. Note, however, that the convergence of the long-range wave function seems a bit less systematic than the convergence of the Coulomb wave function, with the difference between the cc-pVQZ and cc-pV5Z basis being about 3 times larger than the difference between the cc-pVTZ and cc-pVQZ basis. This may hint to the fact the Dunning basis sets have been optimized for the Coulomb interaction and are not optimal for the long-range interaction.

Finally, we note that we have found the same convergence behavior with the short-range exchange-correlation LDA density functional of Ref. 83.

TABLE I. Valence MP2 correlation energies and their errors (in mhartree) for the Coulomb interac-

tion (E c and ∆E c ) and the long-range interaction at µ = 0.5 bohr -1 (E lr,µ c and ∆E lr,µ c ) calculated with Dunning basis sets of increasing sizes for He, Ne, N 2 and H 2 O. The errors are calculated with respect to the estimated CBS limit for the Coulomb interaction and with respect to the cc-pV6Z values for the long-range interaction.

Coulomb interaction

He Ne N 2 H 2 O Basis E c ∆E c E c ∆E c E c ∆E c E c ∆E c cc-pVDZ -

B. Convergence of the correlation energy

We also study the basis convergence of the long-range MP2 correlation energy, given in Eq. ( 7), calculated with RSH orbitals for He, Ne, N 2 and H 2 O.

In Table I we show the valence MP2 correlation energies and their errors as a function of the cardinal number X of the cc-pVXZ basis sets for X ≤ 6. We compare the longrange MP2 correlation energies E lr,µ c,X at µ = 0.5 and the standard Coulomb MP2 correlation energies E c,X corresponding to µ → ∞. For the case of the Coulomb interaction, the error is calculated as ∆E c,X = E c,X -E c,∞ where E c,∞ is the MP2 correlation energy in the estimated CBS limit taken from Refs. 44 and 84. For the range-separated case we do not have an independent estimate of the CBS limit of the long-range MP2 correlation energy for a given value of µ. Observing that the difference between the long-range MP2 correlation energies for X = 5 and X = 6 is below 0.1 mhartree for µ = 0.5, we choose the cc-pV6Z result as a good estimate of the CBS limit. Of course, the accuracy of this CBS estimate will deteriorate for larger values of µ, but in practice this is a good estimate for the range of values of µ in which we are interested, i.e. 0 ≤ µ ≤ 1 [85]. The error on the long-range correlation energy is thus calculated as ∆E lr,µ c,X = E lr,µ c,X -E lr,µ c,6 . The first observation to be made is that the long-range MP2 correlation energies only represent about 1 to 5 % of the Coulomb MP2 correlation energies. Although the long-range correlation energy may appear small, it is nevertheless essential for the description of dispersion interactions for instance. The errors on the long-range MP2 correlation energies are also about two orders of magnitude smaller than the errors on the Coulomb MP2 correlation energies.

Inspired by Ref. 60, we compare two possible forms of convergence for the correlation energy: a power-law form

E lr,µ c,X = E lr,µ c,∞ + AX -α , (26)
and an exponential-law form

E lr,µ c,X = E lr,µ c,∞ + B exp(-βX), (27)
where E lr,µ c,∞ is the CBS limit of the long-range correlation energy and A, B, α, and β are parameters depending on µ, as in Section III. In practice, we actually make linear fits of the logarithm of the error ln(∆E lr,µ c,X ) for the two forms:

ln(∆E lr,µ c,X ) = ln A -α ln X, (28) ln
(∆E lr,µ c,X ) = ln B -βX. (29)
In Table II, we show the results of the fits for the Coulomb interaction and the long-range interaction at µ = 0.5 using either the complete range of X or excluding the value for X = 2.

TABLE II. Results of the fits to the power and exponential laws of the Coulomb valence MP2

correlation energy error ∆E c,X and long-range valence MP2 correlation energy error ∆E lr,µ c,X for µ = 0.5 bohr -1 . Different ranges, X min ≤ X ≤ X max , for the cardinal number X of the Dunning basis sets are tested. The parameters A and B are in mhartree. The squared Pearson correlation coefficients r 2 of the fits are indicated in %. For each line, the largest value of r 2 is indicated in boldface.

Coulomb interaction

Power law

Exponential law We use the squared Pearson correlation coefficient r 2 as a measure of the quality of the fit.

X min X max α A r 2 β B
For the Coulomb interaction and for all the systems studied, the best fit is achieved with the power law AX -α with α ≈ 2.5 -3, which is roughly what is expected [44]. We note however that the fits to the exponential law are also very good with r 2 > 99% when the X = 2 value is excluded. This explains why extrapolations of the total energy based on an exponential formula have also been used for the case of the Coulomb interaction [45,48].

For the case of the long-range interaction, the difference between the fits to the power law and to the exponential law is much bigger. The best fit is by far obtained for the exponential law with r 2 > 99% for all systems and ranges of X considered. This exponential convergence of the long-range correlation energy with respect to X is in accordance with the exponential convergence of the partial-wave expansion of the long-range wave-function observed in Section III.

We have also performed fits for several other values of µ between 0.1 and 1 and always obtained an exponential convergence of the long-range valence MP2 correlation energy with respect to X. However, contrary to what was observed for the partial-wave expansion, we found that when µ decreases β also decreases a bit for the four systems considered. In other words, when the interaction becomes more long range, the convergence of the long-range correlation energy becomes slower. This surprising result may be due to the fact that the cc-pV6Z result may not be as good an estimate of the CBS limit when µ increases. When µ decreases, the prefactor B decreases and goes to zero for µ = 0, as expected. Moreover, we have also checked that we obtain very similar results for the long-range all-electron MP2 correlation energy (including core excitations) with cc-pCVXZ basis sets.

We note that Prendergast et al. [86] have argued that the removal of the electron-electron cusp in a small region around the coalescence point does not significantly improve the convergence of the energy in the millihartree level of accuracy. At first sight, their conclusion might appear to be in contradiction with our observation of the exponential convergence of the long-range correlation energy with X. There are however important differences between the two studies: (1) their form of long-range interaction is different from ours, (2) they consider interelectronic distance "cutoffs" of r c 0.8 bohr whereas we consider larger "cutoffs" r c = 1/µ ≥ 1 bohr, (3) they do not investigate exponential-law versus power-law convergence.

Finally, in the Appendix we provide a complement analysis of the basis convergence of the correlation energy of the He atom for truncated configuration-interaction (CI) calculations in natural orbitals. The analysis shows that, contrary to the case of the Coulomb interaction, the convergence of the long-range correlation energy is no longer limited by the truncation rank the CI wave function but by the basis convergence of the natural orbitals themselves.

This result is consistent with an exponential basis convergence of the long-range correlation energy.

C. Extrapolation scheme

For the long-range interaction case, since both the RSH energy and the long-range correlation energy have an exponential convergence with respect to the cardinal number X, we propose to extrapolate the total energy to the CBS limit by using a three-point extrapolation scheme based on an exponential formula. Suppose that we have calculated three total

energies E X , E Y , E Z for three consecutive cardinal numbers X, Y = X + 1, Z = X + 2. If
we write

E X = E ∞ + B exp(-βX), (30)
E Y = E ∞ + B exp(-βY ), (31)
E Z = E ∞ + B exp(-βZ), (32)
and eliminate the unknown parameters B and β, we obtain the following estimate of the

CBS-limit total energy E ∞ E ∞ = E XY Z = E 2 Y -E X E Z 2E Y -E X -E Z . (33)
In Table III, we report the errors on the RSH+lrMP2 total energy, E µ = E µ RSH + E lr,µ c,MP2 , obtained with the three-point extrapolation formula using either

X = 2, Y = 3, Z = 4 (∆E µ DTQ ) or X = 3, Y = 4, Z = 5 (∆E µ TQ5
), and we compare with the errors obtained with each cc-pVXZ basis set from X = 2 to X = 5 (∆E µ X ). Here again the errors are calculated with respect to the cc-pV6Z total energy, for several values of the range-separation parameter µ = 0.1, 0.5, 1.0, and only valence excitations are included in the MP2 calculations. For all the systems studied the errors ∆E µ DTQ are less than 1.5 mhartree. For Ne, N

= E µ RSH + E lr,µ c,MP2 , obtained with cc-pVXZ basis sets from X = 2 to X = 5 (∆E µ X = E µ X -E µ 6
) and with the three-point extrapolation formula of Eq. ( 33) using

X = 2, Y = 3, Z = 4 (∆E µ DTQ = E µ DTQ -E µ 6 ) or X = 3, Y = 4, Z = 5 (∆E µ TQ5 = E µ TQ5 -E µ 6
). The errors are calculated with respect to the cc-pV6Z total energy for several values of the range-separation parameter µ (in bohr -1 ). Only valence excitations are included in the MP2 calculations. with the errors ∆E µ 5 . Thus, the three-point extrapolation formula with X = 2, Y = 3, Z = 4 provides a useful CBS extrapolation scheme for range-separated DFT. Except for He, the errors ∆E µ TQ5 are negative (i.e. , the extrapolation overshoots the CBS limit) and larger than the errors ∆E µ 5 . Thus, the three-point extrapolation scheme with X = 3, Y = 4, Z = 5 does not seem to be useful. These conclusions extend to calculations including core excitations with cc-pCVXZ basis sets, which are presented in Table IV. All the errors are smaller than for the valence-only calculations. The errors ∆E µ DTQ are now less than 0.9 mhartree, and are smaller or comparable to the errors ∆E µ 5 . The errors ∆E µ TQ5 are always negative and are overall larger than the errors ∆E µ DTQ .

µ ∆E µ D ∆E µ T ∆E µ Q ∆E µ 5 ∆E µ DTQ ∆E µ TQ5 He 0.
= E µ RSH + E lr,µ c,MP2 , obtained with cc-pCVXZ basis sets from X = 2 to X = 5 (∆E µ X = E µ X -E µ 6
) and with the three-point extrapolation formula of Eq. ( 33) using

X = 2, Y = 3, Z = 4 (∆E µ DTQ = E µ DTQ -E µ 6 ) or X = 3, Y = 4, Z = 5 (∆E µ TQ5 = E µ TQ5 -E µ 6
). The errors are calculated with respect to the cc-pCV6Z total energy for several values of the range-separation parameter µ (in bohr -1 ). Core excitations are included in the MP2 calculations. We have also tested a more flexible extrapolation scheme where the RSH energy and the long-range MP2 correlation energy are exponentially extrapolated independently but we have not found significant differences. On the contrary, one may want to use a less flexible two-point extrapolation formula using a predetermined value for β. The difficulty with such an approach is to choose the value of β, which in principle should depend on the system, on the range-separated parameter µ, and on the long-range wave-function method used. For this reason, we do not consider two-point extrapolation schemes.

µ ∆E µ D ∆E µ T ∆E µ Q ∆E µ 5 ∆E µ DTQ ∆E µ TQ5 Ne 0.

V. CONCLUSIONS

We have studied in detail the basis convergence of range-separated DFT. We have shown that the partial-wave expansion of the long-range wave function near the electron-electron coalescence converges exponentially with the maximal angular momentum L. We have also demonstrated on four systems (He, Ne, N 2 , and H 2 O) that the long-range MP2 correla-tion energy converges exponentially with the cardinal number X of the Dunning basis sets cc-p(C)VXZ. This contrasts with the slow X -3 convergence of the correlation energy for the standard case of the Coulomb interaction. Due to this exponential convergence, the extrapolation to the CBS limit is less necessary for range-separated DFT than for standard correlated wave function methods. Nevertheless, we have proposed a CBS extrapolation scheme for the total energy in range-separated DFT based on an exponential formula using calculations from three cardinal numbers X. For the systems studied, the extrapolation using X = 2, 3, 4 gives an error on the total energy with respect to the estimated CBS limit which is always smaller than the error obtained with a single calculation at X = 4, and which is often comparable or smaller than the error obtained with a calculation at X = 5.

We expect the same convergence behavior for range-separated DFT methods in which the long-range part is treated by configuration interaction, coupled-cluster theory, randomphase approximations, or density-matrix functional theory. Finally, it should be pointed out that this rapid convergence is obtained in spite of the fact that the Dunning basis sets have been optimized for the case of the standard Coulomb interaction. The construction of basis sets specially optimized for the case of the long-range interaction may give yet a faster and more systematic convergence. 1s2s2p3s3p3d. Table V shows the Coulomb and long-range correlation energies for the different basis sets and truncation ranks. For the Coulomb interaction, the correlation energy for a fixed rank converges rapidly with the basis size, while the convergence with respect to the rank is much slower. For the long-range interaction, the correlation energy jumps by one order of magnitude when including the 2p natural orbital, which is consistent with the fact that the long-range interaction brings in first angular correlation effects [87].

The long-range correlation energy is essentially converged at rank 1s2s2p, and the overall convergence is now determined by the basis convergence of the natural orbitals.

Finally, we compare two possible forms for the convergence of the truncated CI correlation energies with the cardinal number X, the power law Eq. ( 26) and the exponential law Eq. ( 27). Using as reference the results obtained with the cc-pV6Z basis set, we have calculated for the different truncation ranks the correlation energy errors for the Coulomb

= E c,X -E c,6
, and the long-range interaction, ∆E lr,µ c,X = E lr,µ c,X -E lr,µ c,6 for µ = 0.5 bohr -1 . The range for the cardinal number X of the Dunning basis sets is 2 ≤ X ≤ 6 for the 1s2s and 1s2s2p ranks, and 3 ≤ X ≤ 6 for all the larger ranks. The parameters A and B are in mhartree. The squared Pearson correlation coefficients r 2 of the fits are indicated in %. For each line, the largest value of r 2 is indicated in boldface. interaction, ∆E c,X = E c,X -E c,6 , and for the long-range interaction, ∆E lr,µ c,X = E lr,µ c,X -E lr,µ c,6 , and performed logarithmic fits as in Section IV B. Table VI shows the results of the fits. For both the Coulomb and long-range interactions, for the rank 1s2s2p and larger, the best fit is achieved with the exponential law B exp(-βX). Thus, in comparison with the Coulomb interaction, the long-range interaction does not significantly change the basis convergence of the correlation energy at a fixed truncation rank. However, for the long-range interaction, this exponential convergence at a fixed truncation rank becomes the dominant limitation to the overall basis convergence.

Coulomb interaction

Chapter 3

Self-consistent double-hybrid density-functional theory using the optimized-effective-potential method

In this chapter we wanted to compare two approaches to double-hybrid approximations, combining a fraction of Hartree-Fock exchange and MP2 correlation with exchange and correlation from density-functional approximations. The results of this study has been accepted in The Journal of Chemical Physics and the manuscript is included in this chapter. The global context of this work is the definition of a self-consistent range-separated double-hybrid approximation.

The aim of this first study was to compare two approaches for the one-parameter double-hybrid (1DH) approximation [1]: (i) the standard approach with the MP2 contribution added a posteriori and (ii) a self-consistent approach based on the optimized-effective-potential method and implemented by the team of Prof. Irek Grabowski from the Nicolaus Copernicus University (Poland).

In this work we chose to compare some chemical properties and especially ionization potentials (IP) and electronic affinities (EA). To calculate those properties with standard 1DH the derivatives of the energy are calculated by finite differences [2]

-IP = ∂E 1DH ∂N N -δ ≈ E 1DH (N ) -E 1DH (N -∆) ∆ , -EA = ∂E 1DH ∂N N +δ ≈ E 1DH (N + ∆) -E 1DH (N ) ∆ .
The numerical implementation of IP and EA was done by Bastien Mussard in the Molpro code [3]. The first part of this work was to perform calculations for 1DH and RSH-MP2 for different chemical properties. We then compare the results with the results for the self-consistent 1DH (OEP-1DH) made by our collaborators. To compare these data a work on derivative BIBLIOGRAPHY discontinuity was needed for the OEP-1DH method. The comparison of the RSH-MP2 with self-consistent RSH-MP2 will be part of a later work.

For the standard approach Brillouin's theorem applies and the single excitations do not contribute to second order. However with self-consistency Brillouin's theorem does no longer apply and the single excitations do contribute in the perturbation theory similarly to Görling-Levy perturbation theory [4]. In this study the results of OEP-1DH compared to the standard case only include double excitations. This choice was made to show modifications due to self-consistency without the interfering of further modifications. An interesting open question would be to study the contributions of the single excitations in the OEP-1DH scheme.

I. INTRODUCTION

Density-functional theory (DFT) [1,2] is a powerful approach for electronic-structure calculations of atoms, molecules, and solids. In the Kohn-Sham (KS) formulation, series of approximations for the exchange-correlation energy have been developed for an ever-increasing accuracy: local-density approximation (LDA), semilocal approximations (generalized-gradient approximations (GGA) and meta-GGA), hybrid approximations introducing a fraction of Hartree-Fock (HF) exchange, and nonlocal correlation approximations using virtual KS orbitals (see, e.g., Ref. 3 for a recent review).

In the latter family of approximations, the double-hybrid (DH) approximations are becoming increasingly popular. Introduced in their current form by Grimme [4], they consist in combining a semilocal exchange density functional with HF exchange and a semilocal correlation density functional with second-order Møller-Plesset (MP2) perturbative correlation. Numerous such DH approximations have been developed in the last decade (see Ref. 5 for a review). In general, DH approximations give thermochemistry properties with nearchemical accuracy for molecular systems without important static correlation effects. In virtually all applications of DH approximations, the orbitals are calculated within the generalized KS (GKS) framework [6] (i.e., with a nonlocal HF exchange potential) and without the presence of the MP2 correlation term. The MP2 contribution is then evaluated using the previously self-consistently calculated orbitals and added a posteriori to the total energy.

Recently, Peverati and Head-Gordon [7] proposed an orbital-optimized DH scheme where the orbitals are self-consistently optimized in the presence of the MP2 correlation term. This is a direct extension of orbital-optimized MP2 schemes in which the MP2 total energy is minimized with respect to occupied-virtual orbital rotation parameters [8,9]. Like for regular orbital-optimized MP2 method, also here the optimization of orbitals leads to substantial improvements in spin-unrestricted calculations for symmetry breaking and open-shell situations. Very recently, an approximate orbital-optimized DH scheme was also proposed which confirmed the utility of optimizing the orbitals in complicated electronic-structure problems [10].

In this work, we propose an alternative orbital-optimized DH scheme using the optimized- § abuk@fizyka.umk.pl ¶ ig@fizyka.umk.pl * * eleonora.luppi@upmc.fr † † julien.toulouse@upmc.fr effective-potential (OEP) method [11,12] (see Refs. 13-15 for reviews on OEP). The idea is to optimize the orbitals in the DH total energy expression by using a fully local potential corresponding to the complete exchange-correlation energy expression including the MP2 contribution. This can be considered as an extension of OEP schemes using a second-order correlation energy expression [16][17][18][19][20].

In comparison to the previously-mentioned orbital-optimized DH schemes, we expect that the proposed OEP self-consistent DH scheme to provide additional advantages. First, there is the appeal of staying within the philosophy of the KS scheme with a local potential.

Second, having a local potential and the associated orbital energies can be useful for interpretative purposes. Third, the OEP approach can be advantageous in calculations of excitation energies and response properties, similarly to the advantages of using OEP exact exchange (EXX) versus regular HF. Indeed, contrary to the HF case, with EXX the unoccupied orbitals feel a local potential asymptotically decaying as -1/r, allowing it to support many unoccupied bound states, which are good starting points for calculating high-lying/Rydberg excitation energies and response properties (see, e.g., Refs. 13, 15, and 21).

The paper is organized as follows. In Section II, we review the theory of the standard DH approximations and formulate the proposed self-consistent OEP DH approach. We also explain how the ionization potential and the electron affinity are obtained in both methods. After providing computational details in Section III, we discuss our results in Section IV on total energies, ionization potentials, electron affinities, exchange-correlation and correlation potentials, and correlated densities obtained for a set of atoms (He, Be, Ne, Ar) and molecules (CO and H 2 O). Finally, Section V contains our conclusions.

Throughout the paper, we use the convention that i and j indices label occupied spin orbitals, a and b label virtual ones, and p and q are used for both occupied and virtual spin orbitals. In all equations Hartree atomic units are assumed.

II. THEORY A. Standard double-hybrid approximations

For simplicity, in this work, we consider the one-parameter double-hybrid (1DH) approximation of Ref. 22 in which the density scaling in the correlation functional is neglected. The extension to more general density-scaled or two-parameter double-hybrid approximations [4] is straightforward. The expression of the total energy is thus written as

E = i ϕ * i (x) - 1 2 ∇ 2 + v ne (r) ϕ i (x) dx +E H + E 1DH xc , (1)
where ϕ i (x) are the occupied spin orbitals with x = (r, σ) indicating space-spin coordinates.

In Eq. ( 1), v ne (r) is the nuclei-electron potential,

E H = (1/2) n(x 1 )n(x 2 )/|r 2 -r 1 |dx 1 dx 2
is the Hartree energy written with the spin densities n(x) = i |ϕ i (x)| 2 , and E 1DH xc is the exchange-correlation energy taken as

E 1DH xc = E 1H xc + λ 2 E MP2 c . (2)
In this expression, E 1H xc is the one-parameter hybrid (1H) part of the exchange-correlation energy

E 1H xc = λE HF x + (1 -λ)E DFA x + (1 -λ 2 )E DFA c , (3)
and λ (0 λ 1) is an empirical scaling parameter. The expression of the HF (or exact) exchange energy is

E HF x = - 1 2 i,j ij|ji , (4)
where pq|rs = dx 1 dx 2 ϕ * p (x 1 )ϕ * q (x 2 )ϕ r (x 1 )ϕ s (x 2 )/|r 2 -r 1 | are the two-electron integrals. The expression for the MP2 correlation energy is

E MP2 c = - 1 4 i,j a,b | ij||ab | 2 ε a + ε b -ε i -ε j , (5)
where ij||ab = ij|ab -ij|ba are the antisymmetrized two-electron integrals, and ε p is the energy of the spin orbital p. Finally, E DFA

x and E DFA c are the semilocal density-functional approximations (DFA) evaluated at the spin densities n(x). For example, choosing the Becke 88 (B) exchange functional [23] and the Lee-Yang-Parr (LYP) correlation functional [24] leads to the 1DH-BLYP double-hybrid approximation [22] which is a one-parameter version of the B2-PLYP approximation [4].

We must stress here that the expression of the correlation energy in Eq. ( 5) has a standard MP2 form. However, except for λ = 1, the orbitals are not the HF ones, so the value of the correlation energy in Eq. ( 5) calculated with these orbitals does not correspond to the standard MP2 correlation energy. In the DFT context, the most usual second-order correlation energy expression is given by second-order Görling-Levy (GL2) perturbation theory [16,17], which in addition to the MP2-like double-excitation term also includes a single-excitation term. In this work, following the standard practice for the double-hybrid approximations, we do not include the single-excitation term, which is usually two orders of magnitude smaller than the double-excitation term [25][26][27].

In the standard DH approximations, the spin orbitals are calculated by disregarding the MP2 term in Eq. ( 2) and considering the HF exchange energy as a functional of the one-particle density matrix n 1 (x , x) = i ϕ * i (x)ϕ i (x ), leading to a GKS equation

- 1 2 ∇ 2 + v ne (r) + v H (r) ϕ p (x) + v 1H xc (x, x )ϕ p (x )dx = ε 1H p ϕ p (x), (6)
where v H (r) = n(x )/|r -r|dx is the Hartree potential and v 1H xc (x, x ) is the functional derivative of the three terms in Eq. ( 3) with respect to n 1 (x , x)

v 1H xc (x, x ) = δE 1H xc δn 1 (x , x) = λv HF x (x, x ) + (1 -λ)v DFA x (x)δ(x -x ) +(1 -λ 2 )v DFA c (x)δ(x -x ). (7)
In this expression, v HF x (x, x ) = -n 1 (x, x )/|rr | is the nonlocal HF potential, while v DFA p are thus used in the MP2 correlation expression of Eq. ( 5). We recall that for λ = 0 the 1DH method reduces to the standard KS scheme, while for λ = 1 it recovers the standard MP2 method with HF orbitals. In practice, optimal values of λ are around 0.6-0.8, depending on the density-functional approximations used [22,28,29].

B. Self-consistent OEP double-hybrid approximations

Here, we propose to fully self-consistently calculate the spin orbitals in the DH approximations by taking into account the MP2 term, and considering the HF exchange energy and MP2 correlation energy as implicit functionals of the density. Thus, Eq. ( 6) is replaced by a KS equation

- 1 2 ∇ 2 + v ne (r) + v H (r) + v OEP-1DH xc (x) ϕ p (x) = ε p ϕ p (x), (8)
where v OEP-1DH

xc is a fully local potential obtained by taking the functional derivative with respect to the density of all terms in Eq. ( 2)

v OEP-1DH xc (x) = δE 1DH xc δn(x) = λv EXX x (x) + (1 -λ)v DFA x (x) +(1 -λ 2 )v DFA c (x) + λ 2 v GL2 c (x), (9)
where v EXX (x) must be done with the OEP method, as done in Refs. 18-20. We note that several alternative methods to OEP have been proposed [30][31][32][33][34][35][36], but we do not consider these alternative methods in this work. We will refer to the present approach as the OEP-1DH method. As in the case of the standard DH approach, for λ = 0 the OEP-1DH method reduces to the standard KS scheme. For λ = 1 it reduces to a correlated OEP scheme with the full MP2-like correlation energy expression (but without the single-excitation term), here referred to as the OEP-GL2 scheme.

The OEP equations for the EXX exchange and GL2 correlation potentials

v EXX x (x ) χ s (x , x) dx = Λ x (x), (10)
and

v GL2 c (x ) χ s (x , x) dx = Λ MP2 c (x), (11)
can be obtained after applying a functional-derivative chain rule (see, e.g., Refs. 14, 18, 20, 27, and 37). In these expressions, χ s (x , x) = δn(x )/δv s (x) is the KS static linear-response function which can be expressed in terms of spin orbitals and spin orbital energies,

χ s (x , x) = - i a ϕ * i (x )ϕ a (x )ϕ * a (x)ϕ i (x) ε a -ε i + c.c., (12)
where c.c. stands for the complex conjugate, and

v s (x) = v ne (r) + v H (r) + v OEP-1DH xc (x) is
the total KS potential. The expressions for Λ x (x) and Λ MP2 c (x) are, respectively,

Λ x (x) = δE HF x δv s (x) = i dx δE HF x δϕ i (x ) δϕ i (x ) δv s (x) + c.c. = i,j a ij|ja ϕ * a (x)ϕ i (x) ε a -ε i + c.c. , (13)
and

Λ MP2 c (x) = δE MP2 c δv s (x) = p dx δE MP2 c δϕ p (x ) δϕ p (x ) δv s (x) + c.c. + ∂E MP2 c ∂ε p δε p δv s (x) = 1 2 i,j a,b q =i ij||ab ab||qj ε a + ε b -ε i -ε j ϕ * q (x)ϕ i (x) ε q -ε i + c.c. + 1 2 i,j a,b q =a ij||qb ab||ij ε a + ε b -ε i -ε j ϕ * q (x)ϕ a (x) ε q -ε a + c.c. + 1 2 i,j a,b | ij||ab | 2 (ε a + ε b -ε i -ε j ) 2 |ϕ a (x)| 2 -|ϕ i (x)| 2 . (14)
In practice, in order to solve the OEP equations [Eqs. (10) and ( 11)], the EXX and GL2 potentials are calculated using expansions in a finite Gaussian basis set [38][39][40][41][42]. The EXX potential is thus expanded over orthonormalized auxiliary Gaussian basis functions {g n (r)} as

v EXX x (rσ) = v Slater (rσ) + n c σ x,n g n (r), (15)
where the Slater potential, v Slater (rσ) = -(1/n(x)) dx |n 1 (x, x )| 2 /|rr |, is added to ensure the correct -1/|r| asymptotic behavior of the potential. Similarly, the GL2 potential is expanded as

v GL2 c (rσ) = n c σ c,n g n (r). ( 16
)
Expanding as well the linear-response function in the same basis

χ s (rσ, r σ) = n,m (X σ ) nm g n (r)g m (r ), (17)
and after using Eq. ( 10), the coefficients in Eq. ( 15) are found as

c σ x,n = m (Λ x,σ ) m (X -1 σ ) mn -v σ Slater,n , (18)
where (Λ x,σ ) m = dr g m (r)Λ x (rσ), v σ Slater,n = dr g n (r)v Slater (rσ), and (X -1 σ ) mn are the elements of the (pseudo-)inverse of the matrix X σ . Similarly, after using Eq. ( 11), the coefficients in Eq. ( 16) are found as

c σ c,n = m (Λ MP2 c,σ ) m (X -1 σ ) mn , (19)
where

(Λ MP2 c,σ ) m = dr g m (r)Λ MP2 c
(rσ). In this work, the same basis set is used for expanding the orbitals and the potentials. In practice, our OEP-1DH calculations employ a truncated singular-value decomposition (TSVD) method for the construction of the pseudo-inverse of the linear-response function [used in Eqs. ( 18) and ( 19)] to ensure that stable and physically sound solutions are obtained in the OEP equations [Eqs. (10) and ( 11)].

In principle, this procedure selects the EXX and GL2 potentials which vanish at |r| → ∞.

We note that, when continuum states are included, the GL2 potential actually diverges at infinity for finite systems [43][44][45][46][47]. Nevertheless, this problem is avoided when using a discrete basis set with functions vanishing at infinity (such as the basis set used in this work) [46][47][48]. In practice, the calculated potentials can still be shifted by a function which vanishes at infinity but which is an arbitrary constant in the physically relevant region of space. To remove this arbitrary constant, as in Ref. 18, we impose the HOMO condition on the EXX potential

v EXX x,HH = v HF x,HH , (20)
where

v EXX x,HH = ϕ * H (x)v EXX x (x)ϕ H (x)dx and v HF x,HH = ϕ * H (x)v HF x (x, x )ϕ H (x )dxdx = -j
where

v GL2 c,HH = ϕ * H (x)v GL2 c (x)ϕ H (x)dx and Σ MP2 c,HH (ε H ) = ϕ * H (x)Σ MP2 c (x, x ; ε H )ϕ H (x )
Σ MP2 c,pp (ω) = - 1 2 j a,b | pj||ab | 2 ε a + ε b -ω -ε j + 1 2 i,j b | ij||pb | 2 ω + ε b -ε i -ε j . (22)
Eq. ( 20) can be obtained either by considering the asymptotic limit of Eq. ( 10) and using the fact that the HOMO spin orbital dominates in this limit over all occupied spin orbitals [50],

or by considering the derivative of the HF exchange energy with respect to the electron number (at fixed potential, i.e. at fixed orbitals) and using the chain rule with either the one-particle density or the one-particle density matrix [50,51]. Similarly, Eq. ( 21) can be obtained by considering the derivative of the MP2 correlation energy with respect to the electron number (at fixed potential) and using the chain rule with either the one-particle density or the one-particle Green function [52]. For systems with degenerate HOMO orbitals, we introduce in Eqs. ( 20) and ( 21) sums over the degenerate HOMOs divided by the number of such HOMOs, (1/n H ) H , as done in Ref. 53.

C. Ionization potential and electronic affinity

The ionization potential (IP) and the electronic affinity (EA) can be defined as derivatives of the total energy with respect to the electron number N . For the self-consistent OEP DH approximations, these derivatives can be expressed in terms of frontier spin orbital energies, like in exact KS DFT, [51,52,54,55] -

IP OEP-1DH = ∂E OEP-1DH ∂N N -δ = ε H , (23)
and

-EA OEP-1DH = ∂E OEP-1DH ∂N N +δ = ε L + ∆ xc , (24)
where δ → 0 + , L refers to the LUMO spin orbital, and ∆ xc is the derivative discontinuity of the exchange-correlation energy. For the OEP-1DH approximation, the derivative discontinuity comes from the λ-scaled EXX and GL2 contributions

∆ xc = λ ∆ EXX x + λ 2 ∆ GL2 c . (25)
The terms ∆ EXX x and ∆ GL2 c are given by [56] ∆

EXX x = v HF x,LL -v EXX x,LL -v HF x,HH -v EXX x,HH , (26)
where

v HF x,LL = ϕ * L (x)v HF x (x, x )ϕ L (x )dxdx = -j Lj|jL and v EXX x,LL = ϕ * L (x)v EXX x (x)ϕ L (x)dx
are the expectation values of the HF and EXX exchange potentials over the LUMO spin orbital, and similarly [31,52]

∆ GL2 c = Σ MP2 c,LL (ε L ) -v GL2 c,LL -Σ MP2 c,HH (ε H ) -v GL2 c,HH , (27)
where

Σ MP2 c,LL (ε L ) = ϕ * L (x)Σ MP2 c (x, x ; ε L )ϕ L (x )dxdx and v GL2 c,LL = ϕ * L (x)v GL2 c (x)ϕ L (x)dx
are the expectation values of the MP2 self-energy and of the GL2 local potential over the LUMO spin orbital. Clearly, if the HOMO condition of Eqs. ( 20) and ( 21) is imposed, then the differences of terms in the second parenthesis in Eqs. ( 26) and ( 27) are in fact zero. Note that Eq. ( 27) can be found from the linearized version of the Sham-Schlüter equation [57]. Again, for degenerate HOMOs and/or LUMOs, we introduce in Eqs. (26) and ( 27), sums over the degenerate HOMOs/LUMOs divided by the number of such degenerate HOMOs/LUMOs, i.e. (1/n H ) H and (1/n L ) L .

For standard DH approximations, following Refs. 58 and 59, we obtain IPs and EAs by calculating derivatives of the total energy by finite differences

-IP 1DH = ∂E 1DH ∂N N -δ ≈ E 1DH (N ) -E 1DH (N -∆) ∆ , (28)
and

-EA 1DH = ∂E 1DH ∂N N +δ ≈ E 1DH (N + ∆) -E 1DH (N ) ∆ , (29)
with ∆ = 0.001. To calculate the energies for fractional electron numbers, E 1DH (N -∆) and and

E 1DH (N + ∆),
-IP 1DH = ∂E 1DH ∂N N -δ ≈ ε 1H H + Σ MP2 c,HH (ε 1H H ), (30)
-EA 1DH = ∂E 1DH ∂N N +δ ≈ ε 1H L + Σ MP2 c,LL (ε 1H L ), (31)
which corresponds to standard second-order perturbative propagator theory (see, e.g., Ref. 49). Even though in practice we calculate IP 1DH and EA 1DH using Eqs. ( 28) and ( 29), the approximate connection with the self-energy in Eqs. (30) and ( 31) is useful for comparison and interpretative purposes. For example, it can be shown that Σ MP2 c,HH (ε 1H H ) contains a term corresponding to orbital relaxation in the (N -1)-electron system, and pair-correlation terms for the N -and (N -1)-electron systems [49,61].

III. COMPUTATIONAL DETAILS

The 1DH calculations have been performed with a development version of MOLPRO 2015 [62], and the OEP-1DH ones with a development version of ACES II [START_REF] Stanton | Bartlett Integral packages included are VMOL[END_REF]. In all calculations, we have used the B exchange [23] and the LYP correlation [24] density func-tionals, for E DFA x and E DFA c , respectively. This choice was motivated by the fact that 1DH-BLYP was found to be among the one-parameter double-hybrid approximations giving the most accurate thermochemistry properties on average [22,28,64]. We expect however that the effect of the OEP self consistency to be similar when using other density functional approximations. The performance of both DH methods has been tested against a few atomic (He, Be, Ne, and Ar) and molecular (CO and H 2 O) systems. For the latter, we considered the following equilibrium geometries: for CO d(C-O) = 1.128 Å, and for H 2 O d(H-O) = 0.959 Å and a(H-O-H) = 103.9 • . In all cases, core excitations were included in the second-order correlation term.

In our OEP calculations, for convenience of implementation, the same basis set is used for expanding both the orbitals and the exchange-correlation potential. To ensure that the basis sets chosen were flexible enough for representation of orbitals and exchange-correlation potentials, all basis sets were constructed by full uncontraction of basis sets originally developed for correlated calculations, as in Refs. 65 and 66. In particular, we employed an even tempered 20s10p2d basis for He, and an uncontracted ROOS-ATZP basis [67] for Be and Ne. For Ar, we used a modified basis set [68] which combines s and p basis functions from the uncontracted ROOS-ATZP [67] with d and f functions coming from the uncontracted aug-cc-pwCVQZ basis set [69]. In the case of both molecular systems, the uncontracted cc-pVTZ basis set of Dunning [70] was employed. For all OEP calculations standard convergence criteria were enforced, corresponding to maximum deviations in density-matrix elements of 10 -8 . In practice, the use of the same basis set for expanding both the orbitals and the exchange-correlation potential leads to the necessity of truncating the auxiliary function space by the TSVD method for constructing the pseudo-inverse of the linear-response function. The convergence of the potentials with respect to the TSVD cutoff was studied. Figure 1 shows the example of the convergence of the exchange-correlation and correlation potentials of the Be atom and the CO molecule. For Be, the potentials obtained with the 10 -4 and 10 -6 cutoffs are essentially identical, while for the 10 -8 cutoff the exchangecorrelation potential has non-physical oscillations and the correlation potential diverges. For CO, the potentials obtained with the 10 -4 cutoff are significantly different from the potentials obtained with the 10 -6 cutoff, while no difference can be seen between the potentials obtained with the 10 -6 and 10 -8 cutoffs. A cutoff of 10 -6 was thus chosen for all systems to achieve a compromise between convergence and numerical stability.

In order to assess the quality of the results obtained with the standard and OEP-based DH methods, we considered several reference data. We used estimated exact total energies extracted from numerical calculations [71] for He, Be, Ne, Ar, and from quadratic configuration-interaction calculations extrapolated to the complete basis set [72] for CO and H 2 O. We also used reference data from coupled-cluster singles-doubles with perturbative triples [CCSD(T)] [ [73][74][75][76] calculations performed with the same basis sets. In particular, these CCSD(T) calculations yielded densities which were used as input for generating reference KS potentials by inversion of the KS equations [77,78], using the computational setup described in Refs. 65 and 66. We also used estimated exact total energies extracted from numerical calculations [71] for He, Be, Ne, Ar, and from quadratic configuration-interaction calculations extrapolated to the complete basis set [72] for CO and H 2 O.

IV. RESULTS AND DISCUSSION

A. Total energies

Figure 2 shows the total energy of each system as a function of λ calculated with the 1DH and OEP-1DH approximations. For comparison, CCSD(T) total energies calculated with the same basis sets and estimated exact total energies taken from Refs. 71 and 72 are also reported. Note that the CCSD(T) total energies are significantly higher than the estimated exact energies, which is mostly due to the incompleteness of the basis sets used. Since the explicit density-functional contribution of the DH calculations does not suffer from this large basis incompleteness error, we prefer to use as reference the estimated exact energies. At λ = 0, both DH methods reduce to standard KS using the BLYP functional, which tends to overestimate the total energy by about 3 to 30 mhartree, except for Be atom where it is underestimated (which may be connected to the presence of an important static correlation contribution in this system). At λ = 1, the 1DH approximation reduces to standard MP2 (with HF orbitals), while the OEP-1DH approximation reduces to OEP-GL2 (i.e., the same MP2 total energy expression but with fully self-consistently optimized OEP orbitals) [18,19].

Standard MP2 systematically underestimates the total energy on magnitude (up to more than 100 mhartree for Ar) which is partly due to the missing correlation contribution beyond second order and to the incompleteness of the basis sets used. On the opposite, OEP-GL2 systematically gives too negative total energies, as already known [19,25]. For example, for CO the OEP-GL2 total energy is more than 150 mhartree too low. Note that for Be the OEP-GL2 calculation is unstable, as already reported [19,20,79]. The fact that OEP-GL2

gives much more negative total energies than standard MP2 should be connected to the fact that the HOMO-LUMO orbital energy gap is much smaller with OEP-GL2 orbitals than with HF orbitals (see results in Sections IV B and IV C).

In between the extreme values λ = 0 and λ = 1, the 1DH and OEP-1DH approximation give smooth total energy curves, which start to visually differ for λ 0.2. Note that for Be the OEP-1DH calculations are stable for λ ≤ 0.86. At λ = 0.65, which is the value recommended in Ref. 22 for 1DH with the BLYP functional, both 1DH and OEP-1DH give more accurate total energies than their respective λ = 1 limits (i.e., MP2 and OEP-GL2), but do not perform necessarily better than the λ = 0 limit (KS BLYP). Depending on the system considered, at λ = 0.65, the 1DH total energy is either more accurate or about equally accurate than the OEP-1DH total energy. Thus, even though OEP-1DH provides an important improvement over OEP-GL2, we conclude that the self-consistent optimization of the orbitals in the 1DH approximation (i.e., going from 1DH to OEP-1DH) does not lead to improved ground-state total energies for the few systems considered here. We expect that a similar conclusion generally holds for ground-state energy differences such as atomization energies, similarly to what has been found for the case of the hybrid approximations [80].

B. HOMO orbital energies and ionization potentials

Figure 3 reports, for each system, the HOMO orbital energy in the 1H approximation [Eq. ( 6)] and minus the IPs in the 1DH [Eq. ( 28)] and OEP-1DH [Eq. ( 23)] approximations as a function of λ. The reference IPs are from CCSD(T) calculations with the same basis sets.

The HOMO orbital energy in the 1H approximation represents the simplest approximation to -IP available when doing a 1DH calculation. At λ = 0, the 1H approximation reduces to standard KS with the BLYP functional, and we recover the well-known fact the HOMO orbital energy is much too high (underestimating the IP by about 4 to 8 eV for the systems considered here) with a semilocal DFA like BLYP, most likely due to the self-interaction error in the exchange density functional. At λ = 1, the 1H approximation reduces to standard HF, and in this case the HOMO orbital energy is a much better estimate of -IP, overestimating the IP by about 1 eV except for Be where it is underestimated by about the same amount.

In between λ = 0 and λ = 1, the 1H HOMO orbital energy varies nearly linearly with λ, which suggests that the λ-dependence is largely dominated by the exchange potential in Eq. (7). At λ = 0.65, the 1H HOMO orbital energy is always higher than the reference -IP, by about 1 or 2 eV depending on the system considered.

The IPs obtained with the 1DH method, i.e. taking into account the MP2 correlation term, are smaller than the 1H IPs for all λ and all systems considered here, with the exception of Be for which it is a bit larger. This is consistent with previous works which found that [Eq. ( 28)] and OEP-1DH [Eq. ( 23)] approximations using the BLYP functional as a function of λ.

The reference values were calculated as CCSD(T) total energy differences with the same basis sets.

The vertical lines correspond to λ = 0.65, i.e. the value recommended for 1DH with the BLYP functional in Ref. 22. For Be, the OEP-1DH calculations are unstable for λ > 0.86.

IPs calculating by taking the derivative of the MP2 total energy with respect the electron number are generally too small [58,81]. At λ = 0.65, 1DH gives IPs that are underestimated by about 2 or 3 eV, which is similar to the average accuracy obtained with the two-parameter B2-PLYP double-hybrid approximation [59]. The effect of self-consistency, i.e. going from 1DH to OEP-1DH, is to further reduce the IPs, except for Be for which it increases it. For He the differences between 1DH and OEP-1DH are very small, which may not be surprising since for such a two-electron system the HF and EXX potentials have the same action on occupied orbitals. At λ = 1, OEP-GL2 gives IPs which are generally not very accurate (see also Ref. 82). In particular, for Ne, CO, and H 2 O, OEP-GL2 underestimates the IP 6)] and OEP-1DH [Eq. ( 8)] approximations and minus EAs in the 1DH [Eq. ( 29)] and OEP-1DH [Eq. (24)] approximations using the BLYP functional as a function of λ. The reference EA values were calculated as CCSD(T) total energy differences with the same basis sets, and the reference HOMO energy values were calculated by KS inversion of the CCSD(T) densities. The vertical lines correspond to λ = 0.65, i.e. the value recommended for 1DH with the BLYP functional in Ref. 22. For Be, the OEP-1DH calculations are unstable for λ > 0.86. by more than 3 eV. As a consequence, for these systems, for λ 0.5 self-consistency only deteriorates the accuracy of the IPs. We note that better IPs could be obtained using modified second-order correlated OEP approximations [82].

C. LUMO orbital energies and electronic affinities

Figure 4 reports, for each system, the LUMO orbital energy in the 1H [Eq. ( 6)] and OEP-1DH [Eq. ( 8)] approximations and minus the EAs in the 1DH [Eq. ( 29)] and OEP-1DH [Eq. (24), i.e. including the derivative discontinuity] approximations as a function of λ. The reference EAs are from CCSD(T) calculations with the same basis sets, whereas the reference KS LUMO orbital energies have been obtained by inversion of the KS equations using CCSD(T) densities as input. The reference -EAs are all positive for the systems considered, meaning that the anions are unstable. These positive values are an artifact of the incompleteness of the basis set. In a complete basis set, the EAs should be either negative (i.e., the anion is more stable than the neutral system) or zero (i.e., the anion dissociates into the neutral system and a free electron). Even though the reported EA values are thus not converged with respect to the basis set, the EAs given by different methods can nevertheless be compared for a fixed basis set. By contrast, the reference KS LUMO orbital energies are all correctly negative with the basis set employed. This is due to the fact that the KS LUMO does not represent a state with an additional electron but a bound excited state of the neutral system, which requires much less diffuse basis functions to describe. We note however that, in the case of the CO and H 2 O molecules, the reported KS LUMO orbital energies are not well converged with respect to the basis set due to the lack of diffuse basis functions. Again, we can nevertheless meaningfully compare them with the reference data obtained with the same basis sets.

The LUMO orbital energy in the 1H approximation represents the simplest approximation to -EA (and not to the KS LUMO orbital energy since it is obtained with the nonlocal HF potential) available when doing a 1DH calculation. At λ = 0, the 1H approximation reduces to standard KS with the BLYP functional, and we recover the fact that LUMO orbital energy with a semilocal DFA like BLYP is roughly half way between the exact KS LUMO orbital energy and the exact -EA (i.e., ε DFA L ≈ ε exact L + ∆ xc /2, see e.g. Ref. 83). At λ = 1, the 1H approximation reduces to standard HF, and in this case the LUMO orbital energy is a quite good estimate of -EA (within the finite basis set) for the systems considered here.

At λ = 0.65, the 1H LUMO orbital energy underestimates -EA by about 0.25 to 1 eV, depending on the system. The 1DH approximation gives -EAs rather close to the 1H ones, which indicates that the MP2 correlation term has only a modest effect on this quantity for the systems considered. Coming now to the OEP-1DH results, the LUMO orbital energy obtained in these calculations has a behavior as a function of λ which is clearly distinct from the other curves.

Starting from the BLYP LUMO orbital energy value at λ = 0, it becomes more negative as λ is increased, and becomes an increasingly accurate approximation to the exact KS LUMO orbital energy (and not to -EA). This is an essential difference between having a local EXX (and GL2) potential instead of a nonlocal HF potential. At λ = 0.65, the OEP-1DH LUMO orbital energies underestimate the reference KS LUMO energies by about 1 to 2 eV.

The estimate of -EA in the OEP-1DH approximation is obtained by adding the derivative discontinuity ∆ xc to the OEP-1DH LUMO orbital energy. For the closed-shell systems considered here, the derivative discontinuity is largely dominated by the exchange contribution.

The derivative discontinuity is systematically overestimated in OEP-GL2, leading to -EAs that are much too high. It turns out that, at the recommended value λ = 0.65, OEP-1DH gives -EAs which agree with the reference values within 0.4 eV for the systems considered.

Thus, the OEP-based self-consistency improves the accuracy of EAs.

D. Exchange-correlation and correlation potentials

Figure 5 shows the exchange-correlation potentials calculated by OEP-1DH at the recommended value of λ = 0.65, as well as the potentials obtained at the extreme values of λ, corresponding to KS BLYP (λ = 0) and OEP-GL2 (λ = 1). The reference potentials have been obtained by employing the KS inversion approach using the CCSD(T) densities.

The BLYP exchange-correlation potentials are not negative enough, they do not describe well the shell structure (core/valence transition), and decay too fast at large distances. The OEP-GL2 exchange-correlation potentials have the correct -1/r asymptotic behavior and are quite accurate for the rare-gas atoms (especially for He and Ar), but have too much structure for CO and H 2 O. For Be, the OEP-GL2 calculation is unstable. The OEP-1DH exchange-correlation potentials do not have quite the correct asymptotic behavior since they decay as -λ/r, but for λ = 0.65 they have reasonable shapes in the physically relevant region of space. Note in particular that the OEP-1DH calculation yields a stable solution for Be. For CO and H 2 O, the OEP-1DH exchange-correlation potentials actually improve over both the BLYP and OEP-GL2 exchange-correlation potentials. Therefore, even though the recommended value of λ = 0.65 was determined based on energetical properties of the standard non-self-consistent 1DH scheme, it appears that this value of λ also gives reasonable exchange-correlation potentials as well.

The correlation part of the potentials are plotted in Figure 6. The correlation potentials in BLYP calculations (i.e., the LYP correlation potential evaluated at the self-consistent BLYP density) are unable to reproduce the complex structure of the reference correlation potentials. Note that this is in spite of the fact that LYP correlation energies are usu-

FIG. 6. Correlation potentials calculated with the OEP-1DH approximation [correlation terms in

Eq. ( 9)] using the BLYP functional at the recommended value λ = 0.65, and at the extreme values λ = 0 (standard BLYP) and λ = 1 (OEP-GL2). The reference potentials were calculated by KS inversion of the CCSD(T) densities. For Be, the OEP-GL2 calculation is unstable. For CO, the potential is plotted along the direction of the bond with the C nucleus at -1.21790 bohr and the O nucleus at 0.91371 bohr. For H 2 O, the potential is plotted along the direction of a OH bond with the O nucleus at 0.0 and the H nucleus at 1.81225 bohr. ally reasonably accurate. On the contrary, the OEP-GL2 correlation potentials tend to be largely overestimated, as previously observed [27,65]. Overall, the OEP-1DH correlation potentials at λ = 0.65 have fairly reasonable shapes, providing a good compromise between the understructured BLYP and the overestimated OEP-GL2 correlation potentials. and λ = 1 (OEP-GL2). The reference correlated densities are calculated with CCSD(T). For comparison, correlated densities calculated with standard MP2 are also shown. For Be, the OEP-GL2 calculation is unstable. For CO, the potential is plotted along the direction of the bond with the C nucleus at -1.21790 bohr and the O nucleus at 0.91371 bohr. For H 2 O, the potential is plotted along the direction of a OH bond with the O nucleus at 0.0 and the H nucleus at 1.81225 bohr.

E. Correlated densities

The analysis of the correlated densities provides a useful tool for the detailed examination of the correlation effects on the electronic density and for the test of exchange-correlation approximations in DFT [65,66,[84][85][86][87]. Thus, in Figure 7 we report correlated densities calculated by OEP-1DH at the recommended value of λ = 0.65, as well as the correlated densities obtained at the extreme values of λ, corresponding to KS BLYP (λ = 0) and OEP-GL2 (λ = 1). The correlated density is defined as ∆n c (r) = n(r)n x-only (r) where n(r) is the total density calculated with the full exchange and correlation terms and n x-only (r) is the density calculated with only the exchange terms (see Refs. 65, 66, 84, and 85 for discussions on different definitions of correlated densities). The reference correlated densities have been calculated with the CCSD(T) method as ∆n c,CCSD(T) (r) = n CCSD(T) (r)n HF (r), while n CCSD(T) (r) was obtained from the CCSD(T) relaxed density matrix [88][89][90] constructed using the Lagrangian approach [91][92][93]. For comparison, correlated densities calculated in standard MP2 (1DH with λ = 1) with the same relaxed density-matrix approach are also shown. Due to our current implementation limitations (lack of the relaxed density-matrix approach for the 1DH approximation), the correlated densities for the non-self-consistent 1DH approximation have not been calculated. At λ = 0, i.e. in KS BLYP calculations, the correlated densities are mostly very much underestimated. At λ = 1, the correlated densities are largely overestimated with OEP-GL2. At λ = 0.65, the OEP-1DH correlated densities tend to be quite accurate, achieving a good balance between the underestimated BLYP correlated densities at λ = 0 and the overestimated OEP-GL2 correlated densities at λ = 1. The OEP-1DH correlated densities are overall similar in accuracy to the MP2 and CCSD(T) correlated densities.

V. CONCLUSION

In this work, we have proposed an OEP-based self-consistent DH scheme in which the orbitals are optimized with a local potential including the MP2 correlation contribution. While staying in the philosophy of the KS scheme with a local potential, this scheme constitutes an alternative to the orbital-optimized DH scheme of Peverati and Head-Gordon [7].

We have implemented a one-parameter version of this OEP-based self-consistent DH scheme using the BLYP density-functional approximation and compared it to the corresponding non-self-consistent DH scheme for calculations on a few closed-shell atoms and molecules. While the OEP-based self-consistency does not provide any improvement for the calculations of ground-state total energies and ionization potentials, it does improve the accuracy of electron affinities and restores the meaning of the LUMO orbital energy as being connected to a neutral excitation energy. Moreover, the OEP-based self-consistent DH scheme provides reasonably accurate exchange-correlation potentials and correlated densities. In comparison to the standard OEP-GL2 method [16][17][18][19][20], our OEP-based self-consistent DH scheme is more stable and removes the large overestimation of correlation effects.

Additional work can be foreseen to exploit the full power of the OEP-based self-consistent DH scheme. For example, our scheme should be tested against more systems, including openshell ones for which we expect the OEP self-consistency to provide advantages similar to the orbital-optimized DH scheme of Ref. 7. It would also be interesting to apply linearresponse time-dependent DFT on our OEP-based self-consistent DH scheme to calculate excitation energies. Finally, the present procedure should be applied to the range-separated DH approach [94] which has the advantage of having a fast basis convergence [95].

Chapter 4

Study of the short-range exchange-correlation kernel

Introduction

Time-dependent density-functional theory (TDDFT) [1] and especially the linear-response formalism is one of the most used approach for the calculation of electronic excitation energies. The key quantity in this context is the exchange-correlation kernel which is space and frequency dependent and needs to be approximated. The frequency dependence is very difficult to treat [2] and most of the time this dependence is neglected which is called the adiabatic approximation. Most of the time a further approximation is made by using semilocal density functionals. Despite this crude approximation adiabatic semilocal exchange-correlation kernels give reasonably accurate results in many cases at low computational cost. However, adiabatic semilocal approximations have some limitations: they give too low charge-transfer and Rydberg excitation energies [3,4], they fail to describe double excitations [5,6] or to describe excitation states along bond dissociations coordinates [7].

A way to overcome the limitations of the semilocal approximations is to extend TDDFT to range separation [8]. The description of Rydberg and charge-transfer excitation energies can be improved by introducing a long-range Hartree-Fock exchange kernel [9]. The long-range correlation can be described using long-range linear-response density-matrix-functional theory (DMFT) [10] or long-range linear-response multiconfigurational self-consistent field [11] that give access to double excitations. In practice, in the range-separated TDDFT approaches the short-range exchange-correlation kernel remains described by an adiabatic approximation.

In this study we considered linear-response TDDFT extended to range separation and we studied in particular the frequency dependence of the short-range exchange-correlation kernel.

We were particularly interested by the behavior of the short-range exchange and correlation kernels as a function of the range-separation parameter. The chapter is organized as follows.

A brief review of TDDFT is proposed in Sec. 4.2. The form of the exchange kernel is then explored: the exact-exchange kernel is first generalized to range separation in Sec. 4.3 and the asymptotic behavior of the short-range exact-exchange kernel with respect to the rangeseparation parameter µ is studied in Sec. 4.4. Finally the form of the short-range correlation kernel is explored in the special case of H 2 in a minimal basis in Sec. 4.5.

Review on time-dependent density-functional theory 4.2.1 Time-dependent Schrödinger equation for many-electron systems

The time evolution of a N -particle system in a time-dependent external scalar potential is described by the time-dependent Schrödinger equation that propagates a given initial state

Ψ(t 0 ) = Ψ 0 i ∂ ∂t |Ψ(t) = Ĥ|Ψ(t) , (4.1)
with the N -particle Hamiltonian

Ĥ(t) = T + V (t) + Ŵee , (4.2)
with the kinetic operator T = -1 2 [∇ 2 r n1 (r, r )] r=r dr, the electron-electron interaction Ŵee = 1 2 n2 (r, r )/|r-r |drdr and the time-dependent potential V (t) = v(r, t)n(r)dr that contains the contributions from the electron-nuclei interaction and the external electric-field perturbation.

In these expressions, n1 (r, r ) is the one-particle density matrix operator, n2 (r, r ) is the pairdensity operator and n(r) is the density operator.

The time-dependent many-body Schrödinger equation cannot be solved exactly for most of molecular systems and approximated methods have to be used. One of the most widely used method is time-dependent density-functional theory (TDDFT).

Time-dependent density-functional theory: the Kohn-Sham formalism

The fundamental statement of TDDFT is that the properties of the interacting system can be deduced from the knowledge of the time-dependent density n(r, t) with a fixed initial state condition. This is possible due to the Runge-Gross theorem [1] that can be considered as an extension at the time-dependent level of the first Hohenberg-Kohn theorem. In the timedependent Kohn-Sham (KS) scheme the real system is replaced by a model system of noninteracting particles

i d dt |ϕ i (t) = T + VKS [n](t) |ϕ i (t) (4.3)
with the same density as in the real system n(r, t) = and v xc (r, t) the time-dependent exchange-correlation potential that needs to be approximated.

The most common approximation is the adiabatic approximation.

The adiabatic approximation: In this approximation a functional defined for ground-state DFT is used with a density that corresponds to the time-dependent density at fixed time. This approximation is local in time

v adiabatic xc (r, t) = δE xc [n] δn n=n(r,t) (4.6)
where E xc [n] is the ground-state exchange-correlation density functional.

Linear-response

If the time-dependent external potential is small enough it is not necessary to solve the timedependent Kohn-Sham equation. Calculating the linear change of the density allows us to calculate linear-response properties such as the optical absorption spectrum.

While t < t 0 the time-dependent external potential is switched off and the system is only subject to the nuclear potential and is in its ground state. At t 0 the time-dependent perturbation v (1) is switched on and will induce a change in the density

n(r, t) = n (0) (r) + n (1) (r, t) + n (2) (r, t) + . . . (4.7)
the first-order term of the density depends linearly on v (1) and is sufficient to describe a weak perturbation

n (1) (r, t) = t t 0 dt dr χ(r, t; r , t )v (1) (r , t ) (4.8)
with χ(r, t; r , t ) the linear-response function of the system. The linear change can also be expressed using the Kohn-Sham system which can be expressed in frequency space

n (1) (r, t) = t t 0 dt dr χ 0 (r, t; r , t )v (
χ -1 (r, r , ω) = χ -1 0 (r, r , ω) -f Hxc (r, r , ω) (4.12)
with the Kohn-Sham linear-response function

χ 0 (r, r , ω) = occ. i unocc. a ϕ * i (r)ϕ a (r)ϕ * a (r )ϕ i (r ) i -a + ω + i0 + + occ. i unocc. a ϕ * a (r)ϕ i (r)ϕ * i (r )ϕ a (r ) i -a -ω -i0 + (4.13)
and f Hxc the Hartree-exchange-correlation kernel that needs to be approximated.

Range-separated TDDFT

By decomposing the electron-electron interaction into long-range and short-range parts, the many-electron system can alternatively be replaced by a long-range interacting effective system

i d dt |Ψ lr (t) = T + Ŵ lr ee + V sr (t) |Ψ lr (t) , (4.14)
with the time-dependent long-range wave function Ψ lr which gives the time-dependent density equal to the density of the true system n(r, t). The short-range potential operator is V sr (t) = drn(r)v sr (r, t) with

v sr (r, t) = v(r, t) + v sr Hxc (r, t) (4.15)
and in analogy with the previous section the linear-response theory can then be expressed as

χ -1 (r, r , ω) = (χ lr ) -1 (r, r , ω) -f sr Hxc (r, r , ω) (4.16)
with χ lr the response function of the system described by the Hamiltonian in Eq. (4.14) and the short-range Hartree-exchange-correlation kernel f sr Hxc (r, r , ω) that remains to approximate. The time-dependent exact-exchange (TDEXX) method was introduced by Görling [12,13]. In this section we extend this method to range separation by applying a perturbation theory along the adiabatic connection following the derivation of Ref. [14]. The derivation will be restrained to closed-shell systems. We first consider the range-separated time-dependent Schrödinger equation

i d dt |Ψ lr,λ (t) = T + V sr,λ (t) + λ Ŵ lr ee |Ψ lr,λ (t) (4.17)
with Ŵ lr ee the long-range electron-electron interaction and V sr,λ the potential operator keeping the time-dependent density constant for every value of the coupling constant λ. The initial condition Ψ lr,λ (t 0 ) is chosen to reproduce the initial density n(r, t 0 ) for all λ. When λ = 0 the system is the non-interacting system and can be described by the time-dependent Kohn-Sham

equation i d dt |Φ KS (t) = T + VKS (t) |Φ KS (t) (4.18)
with the initial condition Φ KS (t 0 ) = Ψ lr,λ=0 (t 0 ) and the potential V sr,λ=0 (t) = VKS (t) = V (t) + VHxc (t). At λ = 1 the system is the partially interacting system described by the range-separated time-dependent Schrödinger equation with the potential V sr,λ=1 (t) = V (t) + V sr Hxc (t). The shortrange potential is Taylor expanded with respect to λ

v sr,λ (r, t) = ∞ k=0 λ k v sr,(k) (r, t) = v KS (r, t) + λv sr,(1) (r, t) + λ 2 v sr,(2) (r, t) + . . . (4.19)
The short-range potential can be written as

v sr,λ (r, t) = v(r, t) + v sr,λ Hxc (r, t) = v KS (r, t) -v Hxc (r, t) -v sr,λ Hxc (r, t) , (4.20)
the first-order contribution to the short-range potential is then

v sr,(1) (r, t) = -(v H (r, t) -v sr H (r, t)) -(v x (r, t) -v sr x (r, t)) = -v lr H (r, t) -v lr x (r, t), (4.21)
with v sr H (r, t) and v sr x (r, t) the short-range Hartree and exchange potentials and v lr H (r, t) and v lr x (r, t) the long-range Hartree and exchange potentials. The dependence of the first-order contribution to the short-range potential on the long-range Hartree and exchange potentials can seem counter-intuitive but this can be explained by the long-range nature of the perturbation introduced in Eq. (4.17). The sum of the remaining terms k ≥ 2 corresponds to the negative of the long-range correlation potential. Considering the Taylor expansion of the potential in Eq. (4.19) the Hamiltonian is written as

T + VKS (t) + λ Ŵ lr ee -V lr H (t) -V lr x (t) -λ 2 V lr,(2) (t) + . . . (4.22)
Then we expand the density with respect to λ,

n λ (r, t) = ∞ k=0 λ k n (k) (r, t), (4.23)
every term higher than the zeroth-order being zero because the density remains constant along the adiabatic connection. To express the first-order contribution of the density we finally need to expand the long-range wave function Ψ lr (t) with respect to λ

Ψ lr,λ (t) = ∞ k=0 λ k Ψ lr,(k) (t) = Φ KS (t) + λΨ lr,(1) (t) + λ 2 Ψ lr,( 2
) (t) + . . . (4.24)
In the time-dependent case every Taylor expansion depends on the choice of the initial condition.

The initial-state wave function is chosen to be the ground state of a range-separated timeindependent Schrödinger equation Ψ lr,λ (t 0 ) = Ψ lr,λ 0

T + Û sr,λ + λ Ŵ lr ee |Ψ lr,λ 0 = E λ 0 |Ψ lr,λ 0 , (4.25)
with Û sr,λ is the time-independent potential operator associated with the potential u sr,λ (r) which keeps the density constant along the adiabatic connection. At this point of the derivation the time-dependent potential and the static potential do not need to be equal initially u sr,λ (r) = v sr,λ (r, t 0 ). The expansion of the wave-function in Eq. (4.24) becomes, for the initial state

Ψ lr,λ 0 = ∞ k=0 λ k Ψ lr,(k) 0 = Φ KS 0 + λΨ lr, (1) 0
+ λ 2 Ψ lr, (2) 0
+ . . . (4.26)
Time-dependent perturbation theory [12] yields the first-order contribution to the wave function Ψ lr,λ (t) due to the first-order perturbation in Eq. (

Ψ lr,(1

) (t) = (-i) ∞ k=0 Φ KS k (t) t t 0 dt Φ KS k (t )| Ŵ lr ee -V lr H (t ) -V lr x (t )|Φ KS 0 (t ) + ∞ l=1 Φ KS l (t) Φ KS l | Ŵ lr ee -Û lr H -Û lr x |Φ KS 0 E KS 0 -E KS l , (4.27)
where Φ KS k (t) are the solutions from the time-dependent Kohn-Sham equation that evolve from Φ KS k , the k-th eigenstate of the time-independent Kohn-Sham equation Φ KS k (t 0 ) = Φ KS k . The first-order contribution to the density is then

n (1) (r, t) = Φ KS 0 (t)|n(r)|Ψ lr,(1) (t) + c.c.= 0, (4.28)
where Φ KS 0 (t) is the zeroth-order term of Ψ lr,λ (t) evolving from the initial condition Φ KS (t 0 ) = Φ KS 0 . Combined to the expression of Ψ lr,(1) (t) of Eq. (4.27) it becomes 0 =(-i)

∞ k=0 Φ KS 0 (t)|n(r)|Φ KS k (t) t t 0 dt Φ KS k (t )| Ŵ lr ee -V lr H (t ) -V lr x (t )|Φ KS 0 (t ) + c.c. + ∞ l=1 Φ KS 0 (t)|n(r)| Φ KS l (t) Φ KS l | Ŵ lr ee -Û lr H -Û lr x |Φ KS 0 E KS 0 -E KS l + c.c.. (4.29)
The time-dependent and time-independent Kohn-Sham wave functions are Slater determinants constructed from orbitals ϕ k and ϕ k (t) satisfying the time-independent and time-dependent one-particle Kohn-Sham equations respectively

T + ÛKS |ϕ k = k |ϕ k (4.30) i d dt |ϕ k (t) = T + VKS (t) |ϕ k (t) . (4.31)
We can then express Eq. (4.29) with the Kohn-Sham orbitals and after some algebra we arrive to the range-separated time-dependent exact-exchange equation

t t 0 dt dr χ 0 (r, t; r , t )v lr x (r , t ) = Λ lr x (r, t) (4.32)
with the Kohn-Sham linear-response function

χ 0 (r, t; r , t ) = 2(-i) occ. i unocc. a ϕ * i (r, t)ϕ a (r, t)ϕ * a (r , t )ϕ i (r , t ) + c.c. (4.33)
and the right-hand side of Eq. (4.32)

Λ lr x (r, t) =2(-i) occ. i unocc. a ϕ * i (r, t)ϕ a (r, t) t t 0 dt ϕ a (t )| V NL,lr x (t )|ϕ i (t ) + c.c. (4.34) + 2 occ. i unocc. a ϕ * i (r, t)ϕ a (r, t) × ϕ a | Û NL,lr x -Û lr x |ϕ i i -a + c.c..
The factor 2 is due to the summation over the spin degree of freedom for closed-shell systems.

The non-local (NL) exchange potential corresponds to the Hartree-Fock exchange potential applied to the Kohn-Sham orbitals

V NL,lr x (t)ϕ k (t) (r) = - occ. j ϕ j (r, t)ϕ * j (r , t)ϕ k (r , t)w lr ee (|r -r |)dr (4.35)
and Û NL,lr x and Û lr x are the time-independent long-range non-local and local exchange potentials, respectively.

The short-range exchange potential can be obtained from the difference of the long-range TDEXX equation and the standard TDEXX equation

t t 0 dt dr χ 0 (r, t; r , t ) v x (r , t ) -v lr x (r , t ) = Λ x (r, t) -Λ lr x (r, t) t t 0 dt dr χ 0 (r, t; r , t )v sr x (r , t ) = Λ sr x (r, t) (4.36)
where Λ sr x (r, t) is simply equivalent to Λ lr x (r, t) given in Eq. (4.34) with the long-range electronelectron interaction replaced by the equivalent short-range interaction.

Short-range exact-exchange kernel

To express the short-range exact-exchange (EXX) kernel we need to take the functional derivative of the TDEXX equation (4.36) with respect to the Kohn-Sham time-dependent potential v KS (r, t). The point at which we do the derivation will be v KS (r, t) = u KS (r) such as in the following expressions Vx (t 0 ) = Ûx = Vx . The detailed derivation for the exact-exchange kernel without range separation is presented in Appendix B. With range separation the short-range EXX kernel is given by dr dr χ 0 (r, r , ω)f sr x (r , r , ω)χ 0 (r , r , ω) = h sr x (r, r , ω).

In analogy to the definition of the exact-exchange kernel proposed by Görling [12] we decompose h sr x (r, r , ω) in three terms: h 1,sr x (r, r , ω), h 2,sr x (r, r , ω) and h 3,sr x (r, r , ω) defined as

h 1,sr x (r, r , ω) = -2 occ. i,j unocc. a,b × ϕ * i (r)ϕ a (r) aj|bi sr ϕ * b (r )ϕ j (r ) ( j -b + ω + i0 + )( i -a + ω + i0 + ) + ϕ * a (r)ϕ i (r) ib|ja sr ϕ * j (r )ϕ b (r ) ( j -b -ω -i0 + )( i -a -ω -i0 + ) , (4.38)
h 2,sr x (r, r , ω) = -2 occ. i,j unocc. a,b × ϕ * i (r)ϕ a (r) ab|ji sr ϕ * j (r )ϕ b (r ) ( j -b -ω -i0 + )( i -a + ω + i0 + ) + ϕ * a (r)ϕ i (r) ij|ba sr ϕ * b (r )ϕ j (r ) ( j -b + ω + i0 + )( i -a -ω -i0 + ) , (4.39)
h 3,sr x (r, r , ω) = -2 occ. i,j unocc. a × ϕ * i (r)ϕ a (r) j| V NL,sr x -V sr x |i ϕ * a (r )ϕ j (r ) ( j -a + ω + i0 + )( i -a + ω + i0 + ) + ϕ * a (r)ϕ i (r) i| V NL,sr x -V sr x |j ϕ * j (r )ϕ a (r ) ( j -a -ω -i0 + )( i -a -ω -i0 + ) +2 occ. i unocc. a,b × ϕ * i (r)ϕ a (r) a| V NL,sr x -V sr x |b ϕ * b (r )ϕ i (r ) ( i -b + ω + i0 + )( i -a + ω + i0 + ) + ϕ * a (r)ϕ i (r) b| V NL,sr x -V sr x |a ϕ * i (r )ϕ b (r ) ( i -b -ω -i0 + )( i -a -ω -i0 + ) . (4.40)
The introduced short-range integrals are given by

ij|kl sr = ϕ * i (r)ϕ * j (r )ϕ k (r)ϕ l (r )w sr ee (|r -r |)drdr , (4.41)
with w sr ee (|rr |) = erfc(µ|rr |)/|rr | the short-range interaction. In the following section we will be interested in exploring the frequency dependence of this short-range exact-exchange kernel f sr x for µ → ∞. To do so we will consider a more convenient and compact form of the short-range EXX kernel based on Ref. [15] where the orbitals are real-valued

dr dr ia,jb ϕ ia (r)λ ia (ω)ϕ ai (r )f sr x (r , r , ω)ϕ jb (r )λ jb (ω)ϕ bj (r ) = ia,jb ϕ ia (r)λ ia (ω)X sr ia,jb (ω)λ jb (ω)ϕ bj (r ) (4.42)
with λ ia (ω) = 4 ia /(ω 2 -2 ia ) where ia = ia and ϕ ia (r) = ϕ i (r)ϕ a (r). The function X sr ia,jb (ω) can be decomposed similarly to h sr x in Eq. (4.37), h 1,sr x and h 2,sr x will be given in this new form by

X 1+2,sr ia,jb (ω) = - 1 4 1 + ω 2 ia jb aj|bi sr + 1 - ω 2 ia jb ab|ji sr . (4.43)
The summations run over pairs of occupied/unoccupied orbitals (ia). The next contribution to

X sr ia,jb (ω) will recover h 3,sr x X 3,sr ia,jb (ω) = 1 4 1 + ω 2 ia jb δ ij ϕ a | V NL,sr x -V sr x |ϕ b -δ ab ϕ i | V NL,sr x -V sr x |ϕ j . (4.44)

Asymptotic expansion with respect to the range-separation parameter of the short-range exchange kernel

To explore the behavior of short-range exact-exchange kernel we expand the short-range electronic interaction with respect to the range-separation parameter µ. In the limit of µ going to ∞ the interaction can be described by [16] w sr ee (r

12 ) = µ→∞ 4 √ π m n=0 (-1) n A n n!(n + 2)µ n+2 δ (n) (r 12 ) + O 1 µ m+3 (4.45)
where r 12 = |r 1r 2 | and A n = Γ((n + 3)/2) with Γ the gamma function. In the following we will study the behavior of the short-range exact-exchange kernel when µ goes to ∞.

Leading-order contribution

We consider the first term of the asymptotic expansion of the short-range interaction w sr,(

(r 12 ) = π µ 2 δ(r 12 ). 0) ee
We introduce this interaction in the X sr ia,jb functions defined in Eq.

ϕ ia (r)λ ia (ω)ϕ ai (r )f sr,(0) x (r , r , ω)ϕ jb (r )λ jb (ω)ϕ bj (r ) = ia,jb ϕ ia (r)λ ia (ω)X 1+2,sr,(0) ia,jb (ω)λ jb (ω)ϕ bj (r ) (4.47)
To evaluate X 1+2,sr,(0) ia,jb

(ω) we introduced the leading-order contribution to the short-range interaction in the definition of the short-range integrals

ij|kl sr,(0) = dr π µ 2 ϕ i (r)ϕ j (r)ϕ k (r)ϕ l (r) (4.48) and X 1+2,sr,(0) ia,jb becomes X 1+2,sr,(0) ia,jb = - 1 2 aj|bi sr,(0) = - 1 2 π µ 2 drϕ a (r)ϕ i (r)ϕ b (r)ϕ j (r). (4.49)
We can then rewrite the equation of the kernel as

dr dr ia,jb ϕ ia (r)λ ia (ω)ϕ ai (r )f sr,(0) x (r , r , ω)ϕ jb (r )λ jb (ω)ϕ bj (r ) = ia,jb ϕ ia (r)λ ia (ω) - 1 2 π µ 2 d 3 r ϕ a (r )ϕ i (r )ϕ b (r )ϕ j (r ) λ jb (ω)ϕ jb (r ), (4.50)
and by identification between the right-hand and the left-hand sides of the equation we can determine the expression of the leading-order contribution to short-range exchange kernel

f sr,(0) x (r , r , ω) = - 1 2 π µ 2 δ(r -r ). (4.51)

Next-order contribution

We now consider the next-order contribution of the asymptotic expansion of the short-range interaction w sr,(

(r 12 ) = - 4 3 √ π µ 3 δ (1) (r 12 ) 1) ee
and introduce it in the expression of X sr ia,jb . The next-order contribution vanishes because of the integration over the angular degrees of freedom

dΩ 1 dr 12 δ (1) (r 12 )f (r 12 ) = -dΩ 1 u 12 • ∇f | r 12 =0 = 0, (4.53)
where Ω 1 stands for the angular coordinates associated to the variable r 1 . We can then express the kernel as

f sr x (r 2 , r 3 , ω) = µ→∞ - 1 2 π µ 2 δ(r 2 -r 3 ) + O 1 µ 4 (4.54)
This result is consistent with results obtained from the study of the asymptotic expansion of the exact ground-state short-range spin-independent exchange density functional for µ → ∞ [17] and shows that the leading order term of the exact-exchange kernel is local in space and does not depend on the frequency.

Examples of H 2 and He

In this section we calculate the projection of the short-range exchange kernel in the molecular orbital basis

(ia|f sr x |jb) = ϕ i (r)ϕ a (r)f sr x (r, r )ϕ j (r )ϕ b (r )drdr (4.55)
for two-electron systems in cc-pVDZ (abbreviated as VDZ) basis sets [18]. For two-electrons systems the exact-exchange kernel can be related to the Hartree kernel as f EXX x = -1/2f H .

CHAPTER 4. STUDY OF SHORT-RANGE EXCHANGE-CORRELATION KERNEL

For two-electron systems, the EXX kernel is thus frequency independent. We will compare the behavior of three kernels: the short-range LDA exchange kernel [17], the short-range exactexchange kernel and the short-range exact-exchange kernel when µ → ∞ defined in Eq. (4.54).

The calculations were performed using the quantum chemistry program MOLPRO [19]. The orbitals are calculated at the RSH level [20]. x |ia) with respect to the range-separation parameter µ for H 2 at equilibrium (1.4 bohr) (top) and in the dissociation limit (distance of 10 bohr) (bottom) calculated with the short-range LDA, shortrange exact-exchange and short-range exact exchange in the µ → ∞ limit. For each distance two transitions are shown: the first transition 1σ g → 1σ u and a higher transition 1σ g → 3π g .

The results for H 2 in VDZ are presented in Fig. 4.1. In this basis set we have one occupied orbital and nine unoccupied orbitals, we only considered the diagonal terms of the kernel matrix.

In the limit of high µ we notice that the behavior of the three kernels are the same and that the leading term of the expansion Eq. (4.54) gives a good description in this region. On the other end of the curve, for µ → 0 we can notice the difference of behavior between srLDA and srEXX: srEXX seems to have a "flatter" start than srLDA at µ → 0. This difference of behavior could be explained if we consider the behavior of both potentials (srLDA [17] and srEXX [21]) for µ → 0

v sr,EXX x (r) = v x (r) + 2µ √ π + µ 3 v sr, (3)
x (r) + O µ 5 (4.56)

v sr,LDA x (r) = v LDA x (r) + µ √ π - αµ 2 r s (r) π + O e -1/µ 2 , (4.57)
with r s = (3/(4πn)) 1/3 the Wigner-Seitz radius and α = (4/(9π)) 1/3 . By taking the second-order derivative of both potentials with respect to the density the corresponding kernels are

f sr,EXX x (r, r ) = f EXX x (r, r ) + µ 3 f sr, (3)
x (r, r ) + O µ 5 (4.58)

f sr,LDA x (r, r ) = f LDA x (r)δ(r -r ) + πα 4 r 4 s (r)µ 2 + O e -1/µ 2 , (4.59) with f LDA x (r) = ∂ 2 e LDA x
/∂n 2 the second-order derivative of the LDA exchange energy density with respect to the density. The srLDA exchange kernel in the limit of small µ should thus be quadratic while the short-range EXX exchange kernel should behave as µ 3 . This seems consistent with the behavior observed for H 2 at equilibrium (however in the dissociation limit the variation of the short-range EXX kernel seems to be faster for the first transition 1σ g → 1σ u ).

It is interesting to remark that for localized valence excitation (1σ g → 1σ u at equilibrium) srLDA gives good results with respect to srEXX even for µ close to 0. For Rydberg excitations (1σ g → 3π g at equilibrium and in the dissociation limit) and delocalized valence excitation (1σ g → 1σ u in the dissociation limit) there is a loss of accuracy. We finally considered a second system, the helium atom in the VDZ basis set, and the results are presented in Fig. 4.2. The results for He are similar to H 2 : the local valence excitations are well described by srLDA.

Study of the exact frequency-dependent correlation kernel

In this section we want to extend our study of the short-range exchange kernel to the short-range correlation kernel. The analytic derivation of the exact correlation kernel in a similar way to Sec. 4.3 is complicated due to the definition of the correlation potential v c = ∞ i=2 -v (i) . An example of exploring the orbital-dependent correlation functionals with many-body perturbation theory at second order is given by Bokhan and Bartlett [22]. A way to study the frequency dependence is to consider the exact correlation kernel for a model system: H 2 in minimal basis set.

FCI of H 2 in minimal basis set

The derivation of the short-range correlation kernel in this section will be based on a fullconfiguration interaction (FCI) calculation in a minimal basis set. For H 2 the minimal basis set will be composed of two orthogonal orbitals ϕ 1 and ϕ 2 localized on the two hydrogen atoms.

The different determinants for H 2 in this minimal basis set are presented in Fig. 4.3. the different terms in the full-range Coulomb case can be taken from Ref. [23]. In the following derivation of the short-range correlation kernel we will consider these states as starting point to evaluate the linear-response function for both the full-range and long-range cases.

|1 1 |1 2 | 12 |12 | 12 |2 2
|0 = c 1 |1 1 + c 2 |2 2 E 0 = 2 1 -J 11 + ∆ -∆ 2 -K 2 12 |1 = c 1 |1 1 + c 2 |2 2 E 1 = 2 1 -J 11 + ∆ + ∆ 2 -K 2 12 |2 = 1 √ 2 (|1 2 + |2 1 ) E 3 = h 11 + h 22 + J 12 + K 12 |3 = |12 E 3 = h 11 + h 22 + J 12 -K 12 |4 = | 12 E 4 = h 11 + h 22 + J 12 -K 12 |5 = 1 √ 2 (|1 2 -|2 1 ) E 5 = h 11 + h 22 + J 12 -K 12

Calculation of the linear-response function

We want to calculate the exact linear-response function

χ(x 1 , x 2 ; x 1 , x 2 ; ω) = K =0 0| Ψ † (x 1 ) Ψ(x 1 )|K K| Ψ † (x 2 ) Ψ(x 2 )|0 ω -(E K -E 0 ) + i0 + (4.61) - K =0 0| Ψ † (x 2 ) Ψ(x 2 )|K K| Ψ † (x 1 ) Ψ(x 1 )|0 ω + (E K -E 0 ) -i0 + (4.62)
where

Ψ(x) = i φ i (x)â i , Ψ † (x) = i φ * i (x)â † i , |0
is the ground-state wave function and |K are the excited-states wave functions. It can be expressed in the basis of molecular orbitals [χ(ω)] pq,rs = dx 1 dx 1 dx 2 dx 2 φ p (x 1 )φ * q (x 1 )χ(x 1 , x 2 ; x 1 , x 2 ; ω)φ * r (x 2 )φ s (x 2 ) (4.63) and using Eq. (4.61) it gives

[χ(ω)] pq,rs = K =0 0|â † p âq |K K|â † s âr |0 ω -∆E K + i0 + - K =0 0|â † s âr |K K|â † p âq |0 ω + ∆E K -i0 + (4.64)
with ∆E K = E K -E 0 . The linear-response function can be spin-adapted [24,25,26] where p, q, r and s refer now to spatial orbitals. We can express the singlet linear-response function matrix as

1 χ(ω) pq,rs = 1 2   K =0 0| Êpq |K K| Êsr |0 ω -∆E K + i0 + - K =0 0| Êsr |K K| Êpq |0 ω + ∆E K -i0 +   (4.67)
with Êpq = â † p↑ âq↑ + â † p↓ âq↓ . The triplet contributions 3,0 χ, 3,1 χ and 3,-1 χ can be obtained by using √ 2 T 1,0 = â † p↑ âq↑â † p↓ âq↓ , T 1,1 = -â † p↑ âq↓ and T 1,-1 = â † p↓ âq↑ , respectively, instead of Êpq . The matrix form of the singlet linear-response function is, in the occupied/virtual and virtual/occupied block,

1 χ(ω) =     c 2 1 ω -∆E 2 + i0 + - c 2 2 ω + ∆E 2 -i0 + c 1 c 2 ω -∆E 2 + i0 + - c 1 c 2 ω + ∆E 2 -i0 + c 1 c 2 ω -∆E 2 + i0 + - c 1 c 2 ω + ∆E 2 -i0 + c 2 2 ω -∆E 2 + i0 + - c 2 1 ω + ∆E 2 -i0 +     . (4.68)
In the non-interacting limit (c 1 → 1, c 2 → 0), the Kohn-Sham response function, Eq. (4.13), is obtained. The inverse singlet linear-response function is

1 χ -1 (ω) =     (ω -∆E 2 )c 2 1 -(ω + ∆E 2 )c 2 2 (c 2 1 -c 2 2 ) 2 2∆E 2 c 1 c 2 (c 2 1 -c 2 2 ) 2 2∆E 2 c 1 c 2 (c 2 1 -c 2 2 ) 2 (ω -∆E 2 )c 2 2 -(ω + ∆E 2 )c 2 1 (c 2 1 -c 2 2 ) 2     . (4.69)
Similarly, the matrix form of the triplet linear-response functions is

3,0 χ(ω) =     c 2 1 ω -∆E 5 + i0 + - c 2 2 ω + ∆E 5 -i0 + - c 1 c 2 ω -∆E 5 + i0 + + c 1 c 2 ω + ∆E 5 -i0 + - c 1 c 2 ω -∆E 5 + i0 + + c 1 c 2 ω + ∆E 5 -i0 + c 2 2 ω -∆E 5 + i0 + - c 2 1 ω + ∆E 5 -i0 +     , (4.70)
and the inverse is

3,0 χ -1 (ω) =     -c 2 2 (ω -∆E 5 ) + c 2 1 (ω + ∆E 5 ) (c 2 1 -c 2 2 ) 2 -2c 1 c 2 ∆E 5 (c 2 1 -c 2 2 ) 2 -2c 1 c 2 ∆E 5 (c 2 1 -c 2 2 ) 2 c 2 2 (ω -∆E 5 ) -c 2 1 (ω + ∆E 5 ) (c 2 1 -c 2 2 ) 2     . (4.71)
The χ p↑q↓,r↓s↑ and χ p↓q↑,r↑s↓ elements are equal to 0, 3,-1 χ and 3,1 χ are then equal and, as the energies of the triplet states are degenerate, we can see that 3,0 χ = 3,1 χ = 3,-1 χ.

Derivation of the exact short-range correlation kernel

The short-range kernel can be deduced from the difference of the inverse of the long-range linear-response function and the inverse of the full-range linear-response function, as presented in Eq. (4.16). The long-range linear-response function has the same expression as the full-range

3 f sr H = 0 0 0 0 . (4.78)
The singlet short-range exchange-correlation kernel is then

1 f sr xc (ω) =  ((c lr 1 ) 2 -(c lr 2 ) 2 ) 2   + ω   -1 (c 2 1 -c 2 2 ) 2 + 1 ((c lr 1 ) 2 -(c lr 2 ) 2 ) 2 0 0 1 (c 2 1 -c 2 2 ) 2 - 1 ((c lr 1 ) 2 -(c lr 2 ) 2 ) 2   . (4.80)
We will consider two ways of calculating the short-range exchange kernel. We will first consider the HF exchange kernel If the exchange contribution is the short-range Hartree-Fock exchange kernel the singlet and triplet short-range correlation kernels are

((c lr 1 ) 2 -(c lr 2 ) 2 ) 2   + ω   -1 (c 2 1 -c 2 2 ) 2 + 1 ((c lr 1 ) 2 -(c lr 2 ) 2 ) 2 0 0 1 (c 2 1 -c 2 2 ) 2 - 1 ((c lr 1 ) 2 -(c lr 2 ) 2 ) 2   . (4.87)
In all cases if we compare both kernels we can see that the frequency dependence is the same for the singlet and the triplet kernels.

Calculations on H 2 in STO-3G basis

We calculate the different terms of the singlet short-range correlation kernel [Eqs. (4.84) and (4.86)] by performing calculations using the MOLPRO software [19] with the STO-3G basis set.

The long-range coefficients and energies are calculated with FCI with long-range Hamiltonian using a fixed RSH density [20] and full-range coefficients and energies are calculated with FCI with the full-range Hamiltonian. The Coulomb and exchange integrals J 12 and K 12 are extracted from the TDDFT module. We decompose the kernel matrix in two matrices: the frequencyindependent term 1 f c (ω = 0) and the frequency-dependent term common to both singlet and triplet kernel ∆ 1 f c (ω) = 1 f c (ω) -1 f c (ω = 0). In the calculations the frequency-dependent kernel is evaluated at ω = 1 (hartree) as representative of the order of magnitude of electronic excitations in atoms and molecules. The results are presented in Fig. 4.4. We first consider the evolution of the kernel elements with respect to the range-separation parameter µ and we can notice that the amplitude of the frequency-dependent term is small compared to the elements of the frequency-independent matrix. We considered the ratio of |∆ 1 f c (ω = 1)/ 1 f c (ω = 0)| and we can notice that for µ ≥ 0.9 bohr -1 for postHF and µ ≥ 1.8 bohr -1 for postEXX the ratio is constant showing the similar asymptotic behavior of the frequency-dependent and frequencyindependent contributions of the singlet kernel. The ratio converges to 0.18 this value is small and the adiabatic approximation is thus good for the short range correlation kernel but it remains an approximation. We plotted in Fig. 4.5 the product f c × µ 2 . The fact that the curves go to some constants for µ → ∞ indicates that the asymptotic behavior of the contributions to the correlation kernel behaves as 1/µ 2 . We then considered the evolution of the kernel elements with respect to the internuclear distance R at µ = 0.4 bohr -1 which is an optimal value in range-separated TDDFT. The result is presented in Fig. 4.6. We notice that all the terms have a similar behavior and diverge for large R. We studied the behavior of the correlation term while stretching the H-H bond (at R = 10 bohr) we first notice that there is no longer a difference between postHF and postEXX terms, the leading term of the kernel is given by the difference of the coefficients c 1 and c 2 (and c lr 1 and c lr 2 ) that are nearly degenerate and make the terms diverge. Thus the amplitude of the kernel is more important than for the equilibrium. The ratio |∆ 1 f c (ω = 1)/ 1 f c (ω = 0)| converges slower and to a higher value than at equilibrium ( 1) which implies that in this case the frequency-dependent term is more important and that the adiabatic approximation is less accurate. Finally we can see in Fig. 4.5 that the asymptotic behavior of the frequency-dependent and frequency-independent terms remain as 1/µ 2 . To extend this work it would be interesting to study the effect of this exact short-range correlation kernel on the bond-breaking 1 Σ g -1 Σ u excitation energy to see if it could compensate the decay to zero of the excitation energy observed for the adiabatic approximation.

Conclusion

In this work, we have studied the short-range exchange and correlation kernels. We first extended the EXX time-dependent density-functional theory to range separation and studied the asymptotic expansion of the short-range EXX kernel as a function of the range-separation parameter.

We showed that the two first terms of this asymptotic expansion are frequency-independent and local in space. We then compared the performance of this short-range EXX kernel to the short-range LDA kernel for two simple systems: H 2 and He in the VDZ basis set. We showed that in the limit of large µ the semilocal approximation becomes exact but when µ goes to 0 the semilocal approximation has limitations with a good behavior for localized valence excitation energies but less accurate results for Rydberg and delocalized valence excitation energies.

In the second part of this chapter we studied the short-range correlation kernel by studying the exact short-range correlation kernel for a model system: H 2 in a minimal basis set. We first derived the short-range correlation kernel and noticed that for this system the frequency dependence is the same for both singlet and triplet kernels and that it does not depend on the nature of the exchange contribution (HF or EXX). We continued this study with practical calculations for H 2 in the STO-3G basis and we showed that even in the limit of large µ the adiabatic approximation is not exact. This frequency dependent term is relatively small at equilibrium but becomes large in the dissociation limit.

In further work the study of the short-range exchange kernel could be extended by deriving the next term of the asymptotic expansion. This term would be in 1/µ 4 and should be frequency

Conclusion

In this thesis we were interested in the study of hybrid approximations between wave-function based methods and density-functional theory considering two ways of decomposing the electronelectron interaction: linearly or with range separation. We performed three separate studies.

The first chapter presented in this work was focused on the basis-set convergence in the range-separated scheme. We performed two studies to compare and characterize the basis-set convergence in full-range and range-separated methods. We have first shown that the partialwave expansion of the long-range wave function near the electron-electron coalescence converges exponentially with the maximal angular momentum L. We then have shown that the long-range correlation energy error evaluated with the second-order Møller-Plesset perturbation theory was converging exponentially with respect to the cardinal number of the Dunning basis sets cc-pVXZ. This is an acceleration with respect to the full-range case where the convergence is polynomial (X 3 ). While studying the basis-set convergence we noticed that the Dunning basis sets may not be optimal to perform calculations with range-separated methods. The results obtained in this chapter are expected to be similar when using other methods to describe the long-range correlation such as configuration interaction, coupled-cluster theory or random-phase approximations.

The second chapter was focused on double-hybrid methods and particularly on a self-consistent way to include the MP2 correlation energy using an optimized-effective-potential (OEP) method.

The aim of this study was to compare the performance of such self-consistent double-hybrid approximation to the standard double-hybrid approximation. A one-parameter self-consistent double-hybrid approximation has first been implemented. We considered a set of atoms and molecules and compared some properties such as ionization potentials and electronic affinities.

We observed no improvement for total energies and ionization potentials but we obtained good estimates for LUMO orbital energies and good accuracy for electronic affinities. An interesting result is that this self-consistent approximation give reasonably good accuracy for the exchangecorrelation and correlation potentials and correlated densities.

In the last chapter we considered range-separated linear-response time-dependent densityfunctional theory and we studied the short-range exchange and correlation kernels. After ex-

A.2 Extrapolation scheme

In Chapter 2, based on the exponential convergence of both the RSH and long-range MP2 energies we proposed an extrapolation formula for the total energy in the complete-basisset (BSE) limit based on three calculations in three successive basis sets of cardinal number X, Y and Z (Y = X + 1 and Z = Y + 1) given by

E ∞ = E XY Z = E 2 Y -E X E Z 2E Y -E X -E Z . (A.1)
In this section we compare two ways to apply the extrapolation scheme. The first one, presented in the chapter 2 is simply to apply the formula of Eq. (A.1) to the total RSH+lrMP2 energy , obtained with cc-pVXZ basis sets form X = 2 to X = 5 (∆E X = E X -E 6 ) and the errors obtained with the three point extrapolation formula of Eq. (A.1) using X = 2, Y = 3, Z = 4 (∆E DTQ = E DTQ -E 6 evaluated for the total energy (E RSH+lrMP2 ) or evaluated for RSH and lrMP2 separately and recombined a posteriori (E RSH + E lrMP2 c ). The errors are calculated with respect to the cc-pV6Z total energy for several values of the range-separation parameter µ (in bohr -1 ). Only valence excitations are included in the MP2 calculations. , obtained with cc-pVXZ basis sets form X = 2 to X = 5 (∆E X = E X -E 6 ) and the errors obtained with the three point extrapolation formula of Eq. (A.1) using X = 2, Y = 3, Z = 4 (∆E DTQ = E DTQ -E 6 ) evaluated for the total energy (E RSH+lrMP2 ) or evaluated for RSH and lrMP2 separately and recombined a posteriori (E RSH + E lrMP2 c ). The errors are calculated with respect to the cc-pV6Z total energy for several values of the range-separation parameter µ (in bohr -1 ). Core and valence excitations are included in the MP2 calculations.

E total XY Z = E RSH+lrMP2 XY Z , (A.2)
µ ∆E µ D ∆E µ T ∆E µ Q ∆E µ 5 ∆E µ DTQ ∆E µ DTQ E RSH+lrMP2 E RSH + E lrMP2
µ ∆E µ D ∆E µ T ∆E µ Q ∆E µ 5 ∆E µ DTQ ∆E µ DTQ E RSH+lrMP2 E RSH + E lrMP2 c Ne 0.1 79

Derivation of the exact-exchange kernel

In this appendix, we present a detailed derivation of the exact-exchange (EXX) kernel [1]. This has to be read together with Chapter 4.

To express the EXX kernel we need to take the functional derivative with respect and the right-hand-side term We now consider the calculation of the right-hand side of Eq. (B.12)

Λ x (r, t) =2(-i) occ. i unocc. a ϕ * i (r, t)ϕ a (r, t) t t 0 dt ϕ a (t )| V NL x (t )|ϕ i (t ) + c.c. (B.3) + 2 occ. i unocc. a ϕ * i (r, t)ϕ a (r, t) × ϕ a | Û lr x -Ûx |ϕ i i -a + c.c.
k (r, t) δv KS (r , ω) = e -iωt e ηt all l ϕ l (r)e -i k (t-t 0 ) ϕ * l (r )ϕ k (r ) k -l + ω + iη (B.
χ 0 (r, r , ω) = 2 occ. i unocc. a ϕ * i (r)ϕ a (r)ϕ * a (r )ϕ i (r ) i -a + ω + iη + 2 occ. i unocc. a ϕ * a (r)ϕ i (r)ϕ * i (r )ϕ a (r ) i -a -ω -iη . (B.
δΛ x (r, t) δv KS (r , ω) B - t t 0 dt dr δχ 0 (r, t; r , t ) δv KS (r , ω) v x (r , t ) C (B.20)
and for every term we first insert the derivatives of the orbitals given in Eq. (B.10), then after integrating upon t and some algebra, we obtain the following equations

C = 2e -iωt e ηt occ. i unocc. a all l × ϕ * l (r)ϕ a (r) a| Vx |i ϕ * i (r )ϕ l (r ) ( i -l -ω -iη)( i -a ) + ϕ * i (r)ϕ l (r) a| Vx |i ϕ * l (r )ϕ a (r ) ( a -l + ω + iη)( i -a ) + ϕ * i (r)ϕ a (r) l| Vx |i ϕ * a (r )ϕ l (r ) ( a -l -ω -iη)( i -a + ω + iη) + ϕ * i (r)ϕ a (r) a| Vx |l ϕ * l (r )ϕ i (r ) ( i -l + ω + iη)( i -a + ω + iη) - ϕ * a (r)ϕ l (r) i| Vx |a ϕ * l (r )ϕ i (r ) ( i -l + ω + iη)( a -i ) + ϕ * l (r)ϕ i (r) i| Vx |a ϕ * a (r )ϕ l (r ) ( a -l -ω -iη)( a -i ) + ϕ * a (r)ϕ i (r) i| Vx |l ϕ * l (r )ϕ a (r ) ( a -l + ω + iη)( a -i + ω + iη) + ϕ * a (r)ϕ i (r) l| Vx |a ϕ * i (r )ϕ l (r ) ( i -l -ω -iη)( a -i + ω + iη) (B.21)
where the summation over i, j runs over the occupied orbitals, a, b the unoccupied orbitals and l all the orbitals. We can decompose B in three contributions as follows

B = 2(-i) occ. i unocc. a ϕ * i (r, t)ϕ a (r, t) δ t t 0 ϕ a (t )| V NL x (t )|ϕ i (t ) dt δv KS (r , ω) + c.c B 1 + 2(-i) occ. i unocc. a δ (ϕ * i (r, t)ϕ a (r, t)) δv KS (r , ω) t t 0 ϕ a (t )| V NL x (t )|ϕ i (t ) dt + c.c B 2 + 2 i occ. unocc. a δ (ϕ * i (r, t)ϕ a (r, t)) δv KS (r , ω) × ϕ a | Û NL x -Ûx |ϕ i i -a B 3 (B.22)
where

B 1 =2e -iωt e ηt occ. ij unocc. a all l × ϕ * i (r)ϕ a (r) l| V NL x |i ϕ * a (r )ϕ l (r ) ( a -l -ω -iη)( i -a + ω + iη) - ϕ * i (r)ϕ a (r) aj|li ϕ * l (r )ϕ j (r ) ( j -l + ω + iη)( i -a + ω + iη) - ϕ * i (r)ϕ a (r) al|ji ϕ * j (r )ϕ l (r ) ( j -l -ω -iη)( i -a + ω + iη) + ϕ * i (r)ϕ a (r) a| V NL x |l ϕ * l (r )ϕ i (r ) ( i -l + ω + iη)( i -a + ω + iη) - ϕ * a (r)ϕ i (r) l| V NL x |a ϕ l (r )ϕ * i (r ) ( i -l -ω -iη)( a -i + ω + iη) - ϕ * a (r)ϕ i (r) ij|la ϕ * l (r )ϕ j (r ) ( j -l + ω + iη)( a -i + ω + iη) - ϕ * a (r)ϕ i (r) il|ja ϕ * j (r )ϕ l (r ) ( j -l -ω -iη)( a -i + ω + iη) + ϕ * a (r)ϕ i (r) i| V NL x |l ϕ * l (r )ϕ a (r ) ( a -l + ω + iη)( a -i + ω + iη) (B.23) B 2 =2 occ. i unocc. a all l e -iωt e ηt × ϕ * l (r)ϕ a (r) a| V NL x |i ϕ * i (r )ϕ l (r ) ( i -l -ω -iη)( i -a ) + ϕ * i (r)ϕ l (r) a| V NL x |i ϕ * l (r )ϕ a (r ) ( a -l + ω + iη)( i -a ) - ϕ * l (r)ϕ i (r) i| V NL x |a ϕ * a (r )ϕ l (r ) ( a -l -ω -iη)( a -i ) + ϕ * a (r)ϕ l (r) i| V NL x |a ϕ * l (r )ϕ i (r ) ( i -l + ω + iη)( a -i ) (B.24) B 3 =2e -iωt e ηt occ. i unocc. a all l × ϕ * l (r)ϕ a (r)ϕ * i (r )ϕ l (r ) (i)( i -l -ω -iη) e -i( a-i )(t-t 0 ) a| Û NL x -Ûx |i i -a + ϕ * i (r)ϕ l (r)ϕ * l (r )ϕ a (r ) (-i)( a -l + ω + iη) e -i( a-i )(t-t 0 ) a| Û NL x -Ûx |i i -a + ϕ * l (r)ϕ i (r)ϕ * a (r )ϕ l (r ) (i)( a -l -ω -iη) e -i( i -a)(t-t0) i| Û NL x -Ûx |a i -a + ϕ * a (r)ϕ l (r)ϕ * l (r )ϕ i (r ) (-i)( i -l + ω + iη) e -i( i -a)(t-t0) i| Û NL x -Ûx |a i -a (B.25)
When decomposing the summation over l into a summation over j and b no term of B 3 remains in the occupied-virtual/virtual-occupied space such as B 3 does not contribute in the expression of the EXX kernel.

We can combine the terms of Eq. (B.21) and Eq. (B.23). We remove the phase factor that is common to all the remaining terms (A, B 1 , B 2 and C) and after some algebra the right-hand-side of Eq. (B.12) becomes

B 1 + B 2 -C = -2 occ. ij unocc. a all l × ϕ * i (r)ϕ a (r) aj|li ϕ * l (r )ϕ j (r ) ( j -l + ω + iη)( i -a + ω + iη) + ϕ * i (r)ϕ a (r) al|ji ϕ * j (r )ϕ l (r ) ( j -l -ω -iη)( i -a + ω + iη) + ϕ * a (r)ϕ i (r) ij|la ϕ * l (r )ϕ j (r ) ( j -l + ω + iη)( i -a -ω -iη) + ϕ * a (r)ϕ i (r) il|ja ϕ * j (r )ϕ l (r ) ( j -l -ω -iη)( i -a -ω -iη) + 2 occ. i unocc. a all l × ϕ * i (r)ϕ a (r) l| V NL x -Vx |i ϕ * a (r )ϕ l (r ) ( a -l -ω -iη)( i -a + ω + iη) + ϕ * i (r)ϕ a (r) a| V NL x -Vx |l ϕ * l (r )ϕ i (r ) ( i -l + ω + iη)( i -a + ω + iη) - ϕ * a (r)ϕ i (r) l| V NL x -Vx |a ϕ l (r )ϕ * i (r ) ( i -l -ω -iη)( a -i + ω + iη) - ϕ * a (r)ϕ i (r) i| V NL x -Vx |l ϕ * l (r )ϕ a (r ) ( a -l + ω + iη)( a -i + ω + iη) ϕ * l (r)ϕ a (r) a| V NL x -Vx |i ϕ * i (r )ϕ l (r ) ( i -l -ω -iη)(-i)( i -a ) + ϕ * i (r)ϕ l (r) a| V NL x -Vx |i ϕ * l (r )ϕ a (r ) ( a -l + ω + iη)(-i)( i -a ) - ϕ * l (r)ϕ i (r) i| V NL x -Vx |a ϕ * a (r )ϕ l (r ) ( a -l -ω -iη)(-i)( a -i ) - ϕ * a (r)ϕ l (r) i| V NL x -Vx |a ϕ * l (r )ϕ i (r ) ( i -l + ω + iη)(-i)( a -i ) . (B.26)
Splitting the summations over l in summations over b and j and considering only the terms in the occupied-virtual and virtual-occupied space, this term becomes

B 1 + B 2 -C = -2 occ. ij ab × ϕ * i (r)ϕ a (r) aj|bi ϕ * b (r )ϕ j (r ) ( j -b + ω + iη)( i -a + ω + iη) + ϕ * i (r)ϕ a (r) ab|ji ϕ * j (r )ϕ b (r ) ( j -b -ω -iη)( i -a + ω + iη) + ϕ * a (r)ϕ i (r) ij|ba ϕ * b (r )ϕ j (r ) ( j -b + ω + iη)( i -a -ω -iη) + ϕ * a (r)ϕ i (r) ib|ja ϕ * j (r )ϕ b (r ) ( j -b -ω -iη)( i -a -ω -iη) + 2 occ. ij unocc. a × ϕ * i (r)ϕ a (r) j| V NL x -Vx |i ϕ * a (r )ϕ j (r ) ( a -j -ω -iη)( i -a + ω + iη) - ϕ * a (r)ϕ i (r) i| V NL x -Vx |j ϕ * j (r )ϕ a (r ) ( a -j + ω + iη)( a -i + ω + iη) + 2 occ. i unocc. ab × ϕ * i (r)ϕ a (r) a| V NL x -Vx |b ϕ * b (r )ϕ i (r ) ( i -b + ω + iη)( i -a + ω + iη) - ϕ * a (r)ϕ i (r) b| V NL x -Vx |a ϕ b (r )ϕ * i (r ) ( i -b -ω -iη)( a -i + ω + iη) (B.27)
We can finally define the exchange kernel by combining Eq. (B.19) and Eq. (B.27)

dr 2 dr 3 χ 0 (r 1 , r 2 , ω) f x (r 2 , r 3 , ω) χ 0 (r 3 , r 4 , ω) = h x (r 1 , r 4 , ω) (B.28)
where the function h x (r, r , ω) can be decomposed in four contributions

h x (r, r , ω) = p=1,4 h p x (r, r , ω). (B.29)
The first and the second contribution recover the two first lines of Eq. (B.27)

h 1 x (r, r , ω) = 2 occ. ij unocc. ab × ϕ * i (r)ϕ a (r) aj|bi ϕ * b (r )ϕ j (r ) ( j -b + ω + iη)( i -a + ω + iη) + ϕ * a (r)ϕ i (r) ib|ja ϕ * j (r )ϕ b (r ) ( j -b -ω -iη)( i -a -ω -iη) (B.30) h 2 x (r, r , ω) = 2 occ. ij unocc. ab × ϕ * i (r)ϕ a (r) ab|ji ϕ * j (r )ϕ b (r ) ( j -b -ω -iη)( i -a + ω + iη) + ϕ * a (r)ϕ i (r) ij|ba ϕ * b (r )ϕ j (r ) ( j -b + ω + iη)( i -a -ω -iη) (B.31)
and the third contribution to h x is given by

h 3 x (r, r , ω) = -2 occ. ij unocc. a × ϕ * i (r)ϕ a (r) j| V NL x -Vx |i ϕ * a (r )ϕ j (r ) ( j -a + ω + iη)( i -a + ω + iη) + ϕ * a (r)ϕ i (r) i| V NL x -Vx |j ϕ * j (r )ϕ a (r ) ( j -a -ω -iη)( i -a -ω -iη) +2 occ. i unocc. ab × ϕ * i (r)ϕ a (r) a| V NL x -Vx |b ϕ * b (r )ϕ i (r ) ( i -b + ω + iη)( i -a + ω + iη) + ϕ * a (r)ϕ i (r) b| V NL x -Vx |a ϕ * i (r )ϕ b (r ) ( i -b -ω -iη)( i -a -ω -iη) . (B.32)
It is interesting to note that there is a difference here with the expression given by Görling [1] where there is a fourth contribution h 4 x (r, r , ω). In our derivation we did not include this term to the expression of h x (r, r , ω) because it was not defined in the occupiedvirtual/virtual-occupied space.

(ou du gradient de la densité électronique pour GGA ou de la dérivée des orbitales pour meta-GGA) en un point de l'espace. Ces approximations semilocales donnent souvent une bonne description de la corrélation dynamique de courte portée mais échouent à décrire la corrélation statique et la corrélation dynamique de longue portée. De plus, ces approximations présentent une erreur de self-interaction qui tend à favoriser une délocalisation des électrons et induit une énergie totale trop basse.

Une façon d'améliorer de manière systématique les fonctionnelles approchées, principalement en s'attachant à réduire l'erreur de self-interaction, est de combiner la théorie de la fonctionnelle de la densité avec les méthodes de fonction d'onde en créant des approximations hybrides. La combinaison de ces deux méthodes peut être faite de différentes manières. La méthode la plus simple est de faire une décomposition linéaire de l'interaction coulombienne inter-électronique en deux contributions, une première fraction qui sera traitée avec une méthode fonction d'onde et le reste qui sera traité avec une fonctionnelle de la densité approchée. aléatoires (RPA) [10]. Pour améliorer la description de la corrélation statique, la théorie 123 de la fonctionnelle de la densité peut être associée ç la méthode du champ autocohérent multi-configurationnel (MCSCF) [11].

Pour combiner la théorie de la fonctionnelle de la densité et les méthodes fonction d'onde on peut aller au-delà de la combinaison linéaire avec la séparation de portée. Introduite dans sa forme actuelle par Savin en 1996 [12] et consiste à décomposer l'interaction coulombienne en une contribution de courte portée et une contribution de longue portée : Cette transition est faite de manière progressive grâce à l'utilisation d'une fonction erreur et peut être modulée grâce au paramètre de séparation de portée µ. Dans la pratique, l'interaction de courte portée va être décrite en utilisant une fonctionnelle de la densité approchée et l'interaction de longue portée va être décrite en utilisant une méthode basée sur la fonction d'onde. Une version limitée à la séparation de portée de l'échange a été proposée par Iikura et al. [13] en introduisant de l'échange de longue portée Hartree-Fock alors que l'échange de courte portée et la corrélation sont traités par une fonctionnelle de la densité. La séparation de portée peut aussi être appliquée à la corrélation, par exemple en utilisant MP2 pour la corrélation de longue portée avec l'échange de longue portée en Hartree Fock et une fonctionnelle d'échange et de corrélation de courte portée [14].

Cette décomposition peut aussi être effectuée en utilisant d'autres méthodes pour décrire la corrélation de longue portée comme la RPA [15] ou la méthode coupled-cluster [16] qui sont bien adaptées pour la description des interactions de dispersion de van der Waals.

Un traitement multi-configurationnel de la corrélation de longue portée peut aussi être utilisé pour améliorer la description de la corrélation statique, par exemple en utilisant une méthode MCSCF ou la théorie de la fonctionnelle de la matrice densité (DMFT) [17].

La théorie de la fonctionnelle de la densité a été étendue pour décrire les états excités grâce à la théorie de la fonctionnelle de la densité dépendante du temps (TDDFT). Dans Dans cette thèse nous nous sommes intéressés à différents aspects des ces méthodes hybrides fonction d'onde/fonctionnelle de la densité. Le premier chapitre sera un rappel des bases de la théorie de la fonctionnelle de la densité et de la séparation de portée.

La suite de la thèse est décomposée en trois chapitres qui correspondront à trois projets indépendants autour de l'hybridation entre théorie de la fonctionnelle de la densité et fonction d'onde avec ou sans séparation de portée.

Convergence en base de la théorie de la fonctionnelle de la densité avec séparation de portée

Conclusion et perspectives

Nous nous sommes donc intéressés tout au long de cette thèse aux approximations hybrides entre les méthodes basées sur la fonction d'onde et la théorie de la fonctionnelle de la densité en considérant deux façons de décomposer l'interaction électronique : de manière linéaire ou par la séparation de portée.

Dans le premier chapitre, au travers de l'étude de la convergence en base de la méthode avec séparation de portée nous avons mis en évidence l'accélération de la convergence qui

CHAPTER 1 .

 1 REVIEW OF DFT AND HYBRID METHODSapproximation, considering the electrons moving and the nuclei fixed Ĥ = -

  Ψ [n] is the ground-state wave function associated with n, T [n] = Ψ [n] | T |Ψ [n] the kinetic energy and E ee = Ψ [n] | Ŵee |Ψ [n] the electron-electron interaction energy. In the special case where the considered density is the ground-state density, the energy becomes the energy of the ground state.

F

  min n [n] + n(r)v ne (r)dr . (1.11) Levy's constrained-search formalism Another way to define the universal functional is to use the constrained-search approach proposed by Levy in 1979 [4] F [n] = min Ψ→n Ψ| T + Ŵee |Ψ = Ψ[n]| T + Ŵee |Ψ[n] (1.12) where the minimization is done over all the normalized wave functions Ψ which yield the fixed density n. The minimizing wave function for a given density is Ψ[n]. The constrained approach allows us to easily connect the wave-function variational principle to the density variational principle. Starting fromE 0 = min Ψ Ψ| T + Vne + Ŵee |Ψ ,(1.13)we decompose the variational principle in two steps, a search over the subset of all the antisymmetric wave functions Ψ → n that yield a given density n and a search over all densities E Ŵee |Ψ + n(r)v ne (r)dr .(1.15)

  FIG. 1.Convergence of the truncated partial-wave expansion 1 + f L (r 12 ) for the Coulomb

FIG. 2 .

 2 FIG.2. Convergence rate of the coefficients f lr,µ of the partial-wave expansion with respect to for ≥ 1, for several values of the range-separation parameter µ (in bohr -1 ) and for the Coulomb case (µ → ∞). On the left: Plot of ln |f lr,µ | vs. ln which is linear for a power-law convergence.On the right: Plot of ln |f lr,µ | vs. which is linear for an exponential-law convergence. The curves for the Coulomb interaction and for µ = 10 are nearly superimposed.

ForFIG. 3 .

 3 FIG.3. FCI wave function of the He atom at the Cartesian electron coordinates r 1 = (0.5, 0., 0.) bohr and r 2 = (0.5 cos θ, 0.5 sin θ, 0.) bohr, calculated with Dunning basis sets ranging from cc-pVDZ to cc-pV6Z (abbreviated as VXZ) and shown as a function of the relative angle θ, for the standard Coulomb interaction case (left) and the long-range interaction case for µ = 0.5 bohr -1

  ), in the short-range Hartree-exchange-correlation potential. To facilitate the extraction of the wave function from the program, we use the Löwdin-Shull diagonal representation of the spatial part of the FCI wave function in terms of the spatial natural orbitals (NO) {ϕ µ i } [76, 77]

x

  (x) = δE DFA x /δn(x) and v DFA c (x) = δE DFA c /δn(x) are the local exchange and correlation DFA potentials, respectively. These 1H orbitals and corresponding orbital energies ε 1H

x

  (x) = δE HF x /δn(x) is the EXX potential and v GL2 c (x) = δE MP2 c /δn(x) is here referred to as the GL2 correlation potential (even though it does not contain the singleexcitation term). Since E HF x and E MP2 c are only implicit functionals of the density through the orbitals and orbital energies, the calculation of v EXX x (x) and v GL2 c

  Hj|jH are the expectation values of the EXX and HF exchange potentials over the HOMO spin orbital referred to as H. Similarly, we impose the HOMO condition on the GL2 potential v GL2 c,HH = Σ MP2 c,HH (ε H ),

  dxdx are the expectation values of the GL2 local potential and of the MP2 self-energy over the HOMO spin orbital. The MP2 self-energy is defined as the functional derivative of E MP2 c with respect to the one-particle Green function G(x , x; ω), i.e. Σ MP2 c (x, x ; ω) = 2πi δE MP2 c/δG(x , x; ω), and its diagonal matrix elements Σ MP2 c,pp (ω) are[49]

8 FIG. 1 .

 81 FIG.1. Exchange-correlation and correlation potentials calculated with the OEP-1DH approximation using the BLYP functional at the recommended value λ = 0.65 for the Be atom and the CO molecule using different TSVD cutoffs 10 -4 , 10 -6 , 10 -8 for the pseudo-inversion of the linear-response function. For Be, the potentials for 10 -4 and 10 -6 are superimposed. For CO, the potentials for 10 -6 and 10 -8 are superimposed.

FIG. 3 .

 3 FIG.3. HOMO orbital energies in the 1H approximation [Eq. (6)] and minus IPs in the 1DH

FIG. 4 .

 4 FIG.4. LUMO orbital energies in the 1H [Eq. (6)] and OEP-1DH [Eq. (8)] approximations and mi-

FIG. 5 .

 5 FIG.5. Exchange-correlation potentials calculated with the OEP-1DH approximation [Eq. (9)] using the BLYP functional at the recommended value λ = 0.65, and at the extreme values λ = 0 (standard BLYP) and λ = 1 (OEP-GL2). The reference potentials were calculated by KS inversion of the CCSD(T) densities. For Be, the OEP-GL2 calculations are unstable. For CO, the potential is plotted along the direction of the bond with the C nucleus at -1.21790 bohr and the O nucleus at 0.91371 bohr. For H 2 O, the potential is plotted along the direction of a OH bond with the O nucleus at 0.0 and the H nucleus at 1.81225 bohr.

FIG. 7 .

 7 FIG. 7. Correlated density calculated with the OEP-1DH approximation using the BLYP functional at the recommended value λ = 0.65, and at the extreme values λ = 0 (standard BLYP)

  occ. i|ϕ i (r, t)| 2 and the Kohn-Sham potential operator VKS [n](t) = n(r)v KS (r, t)dr. This density is imposed by the time-dependent Kohn-Sham potentialv KS (r, t) = v(r, t) + v H (r, t) + v xc (r, t)(4.4)with the time-dependent Hartree potential v H (r, t) = dr n(r , t) |rr | (4.5)

4. 3

 3 Study of the short-range exchange kernel 4.3.1 Range-separated time-dependent exact-exchange method

  (4.43) and Eq.(4.44). Due to the Dirac function the non-local short-range exchange potential V NL,srx behaves as a local potential such that X 3,sr,(0) ia,jb (ω) does not contribute and the kernel only depends on X 1+2,sr,(0)

Figure 4 . 1 :

 41 Figure 4.1: Projection of the short-range exchange kernel in the VDZ basis set (ia|f srx |ia) with respect to the range-separation parameter µ for H 2 at equilibrium (1.4 bohr) (top) and in the dissociation limit (distance of 10 bohr) (bottom) calculated with the short-range LDA, shortrange exact-exchange and short-range exact exchange in the µ → ∞ limit. For each distance two transitions are shown: the first transition 1σ g → 1σ u and a higher transition 1σ g → 3π g .

Figure 4 . 2 :

 42 Figure 4.2: Projection of the short-range exchange kernel in VDZ basis set for He calculated with the short-range LDA, short-range exact-exchange and short-range exact exchange in the µ → ∞ limit. Two transitions are shown: the first transition 1s → 2s and a higher transition 1s → 2p.

Figure 4 . 3 :

 43 Figure 4.3: Determinants for H 2 in a minimal basis set.

(4. 60 )

 60 with the coefficients for |0 defined by c 2 1 + c 2 2 = 1 and c 2 = (∆ -∆ 2 + K 2 12 )/K 12 and the coefficients for |1 defined by (c 1 ) 2 + (c 2 ) 2 = 1 andc 2 = (∆ + ∆ 2 + K 2 12 )/K 12 with 2∆ = 2 2 -2 1 +J 22 -4J 12+2K 12 +J 11 . The orbital energies are 1 = h 11 +J 11 and 2 = h 22 +2J 12 -K 12 , 4.5. STUDY OF THE EXACT FREQUENCY-DEPENDENT CORRELATION KERNEL 93 h pp = p| ĥ|p the one-electron integral and J pq = pq|pq and K pq = pq|qp the Coulomb and exchange two-electrons integrals. The wave function |0 is the ground state, |1 and |2 are singlet excited states and |3 , |4 and |5 are triplet excited states. The integrals used to define

+ J lr 12 -K 12 + K lr 12 -K 12 + K lr 12 -J 12 + J lr 12 = 3 1K lr 12 -K 12 + K lr 12 -K 12 + K lr 12 -K 12 + K lr 12 =

 12121212312121212 spin-independent) exact-exchange presented in Sec. 4.

Figure 4 . 4 :

 44 Figure 4.4: Contributions to the singlet short-range correlation kernel for H 2 in STO-3G. We plotted the diagonal element of the frequency-dependent matrix [∆ 1 f c (ω = 1)] 12,12 , the diagonal elements of the frequency-independent matrix [ 1 f c (ω = 0)] 12,12 and the off-diagonal elements of the frequency-independent matrix [ 1 f c (ω = 0)] 12,21 with respect to µ (left) and the ratio |∆ 1 f c (ω)/ 1 f c (ω = 0)| for both HF and EXX exchange with respect to µ (right). These results are calculated at equilibrium distance (1.4 bohr) (top) and when stretching the H-H bond (10 bohr) (bottom). On the bottom plot the curves for the 1 f sr c (ω = 0) are superimposed.

  dependent. An implementation work could be also done to describe larger systems and study the adiabatic semilocal approximation in more complicated examples implying double excitations. Chattaraj, editors, Electronic Structure and Reactivity, Concepts and Methods in Modern Theoretical Chemistry Vol. 1, pages 367-390. CRC Press, 2013.

  where V NL x (t) and Û NL x are the time-dependent and time-independent non-local exchange potential operators and k the orbital energies associated to the time-independent Kohn-Sham orbital ϕ k . We consider the TDEXX equation with the set of Kohn-Sham orbitals satisfying the time-dependent Kohn-Sham equation

  δϕ

17

 17

)

  Finally we apply a chain rule to the derivative of the exchange potential with respect to the Kohn-Sham potentialδv x (r , ω) δv KS (r , ω) = dr δv x (r , ω) δn(r , ω) δn(r , ω) δv KS (r , ω) = dr f x (r , r , ω)χ 0 (r , r , ω) (B.18)and A becomes A = e -iωt e ηt dr dr χ 0 (r, r , ω)f x (r , r , ω)χ 0 (r , r , ω) (B.19)

.

  Cette décomposition a été réalisée pour la première fois en 1993 par Becke qui a proposé une décomposition de l'énergie d'échange half-and-half[7]: une moitié de l'échange est traité par Hartree-Fock alors que la seconde moitié est traitée par une fonctionnelle d'échange approchée, la corrélation étant traitée dans son intégralité en DFT. Cette approximation a été suivie par une nouvelle approximation dont la paramétrisation est empirique et qui contient une part moins importante d'échange Hartree Fock[8]. Ces approximations consistant à introduire une fraction d'échange Hartree-Fock sont dites approximations hybrides. Les versions modernes de ces approximations utilisent une fraction plus réduite d'échange Hartee-Fock (20-25%). L'extension naturelle de cette approximation hybride est de décomposer aussi le terme de corrélation, ceci étant fait en introduisant une fraction d'énergie de corrélation calculée en utilisant une théorie de perturbation Møller-Plesset au second ordre (MP2). Cette décomposition à été introduite à l'origine par Grimme [9] et est qualifiée d'approximation double-hybride. Cette approximation double-hybride permet d'utiliser une fraction plus importante d'échange Hartree-Fock (50-70%) que dans les approximations hybrides sans trop perdre le bénéfice de la compensation d'erreur entre les fonctionnelles d'échange et de corrélation. Cependant la méthode souffre des même limitations que MP2 dans la description de certains phénomènes comme la corrélation statique. La fraction d'énergie de corrélation calculée à l'aide d'une méthode fonction d'onde peut être traitée en utilisant d'autres méthodes comme l'approximation de phases

+ 1 -

 1 erf(µr) r fonctionnelle de la densité .

  son formalisme de réponse linéaire, le terme clé de la TDDFT est le noyau d'échangecorrélation qui comme le potentiel d'échange-corrélation en DFT doit être approché. Dans le cas de la TDDFT, il faut apporter une approximation à la dépendance spatiale et à la dépendance en fréquence du noyau d'échange-corrélation. L'approximation la plus simple est l'approximation adiabatique semilocale, cette approximation est à la fois locale dans le temps (ne dépendant pas de la fréquence) et dans l'espace. Elle offre cependant des résultats raisonnables pour les énergies d'excitations de valence mais échoue dans la description de certains phénomènes comme les excitations multiples, les énergies d'excitation de transfert de charge et les énergies d'excitation de Rydberg. Afin de contourner les lim-itations de l'approximation adiabatique semilocale, une stratégie consiste à étendre la TDDFT à la séparation de portée. La décomposition du noyau d'échange en un noyau d'échange Hartree-Fock de longue portée et un noyau d'échange de courte portée traité avec une fonctionnelle de la densité approchée a été introduite par Tawada et al.[18] et permet de corriger certains problèmes engendrés par l'approximation semilocale comme la description des énergies d'excitation de Rydberg et des énergies d'excitation de transfert de charge. La séparation de portée a ensuite été étendue au noyau de corrélation en combinant un noyau de corrélation de courte portée traité par une fonctionnelle de la densité avec une méthode de réponse linéaire MCSCF[25] ou DMFT[26]. L'introduction de méthodes multi-configurationnelles pour traiter la longue portée permet d'améliorer la description de la corrélation statique et la description des excitations doubles. On peut aussi combiner le noyau de corrélation de courte portée avec un noyau de corrélation de longue portée calculé en utilisant le formalisme des fonctions de Green à N corps[27] qui permet d'introduire une dépendance en fréquence.

  Dans ce chapitre nous nous sommes intéressés à l'étude de la convergence en base dans le cas de la séparation de portée. La convergence en base a été étudiée dans des travaux précédents pour différentes méthodes de fonction d'onde comme MP2, cette convergence dans le cas des bases de Dunning étant polynomiale par rapport au nombre cardinal de la base X[19]. Les méthodes à séparation de portée, qui combinent un traitement DFT de la courte portée et un traitement fonction d'onde de la longue portée, ont montré une convergence en base plus rapide que ce qui avait été observé précédemment pour les méthodes sans séparation de portée. Ces résultats ont pu être observés pour différentes méthodes fonction d'onde par exemple pour MP2[20]. Le chapitre est décomposé en deux études, tout d'abord une étude théorique en développement en ondes partielles et ensuite l'extension de l'étude à la convergence dans les bases mono-atomiques de Dunning. les performances de ces deux approximations double-hybrides (standard et auto-cohérent) pour le calcul de propriétés atomiques et moléculaires tels que le potentiel d'ionisation ou l'affinité électronique. Une première partie de ce chapitre est consacrée à la présentation des méthodes utilisées dans ce chapitre et particulièrement sur le calcul de 1DH autocohérent utilisant la méthode OEP. Comme les propriétés auxquelles nous nous sommes intéressés dans ce chapitre ne peuvent pas forcément être obtenues directement à l'aide de calculs 1DH nous avons eu recours à l'utilisation d'un nombre d'occupation fractionnaire dans le cas de l'approximation 1DH standard. Dans la seconde partie du chapitre nous avons pu comparer différentes propriétés sur quelques systèmes atomiques et moléculaires. Nous avons pu constater que pour les énergies totales et les potentiels d'ionisation l'auto-cohérence n'apporte pas d'améliorations. Nous avons même pu constater que dans le cas du potentiel d'ionisation une bonne première approximation pouvait être obtenue en ne considérant que l'énergie de l'orbitale la plus haute occupée (HOMO) calculée sans le terme MP2. Dans le cas de l'affinité électronique, il est intéressant de noter que la double hybride auto-cohérente permet d'avoir une approximation raisonnable pour l'énergie de l'orbitale LUMO même si elle a une tendance à surestimer l'affinité électronique. Nous avons ensuite étudié les potentiels d'échangecorrélation et de corrélation obtenus avec la 1DH auto-cohérente. Les potentiels obtenus donnent de bons résultats et sont un compromis entre la fonctionnelle BLYP seule qui a un mauvais comportement dans la limite asymptotique et la théorie de la perturbation au second-ordre auto-cohérent (OEP-GL2) qui montre une tendance à surestimer le potentiel. Enfin, nous avons considéré les densités de corrélation qui de la même manière que les potentiels donnent de bons résultats et sont un bon compromis entre BLYP et le OEP-GL2. Étude du noyau d'échange-corrélation de courte portée La théorie de la fonctionnelle dépendante du temps, dans son formalisme de réponse linéaire dépend du noyau d'échange corrélation. L'approximation majoritairement utilisée dans la pratique est l'approximation adiabatique semilocale. Malgré de bon résultats pour les énergies d'excitation électronique de basse énergie, cette approximation connait d'importantes limitations pour décrire les énergies d'excitation de transfert de charge, de Rydberg et les excitations doubles. Une manière de contourner ces limitations de l'approximation adiabatique semilocale consiste à étendre la TDDFT à la séparation de portée. En pratique le noyau d'échange-corrélation de courte portée reste cependant traité à l'aide d'une approximation adiabatique semilocale. Dans ce chapitre, nous nous sommes intéressé à la dépendance en fréquence du noyau d'échange-corrélation de courte portée pour voir si l'utilisation de l'approximation adiabatique à courte portée est suffisante. Pour étudier le noyau d'échange de courte portée nous avons dans un premier temps étendu la méthode de TDDFT d'échange exact à la séparation de portée et défini le noyau d'échange exact de courte portée. Nous avons ensuite étudié le développement asymptotique de noyau dans la limite où µ → ∞. Nous avons pu voir que dans cette limite le premier terme du développement asymptotique n'est pas dépendant de la fréquence et est local. Nous avons finalement considéré deux systèmes modèles : H 2 et He dans la base cc-pVDZ pour comparer le comportement des noyau d'échange exact de courte portée, LDA de courte portée et le premier terme du développement à µ → ∞. Dans un second temps nous nous sommes intéressés au noyau de corrélation de courte portée. Dans ce cas le noyau de corrélation exact n'est pas connu de manière analytique dans le cas général et nous avons choisi d'étudier un système modèle : H 2 en base minimale. Nous avons dans un premier temps calculé de manière analytique le noyau de corrélation de courte portée exact à partir des résultats d'un calcul d'interaction de configurations complète (FCI) dans la base minimale. Nous avons ensuite évalué les résultats obtenus dans la base STO-3G.

  devient exponentielle. Cette étude nous a permis de proposer une méthode d'extrapolation de l'énergie en base complète. Tous les résultats présentés ont été obtenus pour l'énergie de corrélation calculée au niveau MP2 mais des résultats similaires sont attendus si d'autres méthodes sont utilisées pour décrire la corrélation de longue portée. De plus dans cette étude nous avons remarqué que les bases de Dunning utilisée n'étaient pas optimisées pour effectuer les calculs avec séparation de portée. Une façon de poursuivre ce travail pourrait être de définir des bases optimisées pour les calculs avec séparation de portée. Dans le second chapitre nous nous sommes intéressés à une méthode alternative pour traitée les approximations double-hybride de manière auto-cohérente basée sur la méthode OEP. Les résultats que nous avons obtenus sont encourageants, particulièrement pour la description des affinités électroniques, de l'énergie de la LUMO et des potentiels (échange corrélation et corrélation) et densités. Plusieurs pistes peuvent être envisagées pour continuer ce travail : dans un premier temps étendre l'étude en intégrant plus de systèmes et en comparant les résultats obtenus avec d'autres méthodes comme celle proposée par Peverati et Head-Gordon. Une autre étape serait d'étendre cette approximation doublehybride auto-cohérente à la séparation de portée. Enfin il serait aussi intéressant d'étendre cette approximation à la description des états excités en appliquant la TDDFT dans son formalisme de réponse linéaire à l'approximation de double-hybride auto-cohérente. Enfin dans le dernier chapitre nous nous sommes intéressés à l'étude des noyaux d'échange et de corrélation de courte portée. Combinant une approche théorique et des calculs sur des systèmes tels que H 2 ou He nous avons montré que l'approximation semilocale pour le noyau d'échange de courte portée devenait exacte dans la limite où µ. Cette étude a ensuite été étendue au noyau de corrélation de courte portée qui nous a permit de voir que même dans le cas du noyau de corrélation de courte portée l'approximation adiabatique n'était pas exacte et que le terme dépendant en fréquence était encore plus important dans la limite de dissociation de H 2 . Pour continuer cette étude il serait intéressant de calculer le terme suivant dans le développement asymptotique du noyau d'échange exact de courte portée pour obtenir le terme en 1/µ 4 qui devrait voir apparaitre la dépendance en fréquence. Un travail d'implémentation permettrait aussi de décrire des systèmes plus grands pour étudier l'approximation adiabatique semilocale dans des exemples plus compliqués, par exemple avec des excitations doubles pour les noyaux d'échange et de corrélation de courte portée.

  r 2

	He	2	6	2.749	81.84 99.82 0.749 44.04 98.55
		3	6	2.902 104.00 99.97 0.669 29.51 99.19
	Ne	2	6	2.543 843.36 99.59 0.696 479.76 99.00
		3	6	2.754 1169.90 99.93 0.636 355.28 99.37
	N 2	2	6	2.513 697.74 99.71 0.686 397.67 98.76
		3	6	2.693 923.76 99.98 0.621 286.90 99.16
	H 2 O	2	6	2.742 717.04 99.52 0.751 391.27 99.09
		3	6	2.988 1051.28 99.92 0.690 288.61 99.39
					Long-range interaction
					Power law		Exponential law
		X min X max	α	A	r 2	β	B	r 2
	He	2	5	3.128	2.35 96.71 0.974	1.68 99.77
		3	5	3.997	8.16 99.11 1.028	2.13 99.95
	Ne	2	5	3.189	22.06 95.04 0.998 15.97 99.18
		3	5	4.257 101.51 98.27 1.098 24.57 99.65
	N 2	2	5	4.257	73.26 98.63 1.313 44.49 99.97
		3	5	4.997 211.13 99.44 1.283 39.08 99.99
	H 2 O	2	5	3.861

60.27 97.30 1.198 39.24 99.72 3 5 4.686 196.20 97.63 1.210 41.50 99.33

TABLE III .

 III Errors (in mhartree) on the total RSH+lrMP2 energy, E µ

TABLE IV .

 IV Errors (in mhartree) on the total RSH+lrMP2 energy, E µ

TABLE V .

 V Truncated CI correlation energies (in mhartree) of the He atom for the Coulomb interaction and for the long-range interaction (at µ = 0.5 bohr -1 ) for different Dunning basis sets cc-pVXZ (abbreviated as VXZ) and truncation ranks.

			Coulomb interaction
	rank	VDZ	VTZ	VQZ	V5Z	V6Z
	1s2s	-14.997 -15.806 -16.087 -16.204 -16.212
	1s2s2p	-32.434 -35.256 -35.664 -35.794 -35.808
	1s2s2p3s		-35.909 -36.425 -36.595 -36.618
	1s2s2p3s3p		-37.448 -38.068 -38.261 -38.291
	1s2s2p3s3p3d		-39.079 -39.807 -39.999 -40.028
			Long-range interaction
	rank	VDZ	VTZ	VQZ	V5Z	V6Z
	1s2s	-0.018 -0.036 -0.038 -0.039 -0.039
	1s2s2p	-0.155 -0.329 -0.415 -0.449 -0.468
	1s2s2p3s		-0.329 -0.415 -0.450 -0.468
	1s2s2p3s3p		-0.329 -0.415 -0.450 -0.469
	1s2s2p3s3p3d		-0.329 -0.416 -0.451 -0.470

TABLE VI .

 VI Results of the fits to the power and exponential laws of the truncated CI correlation energy error for the Coulomb interaction, ∆E c,X

  we use the extension of the DH total energy expression, including the MP2 correlation term, to fractional orbital occupation numbers, as given in Refs. 58 and 59 (for the details of our implementation, see Ref. 60). As pointed out in Ref. 58, if the variation of the orbitals and orbital energies in the MP2 correlation energy is neglected when taking the derivative of the 1DH total energy with the respect to N , then Eqs. (28) and (29) simplify to

  (r, t; r , t ) the Kohn-Sham linear-response function. The potential v(r, t) can be expressed with respect to the Kohn-Sham potential

	with χ 0 v(r, t) = v KS (r, t) -v Hxc (r, t)	(4.10)
	and can be differentiated with respect to the density
	δv(r, t) δn(r , t )	=	δv KS (r, t) δn(r , t )	-	δv Hxc (r, t) δn(r , t )	(4.11)
					1) KS (r , t )	(4.9)

Table A .

 A 1: MP2 correlation energies and their errors (in mhartree) for the long-range interaction at µ = 0.5 bohr -1 (E lr c and ∆E lr c ) calculated with core Dunning basis sets of increasing sizes for Ne, N 2 and H 2 O. The errors are calculated with respect to the cc-pCV6Z values.

				Long-range interaction
					µ = 0.5
	Basis set	E lr c	Ne	∆E lr c	E lr c	N 2	∆E lr c	E lr c	H 2 O	∆E lr c
	cc-pCVDZ	-0.7003		1.9640	-20.3491	3.2790	-6.4989	3.5277
	cc-pCVTZ	-1.7941		0.8702	-22.8333	0.7949	-8.9992	1.0274
	cc-pCVQZ	-2.3358		0.3285	-23.4009	0.2272	-9.6617	0.3649
	cc-pCV5Z	-2.5675		0.0967	-23.5651	0.0631	-9.9344	0.0922
	cc-pCV6Z	-2.6642			-23.6282			-10.0266

Table A .

 A 2: Results of the fits to the power and exponential laws of the long-range MP2 correlation energy error ∆E lr c,X for µ = 0.5 bohr -1 . Different ranges X min ≤ X ≤ X max , for the cardinal number X of the core Dunning basis sets are tested. The parameters A and B are in mhartree. The squared Pearson correlation coefficients r 2 of the fits are indicated in %.

					Long-range interaction
					Power law			Exponential law
		X min X max	α	A	r 2	β	B	r 2
	Ne	2	5	-3.1969	22.1531	95.04	-1.0007	16.0226	99.17
		3	5	-4.2572	100.9000	98.11	-1.0984	24.4685	99.58
	N 2	2	5	-4.2522	71.3488	98.79	-1.3105	43.1576	99.92
		3	5	-4.9313	188.4132	99.37	-1.2668	35.7243	100.00
	H 2 O	2	5	-3.8597	59.9484	97.39	-1.1968	38.9700	99.73
		3	5	-4.6663	189.9653	97.64	-1.2053	40.4372	99.34

Table A .

 A 3: Errors (in mhartree) on the total RSH+lrMP2 energy, E = E RSH + E lr,MP2 c

Table A .

 A

								c
	He	0.1	8.508	0.772	0.261	0.089	0.224	0.224
		0.5	8.488	0.781	0.245	0.078	0.205	0.161
		1.0	8.258	0.924	0.259	0.078	0.192	0.144
	Ne	0.1	72.999	20.215	5.842	0.716	0.464	0.462
		0.5	74.523	20.337	5.763	0.751	0.401	0.163
		1.0	79.311	20.962	5.726	0.803	0.342	0.342
	N 2	0.1	47.061	13.026	4.136	0.853	0.993	0.992
		0.5	51.581	13.406	4.090	0.810	1.083	1.083
		1.0	61.053	15.108	4.513	0.868	1.337	1.303
	H 2 O	0.1	54.861	15.229	5.005	0.857	1.451	1.451
		0.5	55.850	14.736	4.499	0.726	1.105	1.102
		1.0	61.013	15.212	4.423	0.724	1.099	1.087

4: Errors (in mhartree) on the total RSH+lrMP2 energy, E = E RSH + E lrMP2 c

  10) and using these derivatives applied to Eq. (B.1) we can deduce after some calculations the EXX kernel. The derivative of Eq. (B.1) with respect to v KS (r , ω) is

		t t 0	dt	dr	δχ 0 (rt, r t ) δv KS (r , ω)	v x (r , t ) + χ 0 (rt, r t )	δv x (r , t ) δv KS (r , ω)	=	δΛ x (r, t) δv KS (r , ω)	(B.11)
	which can be rearranged as
	t t 0	dt	dr χ 0 (r, t; r , t )	δv x (r , t ) δv KS (r , ω)	=	δΛ x (r, t) δv KS (r , ω)	-	t t 0	dt	dr	δχ 0 (r, t; r , t ) δv KS (r , ω)	v x (r , t )
					A							B
	δv x (r, t) δv KS (r , ω)	= dr	δv x (r, t) δv x (r , ω)	δv x (r , ω) δv KS (r , ω)	= dr	δv x (r , ω) δv KS (r , ω)	e -iωt e ηt δ(r -r ). (B.14)
	The left-hand-side term becomes
					A = dr	t t 0	dt χ 0 (r, t; r , t ) ×	δv x (r , ω) δv KS (r , ω)	e -iωt e ηt	(B.15)

C (B.12) We first express the derivative of the left-hand-side term A of Eq. (B.12) in the frequency domain by expressing the exchange potential in a similar way as the Kohn-Sham potential [Eq. (B.8)] δv x (r, t) δv x (r , ω) = e -iωt e ηt δ(rr ) (B.13)

and we can then simply express the derivative in A by applying a chain rule based on the fact that in the linear-response regime variations of Fourier components of potentials or densities are only coupled if they have the same frequency

[2]

and integrating over t with the initial time t 0 → -∞ it becomes A = e -iωt e ηt dr χ 0 (r, r , ω) δv x (r , ω) δv KS (r , ω) , (B.

16)

with the Kohn-Sham linear-response function in the frequency domain

Acknowledgments

ACKNOWLEDGEMENTS We thank Andreas Savin, Trygve Helgaker, Éric Cancès, and Gabriel Stoltz for insightful discussions. This work was supported by French state funds managed by CALSIMLAB and the ANR within the Investissements d'Avenir program under reference ANR-11-IDEX-0004-02. ACKNOWLEDGEMENTS We thank Andreas Savin for discussions, as well as Andrew Teale for providing the reference data for the exchange and exchange-correlation potentials. This work was supported by French state funds managed by CALSIMLAB, within the Investissements d'Avenir program under reference ANR-11-IDEX-0004-02, by the Polish-French bilateral programme POLONIUM, and by the Polish National Science Center under Grant No. DEC-2013/11/B/ST4/00771.

Appendix A: Convergence of the correlation energy for truncated CI calculations

In this Appendix, we explore the basis convergence of the correlation energy of the He atom for truncated CI calculations for both the Coulomb and long-range interactions. For a given basis set and interaction, we start by performing a FCI calculation and generating the corresponding natural orbitals. We then use these natural orbitals in truncated CI calculations for increasing orbital active spaces 1s2s, 1s2s2p, 1s2s2p3s, 1s2s2p3s3p, and linear-response function calculated in Sec. 4.5.2 but evaluated with the long-range Hamiltonian Ĥ = T + Ŵee + Vne → Ĥlr = T + Ŵ lr ee + Vne + V sr Hxc (4.72) leading to the long-range orbital energies

with V sr xc,ij = drϕ i (r)v sr xc (r)ϕ j (r) the integral of the short-range exchange-correlation potential and the long-range Coulomb and exchange two-electrons integrals J pq → J lr pq = pq|w lr ee |pq (4.74a)

The coefficients c lr 1 and c lr 2 and the energy differences ∆E lr K have the same form as for the fullrange interaction but are evaluated using the long-range orbital energies and the long-range two-electron integrals.

Introducing the linear-response function of Sec. 4.5.2 into Eq. (4.16) we obtain the singlet and triplet Hartree-exchange-correlation kernels

)

The singlet and triplet short-range Hartree kernels are

CONCLUSION

tending the exact-exchange TDDFT to range-separation we proposed an asymptotic expansion of the short-range exact-exchange kernel with respect to the range-separation parameter µ. The first terms of the asymptotic expansion of the short-range exact-exchange kernel are both adiabatic and local in space. The practical calculations performed for He and H 2 in the VDZ basis set showed an exact behavior for short-range LDA in the limit of large µ. We extended our work on the short-range exact correlation kernel for H 2 in a minimal basis set. Calculations in STO-3G basis set showed that even for the short-range correlation kernel at large value of µ the adiabatic approximation is not exact. The error coming from adiabatic approximation is even more important in the dissociation limit of H 2 .

These studies can be extended and the perspectives are the following:

• A first point which could be of interest would be to define basis sets that would be optimal for performing range-separated calculations.

• To extend the study on the OEP-based double-hybrid approximation a first step would be to consider a comparison based on a larger set of systems and including other orbitaloptimized double-hybrid approximations such as the one proposed by Peverati and Head-Gordon. A second step would be to extend this self-consistent double-hybrid approximation to range-separation or to the calculation of excited states by applying linear-response time-dependent density-functional theory on this self-consistent double-hybrid approximation.

• The short-range exchange kernel could be extended by deriving the next term of the asymptotic expansion. This term would be in 1/µ 4 and should be frequency dependent.

An implementation work could be done also to describe greater systems and study the adiabatic semilocal approximation in more complicated examples implying double excitations for both short-range exchange and correlation kernels.

Appendix A Additional results for the basis-set convergence of the long-range correlation energy

In this appendix, complementary results to the study in Chapter 2 are presented.

A.1 Convergence of the correlation energy including core electrons

The results presented in the article are given for valence-only calculations. We extended the study to calculations including the core excitations. The long-range MP2 correlation energies including core excitations and the errors with respect to the energy calculated in cc-pCV6Z are presented in Tab. A.1. In this case the approximation of taking the calculations at X = 6 as a converged reference is still valid. We then compared the fits for the polynomial law E lr,µ c,X = E lr,µ c,∞ + AX -α and the exponential law E lr,µ c,X = E lr,µ c,∞ + B exp(-βX). The results of the fits for the long-range MP2 correlation energy including core excitations are presented in Tab. A.2. There is an important difference between the results of the two fits. Similarly to the results for valence-only calculations the best fit is obtained for the exponential law by far with r 2 > 99%. and the second way to do this extrapolation is to extrapolate separately the RSH and the long-range MP2 correlation energies

We compared these two ways of performing the three-point extrapolation. In this case we obtain results that have negative errors and generally, in a similar way to the valence-only calculations, the separated extrapolation shows an improvement with respect to the calculation at X = 4. If we compare the global and separated way to do the three-point extrapolations for the calculations with the core excitations we can see that the results are close. We can conclude that there is no notable difference between those two extrapolations at a chemical accuracy level and that the global extrapolation is sufficient.

at a given potential v KS (r, t). Time-dependent perturbation theory yields the derivative of these orbitals with respect to the Kohn-Sham potential

We need to choose the point at which we will do the derivation of the TDEXX equation with respect to the time-dependent KS potential. The equation we need to derive will be composed of KS orbitals corresponding to the chosen KS potential. The point we consider to do the derivative is

where u KS (r) is the time-independent Kohn-Sham potential. In this special case v KS (r, t)

is no longer time-dependent and the TDEXX formalism simplifies to the static EXX formalism. In this special case the time-dependent Kohn-Sham orbitals only depend on time by a phase factor

To be used in time-dependent linear-response density-functional theory, the kernel needs to be expressed as a frequency-dependent quantity. We need to move from the time domain to the frequency domain. The variation of the potential can be expressed as a function of the frequency as rajouter des informations sur l'inhomogénéité de la densité électronique, ceci étant rendu possible en utilisant le gradient de la densité en plus de la densité pour définir la fonctionnelle. Cette nouvelle famille de fonctionnelles est appelée approximation du gradient généralisée (GGA). Un exemple de fonctionnelle GGA est B88 [4] pour l'échange et LYP [5] pour la corrélation. Cette approximation peut encore être améliorée en prenant en compte le Laplacien de la densité électronique et/ou la densité d'énergie cinétique qui nous permet de définir une nouvelle famille de fonctionnelles : meta-GGA. Une des approximations meta-GGA les plus utilisées est TPSS [6]. LDA, GGA et meta-GGA forment un ensemble d'approximations semi-locales car elles ne dépendent que de la densité électronique