en couplant des simulations de dynamique moléculaire à des algorithmes d'apprentissage automatique.

4 Numerical applications to molecules and solids 4.1 Applications to weakly bound molecular dimers . . . . . . . . . 4.

Despite the high computational cost the adiabatic connection fluctuation dissipation theorem (ACFDT) represents a promising approach to improve the description of the electronic correlation within density functional theory. The simplest approximation that can be applied in the context of the ACFDT is the random phase approximation (RPA). First, we show how the RPA can be improved by introducing a kernel containing an approximate electron-hole exchange term that leads to two different beyond-RPA methods. Then, we show that the RPA and beyond-RPA approaches can be efficiently computed within a plane-wave basis set implementation by using dielectric eigenpotentials as a compact auxiliary basis set and the Lanczos algorithm. A series of applications to molecules and solids are presented to demonstrate the efficiency and accuracy of these approximations. Importantly, it will be shown that the highly accurate beyond-RPA methods can be scaled to treat molecular systems with one hundred electrons requiring a basis set with hundreds of thousands of plane-waves. Finally, it is shown how the sophisticated and computationally expensive ACFDT methods can be used to compute finite-temperature properties of realistic materials (adsorption enthalpies of molecules in zeolites) by coupling molecular-dynamics simulations with machine learning algorithms.

Résumé

Malgré un coût de calcul élevé, le théorème de fluctuation-dissipation avec connexion adiabatique (ACFDT) représente une approche prometteuse pour améliorer la description de la corrélation électronique en théorie de la fonctionnelle de la densité. L'approximation pratique la plus simple dans le contexte de l'ACFDT est l'approximation de la phase aléatoire (RPA). D'abord, nous montrons comment la RPA peut être améliorée en introduisant un noyau contenant un terme approximatif d'échange "électron-trou" qui conduit à deux méthodes différentes au-delà de la RPA. Ensuite, nous montrons que les approches RPA et au-delà de la RPA peuvent être calculées d'une façon efficace dans le cadre d'une implémentation basée sur une base d'ondes planes en utilisant les vecteurs propres de la matrice diélectrique comme une base auxiliaire compacte et l'algorithme de Lanczos. Une série d'applications à des molécules et à des solides sont présentées pour démontrer l'efficacité et la précision de ces approximations. De plus, il sera démontré que les méthodes très précises au-delà de la RPA peuvent traiter des systèmes moléculaires avec des centaines d'électrons de valence et nécessitant une base avec des centaines de milliers d'ondes planes. Enfin, nous montrons comment les méthodes ACFDT, très couteuse en termes de temps de calcul, peuvent être utilisées pour calculer les propriétés à température finie de matériaux réalistes (spécifiquement les enthalpies d'adsorption de molécules dans des zéolites)

Chapter 1 Introduction

This manuscript describes my main (but not exclusive) research activities developed starting from September 2012, when I was hired as a Maître de Conférences at the Université de Lorraine. My work has largely focused on methodological development to compute correlation energies based on the random phase approximation (RPA) and its variants.

Because of its reasonable compromise between accuracy and efficiency density functional theory (DFT) [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] has become the most widespread computational method in first principles materials science. DFT is applicable to relatively large systems of up to thousands of atoms, mainly owing its original success to rather simple approximations for the exchange-correlation (xc) functional, such as the local-density (LDA) [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF][START_REF] Perdew | Self-interaction correction to densityfunctional approximations for many-electron systems[END_REF] and generalized-gradient approximations (GGA) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. More recent classes of functionals include meta-generalized gradient approximations [START_REF] Tao | Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids[END_REF][START_REF] Sun | Strongly constrained and appropriately normed semilocal density functional[END_REF] and hybrid functionals, which mix a certain percentage of Hartree-Fock exchange with semi-local contributions [7]. All of these approximations, involving an increasing level of accuracy but also of computational complexity, can be pictorically depicted by the Jacob's ladder scheme of DFT approximations proposed by John Perdew (see Fig. 1.1) [START_REF] Perdew | s ladder of density functional approximations for the exchange-correlation energy[END_REF]. While several approximations for the the unknown xc functional have been proposed, research in this field is still quite active, with the ultimate goal of obtaining a "final theory". Challenging classes of problems for DFT approximations include materials where strong electronic interactions or van der Waals forces play an important role. Additionally, traditional approximate xc functionals do not systematically reach the chemical accuracy threshold (1 kcal/mol). This manuscript will focus on the fifth and highest rung of the ladder of DFT functionals in Fig. 1.1. Specifically, the ground-state electronic correlation energy within the random phase approximation (RPA) will be considered. The RPA was already introduced in the 1950's to describe the correlation energy in the degenerate electron gas [14] and its diagrammatic structure has long been known (see Fig. 1.2) [45]. In the DFT context the theoretical basis for the RPA correlation energy is provided by the adiabatic-connection fluctuationdissipation theorem (ACFDT) [START_REF] Langreth | The exchange-correlation energy of a metallic surface[END_REF][START_REF] Langreth | Exchange-correlation energy of a metallic surface: Wave-vector analysis[END_REF], that will be introduced in Chapter 2. The ACFDT provides in principle an exact equation for the ground-state correlation energy in terms of linear response functions (i.e. the polarizability or the dielectric matrix). Within the ACFDT the random phase approximation (RPA) has provided a first practical scheme for applications to realistic systems, such as molecules [39,[START_REF] Toulouse | Adiabatic-connection fluctuation-dissipation density-functional theory based on range separation[END_REF][START_REF] Nguyen | Efficient calculation of exact exchange and rpa correlation energies in the adiabatic-connection fluctuation-dissipation theory[END_REF][START_REF] Zhu | Rangeseparated density-functional theory with random phase approximation applied to noncovalent intermolecular interactions[END_REF][START_REF] Nguyen | A first-principles study of weakly bound molecules using exact exchange and the random phase approximation[END_REF]35,[START_REF] Ren | Resolution-ofidentity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF][START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF], solids [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Lu | Ab initio calculation of van der waals bonded molecular crystals[END_REF][START_REF] Wilson | Iterative calculations of dielectric eigenvalue spectra[END_REF][START_REF] Harl | Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids[END_REF][START_REF] Kaoui | Random phase approximation correlation energy using a compact representation for linear response functions: application to solids[END_REF], layered materials [START_REF] Marini | First-principles description of correlation effects in layered materials[END_REF][START_REF] Lebègue | Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase approximation[END_REF], and molecules adsorbed on surfaces [START_REF] Ren | Exploring the random phase approximation: Application to co adsorbed on cu(111)[END_REF][START_REF] Schimka | Accurate surface and adsorption energies from many-body perturbation theory[END_REF]. In particular, for condensed matter/solid state applications the traditional quantumchemical methods (e.g. the coupled-cluster theory) can be hardly applied and the RPA is probably the most sophisticated correlated approximation of practical interest. The rise of interest in the RPA is mostly, but not exclusively, related to the capability of the RPA to accurately include van der Waals (vdW) forces in a seamless way [32,41]. Despite its growing success, the RPA is still a niche approach, especially if compared to other DFT approximations routinely applied in materials simulations. This is related to different reasons: (1) While being a highly sophisticated theory, shortcomings remain within the RPA: Binding curves of certain diatomic dimers present an unphysical "bump" at intermediate distances [39], total correlation energies are usually overestimated [START_REF] Kurth | Density-functional correction of randomphase-approximation correlation with results for jellium surface energies[END_REF], and binding energies are often underestimated [41,[START_REF] Harl | Accurate bulk properties from approximate many-body techniques[END_REF]. (2) Implementations to compute the RPA correlation energies are rather complex and practical calculations involve several numer-Figure 1.2: The effective interaction within the RPA is obtained as the sum to infinite order of ring diagrams (figure reproduced from Ref. [START_REF] Ke | Dispersion relation of excitation mode in spinpolarized fermi gas[END_REF].

ical bottlenecks, such as the full diagonalization of the starting-point Hamiltonian to obtain all the virtual (conduction) states, the manipulation of large response function matrices, and the evaluation of improper integrals over the frequency.

(3) The computational time involved in RPA calculations is often several orders of magnitude higher than the time required by basic (semi-)local approximations; additionally, implementations to compute analytic interatomic forces are not trivial. For these reasons applications to realistic materials are highly demanding in terms of computational resources and typically based on frozen geometries optimized using different theories.

In this manuscript I will address these different issues. While a series of methods have already been proposed to improve the RPA in the literature [START_REF] Olsen | Extending the random-phase approximation for electronic correlation energies: The renormalized adiabatic local density approximation[END_REF][START_REF] Lu | Evaluation of model exchange-correlation kernels in the adiabatic connection fluctuation-dissipation theorem for inhomogeneous systems[END_REF]42,[START_REF] Ren | Beyond the random-phase approximation for the electron correlation energy: The importance of single excitations[END_REF][START_REF] Hesselmann | Correct description of the bond dissociation limit without breaking spin symmetry by a random-phaseapproximation correlation functional[END_REF]25,[START_REF] Hellgren | Beyond the random phase approximation with a local exchange vertex[END_REF]6], in Chapter 2 I will discuss my work on the development of new approximations that avoid the self-correlation problem of the RPA by using a non-local exchange kernel. In Chapter 3 I will show how the previously introduced RPA and beyond-RPA methods can be efficiently implemented in an ab initio solid state code based on plane-waves by using a compact auxiliary basis set, the Lanczos algorithm, and density functional perturbation theory techniques. In Chapter 4 this numerical implementation is applied to a series of benchmark molecular test sets and to bulk solids. Chapter 5 will discuss in detail the (slow) numerical convergence of the (beyond-)RPA correlation energy with respect to the basis set size; this is a well known issue for quantum chemical methods that instead does not manifest itself when DFT is employed with more traditional (semi-)local or hybrid functionals. In Chapter 6 I address the issue of the applicability of RPA-based methods to realistic materials including temperature effects. In this context an impressive level of efficiency and accuracy is achieved by coupling machine learning algorithms with thermodynamic perturbation theory to provide the MLPT method. Applications of MLPT to compute molecular enthalpies of adsorption in zeolites shows the feasibility and potential impact of this approach. The future perspectives of this work are presented in Chapter 7. In particular, future research will further address the development of MLPT and, more in general, the development of techniques to bridge highly sophisticated methods for the electronic correlation with applications to realistic materials.

Chapter 2 Correlation energy within the random phase approximation and beyond

This chapter is based on original work previously published in Refs. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]28,31].

In this chapter we will discuss the theoretical formalism which is at the base of this manuscript. After a short presentation of the main ideas behind density functional theory (DFT), the adiabatic connection fluctuation dissipation theorem (ACFDT) is introduced. In principle, this approach provides an exact equation for the exchange-correlation (xc) energy functional of DFT in terms of response functions (i.e. the polarizability for molecules and the dielectric matrix for solids). In practice, approximations for the response functions are necessary and, in this context, the random phase approximation (RPA) is probably the simplest practical approach. The RPA provides a fully non-local xc functional which accurately includes van der Waals forces in a seamless way [32,41]. Still, the RPA is characterized by some issues: Binding curves of certain diatomic dimers present an unphysical "bump" at intermediate distances [39], total correlation energies are usually overestimated [START_REF] Kurth | Density-functional correction of randomphase-approximation correlation with results for jellium surface energies[END_REF], and binding energies are often underestimated [41,[START_REF] Harl | Accurate bulk properties from approximate many-body techniques[END_REF]. With the aim of decreasing the self-correlation of the RPA and alleviating these issues we introduce two new approximations that include exact exchange effects: The electron-hole time-dependent Hartree-Fock (eh-TDHF) and the adiabatic connection second order screened exchange approximation (AC-SOSEX). The eh-TDHF represents a fully original contribution of my work and the AC-SOSEX was not previously used in a dielectric matrix formulation. The accuracy of these RPA and beyond-RPA approaches will be established by the numerical applications in Chapter 4.

Before starting the discussion a few technical notes are useful to follow the development of the formalism:

• Atomic units will be used throughout this manuscript.

• The formalism will be developed in a spin-restricted way and the orbitals φ are considered as space orbitals.

• For the sake of simplicity, we will suppose that all the orbitals φ are real; this is a general assumption in the case of molecular systems with timereversal symmetry. For the RPA the formalism has also been implemented for periodic systems by including the integration over the first Brillouin zone [START_REF] Kaoui | Random phase approximation correlation energy using a compact representation for linear response functions: application to solids[END_REF].

• The subscripts v will be used to denote valence (namely occupied) states and the subscripts c will be used to denote conduction (namely virtual or unoccupied) states.

Density Functional Theory

Density Functional Theory has become impressively popular in the condensed matter physics and quantum chemistry communities because of its good compromise between accuracy and computational cost. The main (formidable) achievement of this approach is the mapping of a many-body electronic problem onto a non-interacting one. DFT was developed in the work of Hohenberg and Kohn [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] and in the work of Kohn and Sham [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]. Hohenberg and Kohn (HK) established the two theorems which constitute the theoretical foundation of DFT:

• Theorem I: For any system of interacting electrons in an external potential V ext (r), there is a one-to-one (apart from an irrelevant additive constant) correspondence between V ext (r) and the ground-state particle density n 0 (r);

• Theorem II: The energy of an interacting-electron system in an external potential V ext (r) can be expressed as a functional of the electronic density n(r) of the form

E HK [n] = F HK [n] + V ext (r)n(r)dr, (2.1)
where F HK [n] is a universal functional of the density which does not depend on V ext (r). Furthermore, the minimum value of E HK [n] is reached for n = n 0 and corresponds to the exact ground-state energy of the system.

The meaning of the first theorem is that the density completely determines all the properties of a given many body system. The second theorem is simply a corollary of the first and of the variational principle of quantum mechanics. The HK theorems provide a quite general theoretical result but they do not contain any practical recipe on how to solve the quantum many-body problem in practice.

A practical formulation of DFT is provided by the Kohn-Sham (KS) ansatz, which results in rewriting Eq. (2.1) as:

E HK [n] = T [n] + E H [n] + E xc [n] + V ext (r)n(r)dr.
(2.2)

In Eq. (2.2) T is the kinetic energy, E H [n] the Hartree energy, and E xc is defined by Eq. (2.2) as the difference

F HK [n] -T [n] -E H [n].
In order to apply DFT in practice, a good approximation for E xc is necessary. In the Introduction we briefly mentioned some popular approaches, such as the local-density approximations (LDA) [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF][START_REF] Perdew | Self-interaction correction to densityfunctional approximations for many-electron systems[END_REF], the generalized-gradient approximation (GGA) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF], the meta-GGA [START_REF] Tao | Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids[END_REF][START_REF] Sun | Strongly constrained and appropriately normed semilocal density functional[END_REF], and the hybrid functionals [7]; in the next section we will introduce an in principle exact formula for E xc that can be used in practice to develop new approximate xc functionals. By applying the variational principle of Theorem II to the HK energy functional in the form of Eq. (2.2) the Euler equation is obtained:

δT [n] δn(r) + V KS (r) = µ, (2.3)
where

V KS (r) = V H (r) + V xc (r) + V ext (r) (2.4) = n(r ) |r -r | dr + δE xc [n] δn(r) + V ext (r), (2.5)
and µ is a Lagrange multiplier that enforces the conservation of the total number of particles. In this way the ground-state density of the interacting-electron system n 0 (r) can be obtained by solving the Schrödinger equation of a fictitious noninteracting system:

- 1 2 ∇ 2 + V KS (r) φ i (r) = i φ i (r). (2.6)
Indeed, the one-particle orbitals φ i can be then used to construct the density

n(r) = Nv i φ * i (r)φ i (r), (2.7)
where N v is the number of occupied states. Eq. (2.6) is named Kohn-Sham (KS) equation. Since V KS is a functional of the exact ground-state density, Eq. (2.6) has to be solved self-consistently. It is worth to note that, strictly speaking, the eigenvalues i and eigenvectors φ i do not have any physical meaning, but they are just mathematical devices used to obtain the ground-state energy and charge density of the full many-body system. Within this approach it is possible to calculate many ground-state properties, such as atomic forces, equilibrium geometries, stress tensors and phonon dispersion curves.

Correlation energy within the adiabatic connection fluctuation-dissipation theorem

In this section we discuss the exact expression of the exchange-correlation energy E xc that can be obtained using the adiabatic connection fluctuation-dissipation theorem (ACFDT) [START_REF] Langreth | The exchange-correlation energy of a metallic surface[END_REF][START_REF] Langreth | Exchange-correlation energy of a metallic surface: Wave-vector analysis[END_REF].

In the context of density functional theory the starting point to develop the ACFDT formalism is the adiabatic connection Hamiltonian for an electronic system

Ĥλ = T + V λ ext + λ Vee , (2.8)
where Vee =

i =j

1 |r i -r j | (2.9)
contains the contribution from the electron-electron Coulomb interaction and the parameter λ can assume values between 0 and 1. When λ = 0 the electrons are not interacting and the corresponding system is the Kohn-Sham system with V λ=0 ext = Vext + VH + Vxc = VKS . When λ = 1 the full electron interaction Vee is turned on and V λ=1 ext = Vext . Because of the particular form of Ĥλ the density n(r) remains constant by switching on the interaction going from λ = 0 to λ = 1; this introduces some simplifications in the the ACFDT formalism.

By using the Hellmann-Feynman theorem and by integrating over the parameter λ it can be shown that

E H + E ACFDT xc = 1 0 dλ d dλ ( ψ λ | Ĥλ |ψ λ -ψ λ | V λ ext |ψ λ ) (2.10) = 1 0 dλ ψ λ | Vee |ψ λ . (2.11)
The expression in Eq. 2.11 can then be recast in terms of density fluctuation operators. By using the fluctuation-dissipation theorem we can arrive to an expression of the correlation energy in terms of response functions Π evaluated on the complex frequency axis iω:

E ACFDT c = - 1 2 1 0 dλ ∞ -∞ dω 2π Tr Π λ (iω) K -Π 0 (iω) K . (2.
12)

The exchange term E ACFDT x is formally equivalent to the Hartree-Fock exact-exchange but is evaluated using KS orbitals.

In the definition of E ACFDT c the matrix K is constituted by four identical blocks K:

K = K K K K , (2.13)
where the matrix elements

K vc,v c = 2 φ v φ c |V |φ v φ c = 2 φ v (r)φ c (r)V (r -r )φ v (r )φ c (r )drdr
(2.14) correspond to the electron-hole representation of the Coulomb potential V ; the independent electron polarizability Π 0 is given by the equation

Π 0 (iω) = -(∆ -iωI) 0 0 -(∆ + iωI) -1 = Π 0 + (iω) 0 0 Π 0 -(iω) , (2.15)
where I denotes the identity matrix and ∆ vc,v c = ( cv )δ vv δ cc contains the KS independent electron excitations, defined as differences between conduction state energies c and valence state energies v . In Eqs. 2.12-2.15 as well as in the rest of this manuscript the bold letters are used to denote matrices. An additional tilde has been added on top of 2N eh × 2N eh matrices (where N eh = N v N c is equal to the product of the number of valence N v and conduction N c states) to distinguish them from the N eh × N eh blocks; this notation will be useful to avoid confusion in the equations of the following sections.

Random phase approximation

In Eq. 2.12, Π λ refers to the polarizability of a fictitious system where the electrons interact with a scaled Coulomb potential λV . This term is unknown, and approximations need to be made to compute the correlation energy expression. The interacting polarizability under the RPA is defined as:

Π λ RPA (iω) = I -λ Π 0 (iω) K -1 Π 0 (iω), (2.16)
which when plugged in Eq. 2.12 gives us the RPA correlation energy after integration over λ and simplification of the trace over the 4 × 4 block matrices:

E RPA c = 1 2 ∞ -∞ dω 2π Tr ln I -χ 0 (iω)K + χ 0 K , (2.17)
where we introduced

χ 0 (iω) = Π 0 + + Π 0 -= 2 1 iω + v -c δ vv δ cc (2.18)
and represents the real part. The Iχ 0 (iω)K term is commonly known as dielectric matrix and will be denoted by ε. The absence of over the matrices in Eq. 2.17 indicates that the dimensionality has been reduced to N eh × N eh . The expression of the correlation energy in Eq. 2.17 is often referred as direct-RPA (or d-RPA) in the quantum chemistry community to emphasize the absence of exchange effects.

Electron-hole time-dependent Hartree-Fock for correlation energies

In this section the response function will be approximated at the time-dependent Hartree-Fock (TDHF) level of theory [START_REF] Mclachlan | Time-dependent hartree-fock theory for molecules[END_REF]16]:

Π λ (iω) = ( I -λ Π 0 (iω) Ξ) -1 Π 0 (iω), (2.19)
where the kernel Ξ is defined as

Ξ = A B B A , (2.20)
with

A vc,v c = 2 φ v φ c |V |φ v φ c -φ v φ v |V |φ c φ c (2.21) B vc,v c = 2 φ v φ c |V |φ v φ c -φ v φ c |V |φ v φ c . (2.22)
In principle, TDHF denotes an approach developed from a Hartree-Fock reference ground-state. In practice, the approximations used in this manuscript will instead use a DFT reference state, which can be efficiently obtained within a PW basis set implementation. A HF ground-state could be also used as a starting point, by rederiving the formalism using the ACFDT based on a HF reference [34].

In order to simplify the formalism we introduce an approximation in the full TDHF response function. This approximation was introduced in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF] and shown to lead to very accurate results. Specifically, the kernel Ξ in Eq. 2.20 can be divided in two different components

Ξ = B + Ξ -B = B B B B + A -B 0 0 A -B. (2.23)
and the identity

I -λ Π 0 (iω) Ξ -1 = I -λ Π 0 (iω) B -1 + λ I -λ( I -λ Π 0 (iω) B -1 Π 0 (iω) Ξ -B -1 × × I -λ Π 0 (iω) B -1 Π 0 (iω) Ξ -B I -λ Π 0 (iω) B -1 (2.24)
can be introduced. By discarding the last term in Eq. 2.24 a TDHF equation with an approximate kernel is obtained:

Π λ (iω) = ( I -λ Π 0 (iω) B) -1 Π 0 (iω); (2.25)
this equation involves exclusively the evaluation of products of valence and conduction states (conduction-conduction and valence-valence products are contained only in the A matrix) and will be denoted as electron-hole time-dependent Hartree-Fock (eh-TDHF). In Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF] it was shown that the eh-TDHF approximation leads to significantly more accurate correlation energies with respect to the RPA. Additional approximations based on the eh kernel (Eq. 2.25) have also been proposed in the literature. In Ref. [START_REF] Hesselmann | Correct description of the bond dissociation limit without breaking spin symmetry by a random-phaseapproximation correlation functional[END_REF], Heßelmann and Görling relied on this approximation to develop a local exact-exchange kernel and showed that the electron-hole kernel is equivalent to the full kernel for a two electron system. In Ref. 6, Bates and Furche used the B kernel to develop the RPA renormalized many-body perturbation theory. In the context of the plane-wave basis set Eq. 2.25 presents important practical advantages with respect to Eq. 2.19 (see Chapter 3 for details).

It is important to stress that the discussion in this section is specifically meant for ground state correlation energies. Indeed, the B kernel and the corresponding response function in Eq. 2.25 are good approximations within the ACFDT framework and are not intended, for example, to compute optical absorption spectra. For this kind of applications the B kernel does not contain the contributions necessary to correctly describe the physics of charge-transfer excitations in molecules [START_REF] Rocca | Ab initio calculations of optical absorption spectra: Solution of the bethe-salpeter equation within density matrix perturbation theory[END_REF] and of excitons in solids [1,[START_REF] Onida | Electronic excitations: density-functional versus many-body green's-function approaches[END_REF]. However, this approximate kernel would satisfy the f-sum rule for any possible reference state used to build the response function [START_REF] Malcio Glu | tur-boTDDFT -A code for the simulation of molecular spectra using the Liouville-Lanczos approach to time-dependent density-functional perturbation theory[END_REF][START_REF] Toulouse | Assessment of range-separated time-dependent density-functional theory for calculating c6 dispersion coefficients[END_REF][START_REF] Rocca | Ab initio optoelectronic properties of silicon nanoparticles: Excitation energies, sum rules, and tamm-dancoff approximation[END_REF]. Since the ACFDT requires an integral between -∞ and ∞ of the response functions (see Eq. 2.12), the fulfillment of the f-sum rule might play a role in making this kernel accurate for correlation energies.

We can now replace the eh-TDHF response function (Eq. 2.25) in the ACFDT equation for the correlation energy (Eq. 2.12). Because of the simplified structure of the electron-hole kernel B, the integration over the parameter λ can be carried out analytically and the trace can be expressed in terms of the N eh × N eh blocks only:

E eh-TDHF c = 1 2 ∞ -∞ dω 2π Tr ln I -χ 0 (iω)B B -1 K + χ 0 (iω)K , (2.26)
where χ 0 (iω) has been defined in Eq. 2.18. If the exchange component of B in Eq. 2.26 is set to 0 (namely B = K) the direct RPA approximation is found. By considering the Taylor expansion of the logarithm in Eq. 2.26 to the second order the MP2 approximation is obtained [30]:

E MP2 c = - 1 4 ∞ -∞ dω 2π Tr χ 0 (iω)Bχ 0 (iω)K = - 1 2 vc,v c K vc,v c B v c ,vc c + c -v -v . (2.27)
As already discussed in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF], the Taylor expansion of Eq. 2.26 is convergent only when all the eigenvalues of χ 0 (iω)B are smaller than 1. Nevertheless, the expression in Eq. 2.27 is mathematically equivalent to MP2 and represents a response function reformulation of this approximation [30].

Adiabatic connection second order screened exchange approximation

The adiabatic connection second order screened exchange (AC-SOSEX) approximation has been first proposed in Ref. 2 and formulated within a dielectric matrix framework in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]. Within the AC-SOSEX approximation the correlation energy is expressed as

E AC-SOSEX c = - 1 2 1 0 dλ ∞ -∞ dω 2π Tr Π λ RP A (iω) Ξ -Π 0 (iω) Ξ . (2.28)
where Ξ is defined by Eqs.2.20-2.22 and Π λ is obtained within the direct RPA:

Π λ RP A (iω) = ( I -λ Π 0 (iω) K) -1 Π 0 (iω). (2.29)
In Eq. 2.28 the antisymmetrization of the correlation energy, which is missing in the direct RPA, is introduced by the exchange part of Ξ. The trace in Eq. 2.28 can be more conveniently expressed as

Tr Π λ RP A Ξ -Π 0 Ξ = λTr ( I -λ Π 0 K) -1 Π 0 K Π 0 B +λTr ( I -λ Π 0 K) -1 Π 0 K Π 0 Ξ -B , (2.30)
where in the last line Ξ has been split as in Eq. 2.23 and the dependence on iω is kept implicit. By exploiting the properties of the trace operator and matrix algebra, the two contributions to the trace in Eq. 2.30 can be rewritten in terms of reduced size N eh × N eh matrices:

λTr ( I -λ Π 0 K) -1 Π 0 K Π 0 B = λTr (I -λ χ 0 K) -1 χ 0 Kχ 0 B (2.31)
and

λTr ( I -λ Π 0 K) -1 Π 0 K Π 0 Ξ -B = λTr (I -λχ 0 K) -1 (Π 0 + KΠ 0 + + Π 0 -KΠ 0 -)(A -B) .
(2.32) In the following we will show that the contribution to the correlation energy coming from Eq. 2.32 is likely to be small. By using the matrix identity

I -λ χ 0 K -1 = I + λ I -λ χ 0 K -1 χ 0 K (2.33)
we notice that the inverse of I -λ χ 0 K is given by the sum of the identity matrix and a smaller contribution. Indeed, since the eigenvalues of -χ 0 K are positive and decay to 0 rapidly [START_REF] Wilson | Efficient iterative method for calculations of dielectric matrices[END_REF], all the eigenvalues of λ(I -λ χ 0 K) -1 χ 0 K are < 1 and most of them are 1. By assuming (I -λ χ 0 K) -1 ≈ I, the contribution of Eq. 2.32 to the correlation energy in Eq. 2.28 can be computed analytically and the final result is exactly 0 [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]. Accordingly, the term in Eq. 2.32 is discarded. The effect of this approximation on the correlation energy was discussed in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF] (see, for example, Figure 2 and Table II) and it was found that neglecting the contribution in Eq. 2.32 leads to more accurate correlation energies. From a practical point of view this approximation allows us to express the correlation energy in terms of electron-hole integrals only. Indeed, by replacing Eq. 2.31 in Eq. 2.28 and carrying out an analytical integration over λ, the following expression for the AC-SOSEX is obtained:

E AC-SOSEX c = 1 2 ∞ -∞ dω 2π Tr ln I -χ 0 (iω)K K -1 B + χ 0 (iω)B . (2.34)
Similarly to Eq. 2.26, also Eq. 2.34 reduces to the RPA in the B = K case and to MP2 by Taylor expanding the logarithm.

Chapter 3

Numerical implementation within a plane-wave basis set using iterative algorithms

This chapter is based on original work previously published in Refs. [START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF]28,31].

The equations of RPA and beyond-RPA methods presented in Chapter 2 are based on an electron-hole representation (only conduction-valence products are present in the formalism). A first "brute-force" implementation of this formalism has been presented in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF] within the Gaussian basis set code Molpro [START_REF] Werner | Molpro: a general-purpose quantum chemistry program package[END_REF]. Results obtained with this implementation will be used in Chapter 4 for benchmark and comparison purposes. However, in this manuscript we will focus on algorithms and applications based on the PW basis set. This numerical approach is more flexible and allows us to treat molecules, bulk solids, and periodic nanostructures within the same framework. However, the size of the basis set can be large, easily reaching several hundreds of thousands for large supercells. In a similar way, the number of electron-hole pairs also becomes rapidly unmanageable. For this reason it is convenient to represent the response functions involved in RPA and beyond-RPA methods using a compact basis set. Here we use the eigenvectors of an approximate polarizability that includes only the kinetic energy contribution [START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF]. This method is an approximation of the projective dielectric eigenpotential technique (PDEP) of Refs. 150, 100. In order to evaluate efficiently the matrix elements of the polarizability on the auxiliary basis set it is useful to resort to additional numerical techniques: The resummation of conduction states to obtain a formalism that depends on only occupied states [3,[START_REF] Walker | Efficient approach to time-dependent density-functional perturbation theory for optical spectroscopy[END_REF][START_REF] Rocca | Turbo charging time-dependent density-functional theory with Lanczos chains[END_REF] and the use of the Lanczos algorithm to compute the resolvent of the Hamitonian [START_REF] Rocca | Turbo charging time-dependent density-functional theory with Lanczos chains[END_REF][START_REF] Umari | Gw quasiparticle spectra from occupied states only[END_REF][START_REF] Nguyen | Improving accuracy and efficiency of calculations of photoemission spectra within the many-body perturbation theory[END_REF]. These methodologies to compute RPA and beyond-RPA correlation energies have been implemented as a separate module in the Quantum Espresso package, which uses plane-waves and pseudopotentials [47,46].

For the sake of simplicity the formalism in this chapter will be discussed exclusively including only the Γ point in the integration over the first Brillouin zone and by considering the orbitals as real. This assumption is valid for calculations involving sufficiently large supercells. The same ideas presented here can be extended to periodic solids. In this case the evaluation of RPA and beyond-RPA correlation energies requires a double integration in the first Brillouin zone (over k-and q-points). Until now this has been done only for the RPA [START_REF] Kaoui | Random phase approximation correlation energy using a compact representation for linear response functions: application to solids[END_REF]. Results obtained with this implementation are presented in Sec. 4.4 of the next chapter.

The sections of this chapter are structured as follows:

• In Sec. 3.1 we recast the equations of Chapter 2 by introducing an auxiliary basis set. At this stage we suppose that an optimally small basis set for response functions exists but we do not define it explicitly.

• In Sec. 3.2 we show how the equations presented in Sec. 3.1 can be reformulated by avoiding any reference to unoccupied states. Within this formalism the response functions are recast as matrix elements of the resolvent of the Hamiltonian.

• In Sec. 3.3 we show how the resolvent of the Hamiltonian can be computed efficiently by using the Lanczos algorithm.

• In Sec. 3.4 we discuss the construction of the auxiliary basis set.

• In Sec. 3.5 we present an additional algorithm, based on the Gram-Schmidt orthogonalization, that is essential to further lower the numerical cost of beyond-RPA methods.

Implementation using an auxiliary basis set

To start the presentation of the formalism we suppose that a basis set exists that allows for a compact representation of the matrices involved in Eq. 2.17, 2.26, and Eq. 2.34 (where compact representation means considerably smaller than the N eh × N eh full size). The details on the generation of this basis set will be discussed later in Sec. 3.4. In the following we will denote the vectors belonging to this auxiliary basis set as Φ and their total number as N aux ; the discussion will be carried out for the eh-TDHF correlation energy but analogous equations can be introduced in the AC-SOSEX case. The equations for the RPA can be simply obtained from the eh-TDHF formalism by setting to 0 the exchange component.

While the focus of this manuscript is mainly related to the PW implementation, the formalism presented in this section is general and could be also employed for a density fitting implementation based on localized basis sets [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]. By using the completeness relation

I = G |Φ G Φ G |
B vc,v c = K vc,v c -K x vc,v c = 2 φ v φ c |V |φ v φ c -φ v φ c |V |φ v φ c , (3.1)
we obtain

K vc,v c = 2 G φ v φ c |V 1 2 |Φ G Φ G |V 1 2 |φ v φ c = G N vc,G N G,v c = N N T vc,v c (3.2) K x vc,v c = G φ v φ c |V 1 2 |Φ G Φ G |V 1 2 |φ v φ c = G M vc,G M G,v c = M M T vc,v c . (3.3)
Since the Coulomb potential is diagonal in a PW representation, the calculation of V

B)B -1 K = Tr ln(I -χ 0 K + χ 0 K x )(K -K x ) -1 K = Tr ln(I -χ 0 K + χ 0 K x )(K -K x ) -1 (χ 0 ) -1 χ 0 K = Tr f χ 0 K -χ 0 K x χ 0 K = Tr f χ 0 K -χ 0 K x χ 0 K(χ 0 K) -1 χ 0 K = Tr f χ 0 N N T -χ 0 M M T χ 0 N N T (χ 0 N N T ) -1 χ 0 N N T = Tr N T f χ 0 N N T -χ 0 M M T χ 0 N N T (χ 0 N N T ) -1 χ 0 N = Tr f N T χ 0 N -N T χ 0 M M T χ 0 N (N T χ 0 N ) -1 N T χ 0 N , (3.4)
where we have defined the function of a matrix f (A) = ln(I -A)A -1 and used the cyclic property of the trace operator; additionally, the result in Eq. 3.4 exploits properties of functions of a matrix that are described in more detail in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]. The dependence of χ 0 on iω is kept implicit in Eq. 3.4.

It is now convenient to introduce the definitions

C QR (iω) = (N T χ 0 N ) QR = Φ Q |V 1 2 χ0 V 1 2 |Φ R , = 4 vc Φ Q |V 1 2 |φ v φ c R 1 iω + v -c φ v φ c |V 1 2 |Φ R , (3.5) J vv QG (iω) = 2 c Φ Q |V 1 2 |φ v φ c R 1 iω + v -c φ v φ c |V 1 2 |Φ G , (3.6)
L vv GR (iω) = 2 c Φ G |V 1 2 |φ v φ c R 1 iω + v -c φ v φ c |V 1 2 |Φ R , (3.7)
and

Z QR (iω) = (N T χ 0 M M T χ 0 N ) QR = 2 vv G J vv QG (iω)L vv GR (iω); (3.8)
the compact notation in terms of the operator χ0 on the second line of Eq 3.5 will be useful later. In Eqs. 3.5-3.8 the subscripts Q, R, and G are used to denote different vectors Φ belonging to the auxiliary basis set. By using the definitions in Eqs. 3.5-3.8 we can finally reformulate the eh-TDHF correlation energy expression in Eq. 2.26 to obtain

E eh-TDHF c = 1 2 ∞ -∞ dω 2π Tr ln I -C(iω) + Z(iω)C -1 (iω) × C(iω) -Z(iω)C -1 (iω) -1 C(iω) + C(iω) . (3.9)
This final expression for the correlation energy is more convenient than the original form Eq. 2.26 since it involves N aux × N aux matrices, where N aux N eh . If the exchange term in the eh-TDHF kernel is set to zero we have that Z(iω) = 0 and the correlation energy reduces to the RPA expression

E RPA c = 1 2 ∞ -∞ dω 2π Tr ln [I -C(iω)] + C(iω) ; (3.10)
where I -C is nothing but the symmetrized dielectric matrix ε represented in the Φ basis set [START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF]. By Taylor expanding the logarithm in Eq. 3.9 up to second order we obtain the MP2 approximation:

E MP2 c = -1 2 ∞ 0 dω 2π tr C 2 (iω) -Z(iω) . (3.11)
This equation reformulates the MP2 correlation energy in a form that could be easily computed using an implementation derived from our eh-TDHF and AC-SOSEX code. However, RPA and beyond-RPA methods work as post-processing to a standard DFT calculation using local or semi-local functionals. While being numerically stable, the MP2 correlation energies computed through Eq. 3.11 using a DFT starting point are incorrect, since the Brillouin theorem is not anymore valid. Indeed the MP2 correlation energy is rather derived in a post-Hartree-Fock framework and is not supposed to use DFT orbitals and energy levels as starting point [START_REF] Szabo | Modern quantum chemistry: introduction to advanced electronic structure theory[END_REF]. We proposed an implementation of Eq. 3.11 in a Hartree-Fock framework in Ref. [30]. However, the calculation of the HF exact exchange is highly expensive within the PW basis set and this also limits the applicability of our MP2 code to systems with few valence electrons. This MP2 implementation will not be further discussed on this manuscript, which will remain focused on RPA-type methods. Also the AC-SOSEX correlation energy (Eq. 2.34) can be rewritten in terms of the compact auxiliary basis set:

E AC-SOSEX c = 1 2 ∞ -∞ dω 2π Tr [ln (I -C(iω)) + C(iω)] × I -C -1 (iω)Z(iω)C -1 (iω) . (3.12)
Similarly to Eq. 3.9, also the AC-SOSEX reduces to MP2 by Taylor expanding the logarithm and to RPA by setting Z(iω) = 0.

To conclude this section it is useful to show explicitly the form of Eq. 3.5 when plane-waves are used as an auxiliary basis set. Indeed, this formulation is used in traditional implementations to compute RPA correlation energies [39,[START_REF] Marini | First-principles description of correlation effects in layered materials[END_REF][START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF]. By denoting the reciprocal-lattice vectors with G and G we have

C G,G (iω) = 4 vc G|V 1 2 |φ v φ c R 1 iω + v -c φ v φ c |V 1 2 |G = 4 vc 1 G φ v |e iG•r |φ c R 1 iω + v -c φ v |e -iG •r |φ c 1 G = 1 G χ 0 G,G 1 G (3.13)
Since in traditional RPA implementations the response functions are represented using the PW basis set, practical calculations involve storage and diagonalization of matrices that become rapidly large by increasing the kinetic energy cut-off and/or the cell size [39,[START_REF] Marini | First-principles description of correlation effects in layered materials[END_REF][START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Shepherd | Convergence of many-body wave-function expansions using a plane-wave basis: From homogeneous electron gas to solid state systems[END_REF]. Additionally, the full diagonalization of the Hamiltonian is required to evaluate RPA and beyond-RPA correlation energies and also this task is influenced by the number of PWs. This last issue will be addressed in the next section.

Elimination of the empty states

Eqs. 3.5-3.8 can be evaluated directly by diagonalizing the ground-state KS Hamiltonian Ĥ to obtain conduction and valence states. However, the large dimension of the PW basis set leads to a very high number of conduction states and fully converged results can be obtained only through extrapolation techniques [START_REF] Marini | First-principles description of correlation effects in layered materials[END_REF] or by diagonalizing the full Hamiltonian. Similarly to recent work on the calculation of RPA correlation energies [START_REF] Nguyen | Efficient calculation of exact exchange and rpa correlation energies in the adiabatic-connection fluctuation-dissipation theory[END_REF][START_REF] Lu | Ab initio calculation of van der waals bonded molecular crystals[END_REF][START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF], our numerical implementation takes instead advantage of DFPT techniques [3]. Indeed, we can notice that Eq. 3.5 and Eq. 3.6 can be rewritten in the form

C QR (iω) = 4 vc Φ Q |V 1 2 |φ v φ c R 1 iω + v -c φ v φ c |V 1 2 |Φ R = 4 vc φ v |Φ Q V 1 2 |φ c φ c | R 1 iω + v - Ĥ |φ c φ c |V 1 2 Φ R |φ v = 4 v φ v |Φ Q V 1 2 P R 1 iω + v - Ĥ P V 1 2 Φ R |φ v (3.14)
and

J vv QG (iω) = 2 c Φ Q |V 1 2 |φ v φ c R 1 iω + v -c φ v φ c |V 1 2 |Φ G = 2 c φ v |Φ Q V 1 2 |φ c φ c | R 1 iω + v - Ĥ |φ c φ c |V 1 2 Φ G |φ v = 2 φ v |Φ Q V 1 2 P R 1 iω + v - Ĥ P V 1 2 Φ G |φ v , (3.15)
respectively; also Eq. 3.7 can be written in an analogous way. In Eqs. 3.14-3.15 P = c |φ c φ c | is the projector onto the conduction state subspace. By exploiting the completeness of the basis set, the projector onto the conduction state subspace can be expressed as P = Îv |φ v φ v |, where in this case the sum is limited to the occupied states only [3]. This approach does not require approximations and automatically includes all the conduction states described by the PW basis set. Once the auxiliary basis set elements Φ Q and Φ R are fixed the matrix elements involved in the last line of Eq. 3.14 and Eq. 3.15 can be efficiently computed by using the Lanczos algorithm for response functions [START_REF] Rocca | Turbo charging time-dependent density-functional theory with Lanczos chains[END_REF][START_REF] Umari | Gw quasiparticle spectra from occupied states only[END_REF][START_REF] Nguyen | Improving accuracy and efficiency of calculations of photoemission spectra within the many-body perturbation theory[END_REF]. This approach will be discussed in detail in the next section.

Lanczos algorithm

Eqs. 3.14-3.15 allow us to compute RPA and beyond-RPA correlation energies without referring explicitly to unoccupied states. However, this formulation requires the computation of the resolvent of the Hamiltonian Ĥ for many different values of iω. Specifically, the mathematical problem that needs to be solved in Eqs. 3.14-3.15 can be formulated in the following general way:

g(iω) = u|(iω -Â) -1 |v , (3.16)
where  is a Hermitian operator. For the specific case of Eqs. 3.14-3.15 we have  = Ĥv . Depending on the number of occupied states and the size of the auxiliary basis set, the matrix element of the resolvent in Eq. 3.16 has to be computed for several different vectors u and v (see also Secs. 3.4 and 3.5). The calculation of g(iω) can be obtained by using standard iterative techniques to solve Hermitian linear systems (e.g., the conjugate gradient algorithm) [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. In this case the linear system (iω -Â) -1 |z = v is solved and g(iω) is computed as u|z . The disadvantage of this approach is that a different linear system has to be solved for each different value of iω.

The Lanczos algorithm is known to be an efficient method to compute the resolvent of Hermitian operators by performing a single iterative recursion independent of ω [START_REF] Haydock | Electronic structure based on the local atomic environment for tight-binding bands[END_REF].The standard Lanczos algorithm is limited to the case u=v; however, for our purpose it is necessary to compute also the off-diagonal elements of g(iω). For this task we use the Lanczos algorithm proposed in Ref. 126 for timedependent density functional theory and later applied to GW calculations [START_REF] Umari | Gw quasiparticle spectra from occupied states only[END_REF][START_REF] Nguyen | Improving accuracy and efficiency of calculations of photoemission spectra within the many-body perturbation theory[END_REF]. This algorithm generates iteratively a series of vectors {q 1 , q 2 , • • • } by using the following procedure:

q 0 = 0, q 1 = v/ v|v , β n+1 |q n+1 = Â|q n -α n |q n -β n |q n-1 (3.17)
where β n+1 is determined in order to impose the normalization condition q n+1 |q n+1 = 1 and α n = q n | Â|q n . In the orthonormal basis set of the vectors {q 1 , q 2 , • • • } the matrix  has tridiagonal form:

T j =         α 1 β 2 0 • • • 0 β 2 α 2 β 3 0 . . . 0 β 3 α 3 . . . 0 . . . 0 . . . . . . β j 0 • • • 0 β j α j         (3.18)
where α and β are the coefficients of the Lanczos recursion (Eq. 3.17) and j is the maximum number of Lanczos iterations performed. The value of g(iω) can then be approximated as

g(iω) ≈ ζ j |(iω -T j ) -1 |e j 1 . (3.19)
where ζ jT is a j-dimensional vector defined as ( u|q 1 , u|q 2 , • • • , u|q j ) and e j 1 is the j-dimensional unit vector (1, 0, • • • , 0). In general, the dimension j (the number of Lanczos iterations) of the matrix T j necessary to obtain an accurate approximation of g(iω) is much smaller than the dimension of the full matrix Â.

There are two important advantages in using Eq. 3.19. First, the matrix T j and the corresponding Lanczos iterative recursion do not depend on iω. Once T j is generated, the value of g(iω) can be computed for many values of iω by simple linear algebra operations in a small j-dimensional space (Eq. 3.19). For all the numerical applications considered in this manuscript 30 Lanczos iterations have been found to be sufficient to achieve a satisfactory level of accuracy. Second, the Lanczos recursion (Eq. 3.17) depends only on the right vector |v and different ζ j vectors can be generated "on the fly" during the Lanczos chain in order to compute matrix elements of g(iω) for different left vectors u|. This feature is particularly convenient to evaluate Eqs. 3.14-3.15, where the calculation of several off-diagonal elements (u = v) is required. One of the main advantages of the present method is the efficient treatment of the dynamical effects in the response function. Once the tridiagonal matrix T l is computed through the Lanczos algorithm the response functions can be efficiently evaluated for several different values of ω (see Eq. 3.19). In this work the integrals over ω in Eqs. 3.9, 3.10, and 3.12 are evaluated by a 20-point Gauss-Legendre quadrature in the interval 0 to 200 Ry. The accuracy of this integration scheme can be estimated by comparison with the simple trapezoidal rule. An example has been presented in Ref. [START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF]: By considering 100 kinetic-PDEP basis vectors for the benzene dimer the difference between the values of ∆E c obtained with the trapezoidal rule (uniform grid with a 0.0001 Ry step) and the Gauss-Legendre quadrature is of the order of 10 -6 kcal/mol. The results are also well converged with respect to the integration range: By integrating between 0 and 100 Ry the change in the value of ∆E c is of the order of 10 -5 kcal/mol.

Construction of a compact auxiliary basis set by diagonalizing an approximate response function

Despite the several advantages of the Lanczos algorithm, in order to make this procedure efficient it is necessary to consider a basis set {Φ I } as small as possible. Indeed, even if the Lanczos procedure requires only a few iterations, it is still necessary to perform a number of Lanczos chains which is proportional to the size of the basis set {Φ I }. In order to discuss the construction of the optimal basis we redefine Eq. 3.5 as:

C QR = 4 vc Φ Q |V 1 2 |φ v φ c R 1 v -c φ v φ c |V 1 2 |Φ R = 4 vc ΦQ |φ v φ c R 1 v -c φ v φ c | ΦR (3.20) = ΦQ (r)χ 0 (r, r ) ΦR (r )drdr , (3.21)
where we focus on the static case (ω = 0). It is possible to switch between the two alternative forms of the auxiliary basis set Φ and Φ simply by multiplying by V 1 2 or its inverse (which are diagonal in reciprocal space). In Eq. 3.21 we used the static polarizability in the real-space representation

χ 0 (r, r ) = 4 cv φ v (r)φ c (r)φ c (r )φ v (r ) v -c ; (3.22)
this definition will be useful below. From Eq. 3.20 it is clear that an optimally small basis set Φ needs to efficiently represent the φ c φ v products. While φ c and φ v are orthonormal orbitals, the products φ c φ v might have a strong linear dependence. This observation suggests the possibility of representing χ 0 on a relatively small basis set (at least much smaller than the number of products φ c φ v ). The resolution of the identity techniques used to speed up localized basis set RPA calculations are based on this idea [35,[START_REF] Ren | Resolution-ofidentity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF]. However, in a PW implementation the number of products φ c φ v could reach several millions and the corresponding subspace would still require a large basis set to be accurately represented. In a PW framework an optimally small basis set for χ 0 (or ˜ ) can be generated by "eliminating" the linear dependence of the products φ c φ v keeping into account at the same time the weighting factor 1/( vc ) (see Eq. The idea of keeping into account the linear dependence of the φ c φ v products together with the weight 1/( vc ) is implicit in projective dielectric eigenpotential (PDEP) method [START_REF] Wilson | Efficient iterative method for calculations of dielectric matrices[END_REF][START_REF] Nguyen | Efficient calculation of exact exchange and rpa correlation energies in the adiabatic-connection fluctuation-dissipation theory[END_REF][START_REF] Lu | Ab initio calculation of van der waals bonded molecular crystals[END_REF][START_REF] Ping | Electronic excitations in light absorbers for photoelectrochemical energy conversion: first principles calculations based on many body perturbation theory[END_REF]. Within this approach the response functions χ 0 or, equivalently, ε are iteratively diagonalized. At each iteration χ 0 is applied to a set of trial potentials { Φtrial I } by using Eqs. 3.21 and 3.14:

χ 0 (r, r ) Φtrial I (r )dr = 4 v φ v (r) r|( v -Ĥ) -1 P | Φtrial I φ v (3.23)
until the iterative algorithm (e.g. conjugate gradient) reaches convergence; Eq. 3.23 can be easily generalized to the dynamical case. The diagonalization of the static ε has been used to build the optimal basis set used in the Lanczos-based GW implementation of Refs. 102, 109. In the PDEP method only a limited number of dominant eigenvalues and eigenpotentials (eigenvectors) are necessary to represent χ 0 or ε. Indeed, the eigenpotentials corresponding to the largest (in absolute value) eigenvalues have a strong superimposition with the products φ c φ v through the term P | Φtrial I φ v in Eq. 3.23, while at the same time ( v -Ĥ) -1 keeps into account the weight of the different transitions. This is done without explicitly computing the empty states φ c of the system. This qualitative argument suggests why the PDEP technique provides a compact representation for χ 0 and ˜ . The PDEP method has been used to compute RPA correlation energies by diagonalizing ˜ for each value of iω included in the numerical integration of Eq. 3.10 [START_REF] Lu | Ab initio calculation of van der waals bonded molecular crystals[END_REF][START_REF] Nguyen | Efficient calculation of exact exchange and rpa correlation energies in the adiabatic-connection fluctuation-dissipation theory[END_REF]. The main drawback of the PDEP procedure is that, each time χ 0 is applied to a vector through Eq. 3.23, the ( v -Ĥ) -1 operator has to be evaluated by solving a linear system.

In this manuscript we discuss a new and more efficient method to construct a basis set to compactly represent dielectric matrices. The main idea consists in noticing that for sufficiently high energies c is mostly determined by the kinetic energy contribution kin c . For this reason we approximate the weighting factor in Eq. 3.22 as 1/( vc ) ≈ 1/( v -kin c ). This approximation describes well high energy excitations but is certainly not valid for the low energy c 's. However, since the main purpose of this procedure is to build a basis set that "discard" the high energy contributions to Eq. 3.22, this approximation is acceptable and works well in practice (see also the numerical results of Chapter 4). Since we want to obtain a compact representation for χ 0 without computing explicitly empty states, similarly to the PDEP procedure, we consider the iterative diagonalization of a χ 0 kin operator containing only the kinetic energy contribution. This operator can be applied to a trial potential in the following way:

χ 0 kin (r, r ) Φtrial I (r )dr = 4 v φ v (r) r| P ( v + ∇ 2 /2) -1 P | Φtrial I φ v , (3.24) -8 -7 -6 -5 -4 -3 -2 -
C 2 H 2 +H 2 O → CH 3 CHO reaction. Specifically, ∆E c = E CH 3 CHO c -E C 2 H 2 c -E H 2 O c .
where the kinetic energy operator -∇ 2 /2 is diagonal in a plane-wave representation (it is proportional to G 2 in reciprocal space). This is one of the main advantages of Eq. 3.24: Since ( v + ∇ 2 /2) is diagonal its inversion is trivial and the implementation to solve Eq. 3.24 is simple and efficient. We will call this procedure kinetic-PDEP. Eq. 3.24 gives an accurate description of high energy excitations while the low energy excitations are approximated in a rather poor way. For example, the term ( v + ∇ 2 /2) does not contain any information on the exact value of the gap of the system. However, it is important to stress that the only purpose of the kinetic-PDEP procedure is to build a compact basis whose size can be systematically increased to achieve convergence. The details of the real electronic energy levels (as well as the dynamical contributions) are included by the Lanczos algorithm. In Eq. 3.24 χ 0 kin is considered exclusively in the static limit (iω = 0); this approximation has been proven reliable to built an effective basis set also for the dynamical case (iω = 0) [START_REF] Nguyen | Improving accuracy and efficiency of calculations of photoemission spectra within the many-body perturbation theory[END_REF][START_REF] Pham | G w calculations using the spectral decomposition of the dielectric matrix: Verification, validation, and comparison of methods[END_REF][START_REF] Govoni | Large scale gw calculations[END_REF]. However, the structure of the matrix elements C QR in Eq. 3.14 is different from the structure of J vv QG (or, similarly, L vv GR ) in Eq. 3.15: J vv QG does not involve a summation over the occupied states. From a rigorous point of view the basis set obtained from the kinetic-PDEP (or PDEP) method is designed to compactly represent C QR in Eq. 3.5, which is the main ingredient for RPA calculations. As a consequence, it is not straightforward to conclude that the kinetic-PDEP basis set would also efficiently represent the matrix elements in Eq. 3.6 and Eq. 3.7. However, as discussed qualitatively in Ref. [START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF], the eigenvectors obtained by diagonalizing the operator in Eq. 3.24 well represent the occupied-virtual orbital products by eliminating at the same time their linear dependence and the excitations too high in energy, that weakly contribute to the response functions. Since Eq. 3.2 and Eq. 3.3 also involve occupied-virtual state products and excitation energies at the denominator, it is reasonable to suppose that the kinetic-PDEP basis set would also well represent these matrix elements.

From numerical calculations we observed that this is indeed the case. In Fig. 3.1, we show the convergence with respect to the auxiliary basis set of the correlation energy contribution to the energy of the C 2 H 2 +H 2 O → CH 3 CHO reaction (this reaction is part of a test set that will be considered in detail in Sec. 4.3). Specifically,

∆E c = E CH 3 CHO c -E C 2 H 2 c -E H 2 O c
. Fig. 3.1 allows us to observe the rapid convergence of ∆E c as a function of the size of the auxiliary basis set. For example, by increasing N aux from 270 to 360 the values ∆E c evaluated at the AC-SOSEX and eh-TDHF are all contained within a 0.04 kcal/mol interval for AC-SOSEX and eh-TDHF and within a 0.1 kcal/mol interval for RPA. In general, we observed that the eh-TDHF and AC-SOSEX correlation energies have a convergence rate with respect to N aux which is similar to the RPA case and sometimes even slightly better. This allows us to use a more compact basis set to solve Eq. 2.26, whose direct solution would be very challenging by expressing the response functions on the eh basis set or the plane-wave basis set. For example, in the case shown in Fig. 3.1 the use of the full plane-wave basis set for the potential to directly represent the response functions in Eq. 2.26 or Eq. 2.34 would involve matrices with a 850000 × 850000 size; alternatively, the use of Eq. 3.9 and Eq. 3.12 with the kinetic-PDEP auxiliary basis set leads to converged results with as little as 360 basis elements corresponding to 360 × 360 matrices.

In order to converge energy differences as done in Fig. 3.1 it is necessary to use a number of basis vectors proportional to the number of electrons for each component. For example, if 360 basis vectors are used to compute the correlation energy of CH 3 CHO, for C 2 H 2 and H 2 O we should use 200 and 160 vectors, respectively. This procedure will be described in detail in Chapter 5.

Gram-Schmidt orthogonalization

Even though the use of the compact auxiliary basis greatly reduces the dimensionality of dielectric matrices, the calculation requires the evaluation (N aux ) 2 matrix elements for C and (N aux N v ) 2 matrix elements for J and L. For the following discussion, we introduce the compact notation

|U v G = P V 1 2 Φ G |φ v and a similar definition for U v G |.
As seen in Eq. 3.10 the calculation of the RPA correlation energy requires the evaluation of matrix elements of the type:

C QR (iω) = 4 v U v Q | iω + v -Ĥ -1 |U v R ; (3.25)
which is just Eq. 3.14 written in a more compact notation. The AC-SOSEX and eh-TDHF require the evaluation of additional terms such as:

J vv QG (iω) = 2 U v Q | iω + v -Ĥ -1 |U v G (3.26)
The matrix elements in Eqs. 3.25-3.26 can be efficiently computed by using the Lanczos algorithm. Since the subscripts v and v are in general different in Eq. 3.26, the numerical complexity of beyond-RPA methods is higher than that of RPA, which requires only the evaluation of Eq. 3.25. This is a major bottleneck for the AC-SOSEX and eh-TDHF methods whose numerical complexity scales as

N P W × N 2 aux × N 2 v ,
where N P W is the number of plane-wave, N aux the number of auxiliary basis vectors, and N v the number of valence states (compared to

N P W × N 2 aux × N v for RPA).
Intuitively, this can be understood by considering that the Lanczos algorithm computes the solution to the N v × N aux linear systems

(iω + v -Ĥ) -1 |U v G ; each solution, represented on N P W plane-waves, is pro- jected over the N v × N aux left vectors U v Q |.
Additionally, in order to store all the left vectors the memory requirements for beyond-RPA calculations scale as N P W × N v × N aux (in RPA it is sufficient to store only one value of v at a time); this problem could be overcome by recomputing these matrix elements on the fly or by storing them on disk but this comes at a big loss in efficiency. Because of these computational and memory issues, the applicability of the AC-SOSEX and eh-TDHF methods would be limited to molecules with a maximum of twenty to thirty electrons. For example, the largest dimers in the S22 test set considered in Sec. 4.1.2 (≈ 100 electrons) would be completely out of reach.

In order to overcome these difficulties in applying the beyond-RPA methods to large systems, we use a procedure inspired by previous work of Umari and Baroni in the context of GW calculations [START_REF] Umari | Gw quasiparticle spectra from occupied states only[END_REF]. Specifically, we reduce the large number of right and left vectors in Eq. 3.26 by applying a Gram-Schmidt orthogonalization procedure to eliminate vectors that are linearly dependent. To do so, we define a new set of orthonormal vectors as,

|O 1 = |U 1 1 /N 1 |O 2 = |U 2 1 -O 1 |U 2 1 |O 1 /N 2 • • • |O Naux×Nv = |U Nv Naux - Naux×Nv i=1 O i |U Nv Naux |O i /N (Naux×Nv) , (3.27)
where the factors N n at the denominator are used to normalize each vector of the basis |O . Up to this point this procedure does not allow for any computational saving but rather introduces additional expensive numerical computations. However, if for a certain |O n we have N n ≈ 0, this vector can be fully represented by the previous ones

(|O 1 , • • • , |O n-1
) and does not need to be included in the orthonormalized basis set. By setting the threshold for the normalizing factors to a reasonable value,this idea can be used to generate a basis set

(|O 1 , • • • , |O N ortho ) where N ortho N aux × N v .
The threshold on the normalizing factors N n has to be chosen with a compromise between computational time and accuracy. In this work we chose a variable threshold, that will be referred to as "default threshold" (DT): 10 -5 for the first 20% of the vectors |U , 10 -4 between 20% and 50%, and 10 -3 for all the remaining vectors. In this way the first vectors, which give a more important contribution to the final result, are more accurately represented. While this specific choice is to a certain extent arbitrary, we will further discuss this point at the end of this section.

The gain in practical calculations can be quantitatively understood by looking at Fig. 3.2, where data for all the dimers and monomers considered in Secs. 4.1.1 and 4.1.2 are presented (92 monomers and 46 dimers). In this figure N aux is fixed to forty times the number of KS valence states N v , which is sufficient to reach convergence for almost all systems (see discussion in Sec. 4.1 and Chapter 5). While by definition N aux × N v grows quadratically with the number of electrons, N ortho has an approximately linear growth. This leads to a significant saving, with the memory necessary to store the vectors reduced to N P W × N ortho and the numerical complexity reduced to

N P W × N 2 ortho .
This linear growth is to a certain extent related to the threshold choice. Indeed, in the case of a sufficiently small threshold N ortho would be comparable to N aux × N v . With the DT this approximately linear behavior has been deduced empirically considering 138 molecular systems. Since this is a rather extended ensemble, we believe that this effect is robust and generalizable to all the molecules. However, in order to be meaningful, this observation of a linear behavior has to be associated with a high level of accuracy when only N ortho vectors are included in the calculations rather than the full N aux × N v vector set. This point is further discussed at the end of this section and in Sec. 4.1.1. As it will be shown in Sec. 4.1.2, this procedure allows for calculations for systems with up to 98 electrons and about 400k plane-waves. In order to compute the correlation energy in Eqs. 3.9 and 3.12 the Lanczos algorithm is used to evaluate the matrix elements

O i | iω + v -Ĥ -1
|O j and the matrix J in Eq. 3.26 is then reconstructed as

J vv QG = 2 N ortho i,j=1 U v Q |O i O i | iω + v -Ĥ -1 |O j O j |U v G . (3.28)
Similar reconstructions can be obtained for C and L. The evaluation of Eq. 3.28 can be performed with a N ortho ×N 2 aux ×N 2 v scaling. Accordingly, this implementation still scales with the fifth power of the system size. However, since N ortho is much smaller than N P W (up to a factor 60 for the largest system considered in this work), computational time is significantly reduced and calculations of much larger systems become feasible. In the following we will use "GSO" to denote the new optimized version of our code and "original" to denote the code that does not use the orthogonalization.

To conclude this section we further discuss the specific choice of the threshold we made. Table 3.1 shows the effect of this parameter on the accuracy of the Table 3.1: Binding energy of the formaldehyde dimer in the A24 test set as computed for different approximations (RPA, eh-TDHF, and AC-SOSEX) with different thresholds for the Gram-Schmidt orthogonalization (GSO). DT denotes the default threshold described in the text. In the last line of the table we report the total time taken to complete the calculation on an Intel Xeon mchine with 14 cores. with the original version of the code; the last line on the table reports the time taken to obtain the results on a node of our cluster (Intel Xeon CPU E5-2680 v4 @ 2.40GHz, 14 cores). Since the RPA, eh-TDHF, and AC-SOSEX methods require the evaluation of the same matrices, the three correlation energies are all provided by a single calculation and, accordingly, a single time is provided in Table 3.1.

Original code GSO DT DT×10 DT

As expected, the GSO implementation significantly speed up calculations and the time decreases by increasing the threshold in the GSO procedure. The difference in the time taken by the DT and DT×10 calculations is rather small; this can be understood by considering that a large threshold significantly decreases N ortho but at the same time the cost of rebuilding the response function (Eq. 3.28) increases. Discarding some small fluctuations in part related to the auxiliary basis set convergence, the original, the GSO DT, and the GSO DT/10 implementations all give results quantitatively in good agreement. Differently, the GSO DT×10 calculation leads to an error of about 0.4 kcal/mol for the RPA result. By considering these results the choice of the DT provides an excellent compromise between accuracy and computer time. In Sec. 4.1.1 we will systematically discuss the influence of the Gram-Schmidt procedure with DT on the accuracy of the binding energies of the full A24 test set.

Chapter 4

Numerical applications to molecules and solids

This chapter is based on original work previously published in Refs. [28,31,[START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF][START_REF] Kaoui | Random phase approximation correlation energy using a compact representation for linear response functions: application to solids[END_REF].

In this chapter we present a series of different applications of the methods introduced in Chapter 2 using the implementation discussed in Chapter 3. Specifically, we will consider weakly bound molecular dimers, molecular reaction energies involving the breaking and formation of covalent bonds, and bulk solids (for the latter the discussion will be limited to the RPA, since an implementation of eh-TDHF and AC-SOSEX including the first Brillouin zone integration is currently not yet available). The accuracy of RPA and beyond-RPA methods will be determined by comparison with high-accuracy reference values obtained from coupled cluster theory (molecules) and experiment (solids). It will be shown that the inclusion of exchange effects beyond the RPA significantly improves the results.

Ideally, ACFDT correlation energies should be computed self-consistently [START_REF] Nguyen | Ab initio selfconsistent total-energy calculations within the exx/rpa formalism[END_REF], as typically done for DFT approaches. In practice, these types of calculations are cumbersome and here we will rely exclusively on a fixed starting point obtained from an affordable DFT approximation. Specifically, all the results for molecules presented in this manuscript are based on response functions built from GGA orbitals and energy levels [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. Instead, GGA-based RPA results for solids will be compared with correlation energies using a LDA starting point.

Applications to weakly bound molecular dimers 4.1.1 A24 Test Set

We present here the results for the A24 test set proposed by Řezáč and Hobza [START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF], which contains interaction energies for 24 molecular dimers. The results are compared with the CCSD(T) and CCSDT(Q) values from the same reference. Our calculations do not include core correlation and relativistic effects, and the reference values have been accordingly chosen. All dimers were placed in a 30 × 30 × 30 a 3 0 box, and the plane-wave cut-off was set to 80 Ry to achieve sufficient convergence. This cut-off corresponds to about 160k plane-waves. Norm-conserving FHI pseudopotentials were used [40]. For the RPA and beyond-RPA calculations, the number of auxiliary basis vectors was set to twenty times the number of electrons. Since we have non-spin polarized calculations, this corresponds to N aux = 40 × N v . The number of elements included in the auxiliary basis is orders of magnitude smaller than the number of plane-waves.

We start the discussion by validating the Gram-Schmidt orthogonalization (GSO) algorithm introduced in Sec. 3.5. Since the dimers in the A24 test set are relatively small, results with and without GSO can be compared to understand the implications of this procedure on the numerical accuracy. For the S22 test set, discussed in the next subsection, the GSO algorithm is indispensable to obtain the correlation energies for the largest dimers. Figure 4.1 shows the differences between the interaction energies calculated with the original and the GSO versions. As can be observed, the differences are small, below 0.085 kcal/mol for all the molecules. As indicated in Table 4.1, the mean absolute error (MAE) of the GSO results with respect to the "exact" results is 0.047 kcal/mol for the RPA, 0.033 kcal/mol for the eh-TDHF, 0.029 kcal/mol for the AC-SOSEX approximations. Therefore, the GSO algorithm can be applied to calculate correlation energies at the RPA and beyond-RPA levels of theory for different types of systems, e.g., weakly interacting systems, which require stricter convergence criteria. From now we will discuss exclusively results for the A24 and S22 (in the next subsection) test sets obtained using the GSO version of the code. To compare our results with CCSD(T) and CCSDT(Q) values, we used slightly more strict convergence parameters for a few systems. Specifically, N aux was increased to sixty times the number of valence states for the borane-methane, ethene-Ar, methane-HF and methane-water dimers to ensure a convergence of the RPA interaction energy to within 0.09 kcal/mol.

The fully converged results are shown in Table 4.2. The maximum absolute deviation (|MAX|) of the RPA and beyond-RPA methods from the CCSDT(Q) reference values is found for the formaldehyde dimer, with almost 1.8 kcal/mol for the RPA, 1.2 kcal/mol for the eh-TDHF and 0.9 kcal/mol for the AC-SOSEX approximations. It could also noticed that in the comparison between the CCSD(T) and CCSDT(Q), the maximum deviation was also found for the formaldehyde dimer. In general, the inclusion of approximate exchange effects in the beyond-RPA methods improve results over the RPA. This is most visible for the hydrogenbonded systems, where the mean absolute errors for the eh-TDHF and AC-SOSEX are 0.356 kcal/mol and 0.339 kcal/mol, respectively, as compared to the value of 0.605 kcal/mol for the RPA, and for the mixed-type systems, where the mean absolute errors for the eh-TDHF and AC-SOSEX are 0.224 kcal/mol and 0.179 kcal/mol, respectively, as compared to the value of 0.442 kcal/mol for the RPA. In the case of dispersion-dominated systems, the improvement is less clear, with the mean absolute error for the eh-TDHF and AC-SOSEX are 0.242 kcal/mol and 0.314 kcal/mol respectively compared to 0.264 kcal/mol for the RPA. While these differences are small compared to numerical accuracy, for these systems only the eh-TDHF method slightly improves over the RPA. The deviations with respect to CCSDT(Q) are plotted in Fig. 4.2. Positive mean errors for the RPA and beyond-RPA methods indicate a general tendency of all these methods to consistently underestimate the absolute binding energies of the various systems. We observe that the beyond-RPA methods improve over the MP2 values reported in the literature [19], where MAEs of 0.389 kcal/mol and 0.364 kcal/mol were found for the counter-poise corrected and uncorrected case, respectively. We also observe a good agreement of our values with earlier works. In the case of RPA, our values have a mean absolute error of 0.409 kcal/mol compared to 0.338 kcal/mol obtained by Bates and Furche [6] and 0.4 kcal/mol obtained by Grundei and Burow [START_REF] Grundei | Random Phase Approximation for Periodic Systems Employing Direct Coulomb Lattice Summation[END_REF] within a localized basis set implementation. For the AC-SOSEX calculations, our results have a mean absolute error of 0.263 kcal/mol compared to 0.208 kcal/mol obtained by Bates and Furche [6]. For the sake of completeness the errors of CCSD(T) [START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF] are also included.

S22 Test Set

We apply here the more efficient algorithm based on the GSO to the S22 test set proposed by Jurečka et al., which consists of interaction energies for 22 molecular dimers [START_REF] Jure Čka | Benchmark database of accurate (mp2 and ccsd (t) complete basis set limit) interaction energies of small model complexes, dna base pairs, and amino acid pairs[END_REF]. This test set has been very valuable for benchmarking new and existing methods for non-covalent complexes, as it contains 22 dimers of various types (H-bonded, dispersion dominated and mixed), and system sizes (from water dimer to adenine-thymine complexes), which has led to a rise of its popularity [START_REF] Hill | Spin-Component Scaling Methods for Weak and Stacking Interactions[END_REF][START_REF] Riley | Assessment of the MP2 Method, along with Several Basis Sets, for the Computation of Interaction Energies of Biologically Relevant Hydrogen Bonded and Dispersion Bound Complexes[END_REF][START_REF] Mcnamara | Semi-empirical molecular orbital methods including dispersion corrections for the accurate prediction of the full range of intermolecular interactions in biomolecules[END_REF][START_REF] Morgado | Can the DFT-D method describe the full range of noncovalent interactions found in large biomolecules?[END_REF][START_REF] Sato | Long-range corrected density functional study on weakly bound systems: Balanced descriptions of various types of molecular interactions[END_REF][START_REF] Marchetti | Accurate calculations of intermolecular interaction energies using explicitly correlated wave functions[END_REF]48,11,[START_REF] Marchetti | Accurate Calculations of Intermolecular Interaction Energies Using Explicitly Correlated Coupled Cluster Wave Functions and a Dispersion-Weighted MP2 Method[END_REF][START_REF] Gráfová | Comparative Study of Selected Wave Function and Density Functional Methods for Noncovalent Interaction Energy Calculations Using the Extended S22 Data Set[END_REF][START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF][START_REF] Hesselmann | Comparison of Intermolecular Interaction Energies from SAPT and DFT Including Empirical Dispersion Contributions[END_REF]34]. It has even been used to parametrize and adjust certain methods so that the relative errors with respect to the S22 values are minimized [START_REF] Schwabe | Double-hybrid density functionals with long-range dispersion corrections: higher accuracy and extended applicability[END_REF]27,[START_REF] Pito Ňák | Scaled MP3 Non-Covalent Interaction Energies Agree Closely with Accurate CCSD(T) Benchmark Data[END_REF]. The kinetic energy cut-off was set to 80 Ry to achieve well-converged PBE orbitals. As with the A24 test set, the number of auxiliary basis vectors was set to forty times the number of valence states to ensure convergence, except for the formic acid (N aux = 60

× N v ), formamide (N aux = 60 × N v ), phenol (N aux = 32 × N v ), uracil (N aux = 32 × N v ) and adenine-thymine dimers (N aux = 32 × N v )
. This ensured that the variations of the RPA and beyond-RPA interaction energies were within 0.1 kcal/mol. The results are tabulated in Table 4. As was observed with the A24 test set, the inclusion of exchange effects generally improve the RPA interaction energy values, with the mean absolute errors for the eh-TDHF and AC-SOSEX methods around 0.43 kcal/mol and 0.42 kcal/mol respectively, compared to 0.81 kcal/mol for the RPA. Specifically, the improvement is more clear in the hydrogen-bonded systems, where the mean absolute errors for the beyond-RPA methods are 0.48 kcal/mol for the eh-TDHF and 0.36 kcal/mol for the AC-SOSEX approximations, compared to the value of 1.26 kcal/mol for the RPA, and in the mixed-type dimers, where the mean absolute errors for eh-TDHF and AC-SOSEX methods is 0.20 kcal/mol and 0.21 kcal/mol respectively, compared to 0.46 kcal/mol for the RPA. In case of dispersion-dominated systems, as was seen with the A24 test set, the improvement is more limited, where the mean absolute errors for the eh-TDHF and AC-SOSEX methods are 0.60 kcal/mol and 0.65 kcal/mol respectively compared to the value of 0.72 kcal/mol for the RPA. A similar behavior was previously noticed for the RPA+SOSEX method within localized basis sets implementations [START_REF] Ren | Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks[END_REF]. Similar to the A24 test set, we observe that all methods have positive mean errors, indicating a general tendency of these methods to consistently underestimate the absolute binding energies when compared with CCSD(T) interaction energy values. We observe that the beyond-RPA methods significantly improve over the MP2 values reported by Takatani et al. [START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF], whose MAE is 0.88 kcal/mol. Overall, we obtain a good agreement of our RPA results with earlier works, with a mean error of 0.81 kcal/mol compared to 0.9 kcal/mol obtained by Ren et al. [START_REF] Ren | Resolution-ofidentity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF][START_REF] Ren | Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks[END_REF] and 0.79 kcal/mol obtained by Eshuis and Furche [34]. The AC-SOSEX values have a mean absolute percentage error of 9.5% in close agreement with the value of 10.5% reported for a localized basis set implementation [START_REF] Ren | Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks[END_REF]. We highlight that our implementation was done within the plane-wave basis set, which can treat molecules and solids on the same footing, and does not use extrapolation techniques to obtain complete basis set-quality results (see Chapter 5). Table 4.3: Interaction energies (in kcal/mol) for the S22 test set as obtained with the RPA, eh-TDHF, and AC-SOSEX methods compared with CCSD(T) reference values [START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF]. The mean absolute errors (MAE), mean errors (ME) and the maximum deviations (|MAX|) are provided with respect to CCSD(T) values [START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF] for the three types of molecular dimers and for the total test set.

System CCSD(T) [START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF] CCSD(T) [START_REF] Jure Čka | Benchmark database of accurate (mp2 and ccsd (t) complete basis set limit) interaction energies of small model complexes, dna base pairs, and amino acid pairs[END_REF] RPA eh-TDHF AC-SOSEX

Binding curves

Up to this point we have only considered weakly bound dimers in their equilibrium geometry. In this subsection we consider the application of RPA, AC- SOSEX, and eh-TDHF to compute the full energy curve of three weakly bound dimers from the S22 test set [START_REF] Jure Čka | Benchmark database of accurate (mp2 and ccsd (t) complete basis set limit) interaction energies of small model complexes, dna base pairs, and amino acid pairs[END_REF]: The ethene dimer, which is categorized as a dispersion-dominated system, the water dimer, which is hydrogen bonded, and the ethene-ethyne dimer, which belongs to the mixed-type complexes. This allows us to estimate the accuracy of the different approximations also for equilibrium positions. In PW calculations a kinetic-energy cut-off of 80 Ry was used and supercell sizes were chosen to guarantee a distance of at least a 10 Å between periodically repeated images. The computed curves for the three dimers are shown in Figs. 4.4-4.6. In the case of the ethene and ethene-ethyne dimers CCSD(T) curves from Ref. 65 are shown as a reference; in the case of the water dimer we included only the CCSD(T) value at equilibrium. A summary of the equilibrium binding energies is also given in Table 4.4. While the RPA reproduces already at least 80% of the binding energy, the methods with exchange systematically improve this result. For example, eh-TDHF gives 92% of the binding energy of the ethene dimer, 99% for the ethene-ethyne dimer, and 93% for the water dimer. A similar result is obtained for the AC-SOSEX method, with a slight overestimation of the binding energy of the ethene-ethyne dimer. The RPA and beyond-RPA methods also perform well in reproducing the minimum of the binding curve, which is exactly reproduced within the 0.1 Å accuracy of the grid used to discretize the intermolecular distance. This is not the case for PBE, that for the ethene dimer predicts a too large binding distance. With exception of the water dimer, whose interaction is dominated by polar forces, PBE also significantly underestimates binding energies.

In Table 4.4 we also compare the results obtained with the present PW implementation and the Gaussian basis set implementation that we introduced in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]. This is also useful to understand the effect of pseudopotentials in PW results. The PBE and HF binding energies were computed by using the aug-cc-pCV5Z basis set; the RPA and beyond-RPA correlation energies were extrapolated to the complete basis set (CBS) limit by using aug-cc-pCVQZ and aug-cc-pCV5Z basis sets. In this case the two different methods lead to close results with differences of the order of about 0.1 kcal/mol at most. This is not surprising, since weak interactions only marginally involve core regions and the effect of pseudopotentials is expected to be weak. While comparing the two different implementations of the RPA, AC-SOSEX, and eh-TDHF methods it is also important to keep into account that the PW results are converged within 0.03-0.04 kcal/mol with respect to the auxiliary basis set and Gaussian basis set results might be affected by errors due to the CBS extrapolation and basis set superimposition. Accordingly, the agreement between the two different implementations can be considered fully satisfactory in the case of the weakly bound dimers.

Application to reaction energies involving covalent bonds

In this section we present the results for the reaction energy test set proposed by Hesselmann [START_REF] Hesselmann | Random-phase-approximation correlation method including exchange interactions[END_REF]. These reactions involve the breaking and formation of covalent bonds. Results obtained within the plane-wave framework will also be compared to values obtained using a previous localized basis set implementation that we developed [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF].

In order to decrease the spurious interaction between periodically repeated images, each molecule from the reaction test set has been placed in a 26 × 26 × 26 a 3 0 supercell. Norm-conserving pseudopotentials were taken from the library of Schlipf and Gygi [START_REF] Schlipf | Optimization algorithm for the generation of oncv pseudopotentials[END_REF], which has been generated within the optimized normconserving Vanderbilt pseudopotential scheme recently introduced by Hamann [START_REF] Hamann | Optimized norm-conserving vanderbilt pseudopotentials[END_REF]. A kinetic-energy cut-off of 80 Ry was used to expand the orbitals.

In the columns 2-6 of Table 4.5 we show the results obtained within the PW implementation. The last three lines of this table contain the mean absolute error (MAE), the mean error (ME), and the absolute value of the largest deviation (|MAX|) in kcal/mol computed with respect to reference localized orbital CCSD(T) results [START_REF] Hesselmann | Random-phase-approximation correlation method including exchange interactions[END_REF]; these results are also pictorially represented in Fig. 4.7. Similarly to the results obtained with the Gaussian basis set of Ref. 99, the RPA significantly improves over PBE (the MAE is decreased from 4.49 kcal/mol to 3.27 kcal/mol) and this result is improved further by including the exchange contribution; in particular, the eh-TDHF is the most accurate approach with a MAE of only 1.19 kcal/mol. This is particularly promising in the context of plane-wave implementations, where access to the most accurate quantum chemical approaches are numerically challenging [15] and applications rely mostly on semi-local DFT.

We now discuss the comparison of the results obtained within the present PW implementation and the localized basis set implementation in Molpro [START_REF] Werner | version 2012.1, a package of ab initio programs[END_REF][START_REF] Werner | Molpro: a general-purpose quantum chemistry program package[END_REF][START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]. For a reliable comparison with the plane-wave basis set, we used localized basis sets larger than in Ref. [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]. Specifically, here the PBE and HF from a PBE starting point (HF@PBE) total energies have been computed by using the aug-cc-pCV6Z basis set. Since HF@PBE results are used to compute RPA and beyond-RPA total energies by adding the correlation term, their inclusion in the discussion is useful to better understand the origin of differences between PW and Gaussian results. The RPA and beyond-RPA correlation energies were extrapolated to the complete basis set (CBS) limit by the usual 1/X 3 formula [START_REF] Kutzelnigg | Rates of convergence of the partial-wave expansions of atomic correlation energies[END_REF] considering augcc-pCVXZ basis sets with X going from Q to 6. In the columns 7-11 of Table 4.5 and in Fig. 4.8 we present the reaction energy values for the localized basis set implementation. The general trend of the errors is similar to the PW case. However, a sizeable difference in the MAE can be seen in the RPA case (which decreases from 3.27 kcal/mol for PWs to 2.42 kcal/mol for the localized basis sets) and differences are present also in the ME and |MAX|. This issue has to be ascribed to the use of norm-conserving pseudopotentials in the plane-wave basis set. Indeed, Hesselmann's test set involves strong covalent interactions and short bond lengths whose energy might be affected by the core regions of pseudopotentials. As shown in Sec. 4.1, in the case of weakly bound systems the interaction energy is mostly determined by long range forces far from the core region and pseudopotentials weakly affect the results.

To better understand the PW-Gaussian basis set comparison in Fig. 4.9 we present the MAE, ME, and |MAX| of the PW basis set results with respect to the corresponding localized basis set results. The first important observation is that errors increase going from PBE, to HF@PBE, and to ACFDT methods. This is not surprising, since the pseudopotentials have been optimized to provide accurate PBE results [START_REF] Schlipf | Optimization algorithm for the generation of oncv pseudopotentials[END_REF]. Another important point is that the MEs for HF@PBE and for ACFDT methods are positive and close in value to the corresponding MAE values; this indicates that the error induced by the pseudopotentials is rather systematic and tends to underestimate the absolute value of reaction energies. This explains why only in the RPA case the PW MAE with respect to CCSD(T) sizeably increases with respect to Gaussian basis set results: Since the RPA approach already presents a systematic underestimation of the absolute value of reaction energies, the error induced by the pseudopotential accumulates over this intrinsic error of the RPA. All the other approaches (PBE, HF@PBE, AC-SOSEX, and eh-TDHF) have more erratic behaviors: The significant difference between the MAE and ME in Figs. 4.7-4.8 implies that these methods sometimes overestimate and sometimes underestimate reaction energies. Accordingly, the effect of the pseudopotential is mostly compensated when computing the MAE. While the effect of pseudopotentials in the PW implementation is relatively large with respect to the level of accuracy of the AC-SOSEX and eh-TDHF methods, by considering MAE and |MAX| in Fig. 4.7 it is possible to clearly appreciate the improvement of these approaches with respect to the RPA.

The pseudopotential approximation affects HF@PBE and ACFDT results mainly in two ways: (1) The lack of valence-core contributions; (2) The use of pseudoorbitals to compute the HF exchange and the response functions. However, it is important to mention that also the accuracy of correlation energies computed within a Gaussian basis set might be affected to a certain extent by the CBS extrapolation.

In principle, it might be possible to improve the comparison of PW basis set results with Gaussian basis set results by using the projector augmented-wave method [13]. However, within this framework the formalism is significantly more involved and approximations are typically introduced to compute response functions.

Applications to solids

In this section we present applications of the RPA within the implementation of Chapter 3 to compute the structural properties of covalent solids (C in the dia- [START_REF] Hesselmann | Random-phase-approximation correlation method including exchange interactions[END_REF]. Errors have been computed for PBE, HF@PBE, RPA, AC-SOSEX, and eh-TDHF results obtained within the plane-wave (PW) basis set implementation with respect to values obtained within the Gaussian basis set implementation of the same methods. mond structure, silicon, and silicon carbide) and of weakly bonded solids (Ne, Ar, and Kr). While this set of bulk materials is not fully comprehensive, it provides a test-bed to prove accuracy and efficiency of the methodology discussed in this work. The discussion will be limited to the RPA, since an implementation of eh-TDHF and AC-SOSEX for solids is currently not yet available. The calculations were performed with norm-conserving pseudopotentials from the Quantum Espresso website [47,46]. Both LDA [START_REF] Perdew | Self-interaction correction to densityfunctional approximations for many-electron systems[END_REF] and GGA [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] have been considered as starting points for the RPA calculations (these approximations will be denoted as RPA@LDA and RPA@GGA). The convergence with respect to the kinetic energy cut-off and the k-point and q-point grids have been tested carefully 1 . All the total energy curves were obtained by interpolating the total energy computed for at least seven different lattice parameters (up to eleven for Ne) with the Murnaghan equation of state (EOS) [START_REF] Murnaghan | The compressibility of media under extreme pressures[END_REF]. For all the materials studied it was found that N aux = 250 is sufficient to obtain accurate results, with the exception of solid Ne, for which the size of the auxiliary basis set was increased up to 480. The use of the auxiliary basis set was instrumental to reach convergence without extrapolation techniques (see also Chapter 5). For example, the dimensionality reduction is particularly important for noble gas solids, with N P W /N aux up to 150 for neon and up to 120 for argon and krypton.

In Tables 4.6 and 4.7 we show the results for lattice constants and bulk moduli, respectively, for the six bulk materials considered in this work. For comparison purposes we also provide experimental (see Refs. 139, 129 and references therein) and other results based on the RPA in the literature [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Harl | Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids[END_REF][START_REF] Lu | Evaluation of model exchange-correlation kernels in the adiabatic connection fluctuation-dissipation theorem for inhomogeneous systems[END_REF].

In general, compared to experimental values, the RPA preserves the good level of accuracy of the LDA and GGA for covalently bonded systems [START_REF] Harl | Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids[END_REF]. In the case of silicon the RPA@LDA lattice constant is slightly less accurate than the LDA lattice constant (the 0.6% underestimation by LDA increases up to 1.5% for RPA@LDA). Previous RPA and RPA+ calculations based on pseudopotentials had already shown a similar trend of RPA@LDA to underestimate the lattice constant of bulk silicon [START_REF] Nguyen | Efficient calculation of exact exchange and rpa correlation energies in the adiabatic-connection fluctuation-dissipation theory[END_REF]43,[START_REF] Lu | Evaluation of model exchange-correlation kernels in the adiabatic connection fluctuation-dissipation theorem for inhomogeneous systems[END_REF]; in particular, as shown in Table 4.6, the RPA result is in excellent agreement with our value [START_REF] Lu | Evaluation of model exchange-correlation kernels in the adiabatic connection fluctuation-dissipation theorem for inhomogeneous systems[END_REF]. In Ref. 60 a RPA@LDA lattice constant was obtained in perfect agreement with experiments; the discrepancy with respect to our results might be due to the details of the implementation, 1 Converged parameters: RPA@LDA for C: 70 Ry (kinetic energy cut-off for wavefunctions), 6 × 6 × 6 (k-point grid), 4 × 4 × 4 (q-point grid shifted from the origin); RPA@GGA for C: 70 since in this previous calculation a code based on the projector augmented wave (PAW) method was used and the RPA correlation energy was obtained by extrapolation techniques. When GGA is used as a starting point for the RPA calculation (RPA@GGA), the lattice constant of silicon is improved with respect to both LDA and GGA. As shown in Table 4.7, the RPA systematically improves the bulk modulus of covalent solids with respect to traditional (semi-)local functionals. The case of noble gas solids is certainly more interesting. In this case, the interaction is dominated by weak dispersion forces and the RPA strongly improves the results of LDA and GGA [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF]. In general, as shown in Table 4.6, the LDA importantly underestimates experimental values of the lattice constant (LDA relative error with respect to experiments: -11.6% for Ne, -5.7% for Ar, and -4.8% for Kr) while the GGA overestimates them (GGA relative error: +5.6% for Ne, +13.6% for Ar, and +14.3 % for Kr). Since the RPA accurately describes weak dispersion forces, the results for lattice constants are significantly improved both using LDA and GGA as a starting point (RPA@LDA relative error: +9.1% for Ne, +3.5% for Ar, and +2.2% for Kr; RPA@GGA relative error: -0.5% for Ne, +1.2% for Ar, and +1.1% for Kr). As shown in Table 4.7, the RPA strongly improves also the bulk moduli of noble gas solids. For example, LDA strongly overestimates the bulk moduli (LDA relative error: +657% for Ne, +182% for Ar, and +70% for Kr) while overall RPA@GGA gives results much closer to experiments (RPA@GGA relative error: +35% for Ne, +17% for Ar, and -17% for Kr). In general, while LDA and GGA lead to very different structural properties, this behavior is less pronounced for the RPA, that gives much closer results to experiment starting from LDA or GGA states. For comparison purposes, in Tables 4.6 and 4.7 we included results from previous RPA calculations available in the literature [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Harl | Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids[END_REF][START_REF] Lu | Evaluation of model exchange-correlation kernels in the adiabatic connection fluctuation-dissipation theorem for inhomogeneous systems[END_REF]. The overall agreement of our results with previous work is good. However, a straightforward comparison is not possible since the majority of previous results were obtained with an implementation based on PAW and extrapolation techniques [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Harl | Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids[END_REF]. Furthermore, bulk moduli of noble gas solids have not been computed before at the RPA level of theory.

By considering the examples shown in this section we can conclude, in agreement with previous work [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Harl | Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids[END_REF], that the RPA preserves the accuracy of (semi-)local functionals for covalently bonded solids, such as C, Si, and SiC, and strongly improves the structural properties of weakly bonded noble gas solids. In particular, the RPA approximated starting from GGA orbitals has been shown to have an impressive accuracy for the structural properties of the systems considered in this work.

Chapter 5 Auxiliary basis set convergence

This chapter is based on original work previously published in Ref. [29].

In this chapter we discuss in detail a point that was only briefly addressed in Chapter 3, the technique used to achieve basis set convergence of correlation energies. Differently from semi-local or hybrid functionals, the ACFDT approaches are characterized by a slow convergence with respect to this parameter. This behavior is analogous to that of traditional quantum chemical methods. Indeed, convergence of the correlation energy is known to be slow in correlated methods such as Møller-Plesset perturbation theory or coupled cluster theory [START_REF] Schwartz | Estimating Convergence Rates of Variational Calculations[END_REF][START_REF] Helgaker | Basis-set convergence of correlated calculations on water[END_REF][START_REF] Halkier | Basis-set convergence in correlated calculations on Ne, N 2 and H 2 O[END_REF]36,[START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF][START_REF] Řezá Č | Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications[END_REF][START_REF] Møller | Note on an approximation treatment for many-electron systems[END_REF]4]. In order to overcome this issue, techniques to extrapolate the results to the complete basis set (CBS) limit are routinely used [START_REF] Kutzelnigg | Rates of convergence of the partial-wave expansions of atomic correlation energies[END_REF][START_REF] Halkier | Basis-set convergence in correlated calculations on Ne, N 2 and H 2 O[END_REF][START_REF] Řezá Č | Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications[END_REF].

In Chapter 4 we did not rely on any complete basis set extrapolation technique for the the application of RPA or beyond-RPA methods, but we rather converged energy differences. This is a peculiar feature of our implementation that takes advantage of a compact auxiliary basis set which is designed to rapidly converge the eigenvalues of the response functions [START_REF] Wilson | Efficient iterative method for calculations of dielectric matrices[END_REF][START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF]. In this chapter we give a detailed description of the methodology that we have been using so far, and compare it to a more traditional approach, the complete basis set extrapolation.

Convergence of energy differences vs. complete basis set extrapolation

To understand how correlation energy differences are converged, a simple observation is necessary, illustrated by a simple example. Let us suppose that we are interested in computing the interaction energy of the argon dimer at the RPA level of theory [START_REF] Nguyen | A first-principles study of weakly bound molecules using exact exchange and the random phase approximation[END_REF][START_REF] Lu | Power series expansion of the random phase approximation correlation energy: The role of the third-and higher-order contributions[END_REF]. We first consider the two atoms at a very large distance so that

E int = E dimer -2E atom ≈ 0, (5.1)
where E int is the interaction energy, E dimer is the total dimer energy, and E atom is the total atom energy. Eq. 5.1 is a consequence of the size consistency of the RPA [START_REF] Gould | Casimir-Polder Size Consistency: A Constraint Violated by Some Dispersion Theories[END_REF]. Clearly, each single contribution to the energy should be exactly compensated for, including the correlation energy:

E int c = E dimer c -2E atom c
≈ 0. By using a single eigenvalue of the polarizability for each frequency λ 1 (iω) (the largest in absolute value), the RPA correlation energy (Eq. 3.10) of each single argon atom is

E atom c = 1 2 ∞ -∞ dω 2π [ln {1 -λ 1 (iω)} + λ 1 (iω)] . (5.2)
In order to obtain E int ≈ 0 it is now clear that the correlation energy of the dimer should be computed including two eigenvalues per frequency in the correlation energy expression:

E dimer c = 1 2 ∞ -∞ dω 2π [ln {1 -λ 1 (iω)} + λ 1 (iω) + ln {1 -λ 2 (iω)} + λ 2 (iω)] .
(5.3) Indeed, since the interatomic distance is large, we have λ 2 (iω) ≈ λ 1 (iω). This error cancellation would have not been possible if we would have naively used the same number of eigenvalues for the atoms and the dimer. The same idea can be also applied to small interatomic distances; there should be a compensation between the eigenvalue contributions of the dimer and the monomers to induce error cancellation. Generally, at small or intermediate distances it is necessary to significantly increase the number of eigenvalues N eig included in the evaluation of the correlation energy to reach convergence, while maintaining the condition N dimer eig = 2N monomer eig . Henceforth, we refer to this approach as the size-consistent energy difference approach.

This idea can also be generalized to systems, whose components are different (e.g. a dimer of two different molecules), or which involve different types of interactions. When the PDEP technique was first introduced it was already understood that the number of eigenvalues necessary to reach a certain level of accuracy in the description of the response function is approximately proportional to the number of electrons [START_REF] Wilson | Efficient iterative method for calculations of dielectric matrices[END_REF]. Accordingly, the number of eigenvalues (namely of auxiliary basis vectors) used to describe the response function should be chosen to be proportional to the number of electrons. For example, if in a certain molecular dimer, the monomer 1 has double the number of electrons than the monomer 2, in the . As shown by the numerical results in Sec. 5.2 this approach works well in practice.

All this discussion has been presented for the RPA. The eh-TDHF and AC-SOSEX correlation energies have a similar structure to the RPA (see Eqs. 3.9 and 3.12), but contain additional exchange terms. As shown in Eqs. 3.6-3.8, the additional exchange contributions have a structure similar to the response function in Eq. 3.5, the main difference being related to the v = v terms. It was already observed that the PDEP basis set, specifically designed for the polarizability (Eq. 3.5), represents well also the exchange terms (Eqs. 3.6-3.8) [28,31]. Because of these affinities between the RPA and the beyond-RPA methods discussed in this work, the same convergence techniques can be used in practical calculations for beyond-RPA methods as well.

By using the size-consistent energy difference technique, the interaction energy converges quite rapidly with respect to the auxiliary basis set size. As an example, we present in Fig. 5.1 the interaction energy of the formaldehyde dimer computed at the RPA and beyond-RPA levels of theory. The convergence with respect to N aux is quite rapid, with the beyond-RPA energies seeming to converge faster than RPA. Since the measurable physical observables are energy differences, convergence of the interaction energy with respect to the auxiliary basis set size is sufficient for practical applications. Empirically we have observed that N aux = 40 × N v is sufficient to achieve convergence for most molecular systems, where N v is the number of Kohn-Sham states (half the number of valence electrons) in the system. This is a very small number compared to the number of plane-waves (N P W = 163000 for the formaldehyde dimer), and convergence of the interaction energy can be achieved without resorting to extrapolation techniques. However, if we consider the correlation energies of the formaldehyde dimer and of the monomer separately, convergence is far from achieved. For example, by increasing the basis set from 35 × N v to 70 × N v , the correlation energy decreases by about 10 kcal/mol for the dimer and by about 5 kcal/mol for the monomers; this leads to the error cancellation shown in Fig. 5.1.

It is important to mention that the size-consistent energy difference technique can in principle be used with different types of auxiliary basis sets. For example, plane-waves are traditionally used to represent dielectric matrices. For a given size of the dielectric matrix, the accuracy of the eigenvalues in the PW representation is lower than for the kinetic-PDEP auxiliary basis sets and, accordingly, also the energy differences less accurate. This implies that the use of PWs or other basis set, while compatible with the energy difference method, might lead to slower convergence.

The approach described in this section to converge energy differences has been previously used for homogeneous [START_REF] Nguyen | A first-principles study of weakly bound molecules using exact exchange and the random phase approximation[END_REF][START_REF] Lu | Power series expansion of the random phase approximation correlation energy: The role of the third-and higher-order contributions[END_REF][START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF] and heterogeneous systems [28, 30,31]. Here we will systematically study its accuracy for weakly bound systems and reaction energies involving the breaking and formation of covalent bonds (see Sec. 5.2). The complete basis set extrapolation scheme for the RPA and beyond-RPA correlation energies used here is,

E Naux c = A + B N aux , (5.4)
where A = E ∞ c is the CBS correlation energy and B is a system-dependent parameter. As shown in Fig. 5.2 for the formaldehyde dimer, the function in Eq. 5.4 provides a reasonable description of the behavior of the correlation energy as a function of N aux . An analogous expression with an additional C/N 2 aux term was previously used to compute the interaction energy of molecular crystals [START_REF] Lu | Ab initio calculation of van der waals bonded molecular crystals[END_REF]. The higher order contribution can be discarded for sufficiently high values of N aux (this is the case for the results of Sec. 5.2). In Ref. [START_REF] Lu | Ab initio calculation of van der waals bonded molecular crystals[END_REF], the PDEP method was used to exactly diagonalize the dielectric matrix for each different value of the frequency. By comparison our approach is approximated, as the basis set generated in the static case is used in the dynamical case as well.

Beyond working well in practice, Eq. 5.4 is similar to the function routinely used in the quantum chemistry community to extrapolate the correlation energy of approximations based on Møller-Plesset perturbation theory and coupled cluster theory. In these approaches the correlation energy is extrapolated to the infinite

N aux = 21×N v to N aux = 70×N v .
The continuous line denotes a linear fit of the data. The dotted and dot-dash lines show partial fits, demonstrating how rapidly the approach converges.

basis set limit as [START_REF] Schwartz | Estimating Convergence Rates of Variational Calculations[END_REF][START_REF] Helgaker | Basis-set convergence of correlated calculations on water[END_REF][START_REF] Halkier | Basis-set convergence in correlated calculations on Ne, N 2 and H 2 O[END_REF]36],

E X c = E ∞ c + B X 3 , (5.5)
with X denoting the cardinal number of the basis set, E ∞ c being its "true" correlation energy evaluated at the complete basis set limit, and B being a system-, and method-dependent constant. The number of the basis set functions is typically

N X = 1 3 (X + 1) X + 3 2 (X + 2) , (5.6)
namely N X ∝ X 3 , and for sufficiently large values of N X we get to leading order

E X c = E ∞ c + B N X , (5.7)
which has a form similar to Eq. 5.4. As with B above, B depends on the system and method.

Alternative extrapolation techniques used in the quantum chemistry community include expressions with variable exponents such as [START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF][START_REF] Řezá Č | Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications[END_REF],

E X c = E ∞ c + B X C , (5.8)
where, typically, small deviations are observed from C = 3 (Eq. 5.5).

In standard implementations using plane-wave basis sets, the RPA correlation energy is often extrapolated to CBS as [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF],

E |G|cut c = E ∞ c + B |G cut | 3 = E ∞ c + B N P W , (5.9)
where |G cut | is the maximum reciprocal lattice vector used to represent the response function and N P W is the number of plane-waves. This is different with respect to Eq. 5.7, where N X is the basis set used to describe the full electronic structure and not just the dielectric matrices (which are not even used in typical coupled-cluster or Møller-Plesset perturbation theory formulations). Our implementation is also based on a plane-wave basis set but the response functions are represented on a system-dependent optimized basis set. Because of this difference, the range of values of N aux used to apply the extrapolation function in Eq. 5.4 might be different from the range of values of N P W used for the extrapolation function in Eq. 5.9. Since the basic numerical framework is the same, the two implementations should give indistinguishable results in the CBS limit.

To complete the literature review, we mention that Björkman et al. [12] and Gulans [START_REF] Gulans | Towards numerically accurate many-body perturbation theory: Short-range correlation effects[END_REF] used a different approach to converge RPA results for layered materials like hexagonal boron nitride and molybdenum di-sulphide. In these works it was observed that the correlation energy differences of layered materials converge as |G cut | -5 , and this observation was used to obtain CBS results. By using this approach it was possible to obtain accurate CBS interaction energies with much lower kinetic energy cut-offs (∼ 100-150 eV). However, it is not clear if this methodology is valid in general for different classes of materials, or is a specific property of spatially infinite planar solids. Though attractive, such behavior of correlation energy differences was not observed in our calculations.

Both approaches based on the CBS extrapolation and on the size-consistent energy differences present advantages and disadvantages. In the case of the CBS extrapolation, the choice of the model is to a certain extent arbitrary (for example, Eq. 5.5 or Eq. 5.8); additionally, the final result might be affected by small fluctuations in the correlation energies used for the fitting. On the plus side this approach can be applied to any kind of system. By contrast, the alternative technique proposed here is not biased by the model chosen, but is based on reasonable assumptions about size-consistency. However, this approach might suffer from slow convergence when energy differences are computed between systems that have a very different dielectric eigenvalue structure.

Let us consider an extreme example: the interaction between a small molecule (monomer1), with only one occupied state, and a larger molecule (monomer2), with nine occupied states. Specifically, we suppose that λ monomer2 400

(iω) > λ monomer1 1 (iω),
where λ i (iω) indicates the i th dielectric eigenvalue for each system. If we compute the RPA correlation energy for the corresponding dimer when the two molecules are separated by a large distance, by using just 400 eigenvalues (40 × N v ), we would have,

E dimer c = 1 2 ∞ -∞ dω 2π 400 i ln 1 -λ dimer i (iω) + λ dimer i (iω) = 1 2 ∞ -∞ dω 2π 400 i ln 1 -λ monomer1 i (iω) + λ monomer1 i (iω) ;
(5.10)

it is clear that in this case, the energy of the dimer is fully dominated by the eigenvalues of the larger system. In other words, there are no eigenvalues of monomer2 present in the eigenvalue listing of the dimer, and no error cancellation can be obtained for N aux = 40 × N v . Accordingly, convergence can be achieved only for very large values of N aux . This slow convergence could also appear in cases where the correlation energy difference is computed between two systems and a corresponding complex containing the same atoms rearranged in a significantly different structure with a significantly different dielectric eigenvalues. In this context the reaction energy test set discussed in Sec. 5.2.C could represent a problematic example, as it contains reactants and products with different molecular structures and possibly significantly different dielectric spectrum; still the size-consistent energy difference method provides accurate estimates of the reaction energies. In general, the size-consistent energy difference approach allows for a systematic study of convergence as a function of N aux , and problematic cases can be identified and potentially dealt with, a significant advantage over the the CBS technique.

In the following section, we present the RPA and beyond-RPA interaction energy results for the A24 test set [START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF], for the methane-formaldehyde complex (MFC), and for reaction energies as computed for different auxiliary basis set sizes (N aux ) with and without the extrapolation technique of Eq. 5.4. The A24 and reaction energy test sets were already considered in Chapter 4; here the discussion will be limited exclusively to the convergence properties with respect to the auxiliary basis set.

Numerical examples

A24 Test Set

To verify the energy differences, and compare the values with those obtained with the CBS extrapolation technique, we start by considering the A24 test set proposed by Řezáč et al., which provides interaction energies for 24 molecular dimers of three types -hydrogen-bonded, dispersion-dominated and mixed-type [START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF]. A comparison of the RPA and beyond-RPA methods with the reference CCSDT(Q) values has been discussed in Chapter 4, and is not repeated here. Here we focus on estimating the inaccuracies due to the use of an incomplete auxiliary basis set; as the molecules included in the A24 are relatively small, a systematic study is possible by increasing N aux to large values. For establishing a reference, the RPA and beyond-RPA calculations were first performed using an auxiliary basis set size of 70 × N v (in Chapter 4 this number was limited to 40 × N v for most systems) and then the absolute values of the correlation energies were extrapolated to CBS using Eq. 5.4 with 30 points from 41 × N v to 70 × N v . Although small numerical fluctuations might still be present, this procedure, by including a large number of points in the fitting and a large N aux , provides reasonably accurate reference values to study the convergence of RPA and beyond-RPA methods.

The choice of the 41N v -70N v CBS value as a reference will also be used for the other numerical examples considered in the following subsections. As we will consider below values from size-consistent energy differences and CBS extrapolations from different intervals of N aux , it is important to keep in mind that our specific choice of the reference introduces a bias in favor of CBS results, which in general will be characterized by a more regular and predictable convergence.

In the following we study the level of convergence of CBS values obtained by extrapolating 20 points in three different ranges: N aux /N v = 11-30, 21-40, and 31-50. We compare these results with size-consistent energy differences computed at N aux /N v = 30, 40, and 50 by using the procedure discussed in Sec. 5.1. For convenience, we denote with E

CBS(n,m) A

the CBS interaction energy computed by extrapolating 20 correlation energies between N aux = n × N v and m × N v at the A level of theory (A = RPA, eh-TDHF, or AC-SOSEX); more specifically, Eq. 5.4 was used to obtain CBS correlation energies for each system (dimer and monomer) and these values were then used to compute correlation energies. Similarly, E Naux=nNv A corresponds to the non-extrapolated (namely based on sizeconsistent energy differences) interaction energy at the A level of theory obtained with

N aux = n × N v .
The top panel of Fig.

green) for the A24 test set interaction energies using CBS extrapolations for E CBS (11,30) , E CBS (21,40) , and E CBS (31,50) , and size-consistent energy differences for E Naux=30 , E Naux=40 , and E Naux=50 . The E CBS (41,70) interaction energy was used as a reference (see text). from energy differences for three different values of N aux , E Naux=30

RPA ,E Naux=40 RPA , and E Naux=50 RPA . In the middle and lower panels of Fig. 5.3, analogous results are shown for the eh-TDHF and AC-SOSEX methods respectively. From the three panels, we observe that the MAEs are rather small in all cases. However, considering the small energy scales involved in the A24 test set (the average interaction energy at the CCSDT(Q) level of theory is about 2 kcal/mol [START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF]), a high level of convergence is also particularly important. By considering Fig. 5.3 we can notice that: (1) RPA presents larger MAEs than the beyond-RPA methods, hinting that its convergence is slower and/or less stable; this empirically observed behavior might be related to a reduced cusp problem in beyond-RPA approximations, although further investigation will be necessary in the future for a better understanding; (2) if the maximum value for N aux is the same, results based on the CBS extrapolation or on the energy differences have similar accuracy; (3) MAEs systematically decrease by increasing the maximum value of N aux ; (4) the value of N aux = 40 × N v used previously in Chapter 4 to compute interaction energies provides a good level of accuracy, especially for the beyond-RPA methods.

According to this discussion for the A24 test set, the CBS results and the size-consistent energy differences have similar convergence rates and the two approaches are equivalent for most practical calculations. However, the sizeconsistent energy differences have the advantage of not being biased by the choice of a particular function to extrapolate to CBS.

Methane-Formaldehyde Complex

To further study the convergence rate of the RPA and beyond-RPA methods with respect to the auxiliary basis set size, we computed six points on the potential energy surface (PES) of the methane-formaldehyde complex (MFC). For these points accurate coupled-cluster benchmark results have been recently published [START_REF] Moore | O Hydrogen Bonding. The Prototypical Methane-Formaldehyde System: A Critical Assessment[END_REF]. The numerical results for the MFC complex are also based on FHI norm-conserving pseudopotentials. This system allows us to verify the accuracy of the RPA and beyond-RPA methods and the reliability of the convergence techniques discussed in Sec. 5.1 also for the case of out-of-equilibrium geometries. As in the case of the A24 test set, we use E CBS (41,70) energies as the reference and we test the convergence of the CBS extrapolated values for three intervals, E CBS (11,30) , E CBS (21,40) , and E CBS (31,50) , and the energy differences at three values of N aux , E Naux=30Nv RPA , E Naux=40Nv

RPA

, and E Naux=50Nv RPA . In Fig. 5.4 we present the corresponding MAEs for the RPA, the eh-TDHF, and the AC-SOSEX. For the MFC, the MAEs decrease by increasing the value of N aux , as was observed with the A24 test set. In contrast though, the convergence of energy differences significantly outperform CBS extrapolations for the beyond-RPA methods, whereas the converse is true for the RPA. To a certain extent this behavior can be interpreted by considering that size-consistent energy differences benefit from error cancellations, whereas the CBS extrapolations are sensitive to small fluctuations in the data. If for example we consider E CBS(41,60) RP A , which includes a subset of the points used for the reference E CBS(41,70) RP A , the MAE actually increases to 0.07 kcal/mol, showing that small fluctuations influence the CBS extrapolation. On the other hand, the MAE on the energy differences continues to decrease with an increase in N aux , reaching a value of 0.12 kcal/mol for E Naux=60Nv RP A . Thus, the convergence of the RPA is slower and more unstable, similar to that observed for the A24 test set. For the beyond-RPA methods, convergence is instead more systematic, with both approaches leading to similar results for sufficiently high values of N aux . The method based on size-consistent energy differences is thus confirmed as a reliable and unbiased approach for the MFC as well, in particular, for the eh-TDHF and AC-SOSEX methods. Fig. 5.5 shows the CBS results for the MFC PES compared to CCSD(T) reference values [START_REF] Moore | O Hydrogen Bonding. The Prototypical Methane-Formaldehyde System: A Critical Assessment[END_REF]. The RPA and beyond-RPA methods all capture the PES correctly, with the correct ordering in energy of the configurations. With respect to the coupled-cluster reference the mean absolute error for the RPA is 0.085 kcal/mol. This value reduces to 0.061 and 0.066 kcal/mol for the eh-TDHF and AC-SOSEX respectively, highlighting the accurate nature of these methods. However, the convergence of RPA results is more problematic and, by considering the small energy scales involved in this PES, the eh-TDHF and AC-SOSEX results are certainly more reliable.

Reaction energies

After considering weakly bound systems, in this section we discuss results for the reaction energy test set by Hesselmann [START_REF] Hesselmann | Random-phase-approximation correlation method including exchange interactions[END_REF], which involves the breaking and formation of covalent bonds. While results for this test set were previously presented in Chapter 4, here we extend the auxiliary basis set up to 70 × N v and compare extrapolated results with energy differences. Studying this test set is particularly interesting, since going from the reactants (e.g. HCOOH and NH 3 ) to the products (e.g. HCONH 2 and H 2 O) involves the breaking and formation of covalent bonds. Accordingly, reactants and products could have a rather different structure of the dielectric eigenspectrum and this could represent a challenge for the size-consistent energy difference method.

In Fig. 5.6 we present the MAEs obtained at the RPA, the eh-TDHF, and the AC-SOSEX levels of theory. As for the weakly bound systems, the MAEs are computed with respect to the E CBS (41,70) extrapolated value. Not surprisingly for this test set the deviation between CBS and size-consistent energy differences is more sizeable than for the weakly-interacting systems. A few observations are AC-SOSEX (green) for the methane-formaldehyde complex interaction energies using CBS extrapolations for E CBS (11,30) , E CBS (21,40) , and E CBS (31,50) , and size-consistent energy differences for E Naux=30 , E Naux=40 , and E Naux=50 . The E CBS (41,70) interaction energy was used as a reference (see text). important for a full understanding of the results in Fig. 5.6: (1) As the average reaction energy at the CCSDT(Q) level of theory is about 29 kcal/mol, the energy scales involved in this test set are usually larger than for weakly interacting systems and, accordingly, higher errors can be tolerated; (2) As mentioned before, the use of an E CBS (41,70) reference is to some extent arbitrary and creates a bias in favor of the CBS results; by using E Naux=70 as a reference instead, the behavior shown in Fig. 5.6 would be swapped, with the size-consistent energy difference method showing a faster convergence; (3) The MAE of E Naux=70 with respect to E CBS (41,70) (or vice versa) is within 0.09 kcal/mol for all the approximations and this observation point to the fact that the two methods studied in this work converge to the same values.

According to this discussion, the size-consistent energy difference method works well also beyond the expected range of applicability. The convergence of correlation energies for applications involving the breaking and formation of covalent bonds is more problematic but can achieve a good level of reliability.

Conclusions

We analyzed the basis set convergence of the RPA, eh-TDHF, and AC-SOSEX methods. Specifically, two techniques were used to converge correlation energies: one based on a complete basis set extrapolation, and the other that converges energy differences by accounting for the size-consistency property. These two methods were tested by considering the A24 test set and six points on the poten-Figure 5.6: Mean absolute errors (MAEs) of RPA (red), eh-TDHF (blue) and AC-SOSEX (green) for the reaction energy test set of Hesselmann [START_REF] Hesselmann | Random-phase-approximation correlation method including exchange interactions[END_REF] using CBS extrapolations for E CBS (11,30) , E CBS (21,40) , and E CBS (31,50) , and size-consistent energy differences for E Naux=30 , E Naux=40 , and E Naux=50 . The E CBS (41,70) interaction energy was used as a reference (see text). tial energy surface of the methane-formaldehyde complex. It was found that both the CBS extrapolation and size-consistent energy differences converge to similar results at similar speeds as a function of the auxiliary basis set size. It was also shown that the RPA is in general characterized by a slower convergence than the beyond-RPA methods. At least for the weakly bound systems considered in this manuscript, the two approaches have been found to be equivalent for practical purposes. If the breaking and formation of covalent bonds is involved, convergence is less trivial but satisfactory levels of accuracy can be reached.

The size-consistent energy difference technique represents a new and alternative scheme to conventional basis set extrapolation for converging results of correlated methods within a dielectric matrix formulation. While possibly less general than the CBS extrapolation, this approach has the advantage of not relying on any assumption on the model used for the fitting.

Chapter 6 Applications to realistic systems by machine learning thermodynamic perturbation theory

This chapter is based on original work previously published in Refs. [START_REF] Rocca | Bridging molecular dynamics and correlated wave-function methods for accurate finite-temperature properties[END_REF]20].

Up to this point we have discussed new approximations for the electronic correlation that hold the promise for a higher level of accuracy than traditional DFT approximations but involve a significantly high computational cost. Additionally, the implementation of analytical gradients, necessary for geometry optimizations or molecular dynamics, is challenging for these methods and further increases the computational requirements. Accordingly, RPA molecular dynamics are by far too expensive to be routinely used in applications that require higher levels of accuracy. In this chapter we will consider examples of computations of adsorption enthalpies in zeolites. Directly performing a RPA ab initio molecular dynamics (AIMD) for these systems is estimated to require several tens of millions of CPU hours. In order to compute finite-temperature properties at the RPA level of theory in a simple and numerically efficient way, here we discuss the MLPT method, that couples machine learning (ML) techniques with thermodynamic perturbation theory (TPT). The numerical applications of this methodology presented below show that a high level of accuracy can be achieved by training a ML model with as few as 10 RPA energies. The MLPT approach represents an important step towards the possibility of applying the RPA and other high-accuracy/high-cost ab initio methodologies to reach a new level of reliability in finite-temperature materials simulations.

Methodological approach: Machine learning thermodynamic perturbation theory

Since the work of Behler and Parrinello (BP) in 2007 [10] the use of ML algorithms in MD simulations is growing in popularity [24,23,[START_REF] Zhang | Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics[END_REF]9]. In the BP approach and its variants, a certain number of configurations of a system is generated by performing a relatively short MD simulation or by properly sampling the configuration space. A ML learning model, usually a neural network (NN), is trained on the energies and/or forces computed for these configurations. A full molecular dynamics (MD) simulation is then performed by using the previously trained neural network. At each step of the MD simulation, instead of directly solving the Schrödinger equation, forces and energies are inexpensively "predicted" by the ML model and used to generate the configuration in the next step. In this way the timescale accessible by a typical MD run is significantly increased while keeping to a large extent the accuracy of an ab initio simulation. However, the training of a ML learning algorithm for a MD simulation typically requires thousands or even tens of thousands of configurations and, accordingly, the BP approach is almost exclusively trained on DFT energies and/or forces.

Similarly to the BP approach, here we use a ML learning regression model to predict the energy surface from a certain number of training examples. Differently from the BP approach, this model is not used to perform a full MD simulation but it is rather coupled with thermodynamic perturbation theory, whose brute force application also requires an impractically high number of RPA calculations. The corresponding MLPT calculation is executed in the following steps: (1) an AIMD simulation is performed using a numerically affordable DFT functional;

(2) a small number N train of configurations generated by the initial MD is selected and RPA calculations are performed on them; (3) a ML model is trained on the selected configurations to predict the RPA energies for all the remaining configurations; (4) based on TPT, the contributions from the configurations generated at step 1 are reweighted using the predicted energies and RPA ensemble averages are computed. For the specific objective of this work, the present MLPT approach has the following advantages with respect to BP MD: (1) The ML regression algorithm is exclusively trained on energies; this is useful since most condensed phase implementations of correlated methods do not support the calculation of analytical gradients. (2) The implementation is easier, since the BP method requires also the derivative of the ML model with respect to atomic positions in order to compute forces. (3) The accuracy of the ML predictions can be enhanced in a straightforward way by using a ∆-machine learning approach [START_REF] Ramakrishnan | Big data meets quantum chemistry approximations: The δ-machine learning approach[END_REF], which trains the model on energy differences instead of total energies (see Sec. 6.1). In the two following subsections we discuss the technical details of thermodynamic perturbation theory and machine learning that are at the base of the MLPT method.

Thermodynamic perturbation theory

Starting from a molecular dynamics driven by a numerically inexpensive approximation, our goal is to use thermodynamic perturbation theory (TPT) to obtain finite-temperature observables at a high level of theory. This methodology is general but for the specific numerical applications of this manuscript the van der Waals corrected PBE+D2 functional [START_REF] Grimme | Semiempirical gga-type density functional constructed with a long-range dispersion correction[END_REF][START_REF] Perdew | Generalized gradient approximation made simple[END_REF] will be used as a (cheap) starting point to compute RPA ensemble averages.

A general review of TPT in the context of classical MD can be found e.g. in Refs. 22, 114. Here we will discuss the main ideas necessary to understand the MLPT method. Molecular dynamics is used to evaluate the finite-temperature expectation values of a certain observable O:

O H = O(q, p)exp {-βH(q, p)} dq dp exp {-βH(q, p)} dq dp , (6.1)
where . . . H is used to compactly denote the average over the canonical ensemble corresponding to the Hamiltonian H, β = 1/k B T , and q and p denote nuclear positions and momenta, respectively. In Eq. 6.1 the nuclear motion is driven by the classical Hamiltonian

H(q, p) = T (p) + V (q), (6.2)
where T (p) is the nuclear kinetic energy and V (q) is the potential acting on the nuclei. Here we will focus on AIMD, where the potential V (q) is determined by the solution of the electronic Schrödinger equation within different approximations, usually based on DFT:

V (q) = Ψ(q)| ĤDFT |Ψ(q) + U II (q), (6.3)
where U II is the repulsive nuclear energy and ĤDFT is the electronic DFT Hamiltonian (not to be confused with the classical Hamiltonian H). As explained above, our main goal is to use thermodynamic perturbation theory to compute observables within a sophisticated quantum chemical approximation (e.g. the RPA) starting from a molecular dynamics driven by a computationally "cheap" DFT approximation (e.g. PBE+D2). Below the Hamiltonian H and the corresponding potential V will be associated with the less expensive quantum-mechanical approximation while primed quantities (H and V ) will be associated with the more expensive approach. By using the identity H (q, p) = H(q, p) + V (q) -V (q) = H(q, p) + ∆V (q), (6.4) the average of the observable O(q, p) over the ensemble corresponding to H can be written as O H = O(q, p)exp {-βH (q, p)} dq dp exp {-βH (q, p)} dq dp = O(q, p)exp {-β∆V (q)} exp {-βH(q, p)} dq dp exp {-β∆V (q)} exp {-βH(q, p)} dq dp

= Oexp {-β∆V } H exp {-β∆V } H , (6.5)
which provides the main result of TPT used in this discussion. This equation has important theoretical and practical implications, showing that the ensemble average O H can be redefined in terms of an ensemble average corresponding to H; this is achieved through the reweighting factor exp {-β∆V }. From a numerical point of view Eq. 6.5 implies that the statistics sampled by a MD driven by H can be used to obtain finite-temperature expectation values corresponding to H . This is a crucial point since directly performing a MD based on H would be extremely costly if feasible at all. From a practical point of view, in order for the perturbative procedure in Eq. 6.5 to be effective, it is necessary that a MD driven by H would sample a configurational space that is similar to that sampled by an hypothetical H MD. If the corresponding high probability ensembles are disjoint or only partially superimposed, the perturbative approach will have low accuracy. A detailed discussion of this issue can be found in Refs. 22, 114. Following Ref. [START_REF] Rocca | Bridging molecular dynamics and correlated wave-function methods for accurate finite-temperature properties[END_REF], we suppose that PBE+D2 represents a reliable starting point for RPA.

In this manuscript the TPT formalism will be applied to the specific problem of computing adsorption energies in zeolites at the RPA level of theory. The adsorption enthalpy can be computed as

∆ ads H(T ) = V S+A H -[ V S H + V A H ] -k B T, (6.6)
where S denotes the clean substrate (a zeolite in this case) and A the adsorbate; the values of V H can be obtained from Eq. 6.5 for the specific case O = V . In order to evaluate numerically the ensemble averages in practical numerical calculations, a MD based on H is first performed and then the ensemble average is computed as

V H ≈ Ntot i=1 V i exp {-β(V i -V i )} Ntot i=1 exp {-β(V i -V i )} = Ntot i=1 V i exp {-β∆V i } Ntot i=1 exp {-β∆V i } , (6.7)
where the index i runs over all the N tot configurations generated by the production MD run. Since in the MD simulations considered below N tot = 190000, Eq. 6.7 highlights the limited practical applicability of TPT alone. Indeed, even by considering a smaller subset of uncorrelated configurations, several thousands of calculations of V (namely of RPA total energies) are expected to be necessary to reach convergence. This issue can be significantly alleviated by using a ML regression algorithm. Within this approach only few V i 's are explicitly computed and used to train a ML model; all the remaining values up to N tot are then predicted at a negligible computational cost.

Machine Learning approach

In order to significantly decrease the number of expensive calculations (i.e. based on the RPA) required to apply thermodynamic perturbation theory, we propose the use of a ML regression model. To this purpose we chose the smooth overlap of atomic positions (SOAP) [5,26] as descriptor for the geometric configurations and the kernel ridge regression (KRR) [START_REF] Rupp | Machine learning for quantum mechanics in a nutshell[END_REF] as ML algorithm to predict total energies. This choice, that allowed us to reach a satisfactory level of accuracy, was partially driven by a recent comparative study that reported a good performance of SOAP in predicting adsorption energies of hydrogen on nanoclusters [START_REF] Jäger | Machine learning hydrogen adsorption on nanoclusters through structural descriptors[END_REF]. Several other descriptors for periodic materials have also been proposed in the literature, such as, for example, the many-body tensor representation (MBTR) [START_REF] Huo | Unified representation of molecules and crystals for machine learning[END_REF] or atom-centered symmetry functions (ACSF) [8]. While a systematic comparison of different ML approaches is beyond the purpose of this discussion, other descriptors and/or ML algorithms might also perform well for the problem we are considering. Below we summarize the main ideas of the SOAP and KRR approaches; a more detailed presentation can be found in Refs. 5, 130, 26. The numerical implementation of the ML models considered here is based on the DScribe libraries [START_REF] Himanen | Library of descriptors for machine learning in materials science[END_REF] 1 .

For a structure A of a given material let us define an environment X A i that includes the atoms surrounding a specific atom i within a certain radius. The density of the atoms ρ in the environment X A i is defined as a sum of Gaussians

ρ X A i (r) = j∈X A i exp - (r -x j ) 2 2σ 2 , ( 6.8)
where the Gaussians of variance σ 2 are centered at the positions x j of all the atoms in the environment X A i (including the central one). For systems containing different atomic species, as those considered here, a different environment has to be defined for each species. The SOAP kernel is then defined as an overlap integral of the densities of different environments, possibly from different structures A and

B: k(X A i , X B j ) = d R ρ X A i (r)ρ X B j ( Rr)dr 2 , (6.9)
where the first integral is performed over the three-dimensional rotations, as required in order to have a rotationally invariant representation. For practical purposes a normalised version of the kernel in Eq. 6.9 is typically used:

k(X A i , X B j ) = k(X A i , X B j ) k(X A i , X A i ) k(X B j , X B j )
. (6.10) This normalized kernel provides a measure of similarity since the maximum value (one) is attained when the two environments X A i and X B j are identical (independently of rotation and translations of the environment). Because of the specific form of the density of atoms, expressed as an expansion of Gaussians, the integral over the rotations in Eq. 6.9 can be carried out analytically. This is achieved by using a basis of spherical harmonics and the corresponding formalism is discussed in details in Refs. 5, 26. The kernels introduced in Eqs. 6.9-6.10 can be used to compare environments surrounding specific atoms and, in this respect, describe only a local similarity. For our specific purpose it is necessary to determine the similarity among different structures (all with the same chemical composition) generated by a MD simulation. This global information is included in the full covariance matrix C(A, B), whose elements

C ij (A, B) = k(X A i , X B j ) (6.11)
contain all the possible pairings of environments belonging to structure A and to structure B. Based on the matrix in Eq. 6.11 different global kernels can be defined. For our purpose the regularized entropy match (REMatch) kernel will be used, that is typically coupled with SOAP descriptors [5,26]. The REMatch kernel is defined as

K γ (A, B) = TrP γ C(A, B), (6.12)
where P is a doubly stochastic matrix that satisfies the condition

P γ = arg min P i,j P ij (1 -C ij + γlnP ij ), (6.13)
and γ is a hyperparameter that controls the entropic contribution. The REMatch kernel in Eq. 6.12 can then be coupled with the kernel ridge regression ML algorithm. Since our purpose is to predict RPA total energies for the configurations generated by a PBE+D2 MD, a certain number of configurations N train is selected to train the ML model (the specific selection criterion used will be discussed later). From these N train configurations the REMatch kernel K γ train is built and used to train the ML model: w = (K γ train + λ1) -1 y train (6.14) where y train in this case contains the RPA energies corresponding to the configurations in the training set and λ is a regularization hyperparameter used to prevent overfitting. The weights w determined by Eq. 6.14 are then used to predict the RPA energies for all the remaining N predict = N tot -N train configurations in the MD:

y predict = K γ predict w, (6.15)
where y predict contains the predicted energies and K γ predict is the N predict ×N train REMatch kernel.

As RPA calculations are highly expensive it is crucial to develop a ML model that can be trained with a number of configurations which is as small as possible. Indeed, N train explicit RPA calculations are required while the remaining N predict are inexpensively predicted by Eq. 6.15. To this purpose it is useful to use the energy difference ∆V i = V i -V i as a variable for training rather than the absolute energy V i . Indeed, the energy difference ∆V i has a much smoother behavior than V i and the error on the prediction can be lowered by one order of magnitude with respect to a ML model that uses total energies. This idea is based on the ∆machine learning approach [START_REF] Ramakrishnan | Big data meets quantum chemistry approximations: The δ-machine learning approach[END_REF]. For example, for CH 4 adsorbed in protonated chabazite (see Sec. 6.2), ∆V i (and, correspondingly, V i ) can be predicted with a root mean square error of about 0.4 kcal/mol by using only N train = 10 and with a systematic decrease of the error as a function of N train .

Statistical error evaluation

The accuracy of a ML model is characterized by a statistical error that can be possibly decreased by increasing the quantity of data used for training. In our specific case the energies V i or, equivalently, the energy differences ∆V i are predicted with a statistical error that has implications on the estimate of the internal energy (Eq. 6.7) and enthalpy (Eq. 6.6).

As a starting point for our discussion it is necessary to determine the error involved in the prediction of ∆V i and V i . To this purpose the root mean square error (RMSE) of the predicted energies with respect to the exact RPA results is computed for a small test set. This set includes configurations that are reasonably independent and uncorrelated with respect to configurations used for training the ML model (see Sec. 6.2 for additional details). Since the prediction is not biased (mean error is always well below 0.1 kcal/mol), the RMSE can be interpreted as the standard deviation of the distribution of the residuals in the regression model and will be denoted by s. The RMSE in the prediction of the RPA energies of individual configurations leads to both a systematic and random error in the internal energy expressed by Eq. 6.7.

To discuss the systematic error we will split the ensemble energy in Eq. 6.7 into two contributions

V H = Ntot i=1 V i exp {-β∆V i } Ntot i=1 exp {-β∆V i } + Ntot i=1 ∆V i exp {-β∆V i } Ntot i=1 exp {-β∆V i } = f 1 + f 2 (6.16)
and we will focus on the last term f 2 that depends only on ∆V i . Because of the statistical error associated with the machine learning model, we suppose that the exact ∆V i is replaced by values E normally distributed with variance s 2 :

P (E) = Kexp - (E -∆V i ) 2 2s 2 , (6.17)
where K is a normalizing factor. The estimate of f 2 that keeps into account this distribution can be written in the following way:

f2 = i dE E exp {-βE} × exp -(E-∆V i ) 2 2s 2 i dE exp {-βE} × exp -(E-∆V i ) 2 2s 2 = i dE E exp {-β∆V i } × exp -(E-∆V i -βs 2 ) 2 2s 2 i dE exp {-β∆V i } × exp -(E-∆V i -βs 2 ) 2 2s 2 = i dE (E + βs 2 )exp {-β∆V i } × exp -(E -∆V i ) 2 2s 2 i dE exp {-β∆V i } × exp -(E -∆V i ) 2 2s 2 = f 2 + βs 2 , (6.18)
where the sum over i is implicitly intended between 1 and N tot . Here f 2 is the reweighted average obtained using the exact energy differences ∆V i while f2 is based on the values of ∆V i predicted by a ML algorithm that introduces a random statistical error. Accordingly, when ML and TPT are coupled, a systematic error of βs 2 in the prediction of ensemble averages has to be expected. In the assumption that the residuals are normally distributed and that the variance s 2 is known, this bias can be easily evaluated using Eq. 6.18. The results that will be discussed below have all been corrected for this systematic error. The value of s and, consequently, of βs 2 decreases by increasing the number of training examples N train .

For example, if we consider CH 4 adsorbed in protonated chabazite (see Sec. 6.2 and in particular Table 6.1) the βs 2 correction is 0.33 kcal/mol for N train =10 and 0.07 kcal/mol for N train =200. By the same logic, it is easy to see that f1 = f 1 and hence the first term on the right hand side of Eq. 6.16 is unaffected by the systematic error arising from s.

We now consider the evaluation of the random error in the ensemble averages associated with the ML model. To this purpose it is important to notice that Eq. 6.7 is a non linear function of V i and, accordingly, an error propagation formula based on first derivatives will be used here. We have that the standard error δ M L in the ensemble average is

δ M L = Ntot i=1 ∂ V H ∂V i 2 s 2 V i = s Ntot i=1 ∂ V H ∂V i 2 (6.19)
where s V i = s is the standard deviation of the V i 's distribution and the derivative is explicitly given by

∂ V H ∂V i = w i (1 -βV i + β V H ) , (6.20)
where we defined the weight

w i = exp {-β(V i -V i )} Ntot i=1 exp {-β(V i -V i )} . ( 6.21)
In all calculations discussed further, δ M L is typically much smaller than s. Finally, we discuss the source of random error that is associated with the finite sample size. This error is not related to the ML approach but is typical of MD simulations, which span a finite time interval. As the term V H in Eq. 6.7 can be considered as a weighted average of the form

V H = Ntot i=1 V i w i (6.22)
with the weight w i defined in Eq. 6.21, the corresponding statistical error can be estimated using the standard formula [44]

δ s = N s (N s -1) Ns i=1 w 2 i (V i -V H ) 2 , (6.23)
where N s N tot is the number of uncorrelated samples. In practice, the latter can be identified e.g. by the block method of Flyvbjerg and Petersen. [37,38]. The value of δ s depends on the length and the quality of the MD trajectory and it can, in principle, be made arbitrarily small by performing a sufficiently long simulation.

In our numerical examples discussed below, the error bars are evaluated as

δ = δ 2 M L + δ 2 s (6.24)
where δ M L is defined in Eq. 6.19 and δ s in Eq. 6.23. The values of δ s are in all cases discussed in Sec. 6.2 significantly larger than those of δ M L and hence the former term represents a dominant contribution to δ.

RPA enthalpies of adsorption in zeolites

In order to establish the accuracy and efficiency of the MLPT methodology we consider the calculation of molecular enthalpies of adsorption in zeolites, which are porous materials used in many applications [START_REF] Speybroeck | Advances in theory and their application within the field of zeolite chemistry[END_REF] such as depollution [21,17,[START_REF] Hessou | Adsorption of no, no 2, co, h 2 o and co 2 over isolated monovalent cations in faujasite zeolite: A periodic dft investigation[END_REF], separation of chemicals [START_REF] Pham | Carbon dioxide and nitrogen adsorption on cation-exchanged ssz-13 zeolites[END_REF][START_REF] Maghsoudi | Adsorption isotherms and ideal selectivities of hydrogen sulfide and carbon dioxide over methane for the si-cha zeolite: comparison of carbon dioxide and methane adsorption with the all-silica dd3r zeolite[END_REF][START_REF] Khalil | Selective capture of phenol from biofuel using protonated faujasite zeolites with different si/al ratios[END_REF], and catalysis [START_REF] Grajciar | Towards operando computational modeling in heterogeneous catalysis[END_REF]18,[START_REF] Losch | Proton mobility, intrinsic acid strength, and acid site location in zeolites revealed by varying temperature infrared spectroscopy and density functional theory studies[END_REF]. Specifically, we consider four applications involving CH 4 and CO 2 in protonated chabazite (denoted as HCHAB) and siliceous chabazite (CHAB). The four corresponding adsorbed systems will be indicated as HCHAB+CH 4 , HCHAB+CO 2 , CHAB+CH 4 , and CHAB+CO 2 . The HCHAB+CH 4 example will be discussed in more detail to describe the numerical properties of the MLPT approach. The corresponding periodic model is shown in Fig. 6.1. Following Ref. [START_REF] Rocca | Bridging molecular dynamics and correlated wave-function methods for accurate finite-temperature properties[END_REF] we used MD trajectories based on the vdW corrected PBE+D2 functional [START_REF] Grimme | Semiempirical gga-type density functional constructed with a long-range dispersion correction[END_REF][START_REF] Perdew | Generalized gradient approximation made simple[END_REF] as a starting point for the MLPT method. MD sim- ulations were performed in the NVT ensemble. Temperature was set to 300 K and the length of trajectories was fixed at 100 ps with an integration step of 0.5 fs. In a production period succeeding the equilibration length of 5 ps, a total of 190,000 configurations were generated. The first step of the MLPT is the selection of a proper training set with N train elements. Since configurations along a MD trajectory are time correlated, we select configurations that are as distanced as possible in time. For example, N train =10 indicates that out of the 190,000 configurations one configuration every 19,000 is selected to obtain a total of 10 training structures. Five training sets are considered here with N train = 10, 25, 50, 100, and 200. Additionally, a test set of 20 configurations is selected to estimate the RMSE in the prediction of the machine learning model. Those configurations are chosen randomly but with the constraint of being at a distance of at least 350 time steps from the geometries contained in the largest training set (N train =200). In this way the test and training sets are reasonably independent. For the HCHAB+CH 4 system the RMSE error (also denoted with s) in the prediction of the RPA energies of configurations is shown in Table 6 Since configurations are selected using time separation as a criterion it is also important to understand to which extent they span the geometric space. This is particularly important for the adsorbed systems, where the molecule can move rather freely in the cage of the zeolite. If certain parts of the configurational space do not have representative structures included in the training set, the ML learning model might have limited predictive power in these regions. Fig. 6.2 helps visualizing the selected configurations for HCHAB+CH 4 using the t-Distributed Stochastic Neighbor Embedding algorithm (t-SNE) [START_REF] Maaten | Visualizing data using t-sne[END_REF]. t-SNE is an unsupervised learning method that allows for a dimensionality reduction to visualize highdimensional data: Points that are closer in the high-dimensional space have higher In this way the analysis is focused on the diffusion of the molecule in the zeolite rather than on the full system, that contains many more environments and degrees of freedom. Fig. 6.2(b) shows the configurations in the N train = 200 training set with respect to the total number of configurations. The sampling is uniform and all the regions of the configuration space are fairly well represented by the training set. For N train = 25 (Fig. 6.2(a)) the sampling is less satisfactory but, as shown below, the model can still be predictive. This is likely to be related to the use of the ∆-ML approach. Indeed, the model is trained on energy differences ∆V and a large contribution to the RPA total energy is already provided by the PBE+D2 total energy, whose exact value is available from the initial MD trajectory.

Detailed MLPT numerical results for CH 4 adsorbed in protonated chabazite are presented in Table 6.2, where values for the RPA ensemble averages of each system (CH 4 , protonated chabazite, and adsorbed system) and enthalpies of adsorption (in kcal/mol) are reported. Results are rather stable with respect to the N train variation and all the values of internal energies are well within a 0.2 kcal/mol interval. As few as 10 configurations can already provide highly accurate estimates of the RPA enthalphy of adsorption. As explained in Sec. 6.1, these values have been corrected for the systematic error βs 2 (see Eq. 6.18). This correction ranges from 0.33 kcal/mol (N train =10) to 0.07 kcal/mol (N train = 200). The table reports also error bars ±δ associated with random errors Eq. 6.24. The values of δ are rather insensitive to the change of N train ; this is mostly due to the fact that the finite sampling component δ s , which is independent of N train , provides a larger contribution than δ M L .

In order to highlight the importance of the ML approach in providing converged values for RPA internal energies and enthalpies, in Table 6.3 we present the results that would be obtained by applying thermodynamic perturbation the- ory only on the small set of training configurations. Namely, in this case no ML model is used and, for example, N conf = 100 means that TPT includes only the 100 configurations computed explicitly. Since configurations along a MD trajectory are time correlated, it is reasonable to expect convergence with a relatively small subset of configurations, as assumed in Ref. 33. However, as shown in Table 6.3 this is not the case for the sets considered in this work. For example, even the largest set N conf = 200 still predicts results for the enthalpy of adsorption that deviate by about 0.5 kcal/mol with respect to the most converged MLPT values in Table 6.2. Importantly, the enthalpies ∆ ads H in Table 6.3 strongly oscillate and no convergence can be inferred from these results. Finally, in Table 6.4 we summarize the RPA enthalpies of adsorption (in kcal/mol) as functions of N train for the all four systems considered in this work. The RPA, which is often considered a reference method for periodic condensed matter applications, provides systematically accurate estimates within chemical accuracy in comparison to experiments. The tendency to slightly underestimate is a well known behavior of the RPA [START_REF] Ren | Beyond the random-phase approximation for the electron correlation energy: The importance of single excitations[END_REF]. Enthalpies obtained from the PBE+D2 MD are instead less consistent, with large deviations for the adsorption in HCHAB and more reliable results for CHAB. Importantly, RPA results with N train = 10 can be obtained with a CPU time comparable to or even shorter than the underlying PBE+D2 MD but can provide a significantly higher level of accuracy and more reliable results for different types of systems. Chapter 7

Perspectives

In Chapter 6 we introduced a general purpose methodology, MLPT, to couple correlated methods with molecular dynamics at a reasonable computational cost. This approach allowed us to compute adsorption enthalpies of molecules in zeolites at the RPA of theory. Considering the high numerical cost of this approximation it would be inconceivable to directly perform a brute-force ab initio molecular dynamics. This recent work on MLPT will be at the base of several future developments and extensions. The ultimate goal will be to develop a highly accurate methodology to describe the properties of materials in realistic (finite-temperature) conditions. Specifically, the following research lines will be pursued:

• Further development of machine learning perturbation theory (MLPT) -The MLPT approach (see Chapter 6) is based on thermodynamic perturbation theory coupled to machine learning (ML) techniques. Within this approach a MD is performed using a numerically inexpensive functional and a small number of structures is chosen to sample the configurational space. The high-level correlated calculations are performed only on these geometries. Part of these high-level results are used to train a ML learning model and part of them are held out to test it (alternatively, a cross-validation approach could be used). The regression model is then used to predict the properties (e.g. the energy) of the full set of (hundreds of thousands) configurations and perturbation theory is applied to compute the finite temperature observables. This procedure can be improved in several aspects: (1) In order to build the training set algorithms could be developed to sample as uniformely as possible the configurational space (2) The performance of different types of descriptors and ML algorithms could be tested to make MLPT more efficient and/or accurate (e.g. the many-body tensor representation (MBTR) [START_REF] Huo | Unified representation of molecules and crystals for machine learning[END_REF] or atom-centered symmetry functions (ACSF) [8].). ( 3)

The dependence of the perturbative approach on the starting point should be better understood and approaches to evaluate the corresponding error should be developed. (4) Using beyond-RPA methods or other more traditional quantum chemistry approaches to further improve the accuracy of thermodynamic calculations. ( 5) Properties other than then the energy could be computed.

The MLPT method is still at an early stage of development but has huge potential to significantly impact the field of materials simulations. To fully establish this approach it will be necessary to understand its range of applicability and numerical properties.

• Machine learning molecular dynamics (MLMD) -This is a methodology that has seen a growing popularity since the seminal work of Behler and Parrinello in 2007 [10]. In this approach a large number of configurations of a certain system is generated by performing a MD simulation or by properly sampling the configurational space. A ML learning model, usually a neural network (NN), is then trained on the energies and/or forces computed for these configurations. A full molecular dynamics is then performed by using this neural network. At each step of the MD simulation, instead of directly solving the Schrödinger equation, forces and energies are inexpensively "predicted" by the ML model and used to generate the next configuration. The accuracy of the model is strictly related to the quality of the sampling used to generate the training examples. One of the advantages of the MLMD approach is that a NN model can be trained by using only the energies and then forces can be obtained by deriving the NN with respects to atomic positions. This avoids the problem of the limited availability of analytic gradient implementations for correlated methods (e.g. RPA) for condensed matter applications. To couple the Behler-Parrinello with RPAtype methods it is crucial to decrease the number of configurations used to train the ML model. Indeed, a high level correlated calculation is required for every training example. Similarly to MLPT, this optimal sampling could be obtained from energetic or geometric arguments. For example, a large number of configurations could be generated at first; then the features corresponding to these geometric configurations could be clustered by using unsupervised learning algorithms; finally, a few representatives of each cluster could be randomly chosen to span the configurational space and to train the NN model. For certain systems it could also be possible to exploit the symmetry to "increase" the amount of information used to train the neural network model [23]. Differently from MLPT, the MLMD approach does not apply perturbation theory but rather performs a full molecular dynamics simulation. Accordingly, the MLMD methods could possibly reach a higher level of accuracy and should not be affected by the choice of the starting point as in perturbative approaches. From a practical point of view some libraries to perform MLMD calculations are already openly available [START_REF] Zhang | Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics[END_REF] or, if necessary, we could develop our own codes by using Tensorflow.

• Direct calculation of gradients -As discussed before, directly preforming a molecular dynamics with correlated approximations (e.g. the RPA) requires impractical computational resources and complex algorithms to compute interatomic forces. An implementation of analytic forces has been recently proposed for the random phase approximation (RPA) within a periodic plane-wave basis set framework [START_REF] Ramberger | Analytic interatomic forces in the random phase approximation[END_REF]; to the best of our knowledge this is the only one based on the PW approach. In the context of this project we will couple the ideas developed in Ref. 116 with the efficient algorithms described in Chapter 3 to obtain an efficient implementation to compute the interatomic forces for the RPA and possibly also for the beyond-RPA and MP2 methods. This part of the work will serve two purposes: (1) At least for the some small systems it will be possible to run a full molecular dynamics to compare and benchmark the MLPT and MLMD approaches; (2) The gradients will be used in the context of MLMD to train machine learning models that specifically require this information [23,[START_REF] Zhang | Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics[END_REF].

These research directions represent different paths heading towards the same goal -the development of a general purpose methodology to compute finite temperature properties using correlated methods, such as the RPA. For the MLPT approach some proof-of-principle results are already available and show great promise. However, most of the ideas proposed here remains unexplored and it is not possible to establish which approach will give the best results in term of accuracy and numerical efficiency. Importantly, the different paths proposed above are not mutually exclusive. For example, the development of different techniques to sample the configurational space could improve both the MLPT and MLDM approaches. Similarly, an implementation to compute interatomic forces could be used to train the models used for MLMD. As the main outcome of our future work, we will establish which methodology is the most general, efficient, and accurate for finite temperature correlated simulations.

Beyond the methodologial development aspect, applications to realistic materials will be considered in the future with the purpose of testing the new methods, establishing them within the scientific community, and shedding new light on important experimental and technological problems. These will include:

• Study of the adsorption of molecules on surfaces or porous materials.

• Study of molecular crystals. These systems, which are important, for example, in pharmaceutical applications, are usually characterized by polymorphs with similar energetic stability. The methods described above will be used to rank the relative stability of different polymorphs.

• Study of few-layer materials and heterostructures of two-dimensional materials. The new methodologies will be used to study the structural parameters of these systems and the interlayer binding energy.

Prospects and Applications, Benasque (Spain), January Awards and Fellowships

Memberships

American Physical Society, 2008-present Appendix C

Résumé étendu en français

Ce manuscrit décrit mes activités de recherche principales (mais non exclusives) développées à partir de septembre 2012, lorsque j'ai été embauché en tant que maître de conférences à l'Université de Lorraine. Mon travail a été principalement concentré sur le développement méthodologique pour le calcul des énergies de corrélation basées sur l'approximation de la phase aléatoire (RPA) et de ses variantes.

En raison de son compromis raisonnable entre précision et efficacité, la théorie de la fonctionnelle de la densité (DFT) [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] est devenue la méthode la plus utilisée pour le calcul ab initio des propriétés des matériaux. La DFT s'applique à des systèmes relativement grands comprenant jusqu'à des milliers d'atomes et son succès initial est fortement lié à des approximations assez simples pour la fonctionnelle d'échange-corrélation (xc), comme les approximations de la densité locale (LDA) [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF][START_REF] Perdew | Self-interaction correction to densityfunctional approximations for many-electron systems[END_REF] et des gradients généralisés (GGA) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. Des classes de fonctionnelles plus récentes incluent la meta-GGA [START_REF] Tao | Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids[END_REF][START_REF] Sun | Strongly constrained and appropriately normed semilocal density functional[END_REF] et les hybrides, qui mélangent un certain pourcentage d'échanges de Hartree-Fock avec des contributions semi-locales [7]. Toutes ces méthodes impliquent un niveau croissant de précision mais également de coût numérique du calcul. Le schéma de l'échelle de Jacob des approximations DFT (voir Dans le contexte de la DFT, la base théorique pour la RPA est fournie par le théorème de fluctuation-dissipation avec connexion adiabatique (ACFDT) [START_REF] Langreth | The exchange-correlation energy of a metallic surface[END_REF][START_REF] Langreth | Exchange-correlation energy of a metallic surface: Wave-vector analysis[END_REF], qui est présenté au chapitre 2. L'ACFDT fournit en principe une expression exacte pour l'énergie de corrélation de l'état fondamental en utilisant des fonctions de réponse linéaires (c'est-à-dire la polarisabilité ou la matrice diélectrique). Dans le contexte de l'ACFDT la RPA est l'approximation pratique la plus simple qui peut être utilisée pour des applications numériques à des systèmes réalistes, tels que des molécules [39,[START_REF] Toulouse | Adiabatic-connection fluctuation-dissipation density-functional theory based on range separation[END_REF][START_REF] Nguyen | Efficient calculation of exact exchange and rpa correlation energies in the adiabatic-connection fluctuation-dissipation theory[END_REF][START_REF] Zhu | Rangeseparated density-functional theory with random phase approximation applied to noncovalent intermolecular interactions[END_REF][START_REF] Nguyen | A first-principles study of weakly bound molecules using exact exchange and the random phase approximation[END_REF]35,[START_REF] Ren | Resolution-ofidentity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF][START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF], des solides [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Lu | Ab initio calculation of van der waals bonded molecular crystals[END_REF][START_REF] Wilson | Iterative calculations of dielectric eigenvalue spectra[END_REF][START_REF] Harl | Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids[END_REF][START_REF] Kaoui | Random phase approximation correlation energy using a compact representation for linear response functions: application to solids[END_REF], des matériaux en couches minces [START_REF] Marini | First-principles description of correlation effects in layered materials[END_REF][START_REF] Lebègue | Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase approximation[END_REF] et des molécules adsorbées sur des surfaces [START_REF] Ren | Exploring the random phase approximation: Application to co adsorbed on cu(111)[END_REF][START_REF] Schimka | Accurate surface and adsorption energies from many-body perturbation theory[END_REF]. En particulier, pour les applications à la matière condensée les méthodes de chimie quantique traditionnelles (par exemple la théorie du cluster couplé) peuvent être difficilement utilisées et la RPA est probablement l'approximation corrélée la plus sophistiquée d'intérêt pratique pour ces systèmes. L'intérêt croissant suscité par la RPA est principalement (mais pas exclusivement) lié à la capacité de la RPA de décrire avec haute précision les forces de van der Waals (vdW) [32,41]. Malgré son succès croissant, la RPA reste une approche de niche, surtout si on la compare à d'autres approximations DFT qui sont appliqués couramment dans les simulations des matériaux. Ceci est lié à différentes raisons : (1) FIGURE C.2 : L'interaction effective au sein de la RPA est obtenue par la somme à l'ordre infini des diagrammes "à anneau" (figure reproduite de Ref. [START_REF] Ke | Dispersion relation of excitation mode in spinpolarized fermi gas[END_REF]).

Bien qu'il s'agisse d'une théorie très sophistiquée, des problèmes subsistent au sein de la RPA : les courbes de liaison de certains dimères diatomiques présentent une "bosse" non physique à des distances intermédiaires [39], les énergies de corrélation totales sont généralement surestimés [START_REF] Kurth | Density-functional correction of randomphase-approximation correlation with results for jellium surface energies[END_REF] et les énergies de liaison sont souvent sous-estimées [41,[START_REF] Harl | Accurate bulk properties from approximate many-body techniques[END_REF]. (2) Les implémentations numériques de la RPA sont plutôt complexes et les calculs pratiques relèvent plusieurs défis, tels que la diagonalisation complète de l'hamiltonien pour obtenir tous les états virtuels (de conduction), la manipulation de matrices de grande taille correspondantes aux fonctions de réponse et l'évaluation d'intégrales impropres sur la fréquence. (3) Le temps de calcul nécessaire pour la RPA est souvent supérieur de plusieurs ordres de grandeur par rapport au temps requis par les approximations (semi-)locales de la DFT ; de plus, le développement d'une implémentation pour le calcul des forces interatomiques analytiques n'est pas banale. Pour ces raisons, les applications aux matériaux réels exigent des ressources de calcul considérables et généralement sont basées sur des géométries fixes optimisées en utilisant des différentes théories.

Dans ce manuscrit, je présente des nouvelles méthodologies pour surmonter ces difficultés. Au chapitre 2, je discute de mes travaux sur le développement de nouvelles approximations qui évitent le problème d'autocorrélation de la RPA en utilisant un noyau d'échange non local. Plus précisément, j'ai développé deux nouvelles approximations pour l'énergie de corrélation de l'état fondamental basées sur l'ACFDT : le "electron-hole time-dependent Hartree-Fock" (eh-TDHF) et le " adiabatic connection second order screened exchange" (AC-SOSEX) [START_REF] Mussard | Dielectric matrix formulation of correlation energies in the random phase approximation: Inclusion of exchange effects[END_REF]28,31]. Au chapitre 3 je présente l'implémentation numérique que j'ai développé pour le calcul des énergies de corrélation pour les méthodes RPA est au-delà en utilisant une base d'ondes planes (PW) et des algorithmes itératifs. La taille de cette base peut atteindre plusieurs centaines de milliers d'éléments pour les grandes mailles. Pour cette raison, il est pratique de représenter les fonctions de réponse impliquées dans les méthodes RPA et au-delà de la RPA à l'aide d'une base auxiliaire compacte. Dans ce manuscrit on utilise les vecteurs propres d'une polarisabilité approximative ne comprenant que la contribution de l'énergie cinétique [START_REF] Rocca | Random-phase approximation correlation energies from lanczos chains and an optimal basis set: Theory and applications to the benzene dimer[END_REF]. La convergence des énergies de corrélation par rapport à la base n'est pas banale et est analysée en détail au chapitre 5. Afin d'évaluer efficacement les éléments de matrice de la polarisabilité sur la base auxiliaire il est utile de recourir à des techniques numériques supplémentaires : l'élimination des états de conduction pour obtenir un formalisme qui ne dépend que des états occupés [3,[START_REF] Walker | Efficient approach to time-dependent density-functional perturbation theory for optical spectroscopy[END_REF][START_REF] Rocca | Turbo charging time-dependent density-functional theory with Lanczos chains[END_REF] et l'utilisation de l'algorithme de Lanczos pour calculer le résolvant de l'Hamiltonien [START_REF] Rocca | Turbo charging time-dependent density-functional theory with Lanczos chains[END_REF][START_REF] Umari | Gw quasiparticle spectra from occupied states only[END_REF][START_REF] Nguyen | Improving accuracy and efficiency of calculations of photoemission spectra within the many-body perturbation theory[END_REF]. Au chapitre 4 cette implémentation numérique est appliquée à une série d'ensembles moléculaires pour lesquels des résultats de haut niveau sont disponibles et à des solides. La le temps de calcul nécessaire pour appliquer la théorie de la perturbation avec la RPA est réduit de plusieurs ordres de grandeur. Les résultats montrant la précision et l'efficacité de la méthode MLPT sont présentés dans le tableau C.1. Plus précisément, l'enthalpie d'adsorption du méthane dans un matériau poreux (la zéolite chabazite) a été calculée en utilisant la RPA. Le modèle périodique correspondant, qui implique 200 électrons de valence, est présenté en figure C.4. Comme le montre le tableau C.1, alors que la fonctionnelle de départ PBE+D2 donne un résultat qui dépasse largement le seuil de précision chimique (1 kcal / mol), la valeur de la RPA est en excellent accord avec la valeur expérimentale. En conséquence, l'utilisation de MLPT associé à la RPA ou à d'autres méthodes corrélées au-delà de la RPA peut atteindre un niveau de précision beaucoup plus élevé tout en limitant dans une large mesure l'investissement en temps de calcul.

  within the random phase approximation and beyond 2.1 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . 2.2 Correlation energy within the adiabatic connection fluctuationdissipation theorem . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Random phase approximation . . . . . . . . . . . . . . . . . . . 2.4 Electron-hole time-dependent Hartree-Fock for correlation energies 2.5 Adiabatic connection second order screened exchange approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Numerical implementation within a plane-wave basis set using iterative algorithms 3.1 Implementation using an auxiliary basis set . . . . . . . . . . . . 3.2 Elimination of the empty states . . . . . . . . . . . . . . . . . . . 3.3 Lanczos algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Construction of a compact auxiliary basis set by diagonalizing an approximate response function . . . . . . . . . . . . . . . . . . . 3.5 Gram-Schmidt orthogonalization . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: The Jacob's ladder of density functional approximations proposed by John Perdew to visualize the different levels of approximation in the exchangecorrelation functionals. By climbing the ladder the accuracy increases at the expense of a higher computational cost (figure reproduced from Ref. 107).

  3.22). For example, for a molecule the importance of the φ LU M O φ HOM O term in the sum in Eq. 3.22 is enhanced by the factor 1/( HOM O -LU M O ); similarly, the contribution of the φ (LU M O+N ) φ HOM O term (where N is a large integer) to the sum in Eq. 3.22 might be much smaller because 1/( HOM O -(LU M O+N ) ) goes to 0 for N → +∞. Accordingly, the optimal basis set should well represent the φ LU M O φ HOM O term rather than φ (LU M O+N ) φ HOM O .

Figure 3 . 1 :

 31 Figure 3.1: Convergence with respect to the auxiliary basis set size of the correlation energy contribution to the energy of the C 2 H 2 +H 2 O → CH 3 CHO reaction. Specifically, ∆E c = E CH 3 CHO

Figure 3 . 2 :

 32 Figure 3.2: This figure shows how the number of vectors involved in the computation of dielectric matrices can be decreased by a Gram-Schmidt orthogonalization (GSO) procedure. The original method requires N aux × N v vectors, which grows quadratically with the number of electrons. The GSO version requires N ortho matrix elements, which grows linearly with the number of electrons. See text for further details.

Figure 4 . 1 :

 41 Figure 4.1: Differences in kcal/mol for the A24 test set introduced by applying the Gram-Schmidt orthogonalization procedure described in the text.

Figure 4 . 2 :

 42 Figure 4.2: Differences in kcal/mol for the A24 test set for the RPA, eh-TDHF, and AC-SOSEX methods computed with respect to CCSDT(Q) reference values.For the sake of completeness the errors of CCSD(T)[START_REF] Řezá Č | Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the "Gold Standard," CCSD(T) at the Complete Basis Set Limit?[END_REF] are also included.

  3, and the errors with respect to the CCSD(T) values of Takatani et al. [142] are shown graphically in Fig. 4.3.

2 )Figure 4 . 3 :

 243 Figure 4.3: Differences in kcal/mol for the S22 test set for the RPA, eh-TDHF, and AC-SOSEX methods computed with respect to CCSD(T) reference values reported by Takatani et al.[START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF]. For the sake of completeness we included the deviations of the CCSD(T) values previously reported by Jurečka et al.[START_REF] Jure Čka | Benchmark database of accurate (mp2 and ccsd (t) complete basis set limit) interaction energies of small model complexes, dna base pairs, and amino acid pairs[END_REF]

Figure 4 . 4 :

 44 Figure 4.4: Binding curve of the ethene dimer.

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Binding curve of the ethene-ethyne dimer.

Figure 4 . 7 :

 47 Figure 4.7: Mean absolute error (MAE), mean error (ME), and maximum absolute deviation (|MAX|) for a test set of 16 reaction energies [65]. Errors have been computed for PBE, HF@PBE, RPA, AC-SOSEX, and eh-TDHF results obtained within the pane-wave (PW) basis set implementation with respect to CCSD(T) reference values.

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Mean absolute error (MAE), mean error (ME), and maximum absolute deviation (|MAX|) for a test set of 16 reaction energies [65]. Errors have been computed for PBE, HF@PBE, RPA, AC-SOSEX, and eh-TDHF results obtained within the Gaussian basis set implementation with respect to CCSD(T) reference values.

Ry, 6 × 6 × 6 , 4 × 4 × 4 ;

 666444 RPA@LDA for Si: 40 Ry, 6 × 6 × 6, and 4 × 4 × 4; RPA@GGA for Si: 40 Ry, 6 × 6 × 6, and 4 × 4 × 4; RPA@LDA for SiC: 70 Ry, 6 × 6 × 6, and 4 × 4 × 4; RPA@GGA for SiC: 70 Ry, 6 × 6 × 6, and 4 × 4 × 4; RPA@LDA for Ne: 160 Ry, 4 × 4 × 4, and 4 × 4 × 4; RPA@GGA for Ne: 160 Ry, 4 × 4 × 4, and 4 × 4 × 4; RPA@LDA for Ar: 60 Ry, 4 × 4 × 4, and 4 × 4 × 4; RPA@GGA for Ar: 80 Ry, 4 × 4 × 4, and 4 × 4 × 4; RPA@LDA for Kr: 70 Ry, 4 × 4 × 4, and 4 × 4 × 4; RPA@GGA for Kr: 80 Ry, 4 × 4 × 4, and 4 × 4 × 4.

  correlation energy equations we should always use N monomer1 aux

Figure 5 . 1 :

 51 Figure 5.1: Interaction energy (kcal/mol) of the formaldehyde dimer at the RPA (red), eh-TDHF (blue) and AC-SOSEX (green) levels of theory as a function of N aux /N v . Dashed lines indicate the CBS interaction energy values computed with Eq. 5.4.

Figure 5 . 2 :

 52 Figure 5.2: Correlation energy (red dots) of the formaldehyde dimer at the RPA level of theory as a function of 1/N aux . The points correspond to the range from N aux = 21×N v to N aux = 70×N v . The continuous line denotes a linear fit of the data. The dotted and dot-dash lines show partial fits, demonstrating how rapidly the approach converges.

5 . 3 Figure 5 . 3 :

 5353 Figure 5.3: Mean absolute errors (MAEs) of RPA (red), eh-TDHF (blue) and AC-SOSEX (green) for the A24 test set interaction energies using CBS extrapolations for E CBS(11,30) , E CBS(21,40) , and E CBS(31,50) , and size-consistent energy differences for E Naux=30 , E Naux=40 , and E Naux=50 . The E CBS(41,70) interaction energy was used as a reference (see text).

Figure 5 . 4 :

 54 Figure 5.4: Mean absolute errors (MAEs) of RPA (red), eh-TDHF (blue) and AC-SOSEX (green) for the methane-formaldehyde complex interaction energies using CBS extrapolations for E CBS(11,30) , E CBS(21,40) , and E CBS(31,50) , and size-consistent energy differences for E Naux=30 , E Naux=40 , and E Naux=50 . The E CBS(41,70) interaction energy was used as a reference (see text).

Figure 5 . 5 :

 55 Figure 5.5: Potential energy surface points (kcal/mol) of the methaneformaldehyde complex at the extrapolated RPA (red), eh-TDHF (blue) and AC-SOSEX (green) levels of theory compared with the CCSD(T) values (black).

Figure 6 . 1 :

 61 Figure 6.1: Periodic model used for CH 4 adsorbed in protonated chabazite. The unit cell used in our calculations is shown by solid lines.

Figure 6 . 2 :

 62 Figure 6.2: For HCHAB+CH 4 the dimensionality reduction t-SNE algorithm is used to visualize the configurations selected to train the machine learning model (in red) with respect to the total set of configurations (in blue). (a) N train = 25; (b) N train = 200.

••

  Prime d'encadrement doctoral et de recherche -PEDR (2017-2021): Salary bonus obtained upon evaluation of research and supervising activities. • Lorraine Université d'Excellence (LUE) Future Leader -ERC track (2018 and 2019) • Chaire d'Excellence CNRS (2012-2017) • Two travel grants from the Materials Computation Center (MCC) at the University of Illinois Urbana-Champaign (2008, 2011) • Undergraduate fellowship of the University of Pavia (1998-2002) • Collegio Borromeo fellowship for undergraduate student (1998-2002) Undergraduate fellowship of the Institute for Advanced Study (IUSS), Pavia (1998-2002).

  FIGURE C.1 : L'"échelle de Jacob" proposée par John Perdew pour visualiser les différents niveaux d'approximation utilisées pour les fonctionnelles d'échangecorrélation en DFT. En grimpant les échelons la précision augmente au détriment du coût de calcul, qui augmente considérablement (figure reproduite de Ref. 107).

2 )FIGURE C. 3 :

 23 FIGURE C.3 : Différences en kcal/mol pour l'ensemble S22 pour les méthodes RPA, eh-TDHF et AC-SOSEX calculées par rapport aux valeurs de référence obtenus pas la théorie du cluster couplé CCSD(T)[START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF]

Table 4 .

 4 1: Mean absolute errors (MAE) and maximum absolute deviations (|MAX|) in kcal/mol for the A24 test set introduced by applying the Gram-Schmidt orthogonalization procedure described in the text.

	System	RPA eh-TDHF AC-SOSEX
	MAE	0.047	0.033	0.029
	|MAX|	0.084	0.075	0.078

Table 4 .

 4 2: Interaction energies (in kcal/mol) for the A24 test set as obtained with the RPA, eh-TDHF, and AC-SOSEX methods compared with CCSDT(Q) and CCSD(T) reference values. The mean absolute errors (MAE), mean errors (ME) and the maximum deviations (|MAX|) are provided with respect to CCSDT(Q) values for the three types of molecular dimers and for the total test set.

	System		CCSDT(Q) [121] CCSD(T) [121]	RPA	eh-TDHF AC-SOSEX
			Hydrogen-bonded systems
	1 water. . . ammonia	C s	-6.492	-6.493	-5.753 -6.037	-6.100
	2 water dimer	C s	-4.994	-5.006	-4.223 -4.603	-4.691
	3 HCN dimer	C s	-4.738	-4.745	-4.506 -4.962	-5.107
	4 HF dimer	C s	-4.564	-4.581	-3.837 -4.244	-4.350
	5 ammonia dimer	C 2h	-3.141	-3.137	-2.584 -2.748	-2.726
	MAE			0.008	0.605	0.356	0.339
	ME			-0.007	0.605	0.267	0.191
	|MAX|			0.017	0.771	0.455	0.415
			Mixed-type systems
	6 methane. . . HF	C 3v	-1.660	-1.654	-1.371 -1.626	-1.651
	7 ammonia. . . methane	C 3v	-0.771	-0.765	-0.511 -0.644	-0.645
	8 methane. . . water	C s	-0.665	-0.663	-0.510 -0.664	-0.660
	9 formaldehyde dimer	C s	-4.479	-4.554	-2.674 -3.255	-3.572
	10 ethene. . . water	C s	-2.564	-2.557	-2.018 -2.364	-2.467
	11 ethene. . . formaldehyde C s	-1.623	-1.621	-1.282 -1.519	-1.562
	12 ethyne dimer	C 2v	-1.529	-1.524	-1.304 -1.526	-1.568
	13 ethene. . . ammonia	C s	-1.382	-1.374	-0.993 -1.153	-1.176
	14 ethene dimer	C 2v	-1.106	-1.090	-0.731 -0.912	-0.898
	15 methane. . . ethene	C s	-0.509	-0.502	-0.324 -0.422	-0.421
	MAE			0.013	0.442	0.224	0.179
	ME			-0.002	0.442	0.224	0.171
	|MAX|			0.075	1.792	1.237	0.927
			Dispersion-dominated systems
	16 borane. . . methane	C s	-1.513	-1.485	-0.636 -0.883	-0.850
	17 methane. . . ethane	C s	-0.836	-0.827	-0.621 -0.717	-0.671
	18 methane. . . ethane	C 3	-0.614	-0.607	-0.455 -0.512	-0.471
	19 methane dimer	D 3d	-0.539	-0.533	-0.416 -0.496	-0.457
	20 methane. . . Ar	C 3v	-0.405	-0.405	-0.263 -0.362	-0.343
	21 ethene. . . Ar	C 2v	-0.365	-0.364	-0.165 -0.226	-0.211
	22 ethene. . . ethyne	C 2v	+0.794	+0.821	+1.007 +1.143	+1.291
	23 ethene dimer	D 2h	+0.909	+0.934	+1.112 +1.263	+1.462
	24 ethyne dimer	D 2h	+1.084	+1.115	+1.311 +1.454	+1.561
	MAE			0.015	0.264	0.242	0.314
	ME			0.015	0.264	0.242	0.314
	|MAX|			0.031	0.887	0.668	0.699
			Total
	MAE			0.013	0.409	0.259	0.263
	ME			0.004	0.409	0.240	0.229
	|MAX|			0.075	1.792	1.237	0.927

  Takatani et al.

			Jurečka et al.
			Hydrogen-bonded systems
	1 ammonia dimer	C 2h	-3.17	-3.17	-2.54	-2.72	-2.71
	2 water dimer	C s	-5.02	-5.02	-4.20	-4.58	-4.65
	3 formic acid dimer	C 2h	-18.80	-18.61 -17.05	-18.24	-18.60
	4 formamide dimer	C 2h	-16.12	-15.96 -14.75	-15.57	-15.71
	5 uracil dimer	C 2h	-20.69	-20.65 -18.78	-20.32	-21.02
	6 2-pyridoxine. . . 2-aminopyridine C 1	-17.00	-16.71 -15.78	-16.46	-16.61
	7 adenine. . . thymine (WC)	C 1	-16.74	-16.37 -15.63	-16.30	-16.37
	MAE			0.15	1.26	0.48	0.36
	ME			0.15	1.26	0.48	0.27
	|MAX|			0.37	1.91	0.56	0.46
		Dispersion-dominated systems
	8 methane dimer	D 3d	-0.53	-0.53	-0.47	-0.52	-0.48
	9 ethene dimer	D 2d	-1.50	-1.51	-1.15	-1.33	-1.31
	10 benzene. . . methane	C 3	-1.45	-1.50	-1.18	-1.32	-1.28
	11 benzene dimer	C 2h	-2.62	-2.73	-1.92	-1.86	-1.70
	12 pyrazine dimer	C s	-4.20	-4.42	-3.40	-3.39	-3.24
	13 uracil dimer	C 2	-9.74	-10.12	-8.61	-9.11	-9.26
	14 indole. . . benzene	C 1	-4.59	-5.22	-3.41	-3.35	-3.19
	15 adenine. . . thymine (stack)	C 1	-11.66	-12.23 -10.36	-10.59	-10.62
	MAE			0.25	0.72	0.60	0.65
	ME			-0.25	0.72	0.60	0.65
	|MAX|			0.63	1.30	1.24	1.40
			Mixed-type systems
	16 ethene. . . ethine	C 2v	-1.51	-1.53	-1.29	-1.50	-1.57
	17 benzene. . . water	C s	-3.29	-3.28	-2.89	-3.10	-3.10
	18 benzene. . . ammonia	C s	-2.32	-2.35	-1.91	-2.07	-2.06
	19 benzene. . . hydrogen cyanide	C s	-4.55	-4.46	-4.03	-4.49	-4.61
	20 benzene dimer	C 2v	-2.71	-2.74	-2.31	-2.45	-2.42
	21 indole. . . benzene (T-shape)	C 1	-5.62	-5.73	-5.18	-5.38	-5.37
	22 phenol dimer	C 1	-7.09	-7.05	-6.31	-6.76	-6.79
	MAE			0.05	0.46	0.20	0.21
	ME			-0.01	0.46	0.20	0.17
	|MAX|			0.11	0.78	0.33	0.33
			Total
	MAE			0.15	0.81	0.43	0.42
	ME			-0.04	0.81	0.43	0.38
	|MAX|			0.63	1.91	1.24	1.40

Table 4 .

 4 4: Binding energies in kcal/mol for three weakly bound dimers obtained by different approximations: PBE, HF@PBE, RPA, AC-SOSEX, and eh-TDHF. Reference CCSD(T) results are provided in the last column. Results have been obtained by using plane-wave (PW) and Gaussian basis sets (see text for discussion and comparison).

	PBE HF@PBE RPA AC-SOSEX eh-TDHF CCSD(T)	Gaussian Gaussian Gaussian Gaussian Gaussian	-1.50 -0.34 1.15 -1.11 -1.27 -1.26	-1.51 -1.18 -0.27 -1.24 -1.47 -1.40	-4.95 -2.64 -4.36 -4.70 -4.63
	RPA AC-SOSEX eh-TDHF	PW PW PW	-1.19 -1.35 -1.37	-1.31 -1.55 -1.49	-4.35 -4.74 -4.67
	HF@PBE	PW	1.12	-0.32	-2.73
	PBE	PW	-0.35	-1.20	-5.08
	Dimer		ethene dimer	ethene-ethyne dimer	water dimer

Table 4 .

 4 5: Reaction energies in kcal/mol for 16 chemical reactions obtained by different approximations: PBE, HF@PBE, RPA, AC-SOSEX, and eh-TDHF. The mean absolute error (MAE), mean error (ME), and maximum absolute deviation (|MAX|) are computed with respect to the CCSD(T) reference values[START_REF] Hesselmann | Random-phase-approximation correlation method including exchange interactions[END_REF] in the last column. Results have been obtained by using plane-wave (PW) and Gaussian basis sets (see text for discussion and comparison).

	AC-SOSEX eh-TDHF CCSD(T)	Gaussian Gaussian Gaussian	-51.38 -49.44 -51.78	-40.36 -39.47 -42.49	-18.32 -18.18 -18.88	-4.60 -5.47 -4.72	-30.74 -29.70 -32.53	-89.37 -87.63 -92.43	-39.51 -38.28 -39.90	-14.56 -14.12 -15.92	-26.42 -25.28 -27.80	-9.28 -10.26 -9.71	-4.08 -6.18 -4.08	-22.00 -20.70 -22.87	-12.10 -13.59 -12.45	-1.56 -1.16 -1.72	-93.45 -93.81 -96.51	-5.43 -3.83 -5.83	1.10 1.97	-0.38 -1.41	4.80
	RPA	Gaussian	-47.63	-37.65	-17.23	-2.63	-26.47	-82.10	-38.37	-13.50	-22.76	-6.55	-5.57	-17.37	-10.03	-0.76	-87.67	-5.30	2.42	2.22	6.14
	HF@PBE	Gaussian	-50.97	-39.60	-21.99	1.55	-30.20	-96.95	-35.39	-10.06	-25.64	2.98	7.06	-19.87	2.75	-1.69	-89.90	-9.22	5.19	2.49	16.34
	PBE	Gaussian	-52.14	-40.59	-18.61	-12.36	-29.60	-81.82	-44.19	-15.57	-23.52	-20.92	-17.55	-18.40	-22.31	-1.75	-99.37	5.04	4.64	-2.28	11.37
	AC-SOSEX eh-TDHF	PW PW	-51.40 -50.55	-42.19 -39.90	-19.13 -18.62	-3.66 -3.41	-31.96 -29.66	-92.07 -88.31	-39.58 -39.07	-15.45 -13.94	-27.27 -25.41	-8.18 -7.48	-2.92 -3.16	-22.33 -20.70	-11.21 -10.71	-1.63 -1.39	-94.99 -91.47	-6.45 -5.70	2.02 1.19	-0.83 0.48	4.44 3.02
	RPA	PW	-46.96	-37.20	-17.39	-1.17	-25.61	-81.05	-38.09	-13.08	-21.95	-4.75	-4.61	-16.68	-8.68	-0.52	-85.66	-6.29	3.27	2.96	8.15
	HF@PBE	PW	-50.49	-39.12	-22.18	2.57	-29.53	-97.51	-34.91	-9.38	-24.96	3.96	7.64	-19.45	3.65	-1.83	-89.87	-9.21	5.53	2.91	17.24
	PBE	PW	-52.08	-40.53	-18.71	-11.97	-29.52	-82.37	-44.05	-15.37	-23.39	-20.89	-16.76	-18.95	-21.84	-2.28	-99.13	5.21	4.49	-2.22	10.63
	Reaction	Basis set	C 2 H 2 +H 2 → C 2 H 4	C 2 H 4 +H 2 → C 2 H 6	C 2 H 6 +H 2 → 2CH 4	CO+H 2 → HCHO	HCHO+H 2 → CH 3 OH	H 2 O 2 +H 2 → 2H 2 O	C 2 H 2 +H 2 O → CH 3 CHO	C 2 H 4 +H 2 O → C 2 H 5 OH	CH 3 CHO+H 2 → C 2 H 5 OH	CO+NH 3 → HCONH 2	CO+H 2 O → CO 2 +H 2	HNCO+NH 3 → NH 2 CONH 2	CO+CH 3 OH → HCOOCH 3	HCOOH+NH 3 → HCONH 2 +H 2 O	CO+H 2 O 2 → CO 2 O 2 +H	4 O+CO C 2 H H 2 CCO+HCHO →	MAE	ME	|MAX|

Table 4 .

 4 6: Lattice parameter as computed within different approximations (LDA, GGA, RPA@LDA, and RPA@GGA) for six different materials. Experimental values and previous RPA results in the literature are also provided for reference.

	Method	C	Si	SiC	Ne	Ar	Kr
	LDA	6.67	10.20	8.21	7.27	9.32	10.09
	RPA@LDA	6.68	10.11	8.17	8.97	10.23	10.83
	GGA	6.73	10.33	8.26	8.69	11.22	12.12
	RPA@GGA	6.73	10.21	8.21	8.18	10.00	10.72
	RPA@LDA Literature 6.75 a	10.27 a , 10.12 b /	8.88 c	10.20 c 10.96 c
	RPA@GGA Literature 6.75 a	10.26 a	8.25 a	8.50 c	10.02 c 10.77 c
	Experiments	6.74 d	10.26 d	8.24 d	8.22 e	9.88 e	10.60 e

a Ref. 60. b Ref. 84. c Ref. 58. d Experimental lattice constants extrapolated at 0 K. See Ref. 139 and references therein. e Experimental lattice constants extrapolated at 0 K. See Ref. 129 and references therein.

Table 4 .

 4 7: Bulk modulus as computed within different approximations (LDA, GGA, RPA@LDA, and RPA@GGA) for six different materials. Experimental values and previous RPA results in the literature are also provided for reference.

	Method	C	Si	SiC	Ne	Ar	Kr
	LDA	449.7	92.10	216.1	8.25	6.71	6.15
	RPA@LDA	441.1	102.1	222.4	0.26	1.58	2.24
	GGA	416.2	85.70	204.4	1.17	0.70	0.61
	RPA@GGA	417.0	97.30	215.0	1.47	2.79	2.98
	RPA@LDA Literature 439 a	98 a , 103 b	/	/	/	/
	RPA@GGA Literature 441	99	223	/	/	/
	Experiments	443.0 c	99.2 c	225.0 c	1.09 d	2.38 d	3.61 d

a Ref. 60. b Ref. 84. c See Ref. 139 and references therein. d See Ref. 129 and references therein.

Table 6 .

 6 1: RMSE (kcal/mol) in the prediction of RPA energies of the configurations of protonated chabazite (HCHAB), CH 4 , and adsorbed system (HCHAB+CH 4 ) as a function of the number of training structures N train . The RMSE is estimated on a test set of 20 configurations (see text).

	N train HCHAB+CH 4 HCHAB CH 4
	10	0.45	0.34	0.03
	25	0.36	0.30	0.01
	50	0.29	0.26	0.01
	100	0.24	0.23	0.01
	200	0.21	0.20	0.01

  .1. The RMSE clearly decreases by increasing the number of training examples N train .

Table 6 .

 6 2: RPA internal energies and enthalpy of adsorption (kcal/mol) of protonated chabazite (HCHAB), CH 4 , and adsorbed system (HCHAB+CH 4 ) as estimated by MLPT for different values of N train . .13 ± 0.21 -20213.93 ± 0.12 -1181.06 ± 0.02 -3.74 ± 0.24 200 -21398.15 ± 0.20 -20213.95 ± 0.12 -1181.06 ± 0.02 -3.74 ± 0.23 probability to be closer in the reduced two-dimensional space used for the plot. In order to generate Fig. 6.2, t-SNE is used to visualize the SOAP descriptors generated by considering only the carbon environment of CH 4 in HCHAB+CH 4 .

	N train	HCHAB+CH 4	HCHAB	CH 4	∆ ads H
	10	-21398.01 ± 0.19 -20213.82 ± 0.13 -1181.05 ± 0.02 -3.74 ± 0.23
	25	-21398.13 ± 0.21 -20213.94 ± 0.13 -1181.07 ± 0.02 -3.72 ± 0.25
	50	-21398.17 ± 0.21 -20213.99 ± 0.13 -1181.06 ± 0.02 -3.72 ± 0.25
	100 -21398

Table 6 .

 6 3: RPA internal energies and enthalpy of adsorption (kcal/mol) of protonated chabazite (HCHAB), CH 4 , and adsorbed system (HCHAB+CH 4 ) as estimated by a brute force application of thermodynamic perturbation theory for a small number of configurations N conf . Differently from Table6.2, the ML approach is not applied in this case, causing thus a dramatic deterioration of convergence.

	N conf HCHAB+CH 4 HCHAB	CH 4	∆ ads H
	10	-21401.99	-20213.75 -1180.97 -7.88
	25	-21398.81	-20214.22 -1181.27 -3.91
	50	-21399.19	-20214.20 -1181.24 -4.34
	100	-21398.06	-20214.56 -1181.32 -2.78
	200	-21398.41	-20214.56 -1181.16 -3.28

Table 6 .

 6 4: PBE+D2, RPA, and experimental enthalpies of adsorption (kcal/mol) of CH 4 and CO 2 adsorbed in protonated (HCHAB) and siliceous (CHAB) chabazite. RPA results are provided as a function of N train .

	∆ ads H	HCHAB+CH 4 HCHAB+CO 2 CHAB+CH 4 CHAB+CO 2
	PBE+D2	-6.09	-9.72	-4.64	-5.14
	RPA (N train = 10)	-3.74 ± 0.23	-8.06 ± 0.44 -3.22 ± 0.19 -4.57 ± 0.25
	RPA (N train = 25)	-3.72 ± 0.25	-7.85 ± 0.37 -3.19 ± 0.19 -4.48 ± 0.22
	RPA (N train = 50)	-3.72 ± 0.25	-7.93 ± 0.38 -3.32 ± 0.19 -4.66 ± 0.21
	RPA (N train = 100) -3.74 ± 0.24	-8.08 ± 0.34 -3.41 ± 0.19 -4.73 ± 0.21
	RPA (N train = 200) -3.74 ± 0.23	-8.01 ± 0.32 -3.36 ± 0.19 -4.65 ± 0.21
	Exp.	-4.06 a	-8.41 b	-4.09 c	-5.02 c
	a Refs. 111
	b Ref. 110
	c Ref. 88
	6.3 Conclusions

in this chapter we discussed the MLPT method that couples machine learning techniques with TPT to efficiently evaluate finite temperature properties using the RPA or other correlated quantum-chemical methods. By computing enthalpies of adsorption of small molecules in zeolites, it is shown that as few as 10 RPA energies used to train a ML model are sufficient to reach a satisfactory level of convergence. This work open the path to a broader use of high-accuracy/high-cost ab initio methodologies to reach a new level of reliability in materials simulations.

  14-17, 2014; Lowscaling and Unconventional Electronic Structure Techniques Conference (LUEST), Telluride (CO, USA), June 1-5, 2014 12. D. Rocca, B. Mussard, G. Jansen, and J. Angyan, Dielectric matrix formulation of correlation energies within the Random Phase Approximation: Inclusion of screened exchange effects, 17th International Workshop on Computational Physics and Materials Science: Total Energy and Force Methods, Trieste, January 15-17, 2015; Green's function methods: The next generation II, Lausanne, May 4-7, 2015 13. D. Rocca, A. Dixit, J. Claudot, and S. Lebègue, A novel implementation to compute MP2 correlation energies without basis set superposition errors and complete basis set extrapolation, APS March Meeting, March 5-9 2018, Los Angeles

TABLE C .

 C 1 : Les energies internes et enthalpie d'adsorption (kcal/mol) au niveau de théorie RPA de la chabazite protonée (HCHAB), de CH 4 et du système adsorbé (HCHAB+CH 4 ) selon les estimations de MLPT pour différentes taille du jeu d'entrainement. .01 ± 0.19 -20213.82 ± 0.13 -1181.05 ± 0.02 -3.74 ± 0.23 25 -21398.13 ± 0.21 -20213.94 ± 0.13 -1181.07 ± 0.02 -3.72 ± 0.25 50 -21398.17 ± 0.21 -20213.99 ± 0.13 -1181.06 ± 0.02 -3.72 ± 0.25 100 -21398.13 ± 0.21 -20213.93 ± 0.12 -1181.06 ± 0.02 -3.74 ± 0.24 200 -21398.15 ± 0.20 -20213.95 ± 0.12 -1181.06 ± 0.02 -3.74 ± 0.23

	N train	HCHAB+CH 4	HCHAB	CH 4	∆ ads H
	10	-21398

The DScribe libraries are available at https://github.com/SINGROUP/dscribe

Teaching Experience

Lecturer at the Université de Lorraine, Nancy, for the classes:

• Mécanique ondulatoire/Wave mechanics (academic year 2012/2013)

Grants

As PI:

• ANR "Jeunes Chercheuses et Jeunes Chercheurs" 2015 (French National Research Agency Young Investigator Grant); Title: Finding Accurate and Scalable Theories for Electronic Correlation (FASTcorrelation). Funding starting January 1st 2016; budget 180000 e (no overhead costs).

• CNRS PEPS "Physique théorique et ses interfaces" 2013; Title: Techniques de théorie quantique des champs appliquées au problème de l'énergie solaire (Quantum field theory techniques applied to the solar energy problem); awarded 4000 e.

• Project Institut Jean Barriol 2013; Title: Techniques de théorie quantique des champs appliquées au problème de l'énergie solaire (Quantum field theory techniques applied to the solar energy problem); awarded 3000 e.

Computer time:

• PRACE 2016; Title: NEWcorrelation; budget 4M core hours.

• GENCI 2015-2019; Title: Modélisation de la structure électronique par calculs ab initio; budget >1M core hours.

Participant in other two projects funded by the ANR (CARAPASS and LIQUI2D).

Professional Services

Publications

Citation metrics available at http://scholar.google.com/citations? user=BzKIlMEAAAAJ