A B S T R A C T

The ability to systematically compute physical properties at chemical accuracy would be of tremendous help for efficient in silico material discovery, drug design or chemical reaction pathways study. While density functional theory has become the workhorse of materials simulations, the quality of results unfortunately often varies depending on the specific choice of the exchange-correlation functional, and this significantly limits the predictive power of this approach. More correlated approaches, such a coupled cluster theory, the random phase approximation, or configuration interaction, provide more reliable results in a systematically improvable way. Their computational cost however render them inaccessible in a lot of applications, and this is particularly true for finite temperature properties, where many samples calculations are necessary to reproduce ensemble averages. Machine learning (ML) on the other, has proven its ability to efficiently reduce computational costs by using knowledge from previous calculations. The aim of this thesis is to contribute to the on-going development of ML assisted methods to obtain chemically accurate calculations. Machine learning perturbation theory (MLPT) is a recent promising approach capable of obtaining chemically accurate finite temperature properties by producing ensemble property estimates for highly expensive computational methods. This is achieved by learning the energy difference between two methods, using few samples from a reference statistics produced from a computationally feasible level of theory, and reweighting the later statistics to the computationally inaccessible target method, with thermodynamic perturbation theory (TPT). In this thesis, the MLPT method is tested against possible limitations of TPT, when a poor overlap between reference and target configurational space produce biased target estimates. A diagnostic test for this problem is proposed, together with a solution in the form of Monte Carlo resampling using the same ML model. MLPT is further applied to compute a coupled cluster estimate of the adsorption enthalpy of carbon dioxide inside a periodic zeolite. Finally, on another topic, a new method is presented, in the framework of configuration interaction (CI), to efficiently iteratively sample new Slater determinants using a generative ML model. This algorithm, called CIgen, is shown to be competitive or outperform other Monte Carlo and ML approaches to the CI problem.

R E M E R C I E M E N T S

Je souhaite ici remercier différentes personnes qui m'ont soutenu, aidé, et appris durant ces dernières années. Dario Rocca tout d'abord qui m'a accueilli sur ce sujet de thèse. Grazie mille per il vostro sostegno, la vostra pazienza, il vostro incoraggiamento e la libertà che mi avete dato. È stato un piacere averla come supervisore. Je remercie Sébastien Lebègue d'avoir accepté de me superviser après le départ de Dario. Merci infiniment pour ta confiance et pour ton aide. À Bastien, Philippe et Aurélia, merci pour tous ces bons moments ici, je suis très heureux de vous avoir rencontré et de pouvoir vous appeler mes amis ; tout comme Maurício muito obrigado por sua ajuda e conhecimento (não sei se posso me considerar um programador, mas agora consigo compilar algumas coisas) e pelos ótimos momentos que passei com você. Elisabeth merci de m'avoir écouté et encouragé, Guillaume pour tes appels et tes conseils. Je tiens à remercier également ma famille (et Léo) pour leur soutien. Léa, mille mercis, tu as tout rendu plus facile.

Figure 1 Jacob's ladder of density functional theory. . . . . . . . . . . Figure 2 Feed forward NN model. . . . . . . . . . . . . . . . . . . . . Figure 3 ∆ -ML: PBE, RPA electronic energies, and their difference for zeolite configurations. . . . . . . . . . . . . . . . . . . . . Figure 4 Restricted Boltzmann Machine. . . . . . . . . . . . . . . . . . Figure 5 Self Consistent construction of the Hartree-Fock wave function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 6 Dissociation curve of the BH molecule using various wave function methods . . . . . . . . . . . . . . . . . . . . . . . . . Figure 7 Self consistent field procedure to find the ground state eigenvectors of the Kohn-Sham Hamiltonian. . . . . . . . . . . . . Figure 8 Illustration of the pseudo-potential method. . . . . . . . . . Figure 9 Dissociation curve of the BH molecule, using various wave function and DFT methods. . . . . . . . . . . . . . . . . . . . Figure 10 Pictorial representation of thermodynamic perturbation theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 11 Smooth Overlap of Atomic Positions descriptor . . . . . . . Figure 12 Kernel Ridge Regression of the sinus function, using a gaussian kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 13 Adsorbed systems studied . . . . . . . . . . . . . . . . . . . . Figure 14 t-SNE representation of the configurational spaces spanned by the different functionals. . . . . . . . . . . . . . . . . . . . Figure 15 Distributions of deviations of SCAN energy from the average value of the corresponding reference calculation . . . . . Figure 16 Radial distribution functions calculated using the SCAN functional for the first and second coordination sphere of O atoms around Si. . . . . . . . . . . . . . . . . . . . . . . . . Figure 17 t-SNE representation of the configurational spaces spanned by the different functionals for the CH 4 @HChab. . . . . . . Figure 18 The unit cell of the system studied in this work, CO 2 in protonated chabazite. . . . . . . . . . . . . . . . . . . . . . . . Figure 19 t-SNE representation of the configurational spaces spanned by the PBE+D2 molecular dynamics (MD) trajectory and the CCSD(T) machine learning Monte Carlo (MLMC) trajectory.

vii

List of Figures

Restricted Boltzmann machine (RBM) as a generative model and its application to the solution of the Schrödinger equation. Figure 22 Pictorial representation of the tower sampling algorithm. . . Figure 23 Effect of the temperature in the generative procedure. . . . . Figure 24 Total energy convergence as a function of the number of iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 25 CIgen dissociation curve of N

chemical accuracy

"A target accuracy must be selected. A model is not likely to be of much value unless it is able to provide clear distinction between possible different modes of molecular behavior.

As the model becomes quantitative, the target should be that data is reproduced and predicted within experimental accuracy. For energies, such as heats of formation or ionization potentials, a global accuracy of 1 kcal/mole would be appropriate." This quote from the 1998 Nobel lecture of John Pople [START_REF] Pople | Nobel Lecture: Quantum chemical models[END_REF] defined the chemical accuracy as an energy resolution goal for computational methods in material science. For example, the ratio between reactants and products of a chemical reaction at equilibrium can be obtained by exponentiation of the ratio of the free energy difference between both states and the thermal energy [START_REF] Bucko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF] (K R→P = exp(-∆A R→P /k B T )), and a variation of 1.4 kcal/mol leads to a variation of ∼ 10.7 in the equilibrium constant K R→P .

While the simplest approximations in quantum chemistry (Hartree-Fock) and in solid state physics (local and semi-local Density Functional Theories (DFT)) are affordable for moderately big systems (up to few thousands atoms, with a numerical cost roughly scaling as O(K 4 ) and O(K 3 ) respectively with basis size [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF]), they cannot be improved in a systematic way, and lack chemical accuracy in a lot of applications [START_REF] Cohen | Challenges for Density Functional Theory[END_REF]. One needs to resort to more involve approximationspost Hartree-Fock methods such as Configuration Interaction, Møller-Plesset Perturbation Theory or Coupled Cluster, or beyond-DFT models such as the Random Phase Approximationbut the later are computationally more demanding and are restricted to smaller systems [START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF]. In particular, while periodic implementations have been recently proposed for condensed-matter applications for second order Møller-Plesset and Coupled Cluster [START_REF] Booth | Towards an exact description of electronic wavefunctions in real solids[END_REF][START_REF] Del Ben | Second-Order Møller-Plesset Perturbation Theory in the Condensed Phase: An Efficient and Massively Parallel Gaussian and Plane Waves Approach[END_REF][START_REF] Dixit | Communication: A novel implementation to compute MP2 correlation energies without basis set superposition errors and complete basis set extrapolation[END_REF][START_REF] Marsman | Second-order Møller-Plesset perturbation theory applied to extended systems. I. Within the projectoraugmented-wave formalism using a plane wave basis set[END_REF][START_REF] Pisani | Periodic local MP2 method for the study of electronic correlation in crystals: Theory and preliminary applications[END_REF], finite temperature results (which require few hundreds of thousands of single point calculations) are still unobtainable for condensed matter systems using those theories. This was conceptualised in 2000 for DFT by Perdew and Schmidt in a famous figure (see the adaptation in Figure 1): the Jacob's ladder of DFT approximations for exchangecorrelation energy [START_REF] Perdew | Jacob's ladder of density functional approximations for the exchange-correlation energy[END_REF]. Perdew and Schmidt. The ladder allows to climb from the mean field description of the Hartree(-Fock) product to chemical accuracy by adding more and more non-locality in the description of the exchange-correlation energy [START_REF] Perdew | Jacob's ladder of density functional approximations for the exchange-correlation energy[END_REF].

In this context Machine Learning methods have seen in the last decade an intensive rise in order to alleviate the computational burden of quantum mechanical methods and make them more affordable [START_REF] Carleo | Machine learning and the physical sciences[END_REF].

prior work

In 1995 the work of Blank et al. introduced the use of a feed-forward neural network (NN) to learn the potential energy surface (PES) of a small molecule (CO) adsorbed on a Nickel surface using data from an empirical force fields, and from density functional theory (DFT) calculations of H 2 desorption from a silicon surface [START_REF] Thomas | Neural network models of potential energy surfaces[END_REF]. A feed-forward NN, pictured in Figure 2, is a graph model with one input layer that encode the data points, and one or more hidden layers that encode high-order correlations between input points. Between two successive layers all vertices are connected. After the hidden layers, the target property is written as a non-linear function of a linear combination of the layers, and is learned using variational optimization algorithm such as gradient-descent [START_REF] David | Information Theory, Inference, and Learning Algorithms[END_REF].

Figure 2: Feed forward NN model consisting on four input neurons, two hidden layers with six neurons each, and a scarlar target y. The hidden neurons take value h j i = f j i (b j i + 6 k=1 a j-1,j k,i h j-1 k ), with biases parameters b j i and weights a j-1,j k,i . The same definition holds for the first hidden layer (fed by the visible input) and the target y. The functions f j i are usually nonlinear such as the sigmoid function.

If reactions coordinate were used as input for the NN, Lorenz et al. performed a similar work in 2004 to study the interaction of H 2 with a Palladium surface, using instead symmetry adapted functions to feed the NN, reducing the need for the model to learn the PES symmetries, and thus reducing the computational cost. With an error comparable to direct DFT computations for a fraction of the computational cost, those models were however reduced to low dimensional systems: the number of weight parameters augment too rapidly with the system size, so as the number of training points required for a proper fit, and the number of iterations necessary in the optimization procedure. Behler and Parinello introduced High-Dimensional Neural Network Potential in a seminal paper in 2007 [START_REF] Behler | Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces[END_REF]. They pursued the aforementioned ideas by defining atomic environments: instead of a single high-dimensional feed-forward NN, they used a architecture where each atom is represented by a single atomic feed-forward NN, and the total energy is given by a sum of all atomic contributions. Thus, the size of the NN model scales linearly with the number of atoms, and the model does not depends on the number and type of atoms, giving it a high flexibility. Each atomic NN describes the configuration of all the neighboring atoms inside a cutoff radius. Symmetry adapted functions representing radial and angular distributions are used in those models. They demonstrated their model by performing molecular dynamics (MD) of periodic Silicon, with 64 atoms in the cell. About ten thousand configurations computed in the local density approximation of DFT were used to train the model, for a negligible validation error (∼ 0.1 kcal/mol). While this approach is able to give large scales molecular dynamics the accuracy of DFT with several order of magnitude of numerical cost avoided, the number of necessary training points remain unattainable for the most accurate quantum chemical methods. Kernel methods, another paradigm of machine learning algorithms, have seen an important parallel development in quantum chemistry and in the solid state. Bartók et al. used Gaussian Processes [START_REF] Edward Rasmussen | Gaussian processes for machine learning[END_REF] to generate inter-atomic potential and learn the PES of bulk Carbon, Silicon and Germanium [START_REF] Bartók | Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons[END_REF]. Rupp et al. predicted atomization energies of organic molecules with Kernel Ridge Regression [START_REF] Rupp | Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning[END_REF]. Databases generation have also been an important step in the process of testing and validating different approaches to learn and predict quantum chemical or solid properties with e.g. the QM9 dataset [START_REF] Ramakrishnan | Quantum chemistry structures and properties of 134 kilo molecules[END_REF] which contain 134 thousand organic molecules with different properties such as atomization energies, dipole moments, polarizability or harmonic frequencies.

Whichever paradigm is used, suitable and efficient representation of atomic configurations are needed. While Behler and coworkers use in their aforementioned model atom-centered symmetry functions [START_REF] Behler | Atom-centered symmetry functions for constructing highdimensional neural network potentials[END_REF], which are sum of radial and angular symmetry functions describing the environment inside a cutoff radius for each atom, other atomic environment descriptions have been proposed. For example Bartók et al. created the Smooth Overlap of Atomic Positions in 2013 [START_REF] Bartók | On representing chemical environments[END_REF], which is a decomposition into spherical harmonics of sum of gaussian centered on atoms inside a cutoff radius. Differently, global representations encode the whole atomic configuration: for instance, the Coulomb matrix depends on pair distances and chemical species [START_REF] Rupp | Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning[END_REF], the Many Body Tensor Representation contain different weighted broadened distance, angle, and dihedral distribution terms [START_REF] Huo | Unified Representation of Molecules and Crystals for Machine Learning[END_REF].

In 2015, Ramakrishnan et al. introduced the idea of ∆-Machine Learning: instead of directly fitting the energy or other property of interest to the configurational space, it can be easier to learn the target value difference as given by two quantum mechanical approximations [START_REF] Ramakrishnan | Big Data Meets Quantum Chemistry Approximations: The ∆-Machine Learning Approach[END_REF]. Considering two different DFT exchange-correlation functional for example, both approximation include the same kinetic term and differ in their treatment of the correlation energy. Learning the difference between the two approximations is easier than directly learn the whole physical description of the target quantum mechanical approximation. In Figure 3 is shown as an example the energies and their difference for several zeolite configurations using the Random Phase Approximation (RPA) and the PBE DFT functional. It is seen how the ∆E curve is smoother and possesses a smaller variation. If several thousands of configurations are necessary to learn the PES of a complex system, this number can be reduced to few tens or one hundred using this technique. Dorner et al. used in 2018 thermodynamics perturbation theory to reweight an ab initio molecular dynamics performed using semi-local DFT to predict the melting temperature of Silicon at RPA level of theory [START_REF] Dorner | Melting Si: Beyond Density Functional Theory[END_REF]. In a similar work, Rocca et al., in 2019, used thermodynamics perturbation theory to reweight ab initio molecular dynamics performed using semi-local DFT onto RPA averages to compute the adsorption enthalpy of carbon dioxide and methane into a zeolite [START_REF] Rocca | Bridging molecular dynamics and correlated wavefunction methods for accurate finite-temperature properties[END_REF]. In thermodynamic perturbation theory, one can expresses the canonical ensemble average of one Hamiltonian as a function of averages in another Hamiltonian, if the difference between both Hamiltonian is known [START_REF] Chipot | Free energy calculations: theory and applications in chemistry and biology[END_REF]. Given the computational cost of RPA calculations, they extracted few tens of configurations on which to compute RPA energy and perform the perturbation theory, and were able to recover chemical accuracy with respect to the experiment. The same year, this work was continued by Chehaibou et al.: a delta machine learning model was devised to learn, from those selected configurations, the energy difference between the DFT and RPA energy using a Kernel Ridge Regression model, thus allowing to perform the thermodynamic perturbation theory on the full trajectory and reducing the error bar [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]. This method was coined Machine Learning Perturbation Theory (MLPT), and has been applied to compute enthalpy of adsorption of small molecules into chabazite [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF], and free energies of activation of the proton exchange reaction in protonated chabazite [START_REF] Bucko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF].

The works mentioned until now assume that quantum mechanical calculations are possible prior to create a model to train on those. It could be advantageous to use ML directly during the process of solving the Schrödinger equation. In 2017, Carleo and Troyer [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF] introduced the idea of Neural Quantum States (NQS). Trial functions have been commonly used in Variational Quantum Monte Carlo (QMC), since their introduction by Jastrow in 1955 [START_REF] Jastrow | Many-Body Problem with Strong Forces[END_REF]: one chooses, in front of a mean-field solution such as Hartree-Fock, a variational ansatz (Jastrow factor) to express the correlation of a quantum mechanical system. The ansatz is subsequently optimized in a variational way, notably using Monte Carlo sampling [START_REF] Toulouse | Chapter Fifteen -Introduction to the Variational and Diffusion Monte Carlo Methods[END_REF]. The Jastrow factor is typically a non-linear function of the pair distances in the problem. While these scheme allows for very accurate calculations (often used as benchmark [START_REF] Jeffrey | Benchmark quantum Monte Carlo calculations[END_REF]), there is a certain degree of arbitrary in the choice of the Jastrow function, and a lack of flexibility in the ansatz. In this context, Carleo and Troyer suggested to use the representative power of neural networks [START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF][START_REF] Roux | 1 Representational Power of Restricted Boltzmann Machines and Deep Belief Networks[END_REF] instead to parametrize the wave function to be optimized. More specifically, they chose a Restricted Boltzmann Machine (RBM) as illustrated in Figure 4.

The associated wave function amplitude for a configuration v of binary input {v i } is written as

Ψ(v) = {h i } e a T v+b T h+h T Wv (1)
and the optimization is performed using a Variational QMC method, the Stochastic Reconfiguration algorithm [START_REF] Sorella | Weak binding between two aromatic rings: Feeling the van der Waals attraction by quantum Monte Carlo methods[END_REF]. Introduced for spins models such as the Antiferromagnetic Heisenberg model were a better convergence compared to Jastrow functions was shown, this work was further brought to electronic structure problems by Choo et al. [START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF]: by mapping the electronic structure Hamiltonian into a spin Hamiltonian (using among other the Jordan-Wigner transformation [START_REF] Jordan | Über das paulische äquivalenzverbot[END_REF]), recovering the Configuration Interaction (CI) framework [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]. They studied total electronic energy calculations of small molecules (up to fourteen electrons) in the STO-3G and 6-31G basis sets. While able to converge most of the calculations to chemical accuracy in the STO-3G basis, this was not feasible in the 6-31G basis set. It appears that it is yet difficult to elaborate efficient schemes to optimize a neural network -moreover with complex weights -on the CI problem. The sampling is complicated by a predominant weight of the most important configurations in the wave function, resulting in a problematic slow-down of the algorithm for the nonminimal basis set 6-31G. Barett et al. addressed that issue using an autoregressive neural network and could obtained chemically accurate calculations for basis set up to twenty-eight spin orbitals [START_REF] Thomas D Barrett | Autoregressive neural-network wavefunctions for ab initio quantum chemistry[END_REF].

Another difficulty for the RBM or other simple models if the highly complex nature of the correlations between occupation of orbitals in the configurations. By simple model is understood a model with few parameters, e.g. a similar number of hidden and input variables, as the representative power of neural networks grows with the number of those hidden variables [START_REF] Roux | 1 Representational Power of Restricted Boltzmann Machines and Deep Belief Networks[END_REF]. One solution is to use then deep neural networks, by, for example, stacking RBM together, one hidden layer being an input for the next hidden layer (deep Restricted Boltzmann Machine [START_REF] Le | Representational power of restricted Boltzmann machines and deep belief networks[END_REF]). Real space wave functions have been encoded as NQS in such deep learning architectures in the last four years [START_REF] Han | Solving many-electron Schrödinger equation using deep neural networks[END_REF][START_REF] Hermann | Deep-neural-network solution of the electronic Schrödinger equation[END_REF][START_REF] Kessler | Artificial neural networks as trial wave functions for quantum monte carlo[END_REF][START_REF] Pfau | Ab initio solution of the many-electron Schrödinger equation with deep neural networks[END_REF], outperforming previous state of the art Variational QMC ansatzes for systems up to thirty electrons. Truncated CI wave functions remains however a highly valuable method to achieve systematic accuracy and calibrate other methods [START_REF] Sherrill | The Configuration Interaction Method: Advances in Highly Correlated Approaches[END_REF]. Since the CI space grows combinatorially with the system size, prior CI algorithms usually select the most important configurations using perturbation theory [START_REF] Huron | Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions[END_REF] or random sampling [START_REF] Greer | Monte Carlo configuration interaction[END_REF] before diagonalization of the resulting Hilbert subspace in an iterative way. Machine learning has been proposed as a mean to replace perturbation theory as selection criterion of the important configurations by Coe et al. in 2018 [36], using a regressive neural network to learn the important configurations.

outline

The aim of this thesis is to contribute to the development of novel machine learning methodologies in order to bring down the computational cost needed in obtaining chemically accurate properties. Machine Learning Perturbation theory is further developed and tested, using various DFT functionals, the Random Phase Approximation, Møller-Plesset and Coupled Cluster theories. On another topic, a generative machine learning model is developed to sample new Slater determinants in the Configuration Interaction framework. The next two chapters are dedicated to understand those diverse quantum mechanical and machine learning methods. The second chapter is a review of the quantum many-electron problem, with an introduction to wavefunction methods and density functional theory. It also introduces ab initio Molecular Dynamics and Thermodynamic perturbation Theory. The third chapter explains the theoretical basics of the machine learning models used in this thesis. The results of this thesis can be found in the fourth, fifth and sixth chapters. In Chapter four, the Machine Learning Perturbation Theory is continued to address the possible configurational mismatch inherent to thermodynamic perturbation theory, by studying total electronic energies and adorption enthalpies with various DFT functionals. A diagnosis test is presented to detect the existence of an overlap problem, in which case a solution is proposed to overcome it, in the form of a Monte Carlo resampling of the biased trajectory, using the original machine learning model. Chapter five carries on with this work by a calculation of the enthalpy of adsorption of carbon dioxide in a zeolite at the coupled cluster single double and perturbative triple (CCSD(T)) level of theory, bringing the first finite temperature simulation of a periodic material at the Coupled Cluster level of theory to our knowledge. Chapter 6 shows how a generative machine learning model can be used to efficiently sample new Slater determinants in an iterative configuration interaction algorithm, improving upon random sampling. The seventh and last chapter discuss the results of this thesis and conclude.

E L E C T R O N I C S T R U C T U R E M E T H O D S

introduction

The aim of this chapter is to introduce the various quantum mechanical methods considered in this thesis. It is organized as follows: the first section is dedicated to wave function-based methods. It begins by an exposition of the many electron problem, and shows its general solution as a linear combination of Slater determinants. The mean field Hartree-Fock Slater determinant is presented as a first level of approximation, together with a definition of the correlation energy. The Configuration Interaction framework is then presented as a straighforward way to recover correlation. Møller-Plesset Perturbation Theory and Coupled Cluster theory are introduced as alternative methodologies, without some of the shortcomings of the Configuration Interaction method. A brief introduction to basis sets in quantum chemistry is given before a comparison of those different methods with a dissociation curve of the Hydrogen-Bore molecule. In the second section, the basis of density functional theory is introduced, with the Hohenberg and Kohn theorems, and the Kohn-Sham self consistent scheme. An exact formulation of the exchange-correlation energy is given, with a presentation of various levels of approximations to the later, ending with the Random Phase Approximation. The previous Hydrogen-Bore dissociation curve is continued with those methods. Finally in the third and last section, finite temperature properties are considered, with an introduction to the adiabatic approximation and the Born-Oppenheimer molecular dynamics, and Monte Carlo sampling. That section ends with a presentation of thermodynamic perturbation theory.

wave function-based methods

2.2.1

The many electron problem

The complete Hamiltonian

Let us consider a non-relativistic system of N electrons of mass m and N α nuclei of masses {M 1 , . . . , M N α } and atomic number {Z 1 , . . . , Z N α } interacting via the Coulomb potential. The corresponding Hamiltonian for this system is

Ĥ = Te + Tn + Vee + Vnn + Ven (2) with • Te = N i=1 -h 2 ∇2 r i 2m • Tn = N α i=1 -h 2 ∇2 R i 2M i • Vee = N i=1 N j>i e 2 4πϵ 0 |r i -r j | • Vnn = N α i=1 N α j>i Z i Z j e 2 4πϵ 0 |R i -R j | • Ven = N i=1 N α j=1 -Z j e 2 4πϵ 0 |r i -R j |
where the {r i }, {R i } are the positions of the electrons and nuclei. The wave function Ψ({r i }, {R i }) of this system is in position basis a function of all electrons and nuclei positions, and a solution of the time-independent Schrödinger equation

Ĥ |Ψ⟩ = E |Ψ⟩ (3)
Additionally, our system is made of half-integer spin particles (fermions), and the wave function must be antisymmetric under the exchange of any two indistinguishable particles.

Electronic Schrödinger equation

In the many electron problem we are interested in the eigenvectors ψ e i and eigenvalues E e i of the electronic Hamiltonian.

Ĥe = Te + Vee + Vext (4)
with Vext = Ven . When solving the electronic Schrödinger equation, the nuclei positions in Vext are fixed parameters and the set of eigenvalues {E e i (R j )} where R j are the nuclei positions is called the potential energy surface of the problem. This is the Born-Oppenheimer approximation, that will be studied more in depth in Section 2.4.1.

General form of the total wave function: Slater determinants

Consider a single electron system, with wave function |ψ⟩ = i c i |ϕ i ⟩, where the set of {|ϕ i ⟩} is a complete orthonormal basis, and i completely specify both space and spin. If we now add another electron, we can describe the new wave function as

|ψ(1, 2)⟩ = i j a ij |ϕ i ⟩ ⊗ ϕ j = i j a ij ϕ i , ϕ j = i j a ij ϕ i (1)ϕ j (2) .
The wave function of fermions being antisymmetric under the exchange of two particles, this translate into |ψ(2, 1)⟩ = i j a ij ϕ j , ϕ i = i j a ij ϕ i (2)ϕ j (1) = -|ψ(1, 2)⟩, so that we must impose a ij = -a ji for all i, j. This yields

|ψ(1, 2)⟩ = i j>i a ij ϕ i (1)ϕ j (2) -ϕ i (2), ϕ j (1) (5)
If we add another electron, the total wave function is

|ψ(1, 2, 3)⟩ = i j k a ijk ϕ i (1)ϕ j (2)ϕ k (3) . ( 6
)
The antisymmetry of |ψ(1, 2, 3)⟩ impose that a ijk change sign under any transposition of two indexes among i, j, k since for example |ψ(1, 2, 3)⟩ = -|ψ(2, 1, 3)⟩ = + |ψ(2, 3, 1)⟩, we will have then

a σ(ijk) =            a ijk if σ(ijk) is an even permutation of i, j, k -a ijk if σ(ijk) is an odd permutation of i, j, k 0 otherwise. (7)
or, using the Levi-Civita symbol,

a σ(ijk) = a ijk ϵ σ(ijk) (8)
and

|ψ(1, 2, 3)⟩ = i j>i k>j a ijk σ(ijk) ϵ σ(ijk) ϕ σ i (1)ϕ σ j (2)ϕ σ k (3) (9)
where σ i gives the corresponding value of i in the permutation σ(ijk). For example, if σ(ijk) = jik, then σ i = j.

For a squared 3 × 3 matrix M with elements M αβ , the determinant can be defined as

det(M) = σ(123) ϵ σ(123) M 1σ 1 M 2σ 2 M 3σ 3 .
We see that we can use the notation of determinants to define

|ij • • • n⟩ = σ(ij•••n) ϵ σ(ij•••n) ϕ σ i (1)ϕ σ j (2) • • • ϕ σ n (N) ≡ ϕ i (1) ϕ j (1) ... ϕ n (1) ϕ i (2) ϕ j (2) ... ϕ n (2) • • • • • • • • • • • • ϕ i (N) ϕ j (N) ... ϕ n (N) (10)
For N electrons, we will obtain the total wave function as

|Ψ(1, 2, • • • , N)⟩ = i j>i • • • n>m a ijk•••lmn |ijk • • • lmn⟩ (11)
The kets |ijk • • • lmn⟩ are called Slater determinants [START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF][START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF] and the total wave function is written as a linear combinations of such Slater determinants. They form a convenient basis of the problem as they directly encode the anti-symmetry requirement.

Second quantization

It is sometime convenient to work in the so-called second quantization formalism where the Hamiltonian reads

Ĥe = ij t ij ĉ † i ĉj + ij v ext ij ĉ † i ĉj + ijkl v ijkl ĉ † j ĉ † i ĉk ĉl ( 12
)
with

t ij = ⟨ϕ i | Te |ϕ j ⟩ (13)
v ext ij = ⟨ϕ i | Vext |ϕ j ⟩ (14)
v ijkl = ⟨ϕ i ϕ j | Vee |ϕ k ϕ l ⟩ (15)
where Te , Vext and Vee are the one-particle kinetic, external and the two-particles Coulomb operators respectively (without sum as in eq. ( 2)). ĉi , ĉ † i are the annihilation and creation operators, obeying (for fermions) the anti-commutation rules

{ ĉi , ĉ † j } = δ ij (16)
{ ĉi , ĉj } = 0 (17)

{ ĉ † i , ĉ † j } = 0 (18)
Those operators act on Slater determinants in the occupation basis

ĉ † i |n 1 , • • • , n i , • • • ,⟩ = (-1) n i-1 k=n 1 n k (1 -n i ) |n 1 , • • • , n i + 1, • • • ,⟩ (19)
ĉi

|n 1 , • • • , n i , • • • ,⟩ = (-1) n i-1 k=n 1 n k n i |n 1 , • • • , n i -1, • • • ,⟩ (20)
where

|n 1 , n 2 , • • • ,⟩ = |1, 0, 0, 1, 1, 0, • • • ,⟩ corresponds to the 3-electron state |1, 4, 5⟩
in eq. [START_REF] Albert P Bartók | On representing chemical environments[END_REF]. We can see then that ĉ1

|1, 4, 5⟩ = |4, 5⟩, ĉ1 |4, 5⟩ = 0, ĉ † 1 |1, 4, 5⟩ = 0 and ĉ † 1 |4, 5⟩ = |1, 4, 5⟩.
It is possible to change the one-electron basis

|ϕ k ⟩ by |x⟩ = k |ϕ k ⟩ ⟨ϕ k |x⟩ = k ϕ * k (x) |ϕ k ⟩ = k ϕ * k (x) ĉ † k |0⟩
, where |x⟩ specifies a new basis, e.g. in position basis |x⟩ = |r, σ⟩ with σ = ± 1 2 the spin projection. This defines the fields creation operators

ψ † (x) = k ϕ * k (x) ĉ † k . ( 21
)
The Hamiltonian can again be rewritten in this basis [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF] with

Ĥe = dx ψ † (x)( Te (x) + Vext (x)) ψ(x) + 1 2 dxdx ′ ψ † (x) ψ † (x ′ ) Vee (x, x ′ ) ψ(x) ψ(x ′ ) (22)

Slater-Condon rules

Using this formalism, we can more easily evaluate the Hamiltonian matrix elements ⟨ijk •

, n 2 , • • •| ĉ † j ĉ † i ĉk ĉl |n 1 , n 2 , • • •⟩ which are only non-zero if (i = k, j = l) or (i = l, j = k).
⟨A| Ĥe |A⟩ = N i=1 (t ii + v ext ii ) + 1 2 N i=1 N j=1 (v ijij -v ijji ) (23)
⟨A| Ĥe |A r a ⟩ = (t ar + v ext ar ) + N i=1 (v airi -v aiir ) (24)
⟨A| Ĥe |A rs ab ⟩ = v abrs -v absr (25)
To summarize, the many-problem consist in solving the eigenvalue problem defined by the time-independent Schrödinger equation, for which there exist no analytical solution. The exponential growth of the Hilbert space size with the number of electrons prevents a direct numerical resolution and we need to resort to approximations. The first of such consist in reducing the one-electron basis size by any kind of numerical cut-off. For a real-space representation (ϕ i (r) = δ(r -r i )σ i , this means introducing a minimal numerical space resolution. Quantum chemistry use atomic orbitals based on the solution of the Schrödinger equation for the Hydrogen atom, and the cut-off here corresponds to the size (the number and type) of the basis set used. An additional type of approximation is to neglect terms in the expansion of the full wave function, retaining only the most important. We will now describe the Hartree-Fock procedure, that express the wave function as a single of such terms.

2.2.2

Hartee-Fock and the correlation energy

Hartree-Fock method

Consider the first element |123 • • • N⟩ of the expansion in eq. [START_REF] Bartók | Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons[END_REF]. Remember that we had defined the single electron basis from |ψ⟩ = i c i |ϕ i ⟩. This basis is not unique and we can obtain all other basis through unitary transformations of the coefficients c i . By doing so, the Hamiltonian matrix element

⟨123 • • • N| Ĥe |123 • • • N⟩ will vary.
By minimizing this energy with respect to the c i , we will obtain the Hartree-Fock Slater determinant. It is the best single Slater determinant approximation to the true ground state N electrons wave function. 1Let construct the following Lagrangian functional with the basis orthonormality constraint:

L[{ϕ i }] = E[{ϕ i }] - ij ϵ ij ⟨ϕ j |ϕ i ⟩ -δ ij . ( 26
)
where E[{ϕ i }] is given by eq. [START_REF] Bucko | Manybody dispersion corrections for periodic systems: an efficient reciprocal space implementation[END_REF]. By cancelling the functional derivatives ∂L/∂ϕ i , one obtains the following eigenvalue problem

Fϕ i (x) =   h(x) + N j=1 (J j (x) -K j (x))   ϕ i (x) = ϵ i ϕ i (x) (27)
with the one-electron Coulomb and exchange operators acting on a state ϕ i (x) as

J j (x)ϕ i (x) = dx ′ |ϕ j (x ′ )| 2 r ij ϕ i (x) (28)
K j (x)ϕ i (x) = dx ′ ϕ * j (x ′ )ϕ i (x ′ ) r ij ϕ j (x) (29)
and the one-electron operator ĥ = Te + Vext . 2 The total operator F is called the Fock operator. This procedure yields a new one-electron basis written as a linear combination of the old one. For molecular systems (by opposition to periodic systems), the one-electron orbitals ϕ i (x) are called molecular orbitals.

The one-electron Coulomb operator acts on each point by averaging the Coulomb repulsion due to the other orbitals. Note the non-locality of the exchange operator whose action on ϕ i at point x depends on its value at all points.

Since both contain the orbitals in their definition, one needs to resort to a self consistent scheme in order to find the solutions of this problem, as illustrated in Figure 5.

The one-electrons energies ϵ i are eigenvalues of the Fock operator:

ϵ i ≡ ⟨ϕ i | F|ϕ i ⟩ = t ii + v ext ii + N j=1 (v ijij -v ijji ). (30)
and the Hartree-Fock state |Φ⟩ = |12 • • • N⟩ has the energy

E HF = N i=1 ϵ i - 1 2 N i=1 N j=1 (v ijij -v ijji ) (31)
The Fock operator can be rewritten as

F = Te + Vext + VHF ( 32
)
where VHF is to be seen as an effective one-particle potential, averaging on the effect of the occupation of the N first orbitals.

In this thesis, we will always consider closed-shell systems, i.e. systems with an equal number of spin up and down electrons (as opposed to open-shell systems).

No Initial orbitals

Fock matrix constuction

Fock matrix diagonalization

New orbitals

Self-consistent ? For such systems, it is possible to consider only spatial atomic orbitals in the minimization yielding the Hartree-Fock state: the same spatial orbital hold two electrons with opposite spins. This procedure is called Restricted Hartree-Fock (RHF). The Unrestricted Hartree-Fock (UHF) method consider spin-orbitals with a different spatial part. Finally, let us note that the molecular spin-orbitals are in practice written as linear combination of atomic spin-orbitals in the form ϕ i (x) = α c αi f α (x), as mentioned earlier. The most computationally intensive step of the Hartree-Fock method is to compute all v ijkl integral terms of eq. ( 15), making it quartic with the atomic basis set size. Those atomic basis sets will be introduced in subsection 2.2.7.

Correlation energy

The correlation energy E c is defined as the difference between the ground state energy E and the Hartree-Fock energy E HF :

E c = E -E HF (33)
Note that in practice, as we are using finite basis sets, the ground state energy E must be understood with respect to that basis set. The correlation energy is commonly decomposed into dynamic correlation near the equilibrium geometry and static correlation near dissociated geometry. The term "dynamic" refers to the lack of relaxation in the Hartree-Fock approximation where the electrons only feel an average potential. The non-dynamic correlation, hence "static", refers to the multi-configurational character of the wave function near dissociated geometries. Considering for example the hydrogen molecule H 2 with only two basis functions, one 1s function attached on each hydrogen. The wave function can be written as a sum of the Hartree-Fock state and a doubly excited determinant. At equilibrium geometry, the Hartree-Fock state represent the largest contribution to the wave function among the two, while in the dissociation limit, the Hartree Fock determinant is degenerated with the doubly excited determinant [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF].

Size-consistency and size-extensivity

A method is said to be size-consistent if, when considering two sub-systems A and B without interaction between each other, the energy of the total system is given by the sum of the energy of the sub-systems, as it should. Moreover, the total wave function should be the antisymmetrized product of the wave functions of A and B. Typically A and B could be the atoms forming a molecule, which become non-interacting in the dissociation limit where the pair distances between atoms goes to infinity. The Fock matrix is constituted of elements F ij = ⟨ϕ i | F|ϕ j ⟩. In the case where A and B are two closed-shell systems, the Fock matrix elements F ij where i (j) are localized on A (B) will be zeros in the dissociation limit of the RHF procedure. Thus the Fock matrix is block diagonal already from the first iteration, and the Hartree-Fock determinant is the product of the two subsystems Hartree-Fock determinants, hence its energy is size-consistent. If A and B are open-shell however, while their sum is not, only the UHF method can recover this block diagonal form. The RHF method is then non size-consistent in general, while the UHF method is.

Size-extensivity is a property obeyed by physical systems such that the energy of M monomers must scales linearly with M when M → ∞ [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]. (This does not mean that the energy is M times the energy of a single monomer). As such, it is necessary for periodic systems to be described by an approximation that respect this constraint. Obviously a non size-consistent method is not size-extensive and truncated CI methods, which lack size-consistency as discussed in the next section, cannot be used for periodic systems.

In the so-called Post Hartree-Fock methods, after the Hartree-Fock Slater determinant has been constructed, we try to recover (part or all) of the correlation energy.

Configuration Interaction

The Configuration Interaction is the most straightforward method to do so, given our previous definitions. Indeed, we can rewrite eq. ( 11) as such substitutions. Thus, FCI scales combinatorially with the basis set size, and only small molecules and small basis sets are attainable. Let consider the CID wave function, made of the Hartree-Fock state plus all doubly substituted states. If we attempt to describe two non-interacting subsystems A and B as defined earlier, we cannot write the CID wave function of the total system as the product of CID wave functions of each subsystem: the former will contain products of singly substituted determinants, while the later will hold quadruply substituted determinants (product of doubly substituted determinants). This is a general feature of the truncated CI method to lack size-consistensy. Of course the FCI wave function, being exact, is size-consistent. 3

|Ψ⟩ = |Φ⟩ + a r c r

Perturbation theory expansion of the correlation energy

We will first describe a formal expansion of the ground state energy in terms of non-degenerate Rayleigh-Schrödinger perturbation theory [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF] (RSPT). In this framework, one expresses the Hamiltonian of interest as

Ĥ = Ĥ0 + λ V (35) Ĥ |Ψ i ⟩ = E i |Ψ i ⟩ (36)
where the chosen unperturbed Hamiltonian Ĥ0 has a known non-degenerate eigenspace: Ĥ0

Φ (0) i = E (0) i Φ (0) i (37)
and the perturbation V is turned on with λ where 0 ⩽ λ ⩽ 1. The eigenstates and eigenvalues are expanded as

|Ψ i ⟩ = ∞ k=0 λ k Φ (k) i (38)
E i = ∞ k=0 λ k E (k) i (39)
Using intermediate normalisation (⟨Φ

i |Ψ i ⟩ = 1), insertion of eq. ( 38) and eq. (39) in eq. [START_REF] Coe | Machine learning configuration interaction[END_REF], identifying power of λ, and insertion of the Ĥ0 eigenstates, one finds for the three first orders

E (1) i = ⟨Φ (0) i | V|Φ (0) i ⟩ (40)
E (2) i = k̸ =i ⟨Φ (0) i | V|Φ (0) k ⟩ 2 E (0) i -E (0) k (41)
E (3) i = k̸ =i l̸ =i ⟨Φ (0) i | V|Φ (0) k ⟩⟨Φ (0) k | V|Φ (0) l ⟩⟨Φ (0) l | V|Φ (0) i ⟩ (E (0) i -E (0) k )(E (0) i -E (0) l ) (42) -E (1) i k̸ =0 |⟨Φ (0) i | V|Φ (0) k ⟩| 2 (E (0) i -E (0) k ) 2
If the first order term E

(1) i is size-consistent, we might suspect the third order term not to be: as E (1) i = e(A) + e(B) for two independent systems A and B, we will find product of A term and B terms in the second sum of E

(3) i . However it turns out that those size inconsistent terms are cancelled by some terms in the first sum. This is a general property and all terms in the perturbative expansion are size-consistent [START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF][START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]. Nonetheless, it is impractical that those size inconsistent terms remains in the formal expansion, if their cancellation in numerical schemes is not totally recovered. An alternative and equivalent derivation of the energy expansion can be found using time-dependent perturbation theory [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Jørgensen | Second Quantization-based Methods in Quantum Chemistry[END_REF]. Using Gell-Mann and Low theorem, the correlation energy is rewritten as [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Gell-Mann | Bound States in Quantum Field Theory[END_REF] E -

E (0) = Φ (0) V(0) Û(0, -∞) Φ (0) Φ (0) Û(0, -∞) Φ (0) (43)
where U I (t, t ′ ) is the evolution operator and interaction picture is used:

U I (t, t ′ ) = T exp -i t t ′ VI (t)dt (44)
T designates the time ordering operator, that reorders terms upon which it act such that the latest is to the left, with a global sign given by the parity of the number of permutations needed. The terms arising from the exponential serie are further decomposed through Wick theorem and finally can be drawn as Goldstone diagrams (time-ordered Feynman diagrams) [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF]. The numerator is then seen to factorize into a product of two sums, one of which cancels the denominator, with only linked diagrams surviving (this is one formulation of the linked-cluster theorem), resulting upon further time-integration in the following formula for the correlation energy, due to Goldstone [START_REF] Goldstone | Derivation of the Brueckner many-body theory[END_REF] E -

E (0) = ⟨Φ (0) | V ∞ n=0 1 E (0) -Ĥ0 V n |Φ (0) ⟩ connected (45)
The subscript connected means that the sum only involves diagrams that are linked, that is, topologically connected diagrams, which correspond to terms that are sizeconsistent. the diagrams are constructed by successively (n times at order n) applying the perturbation, which creates couples of particle-holes pairs for the Coulomb interaction, before the last V returns the pairs to the ground state [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF].

Møller-Plesset Perturbation Theory

Møller-Plesset Perturbation Theory (MPPT) is a RSPT (or equivalently the Goldston expansion of eq. ( 45)) where H 0 is chosen as the Hartree-Fock operator F (which is the sum over electrons of one-particle Fock operators), and the perturbation is then Ĥe -F = Vee -VHF . Thus the corresponding eigenstates are Slater determinants. The zeroth and first orders of the ground state energy are

E (0) = ⟨Φ| F|Φ⟩ = N i=1 ϵ i (46)
E (1) = ⟨Φ| Vee -VHF |Φ⟩ = - 1 2 N i=1 N j=1 (v ijij -v ijji ) (47)
yielding together the Hartree-Fock energy. In second order, only the doubly substituted Slater determinants contribute due to Brillouin theorem. 4 If a, b label the occupied molecular orbitals in the Hartree-Fock state, and r, s the unoccupied ones, we obtain

E (2) = 1 4
a,b,r,s

|v abrs -v absr | 2 ϵ a + ϵ b -ϵ r -ϵ s (48)
or

E (2) = 1 2
a,b,r,s

|v abrs | 2 ϵ a + ϵ b -ϵ r -ϵ s - 1 2 a,b,r,s v * abrs v absr ϵ a + ϵ b -ϵ r -ϵ s ( 49
)
where the correction has been separated into direct and exchange contributions, and the second order MPPT energy, or MP2 energy, reads

E MP2 = E HF + E (2) (50)
In higher order, it is easier to use the diagrammatic expansion of eq. ( 45). The successive diagrams of the expansion can be constructed using the following rules [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]:

1. The Coulomb interaction is represented by an interacting wiggly line between two vertices, with a pair of directed lines (one in and one out) at each vertex 2. At order n, stack n such pieces vertically and connect them with the straight line in all possible ways upon deformation of the lines and horizontal reflection(s) of the interaction(s), allowing only connected diagrams, and forbidding self interacting vertices (Brillouin theorem).

Label each directed line with an index

The associated value of the diagram is computed by

• A factor v klmn = ⟨left in; right in| Vee |left out; right out⟩ per interaction line,

• Count the lines up (particle) and down (hole) between each interaction line and associate a factor

( h ϵ h -p ϵ p ) -1
• Count the number of holes lines h and the number of closed loops l and associate a factor (-1) h+p

• A factor 1 2 for diagrams with horizontal symmetry

• Sum the resulting term over each index 4 ⟨Φ|H|Φ r a ⟩ = h ar + i (v airi -v aiir ) constitutes an off-diagonal element of the Fock matrix.

For example the two second order diagrams (direct and exchange terms above) are

r a s b = (-1) 2+2 1 2 a,b,r,s v * abrs v abrs ϵ a + ϵ b -ϵ r -ϵ s (51) r a s b = (-1) 2+1 1 2 a,b,r,s v * absr v abrs ϵ a + ϵ b -ϵ r -ϵ s ( 52
)
And one of the twelve third order diagrams:

r a s b t c = (-1) 3+3 a,b,c,r,s,t v acrt v btsc v rsab (ϵ a + ϵ c -ϵ r -ϵ t )(ϵ a + ϵ b -ϵ r -ϵ s ) (53)
MPPT is thus size-consistent at all order for closed-shell fragment systems if RHF is used as unperturbed wave function. Convergence of the serie should be considered carefully: first of all the difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied one (LUMO) should be big enough to avoid the divergence of all terms (thus homonuclear molecules will have a divergent MPPT energy near dissociation, see also ref. [START_REF] Grüneis | Second-order Møller-Plesset perturbation theory applied to extended systems. II. Structural and energetic properties[END_REF] for the MP2 divergence in metals). Moreover, the radius of convergence of the serie with λ can be inferior to one. The convergence is usually assured near equilibrium positions but not near dissociation. The MP2 method represent a net improvement over Hartree-Fock near equilibrium, for a cost scaling in M 5 with the number of basis functions [START_REF] Cremer | Møller-Plesset perturbation theory: from small molecule methods to methods for thousands of atoms[END_REF]. The higher orders are less used as the ratio of correlation energy recovered over computational cost is less appealing [START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF]. As a final note, let us mention that since perturbation theory is not variational, the MPPT energy will oscillate with the order of perturbation and can occasionnaly be lower than the FCI energy.

Coupled Cluster

We will finish our presentation of wave function methods with Coupled Cluster (CC) theory. The starting point is the CC exponential ansatz for the ground state wave function:

|Ψ CC ⟩ = e T |Φ⟩ (54)
with |Φ⟩ the Hartree-Fock wave function and the CC operator defined by

T = n Tn , (55)
the Tn being n-substitutions operators. For example we have

T1 = a r t r a ĉ † r ĉa (56) T2 = 1 4 ab rs t rs ab ĉ † r ĉ † s ĉb ĉa ( 57
)
that generate respectively all singly and doubly excited Slater determinant from the reference |Φ⟩, t r a and t rs ab being the single and double CC amplitudes. Expansion of the exponential in eq. ( 54) yields:

e T = (1 + T + 1 2 T 2 + T 3 + • • • ) = (1 + T1 + T2 + • • • + 1 2 T 2 1 + 1 2 T1 T2 + 1 2 T2 T1 + • • • ) (58)
that can be collected into n-fold excitation operators Ĉn :

Ĉ1 = T1 (59) Ĉ2 = 1 2 T 2 1 + T2 (60)
Ĉ3 = 1 3! T 3 1 + 1 2 T1 T2 + T3 (61)
and ( 54) can be rewritten as

|Ψ CC ⟩ = |Φ⟩ + a r c r a |Φ r a ⟩ + a<b r<s c rs ab |Φ rs ab ⟩ + a<b<c r<s<t c rst abc Φ rst abc + • • • (62)
which is exactly the FCI wave function [START_REF] Chmiela | Machine learning of accurate energy-conserving molecular force fields[END_REF]. The coefficients are here given by

c r a = t r a ( 63
)
c rs ab = t rs ab + t r a * t s b ( 64
)
c rst abc = t rst abc + t r a * t st bc + t r a * t s b * t t c ( 65
)
where * is an antisymmetric product with respect to exchange of occupied orbitals and exchange of unoccupied orbitals. (e.g. t r a * t s b = t r a t s b -t r b t s a ). Since the exponential ansatz gives back the FCI wave function with the complication of an additional non-linear parametrization of the amplitudes, it has little interest as it is. Truncation of the CC ansatz is however an improvement over truncation of the CI wave function [START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF]. The truncation is done on the CC operator, and for example the Coupled Cluster Single Double (CCSD) ansatz is given by retaining only T1 and T2 in T . Contrary to the CISD case, triple and higher excitations are still present through product of the single and double substitution operators, while the coefficients of those are expressed as products of the single and double CC amplitudes t r a , t st ab only. The correlation energy is readily found by projecting the Schrödinger equation for |Ψ CC ⟩ on |Φ 0 ⟩ :

E⟨Φ 0 |Ψ CC ⟩ = ⟨Φ 0 |H|Ψ CC ⟩ (66)
E = ⟨Φ 0 |H|Φ 0 ⟩ + a<b r<s c rs ab ⟨Φ 0 |H|Φ rs ab ⟩ (67)
so that

E c = 1 4
a,b,r,s

(t rs ab + 2t r a t s b )(v abrs -v absr ). (68)
Note that this equation is valid for any truncation order of the CC ansatz, the correlation energy is expressed only as function of the single and double CC amplitudes. Those are obtained by projecting the CC Schrödinger equation on single, double, triple, ..., n-substituted Slater determinants for a n truncated CC ansatz. Note that due to the exponential nature of the ansatz, the excitations degrees that were lost in the size-constistency check of the CID wave function are recovered in the CCD one, with additional terms involving single excitations of each fragment A and B. Those last terms however cancels when computing the amplitudes, and the CC method is size-consistent when constructed from a size-consistent reference.

The CC ansatz can also be treated by means of Rayleigh-Schrödinger perturbation theory. The wave function is expanded as in eq. [START_REF] Colonna | Correlation energy within exact-exchange adiabatic connection fluctuation-dissipation theory: Systematic development and simple approximations[END_REF] with

Φ (k) i = µ c (k) µ Ĉ(k) µ Φ (0) i (69)
where c µ design the coefficient associated to the excitation operator Ĉµ . To first order in the ground state we have

E (1) |Φ⟩ = ( F -E (0) )|Φ (1) 0 ⟩ + ( Vee -VKS ) |Φ⟩ (70)
so that by projecting onto a doubly excited determinant we find

⟨Φ rs ab |( F -E (0) |Φ (1) 0 ⟩ = -⟨Φ rs ab |( Vee -VKS )|Φ⟩ (71)
⟨Φ rs ab | Fc (1) 2 Ĉ(1) 2 |Φ⟩ -E (0) ⟨Φ rs ab |c (1) 2 Ĉ(1) 2 |Φ⟩ = -(v abrs -v absr ) * (72) c
(1)

abrs = (v abrs -v absr ) * ϵ a + ϵ a -ϵ r -ϵ s (73)
If we consider the Coupled Cluster Double (CCD) ansatz we recover by inserting (73) into (68) the MP2 correlation energy. The CCSD(T) method consist in solving for the CC amplitudes up to second order, then computing the perturbative triple amplitude to second order by using the CCSD amplitudes (rather than the perturbations to second order), and finally correcting the CCSD energy by the third order perturbative energy term.

Basis sets

For electronic molecular calculations, it is customary to use atomic orbitals basis sets. Pople basis sets were constructed to approximate Slater-type orbitals (STO) [START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF][START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]. Those are nodeless hydrogen-like orbitals that possess an exponential decrease with the distance from the nuclei and thus satisfy the nuclear cusp condition. 5 To facilitate the numerical evaluation of the electrons integrals in [START_REF] Bucko | Manybody dispersion corrections for periodic systems: an efficient reciprocal space implementation[END_REF][START_REF] Carleo | Machine learning and the physical sciences[END_REF][START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF], these orbitals can be approximated by linear combinations of Gaussian functions, which possess a e -αr k are made ofe one Gaussian function (3-21G). To have more flexibility in describing the effect on atoms of the non-uniform electric field due to the others atoms, polarizations functions can be added. Those are orbitals with higher angular momentum, p-like for Hydrogen, d-like for Lithium, etc. The 3-21G* add thus five d k -like orbitals made of one Gaussian function each, for a total of fourteen orbitals. The 3-21G** add again p-like orbitals to the core shell. In order to have more reliable post Hartree-Fock calculations, it is necessary to augment the number of virtual orbitals, while accounting for correlation [START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF][START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]. A first treatment of the sort it to diagonalize the one-particle reduced density matrix ρ(x, x ′ ) = ⟨ ψ † (x) ψ(x ′ )⟩. The obtained eigenvectors are called natural orbitals (NO), and the associated eigenvalues are occupation numbers. They permit a faster convergence of the post-Hartree Fock methods [START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF][START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]: one perform for example a truncated CI calculation on top of the Hartree-Fock canonical orbitals, after what the natural orbitals are computed. Only the ones with the highest occupation numbers are kepts, and the procedure is repeated. Similarily, Atomic Natural Orbitals (ANO) can be computed from post-Hartree Fock methods on atoms alone, and subsequent molecular calculations are performed while retaining only a subset of such ANO. 6The correlation-consistent basis sets where constructed by Dunning [START_REF] Thom | Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen[END_REF] in a similar way than that of ANO, but the orbitals retained in each stage of their hierarchy are those with similar contributions to the correlation energy. The correlationconsistent polarized valence basis sets are denoted by cc-pVXZ, where "cc" stands for correlation consistent, "p" for polarization, "V" for valence (only valence electron correlation has been considered in their design), and "XZ" corresponds to double zeta (DZ), triple zeta (TZ), etc. In Figure 6 is shown a dissociation curve of the closed-shell Hydrogen Bore molecule, using different methods presented supra. The FCI curve, is the exact one in the 6-31G basis set considered here. The RHF solution is qualitatively correct near equilibrium but obviously lack size-consistency in the dissociation limit as both Bore and Hydrogen are open-shell atoms. The CISD wave function, is a net improvement toward the FCI reference, while it overestimates the atomization energy. The MP2 energy is also a good description. The atomization energy overestimation is bigger in MP2, together with a too short equilibrium bond, as it is constructed on Hartree Fock orbitals only, while the orbitals are relaxed in CISD. The coupled cluster energies are the best approximation on this small system, with CCSD(T) giving the most precise atomization energy. Note that both MP2 and CCSD(T) are not qualitatively well behaved near the dissociation limit: neither is variational, the Møller-Plesset perturbative potential (also present in the triple (T) perturbation of CCSD(T) cannot be considered a small perturbation near the dissociation. Note that we only considered single-reference methods in this thesis. Multi-reference methods are extension of the former where additional determinants serve as reference upon which the correlation is restored. This is particularly efficient for static correlation and the correct description of bond breaking.

Methods comparison

density functional theory

The previous methods are commonly called wave-function methods. They rely on constructing the missing correlation energy from the Hartree-Fock single Slater determinant. While they can produce very accurate results, they still possess an unfavorable numerical scaling with the number of electrons. This renders them impracticable to use for moderate size systems (more than hundreds of electrons). In Density functional theory (DFT), the interest is shifted from the wave function to the density: the many electron problem is reformulated as a variational problem over the density of electrons. The later being a function of only three spatial coordinates, while the wave function is a N-fold tensor product of one-electron functions of three spatial coordinates (without considering spin), this is an enormous simplification. In practice however, the exact formulation of DFT involves an unknown universal functional of the density. A self-consistent method has been devised to obtain quantitative results from this theoretical framework, in which the many electron problem is mapped onto a non interacting one with the same density. This non interacting system must be solved exactly, leading to a cubic scaling with the basis set size, still more efficient than the quartic scaling of Hartree-Fock. 7 In this context, DFT has become the very basis of material simulations and systems with thousands of atoms can be studied routinely [START_REF] Maurer | Advances in Density-Functional Calculations for Materials Modeling[END_REF].

Hohenberg-Kohn theorems

Let us rewrite here the electronic Hamiltonian of eq. ( 4) by dropping the e subscript:

Ĥ = T + Vee + Vext ( 74
)
and let Vext = N i=1 v(r i ) be a local (i.e. multiplicative) potential. Since we are not interested in spin polarized systems or non-collinear spin systems, we will drop the spin variable σ i to alleviate the notation, but we shall keep in mind a missing sum over this variable. Moreover, the 4πϵ 0 factor in the Coulomb potential will be dropped from now on. For any N-electron wave function

Ψ(r 1 , r 2 , • • • , r N ) the electron density is given by 8 n(r) = N dr 2 • • • dr N |Ψ(r, r 2 , • • • , r N )| 2 (75)

First Hohenberg-Kohn theorem

The first of Hohenberg-Kohn theorems reads [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF]:

Theorem 1 There exist a one-to-one mapping between the external potential Vext , the non-degenerate ground state, and its corresponding density.

The proof proceeds by assuming that two external potentials V1 ext and V2 ext that differ by more than a constant shift lead to the same density n 1 (r) = n 2 (r) = n(r). Two different external potentials corresponds to two different Hamiltonians Ĥ1 and Ĥ2 , which themselves correspond to two distinct ground states |Ψ 1 ⟩ and |Ψ 2 ⟩. 9We have then

⟨Ψ 1 | Ĥ1 |Ψ 1 ⟩ < ⟨Ψ 2 | Ĥ1 |Ψ 2 ⟩ E 1 < ⟨Ψ 2 |H 2 + V1 ext -V2 ext |Ψ 2 ⟩ E 1 < E 2 + n(r)(v 1 (r) -v 2 (r))
and reciprocally

E 2 < E 1 + n(r)(v 2 (r) -v 1 (r))
so that

E 1 + E 2 < E 1 + E 2 .
Since this is impossible, we must have

n 1 (r) ̸ = n 2 (r).
This unique correspondence between the Hamiltonian (through the external potential) and the ground state density allows one to write either the ground state or the ground state energy as a functional of the ground state density, i.e.

|Ψ 0 ⟩ = |Ψ[n 0 ]⟩ (76)
E 0 = E[n 0 ] (77)
It is important to further note that this correspondence also says that in principle all properties of the system are determined by the ground state density [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF].

Second Hohenberg-Kohn theorem

Theorem 2 The ground state energy associated to a given external potential is the minimum value of the energy functional of the density at the point where the later coincides with the ground state density.

Let us first define the functional

F[n(r)] = ⟨Ψ[n]| T + Vee |Ψ[n]⟩. (78)
Note that this functional is universal, as it does only depend on the electron number, all the specifics of the system under consideration being in the external potential. The energy functional E[n(r)] is now defined as

E[n] = F[n] + drv(r)n(r) (79)
and we have that

⟨Ψ[n]| Ĥ[v]|Ψ[n]⟩ ⩾ E 0 ⟨Ψ[n]| T + Vee |Ψ[n]⟩ + ⟨Ψ[n]| Vext |Ψ[n]⟩ ⩾ E 0 ⟨Ψ[n]| T + Vee |Ψ[n]⟩ + drv(r)n(r) ⩾ E 0 F[n] + drv(r)n(r) ⩾ E 0 E[n] ⩾ E 0
which gives the second theorem according to which

min n E[n] = E[n 0 ] = E 0 (80)
Thus, in principle, to solve the Schrödinger equation for the electronic ground state, one simply needs to minimize the functional E[n] with respect to the density n(r). The functional F[n] is however not known analytically and this remains an open problem [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]. To face this problem, a practical scheme was devised by Kohn and Sham in 1965 [START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF], one year after Hohenberg and Kohn's paper.

Kohn-Sham Density functional theory

The Kohn-Sham scheme for DFT relies on the following assumption:

It is possible to associate to the exact ground state density a non-interacting auxiliary system with the same ground state density

Although this statement remains unproven for general purposes [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF] we will assume its validity.

The HK theorems also apply to non-interacting system. Let |Φ⟩ belong to the space of ground-states of such non-interacting systems in correspondence with external potentials.

Since there are no electrons-electrons interactions, those kets |Φ⟩ are single Slater determinants in some basis {ϕ i }, in accordance with Section 2.2:

|Φ⟩ = 1 √ N! Â |ϕ 1 ϕ 2 • • • ϕ N ⟩ . ( 81
)
Their density reads

n(r) = i |ϕ i (r)| 2 (82)
One can now redefine the universal functional F[n] as

F[n] = T s [n] + E Hxc [n] (83)
where T s [n] is the kinetic energy functional of the auxiliary system:

T s [n] = Φ[n] -h 2 ∇2 r 2m Φ[n] (84)
E Hxc is further decomposed as

E Hxc [n] = E H [n] + E xc [n] (85)
with

E H [n] = 1 2 drdr ′ n(r ′ )n(r) |r -r| ′ . ( 86
) E H [n]
is the so-called Hartree functional and corresponds to the mean-field electrostatic repulsion energy. In the Hartree-Fock formalism we had E HF = ⟨ T ⟩ HF + ⟨ Vee ⟩ HF where ⟨ Vee ⟩ HF decompose into the substraction of a direct and exchange term. E H [n] reproduce the direct term but includes a spurious self-energy term by missing the exchange term. E xc [n] is called the exchange correlation functional and reads explicitly

E xc [n] = T [n] -T s [n] + V ee [n] -E H [n]. (87)
It contains the difference between the kinetic energy of the interacting and noninteracting system, as well as the difference between the Coulomb interaction and the mean-field Hartree term. As such, all the intricacies of the interacting manybody problem are contained in this functional.

In order to minimize E[n] with the universal functional of eq. ( 83), let {ϕ i } be a set of orbitals subject to the constraint ⟨ϕ j |ϕ i ⟩ = δ ij . One can construct the following Lagrangian functional:

L[{ϕ i }] = E[{ϕ i }] - ij ϵ ij ⟨ϕ j |ϕ i ⟩ -δ ij . ( 88
)
By cancelling the functional derivatives ∂L/∂ϕ i , one obtains the Kohn-Sham equations:

- h 2 ∇ 2 2m + v KS (r) ϕ i = ϵ i ϕ i ( 89
)
with the Kohn-Sham potential v KS (r) being defined as

v KS (r) = v ext (r) + ∂E Hxc [n] ∂n(r) (90)
= v ext (r) + dr ′ n(r ′ ) |r -r ′ | + ∂E xc [n] ∂n(r) (91)
= v ext (r) + v H (r) + v xc (r) (92)
Thus according to the Kohn-Sham assumption, there exist a correspondence between the exact interacting many-electron problem and the Kohn-Sham auxiliary non interacting system with Hamiltonian ĤKS = -

h 2 ∇2 2m + VKS (93)
If the exchange-correlation functional was known, it would suffice to solve the later to obtain the exact ground-state density of the interacting problem. Unfortunately the exchange-correlation functional, so as the universal F[n] functional, remains unknown.

In the next section we will see some approximations to E xc [n], but before to go on, let mention a few practical points:

Self consistent field procedure

In the Kohn-Sham Hamiltonian, the potential depends upon the basis functions ϕ i , so that in practice, once an approximation has been chosen for E xc [n], the Kohn-Sham equations needs to be solved in a self-consistent fashion, as illustrated in Figure 7 2.

Basis sets

Until know we have not discussed periodic potentials that appears in solids. We will restrict ourselves to one dimension for simplicity. To describe such systems, the crystal is approximated as the periodic repetition of N unit cells of length a, with periodic boundary conditions. The one-electron wave function is thus such that ψ(x + Na) = ψ(x). The translation operator is T (na) = e na d dx = e -ina p/ h with eigenvalues e ikna . It commutes with the Hamiltonian so that it exist a simultaneous basis ϕ(x) with H e T (na)ϕ(x) = Eϕ(x + a) = Ee ikna ϕ(x).

No

Moreover the boundary condition imposes k = 2πm

Na where m is an integer. The function

u k (x) = e -ikx ϕ(x) is such that T (na)u k (x) = e -ik(x+a) ϕ(x + a) = u k (x).
The one-particle solutions of the Schrödinger equation in a periodic potential can thus be written as (Bloch theorem):

ϕ k (x) = e ikx u k (x) (94)
where u k (x) possesses the periodicity of the potential. k can be restricted to the values 2πm Na , -N 2 < m < N 2 , if we add another quantum number n to specify the congruence. The reduced wavenumber values defines the Brillouin zone, and the n number defines the band index. Since the function u nk (x) is periodic, it can be further decomposed into plane waves through Fourier decomposition. Close to the nuclei, a high number of plane-waves is necessary to reproduce the correct one-electron wave functions, and pseudo-potentials were devised to alleviate the associated computational burden, thus allowing accurate calculations with lower energy cutoff.

Pseudo-potentials

The idea of pseudo-potential (PP) is to use an effective potential for the nuclei and the core electrons, only to consider the valence electrons left, interacting with the effective core potential (see Figure 8). They are constructed by solving an approximate all-electron problem for the atoms alone, and come in different flavors: norm-conserving PP impose the equality of the charge integral inside the core radius R c of the all-electron valence wave functions and the corresponding wave functions obtained with the PP (pseudo-wave function) [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]. Ultrasoft PP release this constraint and report it into the Kohn-Sham eigenvalue problem [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]. Projector augmented waves (PAW), employ a similar strategy than the ultrasoft PP but reproduce the true oscillatory behavior of the valence all-electron wave function. They are implemented in the code VASP [START_REF] Kresse | Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set[END_REF] which was used for the results of chapters 4 and 5

Exchange correlation

While the true exchange functional is not known, we can obtain for it a formal expression that will be convenient to discuss beyond-DFT approaches. We will use the adiabatic connection concept [START_REF] Harris | Adiabatic-connection approach to Kohn-Sham theory[END_REF]. The idea is to define a Hamiltonian parametrized by λ, decomposed into a non-interacting (mean-field) one-body Hamiltonian Ĥ0 and a many-body Hamiltonian Ĥ1 (λ). By slowly (adiabatically) varying λ, one can interpolate between the mean-field and the fully interacting problem. It is moreover assumed that the λ-parametrized external potential is such that the density is kept constant during the adiabatic switch [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF].

Ĥ(λ) = Ĥ0 + λ Ĥ1 (λ) = T + Vext λ + λ Vee , (95)
Ĥ1 = Vee + 1 λ Vext λ -Vext λ=0 (96)
Vext λ=0 = V ext + VMF (98) Vext λ=1 = V ext (99) (97)
From the Hellmann-Feynman theorem we have that

dE λ dλ = ⟨ψ λ | d Ĥλ dλ |ψ λ ⟩ (101)
so that

E(λ = 1) = ⟨ψ 0 |H 0 |ψ 0 ⟩ + 1 0 dλ⟨ψ λ | Ĥ1 (λ) + λ d Ĥ1 (λ) dλ |ψ λ ⟩ (102)
Hence the true energy E(λ = 1) is recovered from the mean-field ground state |ψ 0 ⟩. This is closely related to the Gell-Mann and Low theorem mentioned in subsection 2.2.4. As in the later, level-crossing (from symmetry breaking for example) must be absent on the adiabatic path for this scheme to work (in particular, one must have ⟨ψ 0 |ψ 1 ⟩ ̸ = 0 [START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF]). Going on, we have that

E = E 0 + 1 0 dλ⟨ψ λ | Vee |ψ λ ⟩ + ⟨ψ λ | d Vext λ dλ |ψ λ ⟩ (103)
= E 0 + 1 0 dλ⟨ψ λ | Vee |ψ λ ⟩ + drn(r)(v ext 1 -v ext 0 ) (104)
where we used the fact that the density is kept constant when going from the first to the second line. From second quantization we can rewrite ⟨ Vee ⟩ as [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF] ⟨

Vee ⟩ = 1 2 drdr ′ 1 |r -r ′ | ⟨ n(r) n(r ′ )⟩ -⟨ n(r)⟩δ(r -r ′ ) = 1 2 drdr ′ 1 |r -r ′ | ⟨δ n(r)δ n(r ′ )⟩ + ⟨ n(r)⟩⟨ n(r ′ )⟩ -⟨ n(r)⟩δ(r -r ′ ) (105)
with the density fluctuation

δ n(r) = n(r) -⟨ n(r)⟩. ( 106
)
From Kohn-Sham DFT we have

E[n] = T s [n] + drn(r)v ext (r) + E H [n] + E xc [n] (107)
As the density is kept constant during the integration, we can identify the middle term in eq. ( 105) with E H [n], so that

E xc [n] = 1 2 1 0 dλ drdr ′ ⟨δ n(r)δ n(r ′ )⟩ λ -n(r)δ(r -r ′ ) |r -r ′ | (108)
Let us now separate the exchange correlation energy in an exchange and a correlation contribution :

E xc = E x + E c (109)
The total exchange energy is defined as [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF][START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]]

E x = ⟨ψ 0 | Vee |ψ 0 ⟩ -E H [n] = ⟨V ee ⟩ 0 -E H [n] (110)
where |ψ 0 ⟩ is a the Kohn-Sham Slater determinant. This definition allows to write the correlation energy as

E c = E -⟨ψ 0 | Ĥ|ψ 0 ⟩ (111)
in complete analogy with the definition of eq. ( 33). Using eq. ( 105), we can rewrite eq. ( 110) as

E x = 1 2 drdr ′ 1 |r -r ′ | ⟨δ n(r)δ n(r ′ )⟩ 0 -n(r)δ(r -r ′ ) (112)
and finally

E c = 1 2 1 0 dλ drdr ′ 1 |r -r ′ | ⟨δ n(r)δ n(r ′ )⟩ λ -⟨δ n(r)δ n(r ′ )⟩ 0 (113)
The correlation energy is thus expressed as the integral of the difference between the density-density correlation function of the interacting and non-interacting system times the Coulomb kernel while the interaction is gradually turned on.

Approximations to the exchange correlation energy

One way to approximate the exchange correlation is to replace the true system by an idealization in eq. ( 108). The homogeneous electron gas (HEG or "jellium") is one such approximate system. It consists of N interacting electrons, with a positive uniform background charge to model the nuclei. The HEG energy is function of the density, or equivalently of the average dimensionless distance r s between electrons (or Wigner-Seitz radius). The later is defined by the relation

4 3 π(r s a 0 ) 3 = 1 n (114)
with a 0 = h 2 m e e 2 the Bohr radius. While being an approximate model of reality, the HEG is not analytically tractable, but numerical results of its correlation energy e HEG c (r s ) have been computed and are tabulated [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF][START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]. As an example, in the high-density limit, one can find the following energy density [START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF] e HEG c (r s ≪ 1) ∼ e 2 n a 0 (0.0622 ln r s -0.096).

The exchange energy density is however known and reads [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF] e

HEG x (n) = - 3(3π 2 ) 1/3 e 2 4π dr n(r) 4/3 . ( 116
)
The local density approximation (LDA) consists in using those values of the exchange and correlation energy of the HEG in

E xc [n(r)]. E LDA xc [n] = dr e HEG x (n) + e HEG c (n) n(r) (117)
That is, for each point r, the XC energy is given by the HEG XC energy for the corresponding density n(r), hence the "local" in the name. This scheme is expected to perform well for approximately free electrons systems such as conduction electrons in metals and poorly for localized electronic system such as atoms [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF].

The next class of approximations are coined generalized gradient approximations (GGA), and consist in adding a dependence on the density gradient through some function f:

E GGA xc [n] = dr n(r)e HEG xc (n)f(n, ∇n) (118)
Different such functionals have been proposed, one of the most widely used [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF][START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF] being the GGA functional of Perdew, Burke and Ernzerhof (PBE) [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF], which will be used in this thesis. This gradient dependence allows a better treatment of inhomogeneous systems and improve upon the local aspect of the LDA (GGA functional are usually referred as semi-local functionals [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF]. There exist also Meta-GGA functionals which add a dependence on the density laplacian ∇ 2 n and on the non-interacting kineticenergy density τ(r)

∝ i |∇ϕ i (r)| 2 .
Finally let us mention the hybrid functionals that combine one functional (either LDA/GGA/Meta-GGA) with the so-called exact Hartree-Fock exchange energy 10

E HF x = - 1 2 i,j σ drdr ′ ϕ * iσ (r)ϕ * jσ (r ′ )ϕ jσ (r)ϕ iσ (r ′ ) |r -r ′ | . ( 119
)
The idea behind this is that 108) so that by turning on the interaction one interpolates between E HF -E H = E HF

E xc = 1 0 dλ⟨ Vee ⟩ λ -E H from eq. (
x and the fully interacting E xc [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF]. The PBE0 functional for example is given by

E PBE0 xc = 1 4 E HF x + 3 4 E PBE x + E PBE c ( 120
)
10 Or sometime simply named exact exchange. It is exact in the sense that it is the same expression as for the Hartree-Fock exchange but the two do not coincide as the orbitals differ. Note that spin has been explicitly re-introduced in this expression.

Note that GGA, Meta-GGA or hybrid functionals all incorporate the LDA approximation in their core.

Random Phase Approximation

An alternative route to the LDA and subsequent approximations is to perform an approximation on the exact form of the exchange-correlation of eq. ( 113). Let us define the time-ordered density correlation function

χ(r, t, r ′ , t ′ ) = -i⟨Ψ|T δ n(r, t)δ n(r ′ , t ′ )|Ψ⟩ (121)
which is related to the retarded density correlation function

χ R (r, t, r ′ , t ′ ) = -iΘ(t -t ′ ))⟨Ψ|[δ n(r, t), δ n(r ′ , t ′ )]|Ψ⟩ (122)
defined in linear response theory as the kernel of the density response to an external perturbation(⟨δ n⟩ = χ R V) [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF][START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF]. The time-ordered density correlation is subject to the following Dyson equation written (in Fourier space, for the homogeneous problem)

χ(q, ω) = Π(q, ω) + Π(q, ω)V ee (q)χ(q, ω) (123)
with Π(q, ω) the irreducible polarization, sum of all processes beginning by the creation of a particle-hole pair and ending by its annihilation, such that their corresponding Feynman diagrams are not disconnected when cutting one interaction line. The polarization equally defines the effective, dressed or screened interaction

V eff by V eff (q) = V ee (q) + V ee (q)Π(q, ω)V eff (q). (124)
In this context, the random phase approximation (RPA) consist in keeping only the zeroth order term in the interaction for the polarization, which correspond to the time-ordered correlation function of the unperturbed system:

Π RPA ≡ χ 0 ≡ (125)
And the unperturbed response function χ 0 is given as function of the unperturbed one-electron orbitals. Eqs. ( 123) and [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF] in the random phase approximation are diagrammatically depicted as [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF][START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF]]

χ RPA = + χ RPA (126) = + (127)
After inserting χ RPA and χ 0 in eq. ( 113), one finds after performing the integration over λ, taking the time Fourier transform of the response function, and restricting the frequency integral to positive value using the response function symmetry [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF][START_REF] Ren | Randomphase approximation and its applications in computational chemistry and materials science[END_REF]:

E RPA c = 1 2π ∞ 0 dω drdr ′ [ln 1 -χ 0 (r, r ′ , ω) 1 |r ′ -r| + χ 0 (r, r ′ , ω) 1 |r ′ -r| ] (128)
The exchange energy E x in eq. ( 112) is similarly written using χ 0 . By explicitly expanding χ RPA in eq. ( 121) or the logarithm in eq. ( 128), we obtain the correlation energy as a sum involving terms ∼ (V ee χ 0 ) n , n ⩾ 2. Those correspond to the so-called ring diagrams, the first of them (n = 2) being the MP2 direct diagram, and the next one was given in eq. ( 53), evaluated with the Kohn-Sham orbitals and eigenvalues [START_REF] Engel | Density Functional Theory: An Advanced Course[END_REF][START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF][START_REF] Ren | Randomphase approximation and its applications in computational chemistry and materials science[END_REF].

In the high density regime of the HEG the RPA correlation energy yields the exact correlation energy [START_REF] Fetter | Quantum Theory of Many-particle Systems[END_REF][START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF] (eq. ( 115)). For inhomogeneous systems it is still a valuable approximation, that systematically accounts for correlation at all orders in perturbation, and stays size-consistent. The RPA is particularly suited to study systems where long range weak interaction, such as van der Waals, are important [START_REF] Ren | Randomphase approximation and its applications in computational chemistry and materials science[END_REF]. By taking the exact exchange and computing E RPA c as in eq. ( 128) we are considering the so-called direct RPA method. Beyond RPA methods can include higher order exchange terms [START_REF] Ren | Randomphase approximation and its applications in computational chemistry and materials science[END_REF]. The numerical cost of direct RPA algorithms typically range between O(N 4 ) and O(N 6 ) [START_REF] Ren | Randomphase approximation and its applications in computational chemistry and materials science[END_REF]. This method is also closely related to the so-called GW approximation to the self-energy, the later being approximated as the product of the interacting one-particle green function and the RPA screened interaction [START_REF] Golze | The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy[END_REF][START_REF] Ren | Randomphase approximation and its applications in computational chemistry and materials science[END_REF], yielding accurate quasiparticle energies 11 for solids and thus furnish reliable band gaps values.

Methods Comparison continued

In Figure 9 is shown again the dissociation curve of the Hydrogen Bore molecule in 6-31G basis with various DFT approximations. The RPA curve was computed using PBE Kohn-Sham orbitals as input. Since DFT methods are not constructed on top of a Hartree-Fock reference, all curves has been vertically shifted so that they share the same minimum energy (as one can do: only energy differences are physically meaningful). The LDA present a better behavior than the Hartree-Fock curve but overestimate considerably the atomization energy. PBE is an improvement over LDA, but the atomization energy is again noticeably overestimated. The RPA behavior is much closer to the CCSD(T) one but the RPA equilibrium bond lenght is overestimated. Calculations were realized with Nwchem [START_REF] Aprà | NWChem: Past, present, and future[END_REF] in the 6-31G vasis set. The RPA correlation energy was computed using PBE orbitals

finite temperature properties

We have seen various approximations to solve the electronic time independent Schrödinger equation. We would like now to be able to obtain finite temperatures properties: that is, we are interested in computing averages of some observable O in a given thermodynamical ensemble to which a density operator ρ is associated:

⟨ Ô⟩ = tr ρ Ô tr ρ (129)
In the canonical ensemble the density matrix is given by ρ = e -β Ĥ [START_REF] Mattuck | A Guide to Feynman Diagrams in the Many-body Problem[END_REF] with the inverse temperature β = (k B T ) -1 , and the trace runs over all accessible configurations, that is, all positions of the electrons and nuclei. By assuming the equality of the time and ensemble averages (ergodic hypothesis [START_REF] Landau | Physique théorique[END_REF]), we can obtain the former by studying the dynamics of the system of interest, whether it is a molecule or a crystal. In order to do so, we need to come back to the full Hamiltonian (2) Ĥ = Ĥe + Tn + Vnn . We will describe an approximation in which the electronic state is assumed to remain the lowest one (adiabatic approximation), while the nuclear and electronic degrees of freedom are separated (Born-Oppenheimer approximation). The nucleus dynamics is further solved by solving classical equations of motion using electronic forces.

Ψ(r, R) = ∞ i=0 ψ i (r; R)χ i (R). (131)
Thus, the {χ i (R)} are wave functions of the nuclei positions. Applying Ĥ and taking the product by ψ j yields :

ϵ j (R) + Vnn + Tn -E χ j (R) = ∞ i=0 Λ ji χ i (R) (132)
with

Λ ji = α -h 2 2M α ⟨ψ j |∇ 2 R α |ψ i ⟩ + 2⟨ψ j |∇ R α |ψ i ⟩∇ R α (133)
The Λ term contains all couplings between the nuclear and electronic degrees of freedom. By setting all the off diagonal elements of Λ to zero, we put ourselves in the adiabatic approximation. Indeed, the electronic Hamiltonian contains all its time dependence inside the nuclei position. Expressing the time derivative as ∂ t = α ∇ R α Ṙα , and taking the dot product by ψ j we get for our general electronic wave function ψ(r; R)

= i a i ψ i (r; R) that ȧj (t) = - i h ϵ i (t)a j (t) - i α ⟨ψ j |∇ R α |ψ i ⟩ Ṙα a i (t) (134)
By neglecting the coupling terms between different instantaneous electronic eigenstates, if the system is at some time t 0 on a given PES ϵ k (R), it will remain in the corresponding instantaneous eigenstate ψ k (r; R) for all future times t > t 0 . The adiabatic theorem states that this is accurate as long as the PES are well separated and the electronic Hamiltonian is not varying too quickly [START_REF] Sakurai | Modern Quantum Mechanics[END_REF] :

⟨ψ j , t| Ḣe |ψ i , t⟩ ϵ j (t) -ϵ i (t) ≪ ϵ j (t) h (135)

The Born-Oppenheimer approximation

Moreover, by also equating to zero the diagonal terms in eq. ( 133), we are now performing the Born-Oppenheimer (BO) approximation. In that case we have to solve b

ϵ j (R) + Vnn + Tn χ j (R) = Eχ j (R) (136)
in which the PES ϵ j (R) is a local external potential felt by the nuclei (hence the name potential energy surface. The nuclei and electronic degree of freedom are now totally separated, resulting in a huge simplification of the problem to solve. This is a good approximation near equilibrium points of the nuclei (where ∇ R is close to zero). This however breaks down in the vicinity of avoided crossing points, where the instantaneous electronic eigenstates abruptly change from one diabatic state to the other, or near conical intersections, where the denominator of eq. ( 135) diverge due to the degeneracy. It should moreover not be trusted in metals, for which groundstate and first excited states are not well separated [START_REF] Pisana | Breakdown of the adiabatic Born-Oppenheimer approximation in graphene[END_REF].

Classical motion of the nuclei

A last approximation which is made in what is generally called ab initio Molecular Dynamics (AIMD) is that the nuclei are moved by solving classical equations of motion. That is, once the electronic ground state density has been found for a given set of nuclei positions, one compute the forces acting on the nuclei by using the Hellmann-Feynman theorem while varying the nuclei positions, and update the later. This classical approximation neglects the so-called nuclear quantum effects. Those can be important for light elements, or in order to account for nuclei delocalization. The later for example is responsible for a bandgap reduction (zero point renormalization) of ∼ 0.4 eV in the diamond at zero temperature [START_REF] Karsai | Electronphonon coupling in semiconductors within the GW approximation[END_REF]. Finally, in order to account for the temperature, one must choose a statistical ensemble within which to work and fix the temperature. For example to maintain temperature in the canonical (N,V,T) ensemble, one can choose the Andersen thermostat [START_REF] Hans | Molecular dynamics simulations at constant pressure and/or temperature[END_REF] procedure, in which new nuclei velocities will be drawn from a Maxwell-Boltzmann statistics at random times. To summarize, what we will call AIMD in the next parts of this thesis consist in the Born-Oppenheimer separation of nuclear and electronic degrees of freedom, the later being solved using DFT, before classical equations of motions are applied to the nuclei using forces derived from the DFT.

Monte Carlo sampling

Instead of computing time averages by performing Molecular Dynamics, one can directly compute ensemble averages by sampling configurations in order to reproduce the target probability density given by eq. ( 130). The Metropolis algorithm is a popular and simple method to achieve this: from any starting configuration R with energy ϵ, a new configuration R ′ with energy ϵ ′ is proposed with proposal probability density g(R ′ |R) and is accepted with probability min(1, e -β(ϵ ′ -ϵ) ). The sequence of configurations sampled in this way is ensured to converge to the canonical Boltzmann probability [START_REF] David | Information Theory, Inference, and Learning Algorithms[END_REF][START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF]. The choice of the proposal distribution is arbitrary and should be tuned in order to obtain a reasonable acceptance ratio. For example a Maxwell-Boltzmann distribution can be used to obtain random displacements of the atoms going from R to R ′ , as done in Section 4.3.2. Since it only demands energy calculations, this methodology can be a cheaper alternative to methods where the forces are computationally demanding.

Thermodynamic perturbation theory

Consider two electronic structure methods X and Y with corresponding Hamiltonians H X and H Y = H X -∆V. Let ⟨O⟩ X and ⟨O⟩ Y be canonical ensemble averages of the observable O evaluated with X and Y Hamiltonians. Thermodynamic perturbation theory relates both averages [START_REF] Chipot | Free energy calculations: theory and applications in chemistry and biology[END_REF][START_REF] Pohorille | Good Practices in Free-Energy Calculations[END_REF]:

⟨O⟩ Y = Oe -βH Y dΓ Y e -βH Y dΓ Y = Oe β∆V e -βH X dΓ X e β∆V e -βH X dΓ X (137) so that ⟨O⟩ Y = ⟨Oe β∆V ⟩ X ⟨e β∆V ⟩ X . ( 138
)
Thus the average in the Y ensemble can be expressed as a reweighting of the average in the X ensemble. The corresponding weigths are the Boltzmann factors of the energy difference associated with the configurations upon which the average is performed.

If one has performed an AIMD simulation using some DFT functional (X method) to obtain some average property O, it is then possible to obtain the corresponding ensemble property within another method Y, by computing the energies using said Y method for each configuration of the AIMD. Note that we used the fact that the configurational spaces are identical when going from one integral to the other. In practice the integrals are replaced by summation over finite sets of configurations around equilibrium configurations. In that case, for formula [START_REF] Péter | Optimization of the Linear-Scaling Local Natural Orbital CCSD(T) Method: Improved Algorithm and Benchmark Applications[END_REF] to be a valid, the probability densities of both methods must sufficiently overlap on the space sampled using the X method (see Figure 10).

introduction

In machine learning (ML), one starts with a set of observations -the training dataset -and a model. The model is a function or a set of functions of the training points and of some parameters. One then use an algorithm (machine) to adjust the parameters of the model (learning) in order to optimize some objective. A lot of different tasks can be tackled by ML, such as regression, classification, clustering, sampling, pattern recognition, etc. ML is generally divided into three classes: supervised learning, unsupervised learning, and reinforcement learning. In supervised learning, the training dataset is composed of the input data points and the targets. The inputs characterize the system studied (for example the geometry and chemical species of a molecule) while the target is the observable one would like to predict, given new unseen input. It can be a continuous variable (regression) such as the atomization energy, or a discrete variable (classification) such as the stability of the molecule. In unsupervised learning, there are only input, without target. The objective is to uncover patterns in the training dataset in order to generalize. For example clustering consists in dividing the data into different clusters (the unsupervised pendant of classification), dimensionality reduction aims at extracting the most important features from a set of datapoints, while generative models learn the probability distribution of the dataset, from which one can sample new input points. Finally in reinforcement learning, the ML model learn to interact with an environment in order to maximize its reward.

In the next sections we will first focus on the regression task of fitting a set of atomic configurations to their total energy. To this purpose we will discuss physical descriptors, which are optimized representations of molecular and periodic physical systems in order to produce useful inputs. Then we will describe how to go from linear regression to kernel ridge regression. Dimensionality reduction will be discussed afterward, with the principal component analysis and t-stochastic neighbor embedding algorithms. Finally to prepare the last chapter of this thesis we will introduce generative models and a description of the restricted Boltzmann machine model.

physical descriptors

Suppose we would like to find a ML model that learn the relationship between a set of atomic configurations (input) and their total energy (target), thus avoiding explicitly solving the Schrödinger equation. The later being invariant under global translations and rotations of the chemical configuration, it is desirable that the input representation, or descriptor, also possesses these symmetries, to circumvent the need for the model to learn the later. This disqualifies the traditional Cartesian coordinates or Z-matrix. Additionally, a good representation should also be invariant to permutations of chemical elements. The last condition is less stringent as one can impose a unique way of labelling the different atoms or chemical species. Different solutions have been proposed [START_REF] Behler | Atom-centered symmetry functions for constructing highdimensional neural network potentials[END_REF][START_REF] Hansen | Machine Learning Predictions of Molecular Properties: Accurate Many-Body Potentials and Nonlocality in Chemical Space[END_REF][START_REF] Huo | Unified Representation of Molecules and Crystals for Machine Learning[END_REF]. For example the Coulomb matrix is defined by [START_REF] Rupp | Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning[END_REF]

M ij =          Z 2.4 i 2 if i = j Z i Z j R ij if i ̸ = j (139)
where Z i are the nuclear charges and R ij the distances between atoms. A more involved description consist of considering local environments: a cutoff radius r c is defined, and for each atom, a description of each atoms inside the sphere of radius r c centered at the given atom will be given. This gives much more flexibility to the representation and allow the final model to better generalize to unseen configurations. Considering the size extensivity of the energy, this amounts to decompose the target property into atomic contributions. The Smooth Overlap of Atomic Positions [START_REF] De | Comparing molecules and solids across structural and alchemical space[END_REF] is one of such descriptors, that will be used in the next two chapters. Let us consider a chemical structure χ made of N χ atoms. First, gaussian atomic densities

ρ(r, χ a ) = i∈χ a exp (x i -r)/2σ 2 (140)
are constructed, centered at the atom considered, where the sum runs over atoms within r c around the atom a. Those densities are then decomposed into spherical harmonics, before to take their rotationally invariant power spectrum p(χ a ). The later vector is the SOAP descriptor associated to the atomic environment χ a . Once a descriptor has been chosen, the ML model should be able to assess the degree of similarity between different configuration in order to predict the energy of an unseen configuration. This can be achieved in kernel methods through the use of kernel functions. The later take as input two vectors and return a number, which is in correspondence to some inner product in another space, and will be defined in the next section. In case of local descriptors such as SOAP, the global similarity between two structures is decomposed in local similarities between atomic environments. The local SOAP kernel k(χ a , χ ′ b ) which measures the similarity between atomic region a of structure χ and atomic region b of structure χ ′ is defined as

k(χ a , χ ′ b ) = p(χ a ) • p(χ ′ b ) (141)
This is in fact equivalent [START_REF] De | Comparing molecules and solids across structural and alchemical space[END_REF] to the overlap integration of the two atomic environment densities upon rotation:

k(χ a , χ ′ b ) = d R drρ( Rr, χ a )ρ(r, χ b ) 2 ( 142
)
where R represent the rotation operator. See also Figure 11.

Finally, in order to globally compare two structures, different approaches exist [START_REF] De | Comparing molecules and solids across structural and alchemical space[END_REF], one of such is to define the following global kernel:

K(χ, χ ′ ) = N χ a=1 N χ ′ b=1 k(χ a , χ ′ b ) k(χ a , χ a )k(χ ′ b , χ ′ b ) . ( 143
)
This kernel defines the similarity between the structures χ and χ ′ by summing the local similarities between each atomic environment of the two. This is divided by a normalization factor that ensure that the local similarity is at most one for identical atomic environments.

kernel ridge regression

In regression, the training dataset is composed of a set of inputs x i ∈ R d and targets y i ∈ R. The objective is to optimize a model to reproduce the target given the input, such that the model is able to generalize to unseen input x and predict the target ỹ for those.

In linear regression (LR), we wish to find the hyperplane defined by β • x = ywhere we assumed that our plane goes through zero1 -from which the sum of the distances of the training points is minimized. This corresponds to the model

f( x) = β • x, (144)
We find the parameters β i of the model by minimizing the squared error on the training set:

L = n i=1 (f(x i ) -y i ) 2 = ||Xβ -y|| 2 (145)
2 where X is a matrix whose i th row is the training input x i . By cancelling ∂L/∂β, we find 2 the notorious normal equations

Structures Gaussian Smoothering

β = (X T X) -1 X T y (146)
We can always write β i as a linear function of the x i , such that β = Xα, and our linear model is rewritten as

f( x) = β • x = d i=1 β i xi = d i=1   n j=1 α j x j   i xi = n j=1 α j x j • x Solving for α gives α = (X T X) -1 y. ( 147
)
2 using ∂u T v ∂x = u T ∂v ∂x + v T ∂u ∂x
That is, linear regression only depends of inner products (hence it only depends on distances and angles between input points), and its solution can be expressed through the Gram matrix X T X, or equivalently the covariance matrix, which measure the similarity between inputs. Say now that we are interested in a phenomenon where the relation between input points is non-linear, and we want to perform a polynomial fit instead of a simple linear fit, with some polynomial feature map ϕ(x) ∈ R D , D ⩾ d. We can still use our linear regression model to fit an hyperplane defined in the feature space

R D3 f( x) = β • ϕ( x) = n j=1 α j ϕ(x j ) • ϕ( x) = n j=1 α j k(x j , x). (148)
In the last line we have defined the kernel function k. A kernel function should be positive definite, symmetric, and continuous [START_REF] Rupp | Machine learning for quantum mechanics in a nutshell[END_REF].

The simplest kernel is the linear kernel k(x i , x j ) = x i • x j , which corresponds to linear regression.

Let us now consider the feature map ϕ(x) = e

-x 2 2σ 2 (1, x σ , 1 2! x 2 σ 2 , 1 3! x 3 σ 3 , . . . ). The inner product ϕ(x) • ϕ(y) is equal to the gaussian kernel k(x, y) = e - (x-y) 2 2σ 2 .
We can perform a linear regression in the infinite dimensional space defined by ϕ (or equivalently, an infinite dimension polynomial fit) by computing directly the inner products through k. This in essence in the kernel trick [START_REF] Edward Rasmussen | Gaussian processes for machine learning[END_REF]. The coefficients of α, solution of the kernel regression, are now given by

α = K -1 y ( 149
)
where K ij = k(x i , x j ). The parameter σ in the gaussian kernel is called an hyperparameter and should be optimized separately. Other popular kernels are the laplacian kernel or the polynomial kernel [START_REF] Rupp | Machine learning for quantum mechanics in a nutshell[END_REF].

If the fit is "too good" on the training points, the model can fail to generalize to unseen data, a phenomenon know as overfitting. This can be checked by splitting the data into a training set and a test set. The model is trained on the training set while the error estimate is computed on the test set. The test error should not be much higher than the training error. To prevent overfitting, one can add a λ||f|| 2 = λα T Kα penalty term to the squared error, λ being a hyperparameter to optimize. This is known as Kernel Ridge Regression (KRR), with solution

α = (K + λI) -1 y (150)
KRR has been used to learn the atomization energy of molecules as well as other properties (HOMO/LUMO levels, polarizability, etc.) inside a dataset. In those works, chemical accuracy was achievable with training sets whose size were of the order of 10

dimensionality reduction

Dimensionality reduction is a method to reduce the number of variables in a dataset (thus lowering the dimension d of the input space) while retaining as much of the original information as possible. This can be helpful for a variety of reasons, such as improving computation efficiency, reducing noise, clustering and simplifying the visualization of data, the two later being used in the next two chapters.

Principal component analysis (PCA) is one of the most commonly used techniques in machine learning. PCA works by identifying the directions of maximum variance in the data and projecting the data onto a lower-dimensional subspace that captures the most important features. In practice, PCA is a low-rank approximation of our input matrix X through singular value decomposition (SVD): any n × d real matrix X can be decomposed as

X = UΣV T ( 151
)
where Σ is a n × d diagonal matrix whose diagonal entries are named singular values, and U, V are n × n and d × d orthogonal matrices. Geometrically we can see that most of the information on the application X is contained in those singular values that have the biggest absolute value. Similarly, SVD amounts to diagonalize the covariance matrix C = XX T or the Gram matrix C = X T X. 4 We wish to approximate the n × d matrix X with a k × d matrix X k , k being lower than d5 .

One possibility is to minimize the Frobenius norm of the distance between X and

X k : ||X -X k || F = tr[(X -X k ) T (X -X k )]
. Now the Young-Eckart Theorem states [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] that this is achieved by

X k = k i=1 σ i u i v T i
, where σ i , u i and v i are the i th singular value, in decreasing order of magnitude, and the i th columns of U and V.

Another route for visualizing high-dimensional data consist in mapping the data onto a low-dimensional space while preserving the pairwise distances between the data points. For example one may want to find a linear projection:

x ′ = 1 √ d Ax where A is a k × d matrix such that ||x ′ i -x ′ j || 2 ≈ ||x ′ i -x ′ j || 2 .
This can be easily achieved with A ∼ N(0, 1) according to Johnson-Lindenstrauss lemma [START_REF] Dasgupta | An elementary proof of a theorem of Johnson and Lindenstrauss[END_REF]. Another popular technique is t-distributed stochastic neighbor embedding (t-SNE). t-SNE uses a gaussian kernel instead of the cartesian norm to measure the distance between two points (hence is able to capture nonlinear relationships). It defines a similarity measure by constructing a gaussian conditional probability of the distance bewteen two points, both in the original d-dimensional space (p ij ), and in the reduced k-dimensional space (q ij ). The mapping is found by minimizing the sum of the Kullback-Leibler divergence between both probabilities:

C = i j p ij log p ij q ij = i KL(P i ||Q i ) (152)
The Kullback-Leibler divergence KL(P||Q) = i P(i) log P i Q i , or relative entropy, between two discrete probability distributions P and Q is commonly used as a measure of the similarity between both [START_REF] David | Information Theory, Inference, and Learning Algorithms[END_REF]. Note that the distance itself between points which appears in the Gaussian probabilities in not necessarily Cartesian, and one can use the SOAP kernel to define a distance between two chemical structures χ a and χ b :

D(χ a , χ b ) = (1 -K(χ a , χ b )) 1/2 . ( 153
)

generative models: restricted boltzmann machine

Generative models are a class of unsupervised ML models that can learn to generate new data that is similar to the training data. Among those, energy based models define an energy, function of the model parameters and input points, and a probability distribution over the inputs and parameters. In correspondence with statistical physics, the most probable inputs across the training dataset are then to be associated with the lowest energy and reciprocally during the learning stage. The Boltzmann Machine model is a complete weighted graph with the inputs v i living on the vertices, and the parameters J ij are the link weights. The weight matrix is symmetric, and there is no self interaction (J ii = 0 ∀ i). An energy function E(v, J) = -v T Jv is defined and associated to a Boltzmann probability P(v|J) = exp(-E(v, J))/ {v} exp(-E(v, J)). For binary input values, this is obviously a classical spin system with arbitrary long-range interactions. One can add a bias vector into the energy (an external field with contribution -a T v). This model however, cannot account for correlations higher than of order two as seen from its energy function. A solution is to introduce a hidden layer of vertices h i , that will model higher order statistics between visible input through their interaction with the hidden layer [START_REF] David | Information Theory, Inference, and Learning Algorithms[END_REF]. For the simple Boltzmann machine, the idea is to perform a Hubbard-Stratonovich transformation of the model to mediate the input interactions through a new degree of freedom h [START_REF] Mehta | A high-bias, lowvariance introduction to Machine Learning for physicists[END_REF]. 6 The restricted Boltzmann machine (RBM), pictorically represented in Figure 4, is a neural network model, consisting of one input layer of D visible binary units v i ∈ {0, 1}, one layer of P binary hidden units h j ∈ {0, 1}, and D × P weights W ij between both layers. Two bias vectors of components a i and b j are added to the visible and hidden layer, respectively. By introducing an energy function of Λ = {a, b, W} for a given configuration {v, h} as

E(v, h, Λ) = -(a T v + b T h + v T Wh), (154)
and an inverse temperature β = 1/T , one can define the probability distribution associated with the RBM over the input configurations v to be

P(v) = {h} e -βE(v,h) {h,v} e -βE(v,h) = {h} e -βE(v,h) Z , ( 155
)
where the partition function Z was defined. After training on a set of vectors v distributed according to some (usually unknown) probability distribution, the RBM can be used to generate new vectors v according to the P(v) distribution learnt by the model. To this purpose the Gibbs sampling can be used [START_REF] David | Information Theory, Inference, and Learning Algorithms[END_REF]: Starting from an initial random trial vector as input, one can obtain a hidden layer vector h using the conditional probability p(h i |v i ) and then obtain a new input-layer vector v using the conditional probability p(v i |h i ); the repetition of this operation a certain number of times forms a Markov chain that generates a vector v according to the probability P(v).

The Gibbs sampling algorithm is important for both the training and generative step of restricted Boltzmann machines (RBMs) [START_REF] David | Information Theory, Inference, and Learning Algorithms[END_REF]. It is based on the following steps:

6 The Hubbard-Stratonovich transformation use the simple Gaussian integral identity

K 2π 1 2 dh exp -1 2 Kh 2 + Ks 2 h = exp 1 2 Ks 2 . If J ij is positive semi-definite then ∃W | J = W T W so that E(v) = -v T W T Wv and by taking K = 1 and s = (Wv) T one obtains p(v) = Z -1 exp(-E(v)) = Z -1 exp 1 2 (Wv) T (Wv) = Z -1 dh exp -1 2 h T h + (Wv) T h = Z -1 dh exp(-E(v, h))
, see [START_REF] Mehta | A high-bias, lowvariance introduction to Machine Learning for physicists[END_REF].

1. The components of the initial input vector v (0) are drawn from a uniform probability distribution U(0, 1). Namely, each v (0) i with i = 1, . . . , D contains a random number in the interval [0, 1[.

The values h

(0) j with j = 1, . . . , P are initialized to 0. For each h (0) j a number d j is drawn from U(0, 1) and the value of h (0) j is then updated to 1 if d j < p(h j = 1|v (0) ).

3.

In order to update the vector v, the new components v

(1) i are first set to 0. For each v

(1) i a number d i is drawn from U(0, 1) and the value of v

(1) i is updated to 1 if d i < p(v i = 1|h (0) ).
4. The steps 2 and 3 are repeated L times.

Restricted Boltzmann machines are usually trained using the contrastive divergence learning procedure [START_REF] Hinton | A practical guide to training restricted Boltzmann machines[END_REF]. The contrastive divergence approach maximizes the log-likelihood (divided by M) of this dataset, which is given by

L = 1 M M i=1 log P(v (i) ). (156)
By using eqs. [START_REF] Pople | Nobel Lecture: Quantum chemical models[END_REF][START_REF] Ramakrishnan | Quantum chemistry structures and properties of 134 kilo molecules[END_REF] this becomes explicitly

L = 1 M M i=1 log {h} exp a T v (i) + b T h + v (i)T Wh -log Z = 1 M M i=1 a T v (i) + 1 M M i=1 log j h∈[0,1] exp h j ( k W kj v (i) k + b j ) -log Z = 1 M M i=1 a T v (i) + 1 M M i=1 j log 1 + exp k W kj v (i) k + b j -log Z. (157)
The derivatives of L with respect to the parameters a k , b k , W jk are

∂L ∂a k = 1 M M i=1 v (i) k - 1 Z {h,v} v k e -E(v,h) = ⟨v k ⟩ data -⟨v k ⟩ model , (158)
∂L ∂b k = 1 M M i=1 exp j W jk v (i) j + h k 1 + exp j W jk v (i) j + h k -⟨h k ⟩ model = 1 M M i=1 p(h k = 1|v (i) ) -⟨h k ⟩ model = ⟨h k ⟩ data -⟨h k ⟩ model (159)
∂L ∂W jk = 1 M M i=1 v (i) j p(h k = 1|v (i) ) -⟨v j h k ⟩ model = ⟨v j h k ⟩ data -⟨v j h k ⟩ model . ( 160
)
Given these equations the learning algorithm is based on the following iterative steps: (1) A minibatch of N b training elements v (i) is obtained from the dataset;

(2) The conditional probabilities p(

h k = 1|v (i) ) = 1 1+exp(-β(v (i)T W) k -b k)
a k = a k + ϵ (⟨v k ⟩ data -⟨v k ⟩ model ) (161)
b k = b k + ϵ (⟨h k ⟩ data -⟨h k ⟩ model ) (162)
W jk = W jk + ϵ ⟨v j h k ⟩ data -⟨v j h k ⟩ model . ( 163
)
With the simpler Boltzmann machine, the weights would have been updated by ⟨v j v k ⟩ data -⟨v j v k ⟩ model , thus modelling only the second order statistics of the inputs. To see how the RBM is instead able to go beyond second order, we can use the marginal probability P(v) to write

P(v) = 1 Z ′ e -E(v) = 1 Z {h} e -E(v,h) = e a T v Z {h} e b T h+v T Wh (164)
so that

E(v) = -e a T v -log h e b T h+v T Wh (165)
Let us now consider the probability distribution q a (h a ) = Z -1 e b a h a and its cumulant expansion K a (t) = log E[e th a ] = log q a (h a )e th a = n κ

(n) a t n
n! where the n-th cumulant is given by κ

(n) a = d n K a (0)
dt n . For t = (v T W) a we obtain

E(v) = -a T v - a,n κ (n) a ((v T W) a ) n n! (166)
Thus, by integrating out the hidden variable degree of freedom, one can recovers in the energy of eq. ( 166) all orders of the input correlations. 7 This explains the representative power of the RBM model, and of hidden variables in general [START_REF] David | Information Theory, Inference, and Learning Algorithms[END_REF][START_REF] Mehta | A high-bias, lowvariance introduction to Machine Learning for physicists[END_REF].

The use of machine learning (ML) techniques could be highly beneficial for MD simulations and, since the seminal work of Behler and Parrinello (BP) in 2007 [START_REF] Behler | Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces[END_REF], has seen an increasing popularity in this field [START_REF] Behler | Perspective: Machine learning potentials for atomistic simulations[END_REF][START_REF] Chmiela | Towards exact molecular dynamics simulations with machinelearned force fields[END_REF][START_REF] Chmiela | Machine learning of accurate energy-conserving molecular force fields[END_REF][START_REF] Zhang | Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics[END_REF]. While keeping a level of accuracy comparable to ab initio calculations, ML approaches can be used to replace most of the expensive quantum mechanical calculations with numerically cheap predictions. This allows for an increase in the system size and timescale that would be normally accessible by traditional ab initio MD simulations. However, ML models typically employed in MD simulations require a significant amount of data to be trained and this represents an issue for the most expensive quantum mechanical approximations.

In a previous work was applied a scheme that couples machine learning techniques with thermodynamic perturbation theory (MLPT) to compute enthalpies of adsorption of molecules in the zeolite chabazite at the RPA level of theory [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF] (the same methodology can be used in the context of free energies of activation [START_REF] Bucko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF][START_REF] Gesvandtnerova | Methanol carbonylation over acid mordenite: Insights from ab initio molecular dynamics and machine learning thermodynamic perturbation theory[END_REF]). Achieving these results with a brute force molecular dynamics simulation would be completely out of reach. Indeed, by considering only the cost of the RPA energy calculations for the zeolite, completing a 200000 steps MD simulation would require 65 million CPU hours and more than 110 years on 64 cores at 2.6 GHz (this estimate does not include the significant additional cost of computing RPA forces [START_REF] Ramberger | Analytic interatomic forces in the random phase approximation[END_REF]). With our MLPT technique, as few as 10 single point RPA energy calculations were sufficient to train a machine learning model that was then used to predict the RPA energy for several other configurations [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF].

In this previous work the van der Waals (vdW) corrected PBE+D2 functional [START_REF] Bucko | Improved Description of the Structure of Molecular and Layered Crystals: Ab Initio DFT Calculations with van der Waals Corrections[END_REF][START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF] was used to generate a production MD trajectory (in general, a numerically inexpensive approximation should be chosen as production method). Based on the ∆-ML method [START_REF] Ramakrishnan | Big Data Meets Quantum Chemistry Approximations: The ∆-Machine Learning Approach[END_REF], a model was thereafter trained to predict the difference between RPA and PBE+D2 energies. Following the thermodynamic perturbation approach [START_REF] Chipot | Free energy calculations: Theory and Applications in Chemistry and Biology[END_REF][START_REF] Pohorille | Good Practices in Free-Energy Calculations[END_REF][START_REF] Rocca | Bridging molecular dynamics and correlated wavefunction methods for accurate finite-temperature properties[END_REF], the energy differences were subsequently used to reweight the statistical weights of the PBE+D2 configurations in order to obtain the RPA canonical distribution and compute the RPA ensemble energies and enthalpies (in this context the RPA is considered the "target" level of theory). Thanks to the inexpensive predictions of the ML model, this technique involves a computational cost that is several orders of magnitude smaller than a full RPA MD simulation.

While the error involved in the ML procedure is rather well controlled [START_REF] Bucko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF][START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF], the application of thermodynamic perturbation theory (TPT) might be at the origin of a bias in the results estimated at the target level of theory. Indeed, as explained in subsection 2.4.3 the configurational space sampled by the production MD might have suboptimal superposition with the target configurational space [START_REF] Chipot | Free energy calculations: Theory and Applications in Chemistry and Biology[END_REF][START_REF] Pohorille | Good Practices in Free-Energy Calculations[END_REF]. In certain cases, this might lead to a strong loss of statistical significance when the contributions of the production configurations are reweighted to obtain the target distribution (namely only few configurations could contribute to the whole target level statistics).

In this work, the accuracy of the MLPT approach is assessed using the adsorption of molecules in zeolites as a test case. Indeed, many separation and adsorption processes require the improvement of the dedicated materials, and DFT calculations often provide valuable suggestions of optimized formulations in this regard [START_REF] Adil | Gas/vapour separation using ultra-microporous metal-organic frameworks: insights into the structure/separation relationship[END_REF][START_REF] Chibani | A DFT investigation of the adsorption of iodine compounds and water in H-, Na-, Ag-, and Cu-mordenite[END_REF][START_REF] Paul Hessou | Ab initio screening of zeolite Y formulations for efficient adsorption of thiophene in presence of benzene[END_REF][START_REF] Khalil | Biofuel purification: Coupling experimental and theoretical investigations for efficient separation of phenol from aromatics by zeolites[END_REF]. However, additional efforts have to be done to find a better compromise between accuracy of prediction and calculation cost. By considering five different DFT functionals, including generalized gradient approximation (GGA), meta-GGA, and non-local vdW corrected functionals, full MD simulations are performed to provide references values for the ensemble energies and enthalpies of adsorption. Each one of these five functionals is then used as a MLPT production method to obtain the energy/enthalpy estimates for the remaining four functionals, which are considered as target approximations.

By comparing MLPT estimates with reference results the level of accuracy of this approach is established. Some anomalous cases are found, for example when GGA functionals are used as starting point for meta-GGA functionals and vice versa (in certain cases the MLPT estimate of the ensemble energy can deviate by more than 8 kcal/mol). By using machine learning dimensionality reduction algorithms, a qualitative visualization of the relative distributions of the configurations from the different MD simulations is shown. This analysis confirms that the most problematic cases can be ascribed to a poor superposition of the configurational spaces. In order to detect possible failure of MLPT a diagnostic test is proposed, the I w index, which is evaluated from ML energy predictions and does not rely on the knowledge of any reference results.

We then propose a scheme to significantly improve results even when the I w coefficient is close to 0 (lowest production-target superposition). This approach, denoted as MLMC, is based on a Monte Carlo resampling of the target configurational space that reuses the machine learning model already trained for MLPT. Without any additional target level calculations, this scheme decreases even the largest deviations within the threshold of chemical accuracy.

Finally, the tools developed in this work are applied to analyze previous results based on the RPA [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF], where reference values cannot be produced. The relatively large values of the I w index hints that the PBE+D2 functional provides a reliable starting point for RPA target properties. The results are stable even if a full MLMC resampling of the configurational space is performed. This shows that, with a proper choice of the production approximation, MLPT allows for a quick and accurate estimate of target level properties. The MLPT approach opens the way to a more systematic application of accurate but expensive DFT and quantum chemical approximations in finite temperature simulations of materials.

In the following Section 4.2, the methodological part of this work and computational details are discussed. In Section 4.3, numerical results are presented and analyzed. Section 4.4 contains the conclusions.

methodology

The enthalpies of adsorption(∆ ads H) of the CH 4 molecule in protonated chabazite (HChab) and CO 2 molecule in siliceous chabazite (SiChab) were investigated by employing ab initio molecular dynamics (AIMD) and machine learning thermodynamic perturbation theory [START_REF] Bucko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF][START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF][START_REF] Gesvandtnerova | Methanol carbonylation over acid mordenite: Insights from ab initio molecular dynamics and machine learning thermodynamic perturbation theory[END_REF] (MLPT). The models of the adsorbed molecules are shown in Figure 13. The enthalpy of adsorption is defined by ∆ ads H = ∆ ads U + ∆ ads (pV). In the low coverage limit (few molecules are adsorbed), the change of pV is negligible between the adsorbed system and the substrat (pV(M@Z) ≈ pV(Z)), so that, assuming an ideal gaz behaviour for the adsorbate, ∆ ads (pV) = -k B T . The enthalpy of adsorption is then given by

∆ ads H(M@Z) = ⟨E(M@Z)⟩ -(⟨E(M)⟩ + ⟨E(Z)⟩) -k B T ( 167
)
where ⟨E⟩ is the internal energy computed as ensemble average of potential energy via the AIMD simulations or MLPT (M denotes the molecules, Z the zeolites, and M@Z the adsorbed system), k B is the Boltzmann constant, and T is the system temperature (equal to 300 K in all our simulations). The goal of this work is to benchmark the accuracy of MLPT by comparing its estimates with reference values obtained from full MD simulations. We chose five DFT functionals in order to cover a range of different characteristics but keeping also into account their numerical cost, as long and stable MD trajectories have to be produced for reference. Specifically, was selected: a generalized gradient approximation (GGA) functional, PBE [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] and its version with corrections for van der Waals interactions, PBE+D2 [START_REF] Bucko | Improved Description of the Structure of Molecular and Layered Crystals: Ab Initio DFT Calculations with van der Waals Corrections[END_REF][START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF]; a meta-GGA functional, SCAN [START_REF] Sun | Strongly Constrained and Appropriately Normed Semilocal Density Functional[END_REF][START_REF] Sun | Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional[END_REF]; a non-local vdW functional, vdW-DF2 [START_REF] Dion | Van der Waals Density Functional for General Geometries[END_REF][START_REF] Román-Pérez | Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes[END_REF]; a version of SCAN with non-local vdW corrections, SCAN+rVV10 [START_REF] Peng | Versatile van der Waals Density Functional Based on a Meta-Generalized Gradient Approximation[END_REF]. We will also discuss results based on the random phase approximation (RPA), as implemented in the VASP code [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | Ab initio molecular dynamics for liquid metals[END_REF]; due to the high numerical cost a direct MD simulation at this level of theory is completely out of reach and the MLPT method becomes instrumental to estimate finite temperature properties.

All AIMD simulations were performed within the NVT ensemble, whereby temperature was maintained by Andersen [START_REF] Hans | Molecular dynamics simulations at constant pressure and/or temperature[END_REF] thermostat with collision probability of 0.05. An integration time step of 0.5 fs was considered and the total simulation time of each MD run was 100 ps, i.e. ∼200•10 3 configurations were generated with each method for each system investigated in this work. The initial part of 10 ps of each trajectory was considered as an equilibration period and the corresponding data were discarded. The remaining part of the trajectory was tested for the absence of any drift via Mann-Kendall tests [START_REF] Schiferl | Statistical errors in molecular dynamics averages[END_REF]. All the standard errors on the ensemble energies of these trajectories are below 0.2 kcal/mol.

In the AIMD simulations a cell fixed at the lattice parameters of the chabazite optimized at the PBE level is used (a = b = c = 9.34 Å and α = β = γ = 95 o ). All AIMD simulations as well as single point calculations were carried out with the VASP electronic structure package [START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | Ab initio molecular dynamics for liquid metals[END_REF]. PAW pseudopotentials [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] for all atoms with the default kinetic energy cutoffs were used to expand the wavefunctions. The "accurate" precision setting was used and the hydrogen atomic mass was increased to 3.0 a.u. The Γ point approximation was applied in all calculations.

The results produced by AIMD are used as a reference for our MLPT method. Within this approach an MD simulation is first performed at a certain production level of theory and then "corrected" to obtain results at one or multiple target levels of theory. One of the advantages of MLPT is that an estimate of a finitetemperature property for the target method can be obtained with only a few tens or few hundreds of single point calculations [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]. Out of about 180•10 3 configurations generated in production runs, a reduced sample of about 19•10 3 evenly separated configurations (X i ) and corresponding potential energies obtained with the production method (E i ) are selected to be used for the predictions of the energies at the new target functional level. We note that this configuration reduction is not strictly necessary but, since MD configurations close in time are correlated, this procedure speeds up the predictions of the ML model without compromising the accuracy of prediction.

In order to compute target level ensemble energies ⟨E ′ ⟩, MLPT applies TPT to reweight the production statistical distribution [START_REF] Pohorille | Good Practices in Free-Energy Calculations[END_REF]

⟨E ′ ⟩ = M i E ′ i exp (-β∆E i ) M i exp (-β∆E i ) (168)
where E ′ i denotes the target level energies of the ith configuration, ∆E i = E ′ i -E i is the difference between target and production energies, β = 1/k B T , and M is the number of configurations, in this case 19•10 3 .

The application of eq. ( 168) still requires a large number (M) of calculations at the target level of theory. In MLPT only a very small number N train of these calculations is actually performed and then used to train an ML model that predicts all the remaining M -N train energies. In order to decrease the number of target calculations required to train the ML algorithm instead of considering directly the target energy E ′ i , we build a ML regression model using ∆E i . This approach, which is just one of the many possible applications of the ∆-ML idea [START_REF] Volker | Origins of structural and electronic transitions in disordered silicon[END_REF][START_REF] Ramakrishnan | Big Data Meets Quantum Chemistry Approximations: The ∆-Machine Learning Approach[END_REF], takes advantage of the much smoother dependence of ∆E i on E i in comparison with E ′ i and, accordingly, is easier to predict. To represent the configurations within our ML model we use the Smooth Overlap of Atomic Positions (SOAP) [START_REF] Bartók | On representing chemical environments[END_REF] descriptor, as implemented in the DScribe library [START_REF] Himanen | DScribe: Library of descriptors for machine learning in materials science[END_REF]. While several other descriptors for periodic materials have been proposed in the literature [START_REF] Behler | Atom-centered symmetry functions for constructing highdimensional neural network potentials[END_REF][START_REF] Casier | Hybrid localized graph kernel for machine learning energy-related properties of molecules and solids[END_REF][START_REF] Huo | Unified Representation of Molecules and Crystals for Machine Learning[END_REF], the choice of the SOAP descriptor provides already a satisfactory level of accuracy for MLPT applications [START_REF] Bucko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF][START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF][START_REF] Gesvandtnerova | Methanol carbonylation over acid mordenite: Insights from ab initio molecular dynamics and machine learning thermodynamic perturbation theory[END_REF].

The SOAP approach leads to a natural definition of a kernel for local atomic environments that can be subsequently used to define the global rematch kernel [START_REF] Stuke | Chemical diversity in molecular orbital energy predictions with kernel ridge regression[END_REF]. This kernel, which can be intuitively seen as a measure of similarity between configurations, is then used in the framework of the kernel ridge regression (KRR) [START_REF] Stuke | Chemical diversity in molecular orbital energy predictions with kernel ridge regression[END_REF] ML algorithm; the KRR implementation in the Scikit-learning package is used in this work [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. To train the ML model, 200 evenly distributed configurations are chosen from the 19•10 3 structures of the production trajectory. An additional independent set of 25 configurations is also selected to test the accuracy of the predictions of the ML model; for all the applications considered in this work the root mean square error (RMSE) in the prediction of the energy of single configurations is at most 0.3 kcal/mol. Single point calculations at the target level are performed exclusively for these 225 configurations. The 18775 (i.e., 19000-225) remaining values of ∆E i (and E ′ i = E i + ∆E i ) necessary to evaluate eq. ( 168) are inexpensively predicted by the ML model. This shows the clear advantage of using MLPT, in particular if it is necessary to consider multiple target level theories or highly expensive approximations (e.g., the RPA).

As discussed in the next sections, MLPT does not provide satisfactory results in certain cases. This happens for the estimate of SCAN and SCAN+rVV10 adsorption enthalpies from PBE production calculations. This issue is related to the unsatisfactory superposition of the target and production configurational spaces. To improve the results in these cases we have reused the ML model from the MLPT application to resample the target level configurational space using a Monte Carlo (MC) algorithm [START_REF] Betancourt | The Convergence of Markov Chain Monte Carlo Methods: From the Metropolis Method to Hamiltonian Monte Carlo[END_REF][START_REF] Kremer | Monte Carlo and molecular dynamics simulations of polymers[END_REF][START_REF] Meimaroglou | Review of Monte Carlo Methods for the Prediction of Distributed Molecular and Morphological Polymer Properties[END_REF][START_REF] Paquet | Molecular Dynamics, Monte Carlo Simulations, and Langevin Dynamics: A Computational Review[END_REF]. We will denote this procedure as MLMC.

In our MC procedure, multiple types of trial moves were used. First type of the move applied to all systems was based on a random displacement of atomic positions (x) along the velocity vector (v) drawn at random from the Maxwell-Boltzmann distribution corresponding to our target T of 300 K realized according to the formula x new = x old + v∆t, with ∆t set to 0.5 fs. In the second type of move, used only for the adsorbate+substrate systems (CO 2 @SiChab and CH 4 @HChab), an extra random translation (maximal magnitude of 0.5 Å) and rotation (by up to 35.0 deg.) of the molecular adsorbate was attempted. Further improvements in efficiency could probably be achieved by implementing more sophisticated sampling techniques. This is, however, beyond the scope of the present work, which primarily focuses on design of strategy to overcome the overlap problem occurring in the MLPT simulations.

The MC approach used to sample the canonical ensemble is equivalent to MD but is simpler to implement and does not require the calculation of forces (for certain methods forces are not available in most solid state physics software implementations). By using our in-house software, the MLMC approach recursively sample new configurations following these steps: (1) A new geometric configuration is randomly generated; (2) A production level ab initio calculation is performed to obtain E i for this new configuration; (3) the ML learning model is used to predict ∆E i and E ′ i ; (4) the new configuration is accepted or rejected according to the Metropolis criterion [START_REF] Betancourt | The Convergence of Markov Chain Monte Carlo Methods: From the Metropolis Method to Hamiltonian Monte Carlo[END_REF][START_REF] Kremer | Monte Carlo and molecular dynamics simulations of polymers[END_REF][START_REF] Meimaroglou | Review of Monte Carlo Methods for the Prediction of Distributed Molecular and Morphological Polymer Properties[END_REF][START_REF] Paquet | Molecular Dynamics, Monte Carlo Simulations, and Langevin Dynamics: A Computational Review[END_REF]. The advantage of MLMC is that with respect to MLPT no additional calculations are required at the target level of theory and this is particularly important when computationally expensive quantum chemical methods are used. The trade-off is that each step of the MC procedure requires a new production level calculation. In the applications of next section we have sampled 400•10 3 steps for each MLMC trajectory to get converged ensemble statistics (for all the systems considered in this work the standard error on the MC ensemble energies is at most 0.3 kcal/mol).

results and discussions

Assessing the accuracy of MLPT for different density functional approximations

The experimentally measured enthalpies of adsorption of the investigated systems are rather small, -4.06 kcal/mol and -5.02 kcal/mol for CH 4 @HChab and CO 2 @SiChab, respectively. Indeed, for both systems, the adsorption process is dominated by weak long-range van der Waals (vdW) interactions, which are, due to the small size of both adsorbate molecules, only modest. Ideally, very accurate methods such as RPA, MP2 or CCSD(T) would be necessary to model with systematic improvements in accuracy such processes in which long-range interactions play an important role. While the use of such methods is prohibitive due to their enormous computational cost, machine learning techniques could play an important role in extending their applicability to realistic systems [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]. This requires ML techniques that reasonably preserve the accuracy of these high-level quantum chemical methods and have sufficient predictive power for configurations beyond the training set. In this section we will discuss the accuracy of MLPT for the case of computationally affordable DFT functionals, for which a reference can be easily obtained. At the end of this section we will consider the RPA enthalpies of adsorption.

In Table 1, we report the values of the enthalpies of adsorption (eq. ( 167)) for CH 4 and CO 2 in zeolite as obtained from the full MD simulations. These results will be used as reference values for assessing the quality of our MLPT methodology. Comparing the calculated and the experimental values for the enthalpies of adsorption in Table 1, it can be noticed that only functionals with some vdW corrections can achieve a reasonable agreement with the experimental values, although the accuracy improvement is not systematic. For the CH 4 @HChab system, SCAN+rVV10 was identified as the best performing functional, with a deviation of -0.71 kcal/mol (17% of deviation), whereas for the CO 2 @SiChab the best functional was PBE+D2 with a deviation of -0.46 kcal/mol (9% of deviation). On the other hand, the uncorrected PBE and SCAN functionals show a sizable underestimation in the calculated enthalpies of adsorption for the CH 4 and CO 2 molecules in the zeolites, with relative errors with respect to experiment higher than 35%. However, the main purpose of this work is not to compare the accuracy of different functionals compared to experiment but to assess the performance and limitations of MLPT to reproduce reference AIMD results obtained using different quantum mechanical approximations.

Starting from each of the 5 trajectories obtained using PBE, PBE+D2, SCAN, vdW-DF2, and SCAN+rVV10 production methods, the MLPT approach is used .02 [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] to predict the results for all remaining approximations. The results are shown in Table 2. In this table, columns and lines represent the production and target methods, respectively. If, for example, we consider the PBE column (below Production -MD simulations), and the SCAN line, we find an estimate of the SCAN enthalpy of adsorption that is obtained by applying MLPT to the PBE production trajectory; this estimate avoids completely the generation of a new SCAN molecular dynamics trajectory and requires only 200 additional SCAN calculations that are used to train the ML model. For the sake of completeness, we also report values on the diagonal of this table, namely MLPT values based on the same production and target levels of theory. This shows that some small fluctuations can be introduced by the ML model even in this particularly simple case, additionally to those regular/general fluctuations due to the AIMD convergence and statistical uncertainty.

The accuracy of the MLPT estimates can be established by comparing the MLPT results with reference values of Table 1. The corresponding deviations are reported in parenthesis in Table 2. Altogether, it can be noticed that the predictions within the PBE/PBE+D2/vdW-DF2 group and the SCAN/SCAN+rVV10 are reasonable with deviations typically within few tenths of kcal/mol and only in one case slightly beyond chemical accuracy (a deviation of 1.01 kcal/mol for the prediction of SCAN+rVV10 from SCAN for CO 2 @SiChab). When the production and target methods belong to two different groups of functional, the MLPT is significantly less reliable and significantly large deviations appear. The most problematic cases involve the prediction of CH 4 @HChab enthalpy of adsorption at the SCAN and SCAN+rVV10 levels of theory from PBE, where the deviation is larger than 8 kcal/mol. As the enthalpies of adsorption are computed from differences of internal energies of interacting and non-interacting systems, fortuitous error cancellations might arise that could mask the real accuracy of MLPT predictions. To shed some light onto this problem, we report in Table 3 MLPT total ensemble energies for the individual components involved in the adsorption process (adsorbed system, zeolite alone, molecule alone). These results confirm the conclusions previously drawn from Table 2, with the predictions made for the cases where the production and target methods are both from the same group of functionals (i.e., either PBE/PBE+D2/vdW-DF2 or SCAN/SCAN+rVV10) being well within the chemical accuracy.

MLPT becomes unreliable when predictions mix these two groups, with deviations that often reach several kcal/mol (these observations do not hold for the standalone CH 4 and CO 2 molecules, whose deviations are always below 0.1 kcal/mol). The worst performance is confirmed to correspond to the prediction of SCAN and SCAN+rVV10 from the PBE production trajectory of CH 4 @HChab. For these cases, we notice that the largest deviations in the enthalpies of adsorption (>8 kcal/mol) are not only due to the low accuracy of the MLPT ensemble energy estimates for CH 4 @HChab and HChab, but also to the opposite sign in deviations of predictions made for these two systems, which leads to an error accumulation when evaluating eq. ( 167).

While MLPT allows, in principle, for a quick evaluation of finite-temperature properties at one or more target levels of theory from a single MD production run (typically based on the most computationally inexpensive approximation), the results in Tables 23show that such a strategy does not always allow one to achieve the required level of accuracy. Accordingly, special care should be taken in choosing a production method suitable for the target approximation(s) of interest. The failure of MLPT in certain cases can be explained by the limitations of the thermodynamic perturbation theory (TPT) itself, rather than by an inaccurate ML model. Indeed, if the configurational spaces visited with high likelihood by the production and target approximations do not overlap sufficiently, the TPT has low predictive power [START_REF] Pohorille | Good Practices in Free-Energy Calculations[END_REF]. In the most anomalous cases it can happen that this overlap is so poor that only one or few individual configurations effectively contribute to the ensemble average reweighted according to eq. ( 168). In order to analyze qualitatively this behavior we consider all the production MD simulations based on the 5 different functionals. By using the t-distributed stochastic neighbor embedding (t-SNE) algorithm [START_REF] De | Comparing molecules and solids across structural and alchemical space[END_REF][START_REF] Van | Visualizing Data using t-SNE[END_REF] we present in Figure 14 two-dimensional visualizations of the high-dimensional configurational space spanned by these trajectories (to improve the readability, figures were created using only 500 uncorrelated structures selected from each trajectory). To be consistent with the ML learning algorithm used in this study for the regression, the t-SNE approach was applied using a definition of distance D between two configurations χ A and χ B based on the normalized SOAP kernel K:

D(χ A , χ B ) = (1 -K(χ A , χ B )) 1 2 ; ( 169
)
this definition also inherits some of the properties of the SOAP kernel, such as the rotational and translational invariance required in materials and molecular modelling.

It can be noticed in Figure 14 that for CH 4 @HChab, CO 2 @SiChab, HChab, and SiChab the trajectories generated by the two groups of functionals (PBE/PBE+D2/vdW-DF2 and SCAN/SCAN+rVV10) form two clusters with very limited overlap. This shows that the configurational spaces spanned by the two groups of approximations are to a large extent different, which is at the origin of the poor performance of MLPT in the cases discussed above. For the gas phase molecules CH 4 and CO 2 , the configurational space has a much simpler structure and the five different approximations produce trajectories which always largely overlap; this is also consistent with the results in Table 3, where MLPT is highly accurate for all molecular energies. The t-SNE visualization provides useful insights into the potential failures of TPT and, consequently, MLPT. However, a diagnostic approach that requires the full MD trajectories is highly impractical in reality, since the goal of MLPT is indeed to completely avoid the MD simulation at the target level of theory. In order to overcome this difficulty and introduce a tool to analyze and predict the failures of MLPT, we follow the approach introduced in ref. [START_REF] Bucko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF] and introduce the I w index. By considering the weight factor (w i ) in eq. ( 168)

w i = exp(-β∆E i ) (170)
we define

I w = (M-N) M
where M is the total number of configurations sampled in the production run and N is the lowest integer which satisfies the condition

N i w i M j w j ⩾ 0.5, (171)
where the weights w i have been sorted in ascending order. In practice, the I w index gives the fraction of configurations that contribute, after the reweighting by eq. ( 168), one half of the total statistical weight.

In the case of an extremely small overlap between configurational spaces sampled by the production and target approximations, we would have N ≈ M (since most weights w i would be close to 0, N has to approach M to satisfy the condition in eq. ( 171)) and I w ≈ 0. This corresponds to a significant loss in the statistical significance since the large number of configurations sampled in the production run is reduced to a few individual structures. In the opposite case of excellent overlap we would have N ≈ M/2 and I w ≈ 0.5; if for example the production and target approximations provides energies that differ only by a constant, we would obtain I w = 0.5.

The I w values for all systems investigated in this work are presented in Table 4; in order to simplify the analysis, we are also reporting in parenthesis the deviations of the MLPT ensemble energies with respect to the MD references (same as in Table 3). As expected the I w index is equal to 0.5 when the production and target methods are the same. The values of I w allow for an interpretation of the deviations previously reported in Table 3. All the largest errors of MLPT (often significantly above chemical accuracy) correspond to I w = 0.00, namely to a small correlation/overlap of production and target methods. In few cases, small deviations are found for systems with I w = 0.00 (for example the MLPT evaluation of the SCAN ensemble energy from PBE+D2 for HChab) but these results must be considered as coincidental, whereby the reweighting of a very small number of configurations led, by chance, to a satisfactory reconstruction of the target statisti-cal distribution.

As discussed above, the MLPT estimates become reliable for the production and target method combinations from within the two groups of functionals PBE/PBE+D2/vdW-DF2 and SCAN/SCAN+rVV10. As shown in Table 4 the I w values within this two groups are greater than 0.03. For the isolated gas phase molecules CH 4 and CO 2 we always find large values of I w (>0.1) and consequently the errors on the energy predictions are small (within 0.18 kcal/mol). These remarks point to the fact that the I w can be used for a qualitative a posteriori estimate of the reliability of MLPT. In the case of an I w close to 0, the MLPT estimates are unreliable and often deviate by several kcal/mol from the reference values. Non-vanishing values of I w , ideally as close as possible to 0.5, correspond to accurate MLPT predictions. In the specific case of the systems considered here, an I w of about 0.03 already provides ensemble energies well within chemical accuracy. vdW-DF2 0.03 (-0.38) 0.05 (0.21) 0.00 (0.28) 0.50 (0.06) 0.00 (-3.80) SCAN+rVV10 0.00 (-5.03) 0.00 (-2.04) 0.34 (-0.22) 0.00 ( 1.01) 0.50 (-0.00) PBE 0.50 (0.01) 0.10 (-0.69) 0.00 (-2.95) 0.04 (-0.43) 0.00 (-3.69) PBE+D2 0.08 (0.24) 0.50 (0.04) 0.00 (-2.27) 0.12 (0.02) 0.00 (-1.67)

CO 2 @SiChab SCAN 0.00 ( 2.25) 0.00 ( 1.24) 0.50 (-0.01) 0.00 ( 2.45) 0.35 (-0.59) vdW-DF2 0.03 (-0.05) 0.13 (-0.62) 0.00 (-4.14) 0.50 (-0.00) 0.00 (-2.12) SCAN+rVV10 0.00 ( 3.39) 0.00 ( 2.29) 0.32 (0.90) 0.00 ( 3.25) 0.50 (0.02) PBE 0.50 (0.03) 0.31 (0.03) 0.00 (-3.71) 0.09 (-0.71) 0.00 (-2.95) PBE+D2 0.32 (0.05) 0.50 (0.00) 0.00 (-4.01) 0.11 (-0.63) 0.00 (-3.48) HChab SCAN 0.00 ( 4.14) 0.00 (0.08) 0.50 (0.02) 0.00 (0.72) 0.39 (-0.03) vdW-DF2 0.08 (0.23) 0.11 (0.45) 0.00 (-4.41) 0.50 (-0.04) 0.00 (-3.68) SCAN+rVV10 0.00 ( 3.83) 0.00 (-0.05) 0.39 (0.05) 0.00 (0.36) 0.50 (-0.03) PBE 0.50 (-0.01) 0.36 (-0.38) 0.00 (-2.63) 0.14 (-0.35) 0.00 (-0.23) PBE+D2 0.35 (0.36) 0.50 (0.00) 0.00 (-2.34) 0.14 (0.14) 0.00 (0.35) SiChab SCAN 0.00 (0.66) 0.00 ( 1.05) 0.50 (0.00) 0.00 ( 3.81) 0.40 (0.14) vdW-DF2 0.13 (0.45) 0.15 (0.01) 0.00 (-2.97) 0.50 (0.00) 0.00 (-0.23) SCAN+rVV10 0.00 (0.25) 0.00 (0.74) 0.40 (-0.16) 0.00 ( 3.36) 0.50 (0.00) PBE 0.50 (0.01) 0.46 (-0.05) 0.29 (-0.02) 0.33 (0.01) 0.28 (0.01) PBE+D2 0.46 (0.05) 0.50 (0.00) 0.26 (0.02) 0.29 (0.05) 0.25 (0.06) CH 4 SCAN 0.29 (0.01) 0.26 (-0.04) 0.50 (-0.00) 0.45 (0.03) 0.49 (0.01) vdW-DF2 0.33 (-0.01) 0.30 (-0.07) 0.45 (-0.03) 0.50 (-0.01) 0.44 (-0.01) SCAN+rVV10 0.29 (0.00) 0.25 (-0.05) 0.49 (-0.01) 0.44 (0.01) 0.50 (-0.00) PBE 0.50 (0.00) 0.50 (-0.02) 0.12 (-0.03) 0.48 (-0.02) 0.14 (0.04) PBE+D2 0.50 (0.02) 0.50 (-0.00) 0.11 (-0.01) 0.48 (-0.01) 0.14 (0.06) CO 2 SCAN 0.11 (0.05) 0.12 (-0.02) 0.50 (-0.01) 0.11 (-0.02) 0.45 (0.04) vdW-DF2 0.48 (0.03) 0.48 (0.01) 0.11 (0.00) 0.50 (-0.00) 0.13 (0.07) SCAN+rVV10 0.14 (0.01) 0.15 (-0.05) 0.45 (-0.04) 0.14 (-0.05) 0.50 (0.00)

Machine learning Monte Carlo resampling

Up to this point we have discussed the origin of MLPT inaccuracies occurring in certain cases. Ideally, the production trajectory should be chosen to suit the target approximations of interest. However, this is not always possible and it is of interest to improve MLPT results also in the most anomalous cases with I w ≈ 0. Recently, Rizzi et al. [START_REF] Rizzi | Targeted Free Energy Perturbation Revisited: Accurate Free Energies from Mapped Reference Potentials[END_REF] have discussed a ML-based approach to overcome the limitations of TPT for free energy calculations in cases where there is none or very limited overlap of the configurational spaces of production and target methodologies. Following the previous approach presented by Wirnsberger et al. [192], Rizzi et al. [START_REF] Rizzi | Targeted Free Energy Perturbation Revisited: Accurate Free Energies from Mapped Reference Potentials[END_REF] reported the use of a configurational space transformation built on a neural network model.

While this methodology is appealing, in this work we decided to pursue an approach based on the Monte Carlo resampling of the trajectory. This machine learning Monte Carlo (MLMC) method has the following advantages: The implementation is straightforward and does not require the calculation of forces, which often are not implemented in solid state physics codes for methods such as RPA, MP2, or CC; the same ML model of MLPT is reused and no additional calculations at the target level of theory are required; the MC resampling can be used also to verify results when a good or reasonable value of I w is found. Finally, let us remark that the kernel methods used in our MLPT and MLMC methods typically require much smaller number of training configurations than neural networks [START_REF] Maik | Big-Data Science in Porous Materials: Materials Genomics and Machine Learning[END_REF] used in the method of Wirnsberger et al. [START_REF] Wirnsberger | Targeted free energy estimation via learned mappings[END_REF].

To test the MLMC approach, we consider the prediction of SCAN and SCAN+rVV10 enthalpies of adsorption from PBE production calculations which showed the largest deviations between MLPT estimates and reference values. In Figure 15 we show the following probability distributions of the SCAN energy determined for each zeolite and adsorbed system: the reference from the SCAN molecular dynamics (in blue labeled as reference); the distribution of the ML predicted SCAN energies for the configurations sampled by the PBE AIMD (denoted as ML in green); the distribution of SCAN energies from the MLMC resampling (in red). The reference and ML distributions show an unsatisfactory overlap, further confirming the unsatisfactory sampling of PBE to predict SCAN properties. Since no reweighting of energies (such as in eq. ( 168)) has been performed to generate the histograms shown in Figure 15, the ML distributions are strongly biased, and therefore deviate significantly from the reference MD results. The MLMC approach largely overcomes this issue, as it can be noticed from the good overlap of the blue and red curves. Computing the partial radial distribution function for the Si-O pairs, Figure 16, it is also evident that the MLMC generates geometries which are in a better agreement with those obtained from straightforward MD simulations with the target functional SCAN/SCAN+rVV10. The improvement in the geometric sampling can be also demonstrated by including the MC data in the t-SNE analysis. Figure 17 shows indeed that the clusters of the SCAN and SCAN+rVV10 configurations sampled by the Monte Carlo have an excellent overlap with those obtained from the reference MD simulation performed at the SCAN level. In this context it is important to emphasize again that the MLMC results are obtained without any additional explicit target level calculation. We have shown insofar that the MLMC can significantly improve the sampling of the geometries at the target level. The MLMC also clearly reduces error in computed internal energies of individual systems, as shown in Table 5. For instance, the deviation in the internal energy of the CH 4 @HChab at the SCAN level was reduced from -4.38 kcal/mol from MLPT with PBE production method to about -0.64 kcal/mol. Similarly, the deviations were also reduced for the predicted SCAN and SCAN+rVV10 enthalpies of adsorptions, see Table 6. For CH 4 @HChab enthalpy, for instance, the MLPT presented a deviation of -8.58 kcal/mol which reduced to only 0.12 kcal/mol when the MLMC resampling was employed.

Once again, the greatest improvements have been achieved for the adsorbate with substrate systems, followed by the clean zeolites, while the effect of MLPT or MLMC is negligible for the isolated molecular systems for which most of the observed deviations are close to zero and probably are related to numerical errors. In Tables 56, all the predictions of SCAN/SCAN+rVV10 energies from PBE were significantly improved by the MLMC method. The general improvement reduced the relative internal energy deviations below 1.0 kcal/mol. Exp. -4.06 [START_REF] Piccini | Accurate Adsorption Thermodynamics of Small Alkanes in Zeolites. Ab initio Theory and Experiment for H-Chabazite[END_REF] SCAN+rVV10 -4.77 -13.63 (-8.86) -5.38 (-0.61) CO 2 @SiChab SCAN -2.95 -1.42 (1.54) -3.17 (-0.21) Exp. -5.02 [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] SCAN+rVV10 -7.11 -3.97 (3.14) -5.82 (1.29)

Assessing the accuracy of MLPT for the random phase approximation

In the previous work referenced earlier, (ref. [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]), MLPT has been applied to compute enthalpies of adsorption at the RPA level of theory. In addition to the two systems studied in this work, also CH 4 @CHAB and CO 2 @HChab were considered in this previous paper. The RPA enthalpies of adsorption showed an excellent agreement with experiment and for the molecules adsorbed in HChab zeolite a significant improvement with respect to the production method (PBE+D2) was found. Because of the high computational cost involved in the RPA calculations, in this case it is not possible to produce reference trajectories and use them to repeat some of the previous analysis (indeed, MLPT is meant to make finite-temperature calculations involving computationally expensive approximations possible).

However, the I w index can be straightforwardly obtained and we present its values in Table 7. This type of analysis, which we did not develop in the original paper, ref. [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF], shows that a good performance of MLPT should be expected, as the I w coefficient takes values above 0.1 for all systems. While not strictly necessary, the accuracy can be further confirmed by applying MLMC, which requires several additional production level calculations but avoids completely the highly expensive RPA calculations. To this purpose, we consider all the systems but CO 2 @SiChab, which is the adsorbed system with the highest I w coefficient. (the two adsorbed systems with the lowest I w and the two zeolites). With respect to the MLPT values reported in ref. [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF], the MLMC results present only small deviations: from -0.05 kcal/mol to 0.66 kcal/mol for the studied systems (see Table 7). Keeping in mind that also statistical errors should be considered (the standard error in the MC sampling is about 0.3 kcal/mol), we conclude that the agreement between MLPT and MLMC results is very good from which we deduce that PBE+D2 is a suitable production method for the RPA target level calculations and that the Monte Carlo resampling of the MLPT results is not strictly necessary.

conclusions

We assessed the accuracy of the machine learning thermodynamic perturbation theory approach by computing enthalpies of adsorption of molecules in zeolites. By considering a set of five DFT functionals with different characteristics, MLPT estimates are compared to reference results produced from full MD simulations. The largest deviations are found when the configurational space accessible to the target level has a small overlap with that sampled by the production method. Even when a reference result is not available, the most problematic cases can be detected by using a diagnostic tool, the I w index. In these cases the results can be significantly improved by coupling the Monte Carlo approach to the MLPT model to resample the configurational space at the target level of theory which is, however, performed at the cost of the production method.

Finally, the analysis is extended to some recently published MLPT results that used the RPA as target approximation [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]; in this case the high computational cost of the RPA makes the generation of reference ensemble energies and enthalpies completely unpractical. The application of the I w diagnostic test and of the MC resampling highlights the full reliability of those previous results, which were based on PBE+D2 production MD.

The main advantage of the MLPT approach stands in the possibility of evaluating finite temperature properties at a certain target level of theory by performing only a minimal number of single-point calculations. Additionally, MLPT is also suitable for approximations where the computation of gradients is not available or adds a significant computational cost [START_REF] Ramberger | Analytic interatomic forces in the random phase approximation[END_REF]. This approach opens the possibility of systematically applying high-accuracy/high-cost ab initio methodologies to achieve a new level of predictive power in materials simulations. The content of this chapter is nearly identical to the following Chemrxiv preprint: "Coupled cluster finite temperature simulations of periodic materials via machine learning", by Basile Herzog, Alejandro Gallo, Felix Hummel, Michael Badawi, Tomas Bucko, Sébastien Lebègue, Andreas Grüneis and Dario Rocca, to be published.

C O U P L E D C L U S T E R F I N I T E T E M P E R AT U R E S I M U L AT I O N S O F P E R I O D I C M AT E R I A L S V I A M A C H I N E L E A R N I N G

introduction

The last chapter was devoted to a methodological assessment of the Machine Learning Perturbation Theory. In this chapter, we will be interested in applying those methods to provide correlated quantum chemical wave function finite temperature results. The later, based on post-Hartree-Fock (post-HF) approximations are indeed systematically improvable and could potentially overcome some of the limitations of DFT for materials simulations. Among those, second-order Møller-Plesset perturbation theory (MP2) [137] and coupled cluster theory [START_REF] Bartlett | Coupled-cluster theory in quantum chemistry[END_REF], have been recently implemented for periodic materials [START_REF] Booth | Towards an exact description of electronic wavefunctions in real solids[END_REF][START_REF] Del Ben | Second-Order Møller-Plesset Perturbation Theory in the Condensed Phase: An Efficient and Massively Parallel Gaussian and Plane Waves Approach[END_REF][START_REF] Dixit | Communication: A novel implementation to compute MP2 correlation energies without basis set superposition errors and complete basis set extrapolation[END_REF][START_REF] Marsman | Second-order Møller-Plesset perturbation theory applied to extended systems. I. Within the projectoraugmented-wave formalism using a plane wave basis set[END_REF][START_REF] Pisani | Periodic local MP2 method for the study of electronic correlation in crystals: Theory and preliminary applications[END_REF]. However, their computational cost is significant for most practical applications in materials science and this issue becomes even more dramatic when finite temperature effects have to be included by performing molecular dynamics (MD) simulations or Monte Carlo sampling. For example, a brute-force computation of the enthalpy of adsorption considered in this work would require billions of CPU hours and hundreds of real time years to be completed.

While machine learning (ML) is nowadays a well established tool in the context of MD simulations, using notably Neural Network Potentials [START_REF] Behler | Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces[END_REF][START_REF] Chmiela | Towards exact molecular dynamics simulations with machinelearned force fields[END_REF][START_REF] Chmiela | Machine learning of accurate energy-conserving molecular force fields[END_REF][START_REF] Zhang | Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics[END_REF], the ML-accelerated MD typically requires large amounts of data and becomes rapidly challenging for the more expensive approximations. In this chapter, we show how finite temperature observables for periodic materials can be evaluated using the 'gold standard' coupled cluster ansatz, including single, double, and perturbative triple particle-hole excitation operators (CCSD(T)) in combination with Machine Learning Perturbation Theory and Monte Carlo sampling.

For molecular systems, the application of ML techniques has already been proven to be effective in enhancing the efficiency of CCSD(T) MD simulations [START_REF] Bogojeski | Quantum chemical accuracy from density functional approximations via machine learning[END_REF][START_REF] Chmiela | Machine learning of accurate energy-conserving molecular force fields[END_REF][START_REF] Huziel E Sauceda | Molecular force fields with gradient-domain machine learning: Construction and application to dynamics of small molecules with coupled cluster forces[END_REF][START_REF] Schran | Converged colored noise path integral molecular dynamics study of the zundel cation down to ultralow temperatures at coupled cluster accuracy[END_REF][START_REF] Smith | The ANI-1ccx and ANI-1x data sets, coupled-cluster and density functional theory properties for molecules[END_REF]. Very recently, applications to molecular condensed phase systems, specifically to liquid water, have also been considered. In ref. [START_REF] Daru | Coupled Cluster Molecular Dynamics of Condensed Phase Systems Enabled by Machine Learning Potentials: Liquid Water Benchmark[END_REF] the ML model for periodic water was trained with data produced for finite water clusters using near-linear scaling coupled cluster theory. In ref. [START_REF] Chen | Data-Efficient Machine Learning Potentials from Transfer Learning of Periodic Correlated Electronic Structure Methods: Liquid Water at AFQMC, CCSD, and CCSD(T) Accuracy[END_REF], CCSD(T) calculations were restricted to very small periodic models based on a box of 16 H 2 O molecules and the ML model was then used to compute radial distribution functions, diffusion coefficients, and vibrational densities of states. To the best of our knowledge, this work represents the first application of CCSD(T) to finite temperature simulations of periodic solid materials. This was challenging for several reasons. First, we dealt with a system containing more than 200 electrons, far more than used in any previous report on ML-assisted MD CCSD(T) simulations. Second, we focused on a measurable thermodynamic quantity, the enthalpy of adsorption, whose prediction imposes high demands on the quality of the ML model. This is because any error in energy of a configuration affects not only the underlying phase space function used in ensemble averaging, but also the statistical weight of that contribution (see eq. ( 138)).

This work is based on the combination of the Machine Learning Perturbation Theory (MLPT), together with an efficient periodic coupled cluster theory implementation. This implementation is based on a plane-wave basis set and novel finite size and basis set correction techniques that accelerate the convergence to the complete-basis-set limit and thermodynamic limit significantly [START_REF] Gruber | Applying the coupled-cluster ansatz to solids and surfaces in the thermodynamic limit[END_REF][START_REF] Irmler | Focal-point approach with pair-specific cusp correction for coupled-cluster theory[END_REF]. Using these techniques, it is possible to obtain well converged correlation energies at the CCSD(T)-level of theory for periodic solids and surfaces containing more than 100 electrons on modern supercomputers [START_REF] Liao | A comparative study using state-of-the-art electronic structure theories on solid hydrogen phases under high pressures[END_REF][START_REF] Liao | Structural and electronic properties of solid molecular hydrogen from manyelectron theories[END_REF][START_REF] Schäfer | Surface science using coupled cluster theory via local Wannier functions and in-RPA-embedding: the case of water on graphitic carbon nitride[END_REF][START_REF] Tsatsoulis | Reaction energetics of hydrogen on Si (100) surface: A periodic many-electron theory study[END_REF].

As a specific application for our approach we consider the calculation of the enthalpy of adsorption of carbon dioxide in protonated chabazite (HChab). The adsorption of molecules in zeolites is fundamental for many applications, including depollution, separation of chemicals, and catalysis [START_REF] Grajciar | Towards operando computational modeling in heterogeneous catalysis[END_REF][START_REF] Maik | Big-data science in porous materials: materials genomics and machine learning[END_REF][START_REF] Van Speybroeck | Advances in theory and their application within the field of zeolite chemistry[END_REF]. In this field, more quantitative and systematically improvable theoretical predictions are instrumental to interpret experimental findings and predict new materials. Although the calculations presented in this work are still significantly more expensive than those based on standard density functional theory, this proof-of-principle work paves the way to a more systematic use of highly accurate post-HF methods in materials simulations.

methods

Periodic boundary conditions coupled cluster

The coupled cluster theory calculations are performed using the Cc4s code [START_REF]Cc4s User Documentation[END_REF], which is interfaced to the Vienna ab initio simulation package (VASP) [START_REF] Kresse | Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements[END_REF]. The calculations are performed in several steps involving Hartree-Fock and MP2 theory to obtain corresponding energies and optimized approximate natural orbitals [START_REF] Grüneis | Natural orbitals for wave function based correlated calculations using a plane wave basis set[END_REF]. Once the natural orbitals have been computed, the Cc4s interface to VASP is used to compute intermediate quantities [START_REF] Hummel | Low rank factorization of the Coulomb integrals for periodic coupled cluster theory[END_REF] that are needed for the subsequent coupled cluster energy calculations including the corresponding finite size [START_REF] Gruber | Applying the coupled-cluster ansatz to solids and surfaces in the thermodynamic limit[END_REF] and basis set corrections [START_REF] Irmler | Focal-point approach with pair-specific cusp correction for coupled-cluster theory[END_REF]. In ref. [START_REF] Schäfer | Surface science using coupled cluster theory via local Wannier functions and in-RPA-embedding: the case of water on graphitic carbon nitride[END_REF], all individual steps are described when combined with an embedding approach, which was not necessary for the present system due to its relatively small unit cell containing up to 40 atoms only. In the present calculation, 10 unoccupied approximate natural orbitals per occupied orbital are used for the CCSD calculations, whereas only 5 unoccupied approximate natural orbitals per occupied orbital are employed to evaluate the (T) contribution. A single CCSD(T) calculation for the given structures containing up to 40 atoms took about 10,000 core hours.

The convergence of the CCSD and (T) correlation energy contributions to the molecular adsorption energy was tested on a single configuration. In particular, the following difference has been considered: ∆E corr. = E corr.

CO 2 -HChab -E corr. CO 2 -E corr. HChab , where E corr.

CO 2 -HChab , E corr. CO 2 , and E corr. HChab are, respectively, the correlation energies of the chabazite containing the CO 2 molecule, the isolated CO 2 molecule and the pristine chabazite. Table 8 lists the ∆E corr. values obtained for different basis set sizes and methods. In particular, the CCSD results with (CCSD-FP) and without (CCSD) a focal point finite basis set error correction and the perturbative triples contributions with (T * ) and without (T) an approximate finite basis set error correction have been compared. We stress that CCSD-FP and (T * ) have been thoroughly tested in ref. [START_REF] Irmler | Focal-point approach with pair-specific cusp correction for coupled-cluster theory[END_REF]. As discussed in ref. [START_REF] Irmler | Focal-point approach with pair-specific cusp correction for coupled-cluster theory[END_REF], 10 virtual natural orbitals per occupied state suffice in combination with the focal point basis set correction to achieve chemical accuracy for a large number of reaction energies on the level of CCSD and (T). Although we can not fully converge our binding energy estimates for the present system, we find satisfactory convergence of the correlation energy contributions with respect to the number of virtual orbitals. From the second-largest to the largest basis set size the CCSD-FP and (T*) estimates change by 0.5 kcal/mol and 0.4 kcal/mol only. From this we conclude that the remaining basis set error when using N v /N o = 10 for CCSD-FP and N v /N o = 5 for (T * ) will not exceed 1 kcal/mol. CCSD refers to the CCSD correlation energy, whereas CCSD-FP includes a focal point finite basis set error correction, as thoroughly tested in ref. [START_REF] Irmler | Focal-point approach with pair-specific cusp correction for coupled-cluster theory[END_REF]. (T) and (T * ) correlation energy estimates are without and with an approximate finite basis set error correction, respectively.

N v /N o CCSD CCSD-FP (T) (T * ) 5

Machine Learning Thermodynamic Perturbation Theory and Monte Carlo resampling

In this work we consider the calculation of the enthalpy of adsorption of carbon dioxide in a porous zeolitic material (see Section 5.3). In practice this quantity is computed as

∆ ads H(M@Z) = ∆ ads U(M@Z) + ∆ ads (pV)(M@Z) = ⟨E(M@Z)⟩ -(⟨E(M)⟩ + ⟨E(Z)⟩) -k B T , (172)
where ∆ ads U is the internal energy of adsorption, ⟨E(i)⟩ denotes the ensemble average of total energy of the system i corresponding to a gas phase molecule (M), clean zeolite (Z), and the adsorbed system (M@Z), and the identity ∆ ads (pV)(M@Z) = -k B T is obtained by assuming an ideal gas behavior of M and a negligible change of pV of the zeolite due to adsorption. The canonical ensemble energy can be evaluated by directly performing an ab initio molecular dynamics (AIMD) simulation but, because of the high computational cost of CCSD(T) and MP2, this approach is unpractical at these levels of theory. Starting from an AIMD trajectory obtained using numerically affordable semi-local DFT with empirical van der Waals corrections (PBE+D2) [START_REF] Bucko | Manybody dispersion corrections for periodic systems: an efficient reciprocal space implementation[END_REF][START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF], the post-HF ensemble energies are estimated using the MLPT approach trained on a small number of single point calculations. This approach was described in Chapter 4 and the two main steps are summarized here:

1. Given a set of configurations {R i } M i=1 from an AIMD trajectory in an NVT ensemble with the PBE+D2 reference Hamiltonian H 0 and potential energy E 0 , the ensemble average energy generated by the target Hamiltonian H 1 with potential energy E 1 can be obtained from thermodynamic perturbation theory by reweighting:

⟨E 1 ⟩ 1 = M i=1 E 1 (R i ) exp(-β∆E(R i )) M i=1 exp(-β∆E(R i )) , (173)
where ∆E(R) denotes the energy difference E 1 (R) -E 0 (R) for a specific atomic configuration R. In this work, E 1 denotes either the MP2 or the CCSD(T) target method potential energy. The trajectory obtained with the reference Hamiltonian is called production trajectory.

2. While the application of eq. ( 173) requires a large number of high-level calculations, in practice those can be largely replaced by inexpensive predictions of a machine learning model. MLPT limits the amount of data required for the training by using efficient algorithms based on the kernel ridge regression with the SOAP kernel [START_REF] Albert P Bartók | On representing chemical environments[END_REF][START_REF] De | Comparing molecules and solids across structural and alchemical space[END_REF] and ∆-ML [START_REF] Ramakrishnan | Big data meets quantum chemistry approximations: The ∆-machine learning approach[END_REF]. E 0 (R) is known and the evaluation of eq. ( 173) requires only the energy difference ∆E(R). In this work the training set is based on 100 uncorrelated configurations evenly spaced along the PBE+D2 trajectories and 10 randomly chosen configurations for the test set. The MP2 and CCSD(T) calculations are performed only for those selected geometries.

In order to correct for the basis set incompleteness, the (T * ) correlation energy was computed for 25 configurations of each system. The standard deviation of the difference between the (T * ) and (T) correlation energies across the 25 configurations of the M@Z system was found to be very small (within 0.19 kcal/mol). Hence, to a very good accuracy, (T * ) can be obtained from (T) by a uniform shift of contribution of each configuration by the average difference between the (T * ) and (T) energies obtained using 25 configurations. We applied this procedure to obtain our final CCSD(T) MLPT result. To validate this approach, we computed the MLPT enthalpy result using only those 25 configurations as training set, both with the (T*) energy and with the (T) energy plus the uniform shift, and found the two results to agree within 0.02 kcal/mol.

Since MLPT is based on thermodynamic perturbation theory, a limited overlap between the production and target configurational spaces can lead to inaccurate results. If a suboptimal overlap is suspected, a Monte Carlo (MC) resampling can be performed. This procedure, described in detail in Chapter 4, uses Metropolis MC [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] to resample the canonical ensemble at the CCSD(T) and MP2 levels of theory. At each MC step, configurational energies are computed with the production approximation (PBE+D2), and subsequently evaluated at the post-HF level level using the same ML model of MLPT. The Metropolis acceptance criterion is applied at the target level of theory and, accordingly, the correct target configurational space is sampled without bias from the starting point.

results and discussion

In this section we present and discuss the adsorption enthalpies of CO 2 in protonated chabazite as computed at the MP2 and CCSD(T) levels of theory. The latter approximation is commonly described as the 'gold standard' of quantum chemical simulations and is routinely used to produce reference test sets to benchmark the accuracy of other methods [START_REF] Mutasem | Estimates of the ab initio limit for π-π interactions: The benzene dimer[END_REF][START_REF] Takatani | Basis set consistent revision of the S22 test set of noncovalent interaction energies[END_REF]. The primitive cell of the model considered here is shown in Figure 18. The experimental value of the enthalpy of adsorption of CO 2 in HChab, -8.41 kcal/mol [START_REF] Trong D Pham | Carbon dioxide and nitrogen adsorption on cation-exchanged SSZ-13 zeolites[END_REF], is used as a reference for the computational results. This experimental estimate is obtained by extrapolating measurements to the zero coverage limit. The errors possibly arising from this procedure are not discussed in ref. [START_REF] Trong D Pham | Carbon dioxide and nitrogen adsorption on cation-exchanged SSZ-13 zeolites[END_REF] and we cannot exactly quantify the uncertainty in the experimental reference.

Al

The computed results are presented in Table 9, where the error bars related to the finite sampling and the ML model are also indicated [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]. The molecular dynamics at PBE+D2 level leads to an estimate for the adsorption energy which is more than 1 kcal/mol below the experimental value, corresponding to a deviation well beyond chemical accuracy. This MD trajectory is used as a starting point for MLPT to obtain post-HF enthalpies. Similarly, the MP2 approximation obtained from MLPT also tends to overbind and leads to results that do not qualitatively differ from PBE+D2. This is not surprising and we believe that this overestimation is caused by the lack of screening of long-ranged correlation effects in MP2 theory. The computational estimate of the enthalpy significantly improves at the CCSD(T) level, which provides a value in excellent agreement with experiment. This result demonstrates the high accuracy and predictive power of the CCSD(T) approximation also for finite temperature simulations of materials.

In a previous work, it was demonstrated that also the RPA provides accurate enthalpies of adsorption of molecules in zeolites [START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]. Specifically, the value for CO 2 in protonated chabazite is -8.01 kcal/mol. Although the RPA has a diagrammatic structure it is not as straightforward to systematically improve its accuracy as for post-HF methods [START_REF] Bates | Communication: Random phase approximation renormalized many-body perturbation theory[END_REF][START_REF] Dixit | Improving the accuracy of ground-state correlation energies within a plane-wave basis set: The electron-hole exchange kernel[END_REF][START_REF] Furche | Fluctuation-dissipation theorem densityfunctional theory[END_REF][START_REF] Grüneis | Making the random phase approximation to electronic correlation accurate[END_REF][START_REF] Hellgren | Beyond the random phase approximation with a local exchange vertex[END_REF][START_REF] Hummel | Screened Exchange Corrections to the Random Phase Approximation from Many-Body Perturbation Theory[END_REF][START_REF] Olsen | Extending the random-phase approximation for electronic correlation energies: The renormalized adiabatic local density approximation[END_REF][START_REF] Xinguo Ren | Beyond the Random-Phase Approximation for the Electron Correlation Energy: The Importance of Single Excitations[END_REF]. In practice, the RPA often provides more realistic results starting from a DFT approximation rather than from HF [START_REF] Xinguo Ren | Beyond the Random-Phase Approximation for the Electron Correlation Energy: The Importance of Single Excitations[END_REF], and this starting point dependence makes this approximation less reliable as a general predictive method. To fully prove the accuracy of the MLPT approach for MP2 and CCSD(T) a crucial point concerns the reliability of the PBE+D2 trajectory used as starting point for thermodynamic perturbation theory. Specifically, if the target (MP2 or CCSD(T)) configurational space has a small overlap with the production (PBE+D2) configurational space, the results of TPT may be affected by a strong systematic error. As discussed thoroughly in Chapter 4 for systems similar to the one considered here, the occurrence of this issue can be identified even if the exact target trajectory is unknown. Thermodynamic perturbation theory is based on the reweighting of the statistics sampled by the production trajectory to obtain the target level statistics (see eq. ( 173)); in case of a poor overlap only few configurations contribute to the total weight, leading to poor ensemble estimates. In practice, this effect can be measured by the I w index, as defined in Chapter 4. This index assumes the value of 0.5 in the optimal configuration overlap case and tends to 0 for decreasing overlaps. For the adsorption of molecules in zeolites it has been shown that even relatively small values of I w around 0.03 -0.05 still allow for reliable MLPT estimates. The reweighting of the trajectories at the MP2 level provides large values for I w (> 0.15) and the corresponding enthalpies in Table 9 should be considered fully reliable. For the CCSD level, a very low I w value for the adsorbed system (0.008) precludes making any reliable predictions of adsorption enthalpy; for this reason this level of theory is not discussed here. For the CCSD(T) level of theory, the I w coefficient is one order of magnitude higher: 0.07 for HChab and 0.05 for the adsorbed system, indicating a better match between with the PBE+D2 equilibrium structure as compared to the CCSD level. While these I w values are likely to be sufficient to confirm the reliability of our results, considering the pioneering nature of our work and the lack of any previous finite temperature benchmark results for periodic CCSD(T), we further investigated the robustness of the MLPT estimate by resampling the CCSD(T) trajectory. This is achieved by performing a Metropolis Monte Carlo (MC) sampling of the canonical ensemble at the CCSD(T) level by replacing the expensive coupled cluster calculations with the predictions of the same machine learning model previously trained for MLPT. Differently from most machine learning-based MD approaches [START_REF] Behler | Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces[END_REF][START_REF] Chmiela | Towards exact molecular dynamics simulations with machinelearned force fields[END_REF][START_REF] Chmiela | Machine learning of accurate energy-conserving molecular force fields[END_REF][START_REF] Zhang | Deep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum Mechanics[END_REF], this MLMC approach avoids training on atomic forces, which are not readily available in the current periodic CCSD(T) implementation and that would require a significant overhead cost. Since thermodynamic perturbation theory is not used and a new trajectory is instead sampled from scratch, MLMC avoids the starting point bias. The corresponding result for the enthalpy of adsorption, shown in Table 9, differ by only 0.2 kcal/mol from the MLPT value. In the MLMC case the error bar is however sizeably larger because of the long auto-correlation length of this trajectory (about a factor 10 longer than for the MD trajectory) but this is sufficient to support our conclusion that the PBE+D2 trajectory provides a reliable starting point to compute CCSD(T) ensemble energies. This is also qualitatively confirmed by visualizing the (high-dimensional) geometries sampled by the MD and MC methods with the t-distributed stochastic neighbor embedding (t-SNE) algorithm [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF]. As shown in Figure 19, the PBE+D2 molecular dynamics and the CCSD(T) Monte Carlo trajectories span configurational spaces that overlap well. This figure also demonstrates that the training set provides a rather uniform sampling of the data, as required for a balanced training of the ML model.

To further analyze the overlap between the configurational space of the PBE+D2 functional and of the post-HF methods we consider the structure of the protonated chabazite cage. For this purpose, the radial distribution function of the Si-O pairs has been computed for the PBE+D2 molecular dynamics trajectory and for MP2 and CCSD(T) approaches using MLPT and MLMC. As previously shown in Chapter 4, the most spectacular failures of MLPT are encountered when the production approximation predicts equilibrium distances of covalent bonds that differ from the target theory; this translates to very different configurational spaces and fully unreliable perturbative estimates. For protonated chabazite, Figure 20 clearly shows that the radial distribution functions computed for the Si-O pairs are similar at different levels of theory and problematic behaviors of MLPT should not be expected.

conclusions

We have presented an application of CCSD(T) to compute the enthalpy of adsorption of carbon dioxide in a periodic model of zeolite. Due to the high computational cost, applications of CCSD(T) to periodic materials are so far limited and direct calculations of finite temperature observables are unpractical in terms of required computational resources and execution time. Here we showed that these challenges can be overcome by coupling machine learning models requiring small training sets with an efficient implementation of periodic coupled cluster theory.

The computed enthalpy of adsorption of carbon dioxide in protonated chabazite was found to be in excellent agreement with experiment. While still significantly more expensive than approaches based on density functional theory, our pioneering work opens the door to more reliable and predictive simulations of materials in finite temperature conditions. Future work will be aimed at demonstrating the accuracy of ML-based CCSD(T) in broader classes of problems, including, for example, the computation of free energies of activation, which play a fundamental role in the modelling of catalytic reactions. The content of this chapter is nearly identical to the following published work: The application of machine learning (ML) to quantum chemistry and computational materials science has experienced an impressive growth in the past few years. However, methods that apply ML to molecular property predictions [START_REF] O Anatole Von Lilienfeld | Exploring chemical compound space with quantum-based machine learning[END_REF][START_REF] Sanchez | Inverse molecular design using machine learning: Generative models for matter engineering[END_REF] or molecular dynamics [START_REF] Behler | Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces[END_REF][START_REF] Oliver T Unke | Machine learning force fields[END_REF] usually imply the availability of a certain amount of data previously produced by approximating the solution of the Schrödinger equation, as seen for instance in the two previous chapters. By considering the exponential numerical complexity involved in the exact solution of this equation it would be highly desirable to take advantage of ML techniques also in this context. For realistic molecules and materials this is a widely open field of research and the full potential of ML has yet to be widely exploited in its full potential.

G E N E R AT I V E M A C H I N E L E A R N I N G T O S O LV E T H E S C H R Ö D I N G E R E Q U AT I O N I N T H E C O N F I G U R AT I O N I N T E R A C T I O N S PA C E
Following earlier work on quantum neural networks [START_REF] Schuld | The quest for a quantum neural network[END_REF], in 2017 Carleo and Troyer showed that neural networks can effectively represent the quantum states of spin models, thus reducing the exponential complexity of the many-body problem [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF]. Specifically, they showed that a restricted Boltzmann machine (RBM) used as an ansatz within variational Monte Carlo can achieve variational energies lower than those obtained with traditional approaches. This wavefunction ansatz and similar variants based on different architectures are generally referred to as neuralnetwork quantum states (NQS). From a theoretical point of view this success relies on universal approximation theorems [START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF][START_REF] Le | Representational power of restricted Boltzmann machines and deep belief networks[END_REF], which imply that neural networks can approximate complex (but "reasonably smooth") high-dimensional functions, including quantum states. This method was demonstrated numerically considering the one-and two-dimensional Ising and Heisenberg models. More recently Carleo and coworkers have applied a similar approach to realistic Hamiltonians of small molecules [START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF]. This was achieved by mapping the electronic structure Hamiltonian into a spin-like Hamiltonian by using quantum information encodings. The most accurate results where obtained with the Jordan-Wigner mapping [START_REF] Jordan | Über das paulische äquivalenzverbot[END_REF], which leads to an approach equivalent to the configuration interaction (CI) [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]. Within the CI method the fully interacting wavefunction is expressed as a linear combination of excited Slater determinants (the "configurations"); while the coefficients of this expansion are typically computed as solution of an eigenvalue problem, in ref. [START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF] they were learnt by the RBM in an unsupervised way using a Monte Carlo sampling and the variational principle. While this approach was achieving chemical accuracy for small basis sets, the Monte Carlo sampling was repeatedly drawing the most dominant state (i.e. the Hartree-Fock determinant) and this was at the origin of a significant slow down of convergence for basis sets beyond the minimal STO-3G. It was recently shown that this issue can be alleviated by using autoregressive neural networks, that allowed calculations with up to 30 spin orbitals [START_REF] Thomas D Barrett | Autoregressive neural-network wavefunctions for ab initio quantum chemistry[END_REF]. Our numerical results presented below show that a simple RBM architecture used as a generative model can easily double this number, performing calculations with up to 56 spin-orbitals. More in general it should be noticed that the representation of the wavefunction in the configuration space is strongly non-smooth (determinants with similar occupations can provide significantly different contributions to the wavefunction) and this could be challenging for ML regression approaches. To overcome this issue alternative methods use instead the high expressive power of deep neural networks to represent electronic structure wavefunctions in real-space [START_REF] Han | Solving many-electron Schrödinger equation using deep neural networks[END_REF][START_REF] Hermann | Deep-neural-network solution of the electronic Schrödinger equation[END_REF][START_REF] Kessler | Artificial neural networks as trial wave functions for quantum monte carlo[END_REF][START_REF] Pfau | Ab initio solution of the many-electron Schrödinger equation with deep neural networks[END_REF]; as an advantage some of these techniques can achieve the complete basis set limit in a rather straightforward way but on the other side they require deep architectures involving the optimization of a large number of parameters and special care to keep into account the antisymmetry of the electronic wavefunction. In a recent work this type of approaches has also been extended to efficiently describe the wavefunction at different molecular geometries: Most of the weights in the neural network are shared among different different structures and only a small percentage of them is actually reoptimized [START_REF] Scherbela | Solving the electronic Schrödinger equation for multiple nuclear geometries with weight-sharing deep neural networks[END_REF].

The antisymmetry is instead naturally included in the CI space which, however, grows unfavorably with the system size. This has lead to the development of a series of methods that select the excited determinants that contribute the most to the wavefunction, either based on perturbation theory [START_REF] Huron | Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions[END_REF] or Monte Carlo sampling [START_REF] Greer | Monte Carlo configuration interaction[END_REF]. More recently machine learning techniques have been coupled to these methods [START_REF] Coe | Machine learning configuration interaction[END_REF][START_REF] Jeong | Active Learning Configuration Interaction for Excited-State Calculations of Polycyclic Aromatic Hydrocarbons[END_REF][START_REF] Flores | Chembot: A Machine Learning Approach to Selective Configuration Interaction[END_REF]. In this context Coe was the first to propose a machine learning approach that explores the configuration space but, instead of explicitly computing the coefficient of each configuration, this is predicted by a regression neural network [START_REF] Coe | Machine learning configuration interaction[END_REF]. While Carleo and coworkers apply a neural network approach to "exactly" fit the wavefunction in the configuration space using the variational principle [START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF], the alternative approaches of refs. [START_REF] Coe | Machine learning configuration interaction[END_REF][START_REF] Jeong | Active Learning Configuration Interaction for Excited-State Calculations of Polycyclic Aromatic Hydrocarbons[END_REF][START_REF] Flores | Chembot: A Machine Learning Approach to Selective Configuration Interaction[END_REF] are more qualitative and use supervised learning techniques. Restricted Boltzmann Machine (RBM) In this work we propose an alternative approach named CIgen that uses generative machine learning to directly generate the configurations that contribute the most to the electronic wavefunction. This is done with the general aim of avoiding the exploration of the overwhelmingly large configuration space to search for the most important determinants. In this respect our approach lies between the selective CI and NQS approaches: Similarly to selective CI, the determinants with the most important contributions are selected and used in an exact diagonalization procedure; this selection procedure is based on a generative NQS.

v i ∈ {0,1} h i ∈ {0,1} W ij E(v,
Notable examples of generative machine learning algorithms include the restricted Boltzmann machines [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF][START_REF] Roger G Melko | Restricted Boltzmann machines in quantum physics[END_REF], the variational autoencoders [START_REF] Diederik | Auto-encoding variational bayes[END_REF], and the generative adversarial networks (GANs) [START_REF] Goodfellow | Generative adversarial nets[END_REF]. CIgen is based on RBMs, a type of generative neural network whose architecture is shown in Figure 21(a). As for other generative models the RBM can learn in an unsupervised way the statistical distribution behind a series of input objects and then be used to generate new ones. A qualitative example is shown in Figure 21(b): An RBM can be trained with a series of pictures, e.g. of handwritten digits, in practice represented as a flattened matrix of pixels and afterward used to generate new realistic images of the same digit.

For a given molecule or material the excited Slater determinants are associated to an underlying probability distribution that could only be exactly evaluated by solving the Schrödinger equation to obtain the wavefunction. By representing the space of excited determinants simply as binary vectors (a sort of onedimensional binary image), the RBM is here iteratively trained using data from the current approximation of the wavefunction and subsequently used to generate new important determinants to improve this approximation (see Figure 21(c-d)).

Within the previous approach of Carleo and coworkers [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF][START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF] the RBM was not used as a generative model but rather as a regression model for the wavefunction (to this purpose the RBM architecture was generalized to include complex weights); here we use instead the RBM as a model to qualitatively represent the probability distribution associated with the wavefunction, which is afterward used to generate the most likely configurations. In this respect the CIgen approach has some similarities with the use of RBMs to learn the statistical distribution of data from experimental measurements of quantum states and to subsequently generate new configurations to compute quantum averages [START_REF] Roger G Melko | Restricted Boltzmann machines in quantum physics[END_REF][START_REF] Torlai | Neural-network quantum state tomography[END_REF]. Differently from selective CI or Monte Carlo CI our approach does not "explore" the huge determinant space to find significant contributions but rather directly generate them with the current approximation of the probability distribution associated with the RBM. Comparison with previous approaches will be further discussed below.

methodological approach

Configuration interaction approach

Considering the non-relativistic N electron problem in the Born-Oppenheimer approximation, the configuration interaction method solves the Schrödinger equation for a fixed basis set by expanding the wavefunction in the following way [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF]:

|Ψ⟩ = c 0 |Φ 0 ⟩ + ra c r a |Φ r a ⟩ + a<b r<s c rs ab |Φ rs ab ⟩ + a<b<c r<s<t c rst abc |Φ rst abc ⟩ + • • • (174)
where |Φ 0 ⟩ is the Hatree-Fock ground state, |Φ r a ⟩ is a singly excited Slater determinant from occupied spin orbital a to unoccupied spin orbital r, and all the other terms correspond to multiple excitations; while |Ψ⟩ could denote any excited state of a given system here the discussion will be focused on the ground-state. When all the possible excited determinants are included in eq. ( 174) the approach is called full configuration interaction (FCI) and the solution of the Schrödinger equation becomes exact for a given basis set. Within the "brute-force" FCI approach the coefficients of the expansion c 0 , c r a , • • • are determined by computing the expectation value of the Hamiltonian ⟨Ψ| Ĥ|Ψ⟩ and applying the variational principle. Since the total number of determinants grows as 2K N , where 2K is the number of spinorbitals, the FCI approach becomes quickly unpractical for most applications. It is however well known that often the wavefunction can be accurately represented by a limited number of excited determinants and here a ML model is proposed to directly generate the important contributions.

Generative model

In this work, it is shown that restricted Boltzmann machines (RBM) can be used to efficiently generate the excited determinants that contribute the most to the wavefunction in eq. ( 174). The RBM model was defined in Chapter 3 and we refer the reader to this section for more details. For the purpose of this work the input vectors v of the visible layer represent the Slater determinants and their associated probability P(v) should ideally be proportional to their contribution to the wavefunction (the square of the c coefficients in eq. ( 174)). The determinants are represented as binary vectors, where 1 denotes occupied and 0 unoccupied states, with size equal to the number of spin-orbitals (see Figure 21). The architecture of the RBM used in this work is rather simple, with 2K (the number of spin orbitals) neurons in the visible and in the hidden layers. Given a certain fixed subspace of determinants, the Hamiltonian is diagonalized to obtain the wavefunction coeffi- Left: Considering at first the spin down states (odd indexes), a random number in the interval between 0 and odd i m v i is drawn (depicted as a red dot) resulting in the population of the corresponding spin-orbital. Right: In a second iteration the occupation probability m v 3 of the third state is now set to zero to prevent a new sampling of this already occupied state and the procedure is repeated, resulting in the population of a new spin-orbital. The procedure is repeated until N/2 spin-orbitals for each spin are occupied. cients; the model is then trained with random configurations generated according to the corresponding probability distribution by updating the RBM parameters with the contrastive divergence algorithm [START_REF] Hinton | A practical guide to training restricted Boltzmann machines[END_REF]. The computational parameters are fixed to M = 50000 for the training set size, N b = 10 for the minibatch size, L = 20 for the number of Gibbs sampling iterations, and ϵ = 0.01 as learning rate (see Chapter 3). In Chapter 3 was given the Gibbs sampling algorithm as applied to the RBM. For our specific application we have a constraint for the values of v i , whose sum has to be equal to the number of electrons i v i = N. Since we exclusively consider closed-shell systems, we have the additional constraint that the number of spin up and down electrons must be equal. While we could generate vectors without restrictions and reject them a posteriori, it is significantly more efficient to intrinsically include this constraint in the procedure. This is achieved by replacing the step 3 of the Gibbs sampling with an approach based on the tower sampling algorithm [START_REF] Krauth | Statistical mechanics: algorithms and computations[END_REF]. Within this approach the N states are randomly populated with a probability proportional to m v i = p(v i = 1|h) (this quantity can indeed be interpreted as the occupation probability of each visible neuron). The basic idea of the application of the tower sampling algorithm to our problem is depicted in Figure 22. It is important to mention that in our approach we represent Slater determinants with 2K component binary vectors whose even and odd indexes label spin up and down, respectively. The Gibbs sampling approach is used both in the training of the model and in the generative phase. While L = 20 iterations are used in the training part, the generative process, which expands the configuration space, uses only one step L = 1. This (underconverged) choice introduces some "noise" in the model that allows to more easily sample determinants that are not already in the training set. As also shown by numerical experiments that we performed, a similar behavior can be obtained by lowering the number of iterations in the training or by increasing the temperature of the RBM during the generative procedure. However, in the former case it is harder to determine the reliability of the model and in the latter the generative procedure loses efficiency as larger values of L are used. The generation of new determinants and the retraining of the model are performed iteratively. Specifically, the CIgen algorithm performs the following steps:

1. Start: The configuration interaction singles and doubles (CISD) is used to generate an initial guess wavefunction.

2. Pruning: The determinants whose squared wavefunction coefficients (eq. ( 174)) are below a certain threshold p min are pruned.

3.

Training: The RBM model is trained only on the non pruned determinants but the Hartree-Fock (HF) determinant, which is typically associated to a very high probability, is not included in the training set (otherwise the ML would mainly generate the Hartree-Fock state in the following step). While discarded at the training stage, the HF determinant is included by default in the diagonalization step below.

Generation:

The trained RBM is used to generate a set of new determinants, whose number is proportional to the number of determinants already included in the current approximation of the wavefunction. During this procedure determinants are automatically discarded if they are already included in the wavefunction, if they do not couple through single or double substitutions with the current determinant set, and if they do not have correct spin and point group symmetries.

Diagonalization:

The Hamiltonian is diagonalized in the subspace that includes the newly generated determinants and new coefficients are obtained for eq. ( 174)

6. Iterate: The procedure is repeated from step 2 until convergence of the energy is achieved.

When sampling from the RBM model, the conditional probabilities p(v i = 1|h) and p(h i = 1|v) are sigmoid functions, as one can deduce from eqs. (154-155) For example, for the visible cells we have

p(v k |h) = exp βv k ( j W kj h j + a k ) v k exp βv k ( j W kj h j + a k ) = exp βv k ( j W kj h j + a k ) 1 + exp βv k ( j W kj h j + a k ) (175) and p(v i = 1|h) = 1 1 + exp(-β((Wh) k -a k )) . ( 176
)
One can easily see that rising the temperature will increase noise in the sampling and for T → ∞ all probabilities of occupation will be equal to one half. In this limit the specific parameters obtained training the model become irrelevant and the RBM generates new configurations completely randomly (with uniform probability). This trend as a function of the temperature can be seen in Figure 23, that shows the performance of CIgen to predict the total energy of H 2 O for three different temperatures. In our problem, once the energy is close to convergence and there are few determinants left to find, we will repeatedly sample the same determinants that are already in the wavefunction, which slows down the procedure. To reduce this computational burden and to obtain more out of equilibrium samples, we add a temperature increment in the generative procedure in the following form: at a given iteration, with N d the number of determinants after pruning, we try to generate f gen • N d new determinants. Let k count each attempt to add a new determinant. The temperature in the sampling is fixed to be T

= T + ⌊k/(50 • f gen • N d )⌋.
This number was fixed by numerical experiments. Differently from previous approaches based on RBMs for quantum states [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF][START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF], CIgen does not require an RBM architecture generalized to include complex weights and, while the model still learns in an unsupervised way, the procedure is not based on the variational principle. It can be further noted that if traditional NQS aim at fitting the wave function |Ψ⟩, our approach is instead a representation of the squared amplitude |Ψ| 2 . While the training of a ML model based on the variational principle is more elegant and physically motivated [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF][START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF], the CIgen approach requires a less strict training as it does not aim to exactly fit the wavefunction.

Direct generation of single and double excitations

The CIgen algorithm is based on a rather straightforward way of using a generative neural network model to enlarge the configuration subspace of interest. are connected to b. A uniform sampling will choose among the four with probability one fourth. If instead we first choose a or b with probability one half, before to randomly choose one of their excited determinants, d 1 and d 4 will be sampled with probability one sixth, and d 2 , d 3 with probability one third. After sampling this way, the resulting Hamiltonian to be diagonalized will be less sparse, resulting in more precise estimations of the FCI coefficients. Hence there is less chance to discard important determinants. To take advantage of this feature of the configuration space we developed a variant of our approach named CIgen-SD, using the following algorithm to generate one determinant:

1. A determinant v is randomly selected (with uniform probability) among those in the current wavefunction.

2. Gibbs sampling generates the probabilities of occupation of each visible neuron m v i = p(v i = 1|h). In this case, instead of directly applying the tower sampling algorithm, these probabilities are used to built a transition matrix whose elements are defined as

S ij = (1 -m v i ) • m v j
, where i is an occupied state and j an unoccupied state of v.

3.

It is randomly decided if a single or a double substitution to v is constructed (the two options having the same probability).

4.

The specific substitution is sampled using the tower sampling algorithm (Figure 22) but, in this case, the transition probability matrix elements S ij replace m v i .

The CIgen, CIgen-SD, and MCCI will be compared in the Numerical Results subsection.

Computational complexity

This section will discuss the numerical complexity of CIgen(-SD) and compare it with the Monte Carlo configuration interaction and the NQS method of Choo et al. [START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF]. We start the discussion by evaluating the cost of a single iteration. Let N det be the number of retained determinants (namely whose squared coefficients are larger than p min ) at a given iteration of the CIgen(-SD)/MCCI methods. Similarly, in the NQS method of Choo et al., N det corresponds to the number of unique determinants inside the sample of configurations. In each of those methods, the cost of a single iteration scales as O(N 2 det K 4 ), where the dependence on K 4 can be improved when the sparsity of the Hamiltonian matrix elements can be exploited. Indeed, selective CI methods such as CIgen(-SD) and MCCI involve a diagonalization procedure that seeks the lowest eigenstate; by using iterative algorithms this can be achieved by a repeated product of the Hamiltonian matrix and a guess vector, which costs O(N 2 det K 4 ) operations. Similarly, the optimization procedure in the NQS-based variational quantum Monte Carlo implies the computation of the expectation value of the Hamiltonian for a given sample of configurations. This is again a O(N 2 det K 4 ) procedure (more precisely, the calculation of the local energy of each sampled configuration requires the evaluation of expectation values with all the other configurations already sampled). It should be noticed that the sampling of new determinants scales better than the diagonalization/optimization step. All of the methods considered here propose new determinants at a cost that is O(K) or O(K 2 ). Since K ≪ N det in practice the most expensive task in this phase is to verify if the proposed determinants are already included in the current wavefunction set, that leads to a cost of O(N det log N det ). This can be explained considering that determinants can be stored as a sorted list of decimal integers corresponding to the binary occupation arrays. Using the binary search algorithm the identification of duplicates scales as O(log N det ) and this task has to be repeated for all the newly generated determinants, whose number is itself proportional to N det . Differently from the other methods, CIgen also requires a check to determine if the proposed determinants are among the single-double substituted determinants included in the current wavefunction. This task does not take advantage of the sorting of the list and, accordingly, requires a brute-force search through the list of the N det determinants in the wavefunction; repeating this task for each one of the newly proposed determinants yields a O(N 2 det ) complexity. For the sake of clarity the computational complexity involved in the sampling of new determinants are summarized in Table 10.

As shown by the numerical applications in Sec. 6.3 the MCCI sampling requires a higher number of iterations to converge than both CIgen and CIgen-SD. While we do not report NQS results in this paper in the original reference of Choo et al. it was shown that achieving chemical accuracy for the energy of H 2 O in the 6-31G basis set could be already problematic and a number of required iterations larger than CIgen(-SD) is expected (our results for this system are discussed in Sec.

6.3). These comparisons show that CIgen(-SD) is competitive with similar methods based on the configuration space sampling and, by decreasing the number of required iterations, avoids several computationally expensive diagonalization/optimization steps, which scale as O(N 2 det K 4 ) and dominate the computational cost per iteration.

A direct comparison of the numerical complexity of CIgen(-SD) with deep neural network methods such as FermiNet [START_REF] Pfau | Ab initio solution of the many-electron Schrödinger equation with deep neural networks[END_REF] or PauliNet [START_REF] Hermann | Deep-neural-network solution of the electronic Schrödinger equation[END_REF] is not straightforward. As a main difference those deep learning methods use a continuous space representation rather than the configuration space of CIgen. While the Monte Carlo wavefunction ansatzes of FermiNet and PauliNet still involve a linear combination of Slater determinants, their number is limited and the electronic correlation is to a large extent included by optimized orbitals and the Jastrow factor. These methodologies achieve the complete basis limit in a seamless way but require a number of parameters that is of magnitude larger than the RBM model used by CIgen.

Method Sampling complexity Diagonalization/optimization complexity

CIgen O(N 2 det ) O(N 2 det K 4 ) CIgen-SD O(N det log N det ) O(N 2 det K 4 ) MCCI O(N det log N det ) O(N 2 det K 4 ) NQS O(N det log N det ) O(N 2 det K 4 )
Table 10: Computational complexity of the different methods discussed in this work. This summary shows that the cost involved in the sampling of new determinants is smaller than the diagonalization/optimization part, which dominates the complexity of all these methods. Accordingly, the CIgen methods, by requiring less iterations to converge, can decrease the number of required diagonalizations and the overall computational cost (see text).

Computational details

The equilibrium geometries of the molecules considered in this work are optimized at the CCSD(T) level of theory in the corresponding basis set (6-31G or cc-pVDZ) and can be obtained from the Computational Chemistry Comparison and Benchmark DataBase [START_REF]NIST Computational Chemistry Comparison and Benchmark Database[END_REF]. The reference results for these molecules at the CCSD(T) and FCI level of theory (Table 11 and Figs. [START_REF] Carleo | Machine learning and the physical sciences[END_REF][START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF] have been obtained from the Molpro code [START_REF] Knowles | MOL-PRO, a package of ab initio programs[END_REF][START_REF] Peter | A determinant based full configuration interaction program[END_REF].

The electronic structure part of our implementation (specifically the diagonalization of the Hamiltonian in the space of the generated determinants) is based on the Quantum Package code [START_REF] Garniron | Quantum Package 2.0: An Open-Source Determinant-Driven Suite of Programs[END_REF]. For the CIgen, CIgen-SD, and MCCI methods the selection of the determinants to be included in the wavefunction in eq. ( 174) is based on an acceptance threshold of 10 -12 on the squared coefficients.

numerical results

We now discuss the efficiency and accuracy of the CIgen approach by considering applications to molecular systems. To this purpose total energy calculations for the C 2 , N 2 , and H 2 O molecules are performed for the 6-31G and cc-pVDZ basis sets considering the active spaces indicated in the third column of Table 11. For the largest system the number of spin orbitals is about double than what previously achieved with neural network quantum states [START_REF] Thomas D Barrett | Autoregressive neural-network wavefunctions for ab initio quantum chemistry[END_REF].

As a preliminary step it is important to discuss the effect of the temperature on the generative power of the RBM (see eq. ( 155)). Figure 23 shows the convergence of the total energy of H 2 O in the cc-pVDZ for three values of the temperature. The convergence of the CIgen approach is optimal at T=1 (chemical accuracy is achieved within 11 iterations) but sizeably slows down when T is increased to 5. In the limit of a very large temperature (T=500) the method performs poorly. Indeed, as explained earlier, the T → ∞ limit leads to a completely random generation of determinants with uniform probability. This shows that for reasonable choices of the temperature the CIgen method actually learns the probability associated with the wavefunction and significantly speeds up the generation of the excited determinants with respect to random sampling.

We now consider the convergence of the total energy in the cc-pVDZ basis for the three molecules C 2 , N 2 , and H 2 O in the cc-pVDZ basis and compare the perfomance of the CIgen and CIgen-SD algorithms with the MCCI algorithm. In all three cases the two CIgen variants considered here outperform the MCCI method. This is particularly true for the CIgen-SD approach that achieves chemical accuracy with a significantly smaller number of iterations. The detailed results for the total energies of the molecules considered here are presented in Table 11, where CIgen-SD values are compared to coupled-cluster with singles, doubles and perturbative triples CCSD(T) and to FCI reference values. By considering a limited number of iterations (from a minimum of 9 for C 2 in 6-31G to a maximum of 17 for water in cc-pVDZ), CIgen-SD converges to FCI values within chemical accuracy. This is achieved by generating an excited determinant subspace that is significantly smaller than the full space. For example, for the molecule N 2 in the cc-pVDZ basis set only 18 millions determinants are generated and included in the wavefunction out of the about 540 millions that would be allowed by spin and space symmetries. This significant reduction in the number of determinants is not achieved by a search in the configuration space but rather by an iterative generation (and model training) of the determinants that contribute the most to the total energy. If for certain applications a lower level of accuracy can be tolerated the number of selected determinants can be further decreased by increasing the value of p min . Using a threshold of 10 -8 instead of 10 As a final result in Figure 25 we present a full binding curve for the N 2 molecule. This represent a much more challenging example since for large interatomic distances the binding is characterized by strong static correlation whose description could become problematic for several traditional quantum chemical approaches. This is the case of CCSD(T) that starting at around 1.8 Å produces a binding curve 0.8 1.0 1.2 1. [START_REF] Hans | Molecular dynamics simulations at constant pressure and/or temperature[END_REF] with an unphysical behavior. The CIgen approach well describes the N 2 binding curve at every distance and produces a curve in excellent agreement with FCI. The comparison with the CISD curve, which provides the starting training data for the generative model, shows clearly that CIgen recovers an increasing amount of electronic correlation energy as a function of the interatomic distance. The increase in the electronic correlation manifests itself in the number of determinants that have to be included in the wavefunction to achieve convergence: At equilibrium distance the CIgen method selects a few hundreds of thousands of determinants while closer to dissociation (2.4 Å), around three millions determinants are necessary. While the number of the required determinants changes at different distances it still remains significantly smaller the the 127 millions of determinants in the full CI solution.

conclusions

In summary we have shown how a generative model can be used to solve the Schrödinger equation by sampling the excited Slater determinants that contribute the most to the wavefunction. Numerical applications show that this approach is already competitive with previous approaches based on Monte Carlo sampling of the excited determinants [START_REF] Greer | Monte Carlo configuration interaction[END_REF] or on machine learning ansatzes to represent the wavefunction in the configuration space [START_REF] Thomas D Barrett | Autoregressive neural-network wavefunctions for ab initio quantum chemistry[END_REF][START_REF] Choo | Fermionic neuralnetwork states for ab-initio electronic structure[END_REF]. An improvement that should be addressed in future work involves the development of a generative model that intrinsically takes into account the symmetry of the determinants and, more in general, the properties of the determinants contributing to the wavefunction. It has already been shown that a sizeable improvement in the convergence speed can be achieved by directly generating the determinants that couple with double and single substitutions to the determinants already included in the wavefunction.

In the current implementation, however, spin and spatial symmetries of the generated determinants are verified and enforced only a posteriori. The development of a generative machine learning model that keeps into account these symmetries would certainly increase the numerical efficiency of CIgen to possibly perform calculations with significantly larger numbers of spin orbitals. This will be considered in future work by developing new symmetry compliant architectures or by directly excluding symmetry forbidden substitutions in the transition probability matrix. This thesis has brought contributions to the computational physics emerging field of machine learning (ML) aided ab initio calculations. More specifically, it has been seen how machine learning can significantly lower the computational cost of correlated quantum methods to ease the reach of chemical accuracy.

In Chapter 4, an assessment of the accuracy and of the limitations of Machine Learning Perturbation Theory (MLPT) was provided. MLPT uses a ML model to learn the energy difference between two quantum mechanical approximations, and use that model to perform thermodynamic perturbation theory (TPT) to a fraction of its cost, allowing to obtain finite temperature estimates in theories that would be inaccessible from a straighforward molecular dynamics. Using different exchange-correlation functionals to study the enthalpy of adsorption of small molecules inside zeolites, it was shown how TPT can be sensible to a potential mismatch between reference and target configurational spaces. This was made apparent from a quantitative deviation of the MLPT estimates with respect to the reference values of total electronic energies and enthalpies of adsoprtion, and from a clustering of the non overlapping configurational spaces using dimensionality reduction techniques. A diagnostic test, in the form of a number named I w , was proposed to detect those biased estimates. This number measures the fraction of configurations necessary to recover half of the sorted statistical weights after the TPT reweighting. A half unit for I w corresponds to a uniform shift of all energies, while a I w close to zero indicates that a negligible fraction of the configurations holds most of the statistical weights, resulting in a biased ensemble average. A solution to this problem was devised in the form of a Monte Carlo resampling of the statistics, in combination with the same ML model (MLMC). The MLMC method is a straightforward Metropolis algorithm, where the energy evaluation is performed in two steps: the energy of the configuration is first computed using the reference theory, and the previous MLPT ∆-ML model is used in a second step to infer the target value of the energy. The correct configurational space can be recovered in this way, at the price of a new ab initio simulation, but without any new target level calculation. MLMC can also be used to further confirm some MLPT estimate with a good a priori overlap between reference and target levels of theory, as was done to study previous results in the Random Phase Approximation.

Those results were used in Chapter 5, to provide an estimate of the ambient temperature enthalpy of adsorption of carbon dioxide in a protonated chabazite with Coupled Cluster Single Double and perturbative Triple (CCSD(T)). Using a periodic implementation of the Coupled Cluster method, in a localized basis set formed by the MP2 natural orbitals, one hundred and ten CCSD(T) calculations were performed for each (molecule, zeolite and adsorbed zeolite) system. The further MLPT estimate was found to be in excellent agreement with the experimental reference, with an error below 0.1 kcal/mol. Because of a rather small value of I w in this work, even if reliable from the results of Chapter 4, a MLMC estimate was computed, confirming the MLPT result, and providing what is to our knowledge the first finite temperature CCSD(T) estimate for a periodic material.

There is however still many roads to further developp and assess MLPT and MLMC. In the work presented here, the training set was chosen from evenly separated points in the reference molecular dynamics. A more informed training set choice could allow to reduce even more the size of the training set or allow for a better learning. Other level of reference theories could be considered, such as semi empirical tight binding models, for example associated to density functional theory (DFT) to provide fast estimates with DFT accuracy. More complex systems could be studied, chemical reactions for instance, to have more insight into the capabilities of MLPT. A perfectionned MLMC algorithm, using, e.g., parallel tempering, could be equally profitable, in order to reduce the autocorrelation length of the samples, which induce a higher error bar compared to the MD.

Finally, in Chapter 6, it was shown, by computing total energies and dissociation curve of N 2 ,C 2 and H 2 O in the 6-31G and cc-PVDZ basis sets, that a generative ML model can be used to efficiently generate new Slater determinants in an iterative selective Configuration Interaction (CI) algorithm. The ML model is a restricted Boltzmann machine (RBM): a two-layers energy based generative neural network, that models a probability distribution, from which one can sample. Contrary to neural quantum states where a ML model is used to model the wave function, this is the wave function squared amplitude that is fitted here. From a starting wave function, the RBM is trained to learn the Born probabilities associated to each Slater determinant in the wave function. In a second step new determinants are sampled from the RBM, providing a bigger Hilbert space to be diagonalized for the ground state. The new ground state coefficients provide new probabilites to be learned, and this is repeated iteratively until convergence. Two versions of this method were proposed: the first one sample new determinants in a direct manner, by populating orbitals according to their RBM occupation probability. The second version exploit the properties of the connections among determinants in the CI space, by directly proposing single and double substitutions to the current wave function. This is achieved by constructing a transition probability matrix from the RBM occupation probabilities. This methodology, coined CIgen, was seen to outperform the uniform sampling of the determinants, and to be already competitive with previous Monte Carlo approaches or ML quantum Monte Carlo Configuration Interaction ansatzes. By reducing the number of iterations necessary to converge, several of the most costly iteration (the obtention of the ground state) are avoided. One future improvement could be to directly learn the symmetries of the system, instead of enforcing them a posteriori. Other ML models on which exact sampling is possible could be tested, as well as other physical systems. Encouraging results were found using the CIgen code to solve spin systems such as the antiferromagnetic J1-J2 model in one and two dimensions. From a more general point of view, the CIgen model is essentially about replacing a Monte Carlo proposal distribution by a more informed one, in a Bayesian inference framework. This idea has a lot of potential to bring more efficient exploration of the configurational space, allowing in the case of CIgen for a faster convergence of the CI energy minimization. It can also be used in a more traditional Monte Carlo setup to reduce the autocorrelation length of the sample, as was done in the context of lattice field theory with normalizing flow models [START_REF] Kanwar | Machine Learning and Variational Algorithms for Lattice Field Theory[END_REF] for example. This could be brought to continuous-time quantum monte carlo as well, or in finite temperature simulation, such as the MLMC algorithm, by learning the configurational space produced by the reference dynamics, and propose Monte Carlo updates in an active or reinforcement learning framework. R É S U M É La précision chimique, définie comme une kilocalorie par mole, est un objectif de résolution énergétique pour les méthodes computationelles en science des matériaux. Alors que les approximations les plus simples en chimie quantique (Hartree-Fock) et en physique du solide (théories fonctionnelles de la densité locales et semi-locales (DFT)) sont numériquement abordables pour des systèmes modérément grands (jusqu'à quelques milliers d'atomes, avec un coût numérique typiquement en O(K 4 ) et O(K 3 ) respectivement avec la taille de la base), elles ne peuvent pas être améliorés de manière systématique et ne parviennent pas à la précision chimique dans beaucoup d'applications. Il faut alors recourir à des approximations plus précisesméthodes post-Hartree-Fock telles que Configuration Interaction, la théorie des perturbations Mollër-Plesset, Coupled Cluster, ou des modèles au-delà de la DFT tels que la Random Phase Approximation (RPA) mais ces dernières sont numériquement plus coûteuses et donc limitées à de plus petits systèmes. En particulier, alors que des implémentations périodiques de Mollër-Plesset ou Coupled Cluster ont été récemment proposées, les résultats à température finie (qui nécessitent quelques centaines de milliers de configurations calculées) sont toujours inatteignables pour les systèmes de matière condensée utilisant ces théories. Dans ce contexte, les méthodes d'apprentissage automatique ont connu au cours de la dernière décennie un important développement afin de réduire le coût numérique des méthodes de mécanique quantique. Parmi ces méthodes de machine learning (ML), nous nous intéresseront tout particulièrement dans cette thèse à la Machine Learning Perturbation Theory (MLPT). Cet algorithme repose sur deux ingrédients : le ∆-Machine Learning et la théorie de la perturbation thermodynamique. L'idée du ∆-Machine Learning est la suivante : plutôt que d'apprendre directement l'énergie (ou une autre propriété) en fonction de la configuration atomique, il est plus facile d'apprendre la différence d'énergie telle que donnée par deux approximations quantiques différentes. En considérant deux fonctionnelles d'échange et corrélation DFT différentes par exemple, les deux approximations incluent le même terme cinétique et diffèrent dans leur traitement de l'énergie de corrélation. Apprendre la différence entre les deux approximations est plus facile que d'apprendre directement toute la description physique contenue dans l'une des deux. Ainsi, si plusieurs milliers de configurations sont nécessaires pour apprendre la surface d'énergie potentielle d'un système complexe, ce nombre peut être réduit à quelques dizaines ou centaines grâce à cette technique. Dans la théorie des perturbations thermodynamiques maintenant, on peut exprimer la moyenne d'ensemble canonique d'un hamiltonien en fonction de moyennes d'un autre hamiltonien, si la différence entre les deux hamiltoniens est connue. Plus explicitement, étant donnée une dynamique moléculaire réalisée avec une approximation X, on peut écrire la moyenne telle qu'elle aurait été obtenue avec une autre approximation Y si, pour chaque configuration de la dynamique moléculaire, la différence d'énergie entre les approximations X et Y est connue. La technique MLPT, introduite en 2019 par Chehaibou et al., utilise la théorie de la perturbation thermodynamique pour repondérer une dynamique moléculaire ab initio, réalisée avec une référence DFT semi locale, pour obtenir des moyennes d'ensemble RPA, et ainsi pouvoir calculer l'enthalpie d'adsorption du dioxyde de carbone et du méthane dans une zéolithe. Compte tenu du coût de calcul RPA, quelques dizaines de configurations ont été extraites des dynamiques moléculaire pour entraîner un modèle de ∆-Machine Learning et ainsi apprendre la différence d'énergie entre l'énergie DFT et RPA. Un travail utlérieur similaire a permis de calculer des énergies libres d'activation pour une réaction d'échange de proton dans une chabazite protonée, également au niveau RPA. Le travail présenté dans cette thèse se propose de continuer le développement de cette méthode MLPT. D'une part, la théorie de la perturbation thermodynamique repose sur l'égalité des espaces configurationnels accessibles aux théories de référence X et cible Y. En pratique, les dynamiques moléculaires numériques produisent des échantillons discrets, proches des positions d'équilibre associées au hamiltonien utilisé. Si ces configurations d'équilibre sont trop distantes entre les théories de référence et cible, l'estimation issue de la repondération risque d'être biaisée. Ce phénomène est étudié dans cette thèse : un test de diagnostic est proposé pour détecter ces cas problématiques, ainsi qu'une solution sous la forme d'un ré-échantillonnage Monte Carlo. Ce dernier utilise l'algorithme de Metropolis, mais le modèle de ∆-Machine Learning est réutilisé ici pour évaluer l'énergie au niveau cible en la calculant d'abord au niveau de référence. Ainsi on s'assure de procéder à un échantillonage de l'espace configurationnel de la théorie cible et l'on peut retrouver des estimations non biaisées. Pour procéder à cette analyse, cinq fonctionnelles DFT différentes sont utilisées, pour lesquelles on dispose à chaque fois d'une dynamique de référence. La méthode MLPT est appliquée entre chaque couple parmi les cinq approximations. Le test de diagnostic permet bien de détecter les cas pathologiques, et le ré-échantillonage Monte Carlo redonne des valeurs correctes dans ces cas là, à la précision chimique. Ce ré-échantillonage est également réalisé pour confirmer de précédents calculs RPA, et le test de diagnostic semble indiquer un recouvrement configurationnel satisfaisant entre la fonctionnelle semi locale de Perdew-Burke-Ernzerhof plus une correction semi-empirique pour les interactions van der Waals (PBE+D2), et la RPA. Dans un second temps, on utilisera cette méthodologie pour calculer l'enthalpie d'adsorption du dioxyde de carbone dans une zeolite, au niveau de la théorie du cluster couplé, avec excitations simples, doubles, et triples perturbatives (CCSD(T)).

Cette dernière est généralement considérée comme la méthode de référence en chimie quantique, permettant d'obtenir la précision chimique dans la plupart des applications. Son coût numérique est cependant bien trop élevé pour des applications à température finie en matière condensée, et l'on présente ici la première estimation CCSD(T) d'une propriété à température finie dans un matériau périodique.

Les travaux mentionnés jusqu'à présent nécessitent dans le cadre de l'entrainement d'un modèle d'apprentissage automatique qu'un certain nombre de calculs coûteux aient été réalisés au préalable. Il pourrait être avantageux d'utiliser le machine learning directement durant le processus de résolution de l'équation de Schrödinger. En 2017, Carleo et Troyer ont introduit l'idée des états quantiques neuronaux (NQS). Dans le cadre du Monte Carlo variationnel quantique, on utilise des ansatzes de fonction d'onde qui sont optimisées variationnellement. Dans ce contexte, Carleo et Troyer ont suggéré d'utiliser plutôt le pouvoir représentatif des réseaux de neurones pour paramétrer la fonction d'onde à optimiser. Plus précisément, ils choisissent une machine de Boltzmann restreinte (RBM) et optimisent ses paramètres au moyen d'un algorithme provenant des méthodes de Monte Carlo variationnel quantique. Cet ansatz neuronnal est appliqué à des modèles de spins mais a ensuite été étendu aux problèmes de structure électronique par Choo et al. : en reformulant l'hamiltonien de structure électronique en un hamiltonien de spin (en utilisant la transformation de Jordan-Wigner par exemple), on peut retrouver la théorie des interactions de configuration (CI). Ils ont effectué des calculs d'énergie électronique sur de petites molécules (jusqu'à quatorze électrons) dans les bases STO-3G et 6-31G. Bien que la convergence des calculs à la précision chimique était observée dans la base STO-3G, ça n'était pas le cas pour la base 6-31G. Il apparaît encore difficile d'élaborer des algorithmes efficaces pour optimiser un réseau de neurones -de surcroît avec des poids complexes -sur le problème CI. L'échantillonnage est compliqué par un poids prédominant des configurations les plus importantes dans la fonction d'onde, entraînant un ralentissement problématique de l'algorithme pour l'ensemble de base non minimal 6-31G. Une autre difficulté pour le RBM ou d'autres modèles simples est la nature très complexe des corrélations entre l'occupation des orbitales dans les configurations. Par modèle simple, on entend un modèle avec peu de paramètres, e.g. un nombre similaire de variables cachées et d'entrée, car le pouvoir représentatif des réseaux de neurones augmente avec le nombre de ces variables cachées. Une solution consiste à utiliser des réseaux de neurones profonds, par exemple en empilant des RBM, une couche cachée étant une entrée pour la couche cachée suivante (machine de Boltzmann restreinte profonde). Des fonctions d'ondes spatiales ont été codées en tant que NQS dans de telles architectures d'apprentissage profond au cours des quatre dernières années, surpassant l'état de l'art précédent des ansatzes varia-tionnels QMC, avec des systèmes jusqu'à trente électrons. Les fonctions d'onde CI tronquées restent cependant une méthode précieuse pour obtenir une précision systématique et calibrer d'autres méthodes. Étant donné que l'espace CI croît de manière combinatoire avec la taille du système, les algorithmes CI antérieurs sélectionnent généralement les configurations les plus importantes en utilisant la théorie des perturbations ou un échantillonnage aléatoire, avant une diagonalisation du sous-espace de Hilbert, dans un processus itératif. L'apprentissage automatique a été proposé comme moyen de remplacer la théorie des perturbations comme critère de sélection des configurations importantes par Coe et al. en 2018, en utilisant un réseau de neurones régressif pour apprendre les configurations importantes. Dans cette thèse, on propose une voie hybride. On se place dans le même type d'algorithmes itératifs, avec un échantillonage aléatoire. Mais on remplace la distribution de probabilité uniforme pour la proposition de nouveaux déterminants par une distribution contenant de l'information quant aux déterminants actuellement dans la fonction d'onde. Ceci est réalisé en entrainant une machine de Boltzmann restreinte à reproduire les probabilités de Born associées aux coefficients des determinants de la fonctions d'onde CI. On échantillonne ensuite de nouveaux déterminants à partir de ce modèle, avant diagonalisation. Contrairement aux états quantiques neuronaux, le réseau de neurone ne sert pas d'ansatz à la fonction d'onde, mais au carré de son amplitude, qui est donc une amplitude de probabilité. Le modèle RBM permettant l'échantillonage, on peut donc proposer de nouvelles configurations plus probables qu'avec un échantillonnage non informé. Ce modèle est testé sur de petites molécules (jusqu'à quatorze éléctrons) dans les bases 6-31G et cc-PVDZ, et l'on montre que la convergence à la précision chimique est plus rapide qu'avec un Monte Carlo traditionnel.

Cette thèse est donc une contribution au développement de nouvelles méthodologies d'apprentissage automatique afin de réduire le coût de calcul nécessaire à l'obtention de propriétés chimiquement précises. Après une introduction aux problématiques considérées ici dans le premier chapitre, le chapitre deux procède à une revue du problème à N éléctrons et des méthodes considérées dans cette thèse, avec une première approximation appelée Hartree-Fock. Les méthodes dites post-Hartree Fock sont étudiées, notamment les interactions de configuration, la théorie du cluster couplé, et la perturbation Mollër-Plesset. La théorie de la fonctionnelle de la densité est alors introduite, ainsi que la RPA. Dans le chapitre trois on s'intéressera aux différens concepts et méthodes d'apprentissage automatique utilisées dans le reste de ce manuscript. Les chapitres quatre, cinq et six constituent les contributions de cette thèse. Le chapitre quatre est une étude de la méthode MLPT et du problème de recouvrement configurationnel mentionné plus haut. On y présente le ré-échantillonage Monte Carlo comme solution. Le chapitre cinq applique ces méthodes (MLPT et ré-échantillonage Monte Carlo) au calcul d'une estimation de l'enthalpie d'adsorption du dioxyde de carbone dans une zeolite au niveau CCSD(T). Le chapitre six est consacré à un modèle génératif pour proposer de nouvelles configurations dans le cadre de la méthode CI. Enfin le septième et dernier chapitre discute les résultats de cette thèse et conclut.

  Figure 20First (a) and second (b) series of peaks of the partial radial distribution function for the Si-O pairs determined at different levels of theory. . . . . . . . . . . . . . . . . . . . . . . Figure21Restricted Boltzmann machine (RBM) as a generative model and its application to the solution of the Schrödinger equation. Figure22Pictorial representation of the tower sampling algorithm. . . Figure23Effect of the temperature in the generative procedure. . . . . Figure24Total energy convergence as a function of the number of iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure25CIgen dissociation curve of N 2 . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Jacob's ladder of density functional theory, proposed by Perdew andSchmidt. The ladder allows to climb from the mean field description of the Hartree(-Fock) product to chemical accuracy by adding more and more non-locality in the description of the exchange-correlation energy[START_REF] Perdew | Jacob's ladder of density functional approximations for the exchange-correlation energy[END_REF].

Figure 3 :

 3 Figure 3: Total electronic energies, centered around their mean, of forty protonated chabazite configurations, within the PBE DFT functional and the random phase approximation. The black curve shows their difference.

Figure 4 :

 4 Figure 4: Restricted Boltzmann Machine: this neural network model is composed of one input (or visible) layer (the v i ), and one unique hidden layer (h i ).The model parameter is made of the weight W ij between both layers, plus two additional bias vectors composed of a i in the visible layer and b i in the hidden one.

Figure 5 :

 5 Figure 5: Self Consistent construction of the Hartree-Fock wave function: Af-ter construction and diagonalization of the Fock matrix, starting from a guess of the orbitals (possibly random), a new set of orbitals ϕ i and eigenvalues ϵ i are founds. This is repeated until some convergence criterion is met on the later.

Figure 6 :

 6 Figure 6: Dissociation curve of the BH molecule using various wave function methods presented in this chapter. Calculations were realized with Nwchem [5] in the 6-31G basis set.

Figure 7 :

 7 Figure 7: Self consistent field procedure to find the ground state eigenvectors of the Kohn-Sham Hamiltonian associated to a given external potential V ext in a given one-electron basis set {ϕ i }.

Figure 8 :

 8 Figure 8: Illustration of the pseudo-potential method: the real diverging Coulomb potential (dashed blue line, bottom), is replaced by a an effective potential (red line, bottom) inside some cutoff radius around the nuclei. The highly oscillatory behavior of the wave function (dashed blue line, above) is thus suppressed in this region (red line, above), and a smaller number of plane waves are necessary.

Figure 9 :

 9 Figure 9: Dissociation curve of the BH molecule using various wave functions methods, the LDA and PBE DFT functionals, and the RPA method.Calculations were realized with Nwchem[START_REF] Aprà | NWChem: Past, present, and future[END_REF] in the 6-31G vasis set. The RPA correlation energy was computed using PBE orbitals

2. 4 . 1

 41 Ab initio Molecular Dynamics 2.4.1.1 The adiabatic approximation Say {ψ i (r|R)} are the eigenstates of the electronic Hamiltonian Ĥe with eigenvalues {ϵ i }, where r (R) is a shorthand for all the electrons (nuclei) positions. Let us call Potential Energy Surfaces (PES) the energy surfaces ϵ i (R). It is possible to decompose a wave function belonging to the full Hilbert space of Ĥ as the following product of wave functions belonging to the Hilbert spaces of the electronic and nuclear Hamiltonian 12 :

Figure 10 :

 10 Figure 10: Pictorial representation of thermodynamic perturbation theory: the black curves are the potential energies of two harmonic Hamiltonian with different equilibrium positions (black dots). The bell curves are the associated density probabilities distributions along the configurational space.

Figure 11 :

 11 Figure 11: Smooth Overlap of Atomic Positions descriptor: two local atomic environments are described as gaussian centered densities inside some cutoff. The similarity between both environments is assessed by the amplitude of the overlap of those densities, upon rotation to preserve rotational invariance of the similarity measure.

  3 to 10 4[START_REF] Bartok | Machine Learning Unifies the Modelling of Materials and Molecules[END_REF][START_REF] Hansen | Machine Learning Predictions of Molecular Properties: Accurate Many-Body Potentials and Nonlocality in Chemical Space[END_REF][START_REF] Huo | Unified Representation of Molecules and Crystals for Machine Learning[END_REF].

Figure 12 :

 12 Figure 12: Kernel Ridge Regression of the sinus function, using a gaussian kernel.

  are computed and used to evaluate the averages ⟨• • • ⟩ data ; (3) N b pairs {v (L) , h (L) } are generated by a L step Gibbs sampling and used to evaluate the averages ⟨. . .⟩ model ; (4) the parameters of the RBM model are updated with a learning rate ϵ:

Figure 13 :

 13 Figure 13: Adsorbed systems studied. (a) CH 4 in protonated Chabazite (HChab) and (b) CO 2 in siliceous chabazite (SiChab).

Figure 14 :

 14 Figure 14: t-SNE representation of the configurational spaces spanned by the different functionals. This visualization is based on 500 selected configurations from each molecular dynamics trajectory. principal component of the SOAP descriptor for each 500 configurations of each functionals.

Figure 15 :Figure 16 :Figure 17 :

 151617 Figure 15: Distributions of deviations of SCAN energy from the average value of the corresponding reference calculation (µ(E Reference )) obtained in three different simulations: straightforward MD with the SCAN functional (blue), ML (green) and MLMC (red) with PBE production method.

Figure 18 :

 18 Figure 18: The unit cell of the system studied in this work, CO 2 in protonated chabazite.

Figure 19 :

 19 Figure 19: t-SNE representation of the configurational spaces spanned by the PBE+D2 molecular dynamics (MD) trajectory and the CCSD(T) machine learning Monte Carlo (MLMC) trajectory. The configurations included in the training set are also shown to demonstrate that they cover essentially whole relevant part of the configurational space sampled at the CCSD(T) target level. The axes represent the two components of the t-SNE projection.

Figure 20 :

 20 Figure 20: First (a) and second (b) series of peaks of the partial radial distribution function for the Si-O pairs determined at different levels of theory.

[ 83 ]

 83 "Solving the Schrödinger Equation in the Configuration Space with Generative Machine Learning", by Basile Herzog, Bastien Casier, Sébastien Lebègue, and Dario Rocca, published in the Journal of Chemical Theory and Computation 2023 19 (9),

2 = 111

 2111 among 2K spinorbitals (SO): factorial number of Slater determinants v i Probability: P(v i )=|c i |

Figure 21 :

 21 Figure 21: Restricted Boltzmann machine (RBM) as a generative model and its application to the solution of the Schrödinger equation. (a) RBM architecture consisting of one visible input layer and one hidden layer of binary values (differently from standard neural networks a proper output layer is not present); for a given configuration (v, h) the parameters (a, b, W) are used to define an energy function E and an associated Boltzmannlike probability density P. (b) As an example, the RBM can be trained on a set of handwritten digits and afterward used to generate new realistic ones; to this purpose the digit's images are flattened to become unidimensional binary vectors v (k) where 1 and 0 correspond to the digit and background pixels, respectively. (c) The configuration interaction (CI) approach expands the wavefunction of a molecule as a a linear combination of excited Slater determinants, which can be represented as a sort of unidimensional binary image. (d) The CIgen algorithm presented in this work iteratively trains an RBM on the distribution of determinants in the current approximation of the wavefunction and subsequently uses it to expand it by generating new important contributions.

5 Figure 22 :

 522 Figure 22: Pictorial representation of the tower sampling algorithm, in which the probability of sampling the i-th orbital is proportional to the height m v i of each rectangle.Left: Considering at first the spin down states (odd indexes), a random number in the interval between 0 and odd i m v i is drawn (depicted as a red dot) resulting in the population of the corresponding spin-orbital. Right: In a second iteration the occupation probability m v 3 of the third state is now set to zero to prevent a new sampling of this already occupied state and the procedure is repeated, resulting in the population of a new spin-orbital. The procedure is repeated until N/2 spin-orbitals for each spin are occupied.

Figure 24 :

 24 Figure 24: Total energy convergence as a function of the number of iterations: The three molecules N 2 , C 2 , and H 2 O in the cc-pVDZ basis set are considered and the plot compares the CIgen, CIgen-SD and MCCI algorithms. The full configuration interaction (FCI) reference values and the chemical accuracy threshold are represented by horizontal lines.

Figure 25 : 2 :

 252 Figure 25: CIgen dissociation curve of N 2 : The results are obtained using the 6-31G basis set and the CCSD(T), FCI, and CISD curves are provided for comparison purposes.

note

  The programs needed to perform CIgen and CIgen-SD calculations are available at https://github.com/bslhrzg/cigen. C O N C L U S I O N A N D P E R S P E C T I V E

  2 . . . . . . . . . . . . . . . . . .Table Experimental and reference theoretical enthalpies of adsorption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Table Enthalpies of adsorption estimates of CH 4 @HChab and CO 2 @SiChab 66 Table Deviations of the MLPT estimates . . . . . . . . . . . . . . . 68 Table Values of the I w index corresponding to each individual MLPT estimate of the ensemble internal energy. . . . . . . . 73 Table Deviations of the SCAN and SCAN+rVV10 internal energies predicted using MLPT and MLMC. . . . . . . . . . . . . 80 Table Experimental enthalpies of adsorption of CH 4 @HChab and CO 2 @SiChab and their calculated reference MLPT, and MLMC values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 Table Values of the I Chemical accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Prior work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

	L I S T O F TA B L E S A C R O N Y M S
	I N T R O D U C T I O N
	ML Machine Learning
	MLPT Machine Learning Perturbation Theory
	1.1
	MLMC Machine Learning Monte Carlo
	TPT Thermodynamic Perturbation Theory
	ix

w index corresponding to the MLPT predictions of RPA energies from PBE+D2 trajectories. . . . . . . . 83 Table Correlation energy contributions to the adsorption energy (kcal/mol) computed using different levels of theory and basis sets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Table Enthalpy of adsorption of CO 2 in protonated chabazite (kcal/mol) computed using different target and sampling methods. 91 Table Computational complexity of the different methods discussed in this work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 Table Total energies (Ha) for the different systems/basis sets . . . 110 DFT Density Functional Theory LDA Local Density Approximation GGA Generalized Gradient Approximation PBE Perdew-Burke-Ernzerhof HF Hartree-Fock CI Configuration Interaction MP Møller-Plesset CC Coupled Cluster RPA Random Phase Approximation MD Molecular Dynamics RBM Restricted Boltzmann Machine KRR Kernel Ridge Regression x 1

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Wave function-based methods . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 The many electron problem . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2 Hartee-Fock and the correlation energy . . . . . . . . . . . . . . . .

2.1 14 2.2.3 Configuration Interaction . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.4 Perturbation theory expansion of the correlation energy . . . . . . 19 2.2.5 Møller-Plesset Perturbation Theory . . . . . . . . . . . . . . . . . . 20 2.2.6 Coupled Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.7 Basis sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.8 Methods comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3 Density functional theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.1 Hohenberg-Kohn theorems . . . . . . . . . . . . . . . . . . . . . . . 27 2.3.2 Kohn-Sham Density functional theory . . . . . . . . . . . . . . . . . 29 2.3.3 Exchange correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3.4 Approximations to the exchange correlation energy . . . . . . . . . 36 2.3.5 Random Phase Approximation . . . . . . . . . . . . . . . . . . . . . 38 2.3.6 Methods Comparison continued . . . . . . . . . . . . . . . . . . . . 39 2.4 Finite temperature properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.4.1 Ab initio Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . 41 2.4.2 Monte Carlo sampling . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.4.3 Thermodynamic perturbation theory . . . . . . . . . . . . . . . . . 43

  The two-body operator Vee involve matrix elements ⟨n 1

• • lmn| Ĥe |ijk • • • lmn⟩. The one-body operators Te and Vext involve matrix elements ⟨n 1 , n 2 , • • •| ĉ † i ĉj |n 1 , n 2 , • • •⟩ which are only non-zero if i = j.

is called Full Configuration Interaction (FCI) and yields the

	a |Φ r a ⟩ +	c rs ab |Φ rs ab ⟩ +	c rst abc Φ rst abc + • • •	(34)
		a<b r<s	a<b<c r<s<t
	where the reference state |Φ⟩ is chosen to be the Hartree-Fock Slater determinant,
	and the wave function possesses the intermediate normalization condition (both
	|Ψ⟩ and |Φ⟩ have unit norm). In that formalism, we can restrict the wave function
	to a given number of orbital substitutions with respect to the Hartree-Fock state
	and solve the Schrödinger equation in that particular truncated Hilbert space. De-
	noting by S,D,T,Q all the singly, doubly, triply and quadruply substituted states,
	there exist different approximated wave functions such as CIS, CISD, CID, CIS-
	DTQ. Considering all possible substitutions (and thus recovering the total wave
	function form) exact
	(in the basis set) energy E 0 .
	For a RHF solution with 2K molecular spin-orbitals, N α = N β spin up and spin
	N α down electrons, there are K	2

STO-LG basis set are minimal basis sets, meaning

  2 radial dependence. The they use only one STO for each atomic orbital, up to the valence shell. The STO are formed by linear combinations of L Gaussian functions. For example, the STO-3G basis for the Lithium is made of five atomic orbitals (ϕ 1s , ϕ 2s , ϕ 2p x , ϕ 2p y , ϕ 2p z , each of which is a sum of three Gaussian functions.

	Double zeta basis sets are constructed by adding one STO to each orbital, usually
	only to the valence ones (split-valence basis sets). The 3-21G basis for the Lithium
	is made of nine orbitals: one ϕ 1s core orbital made of three gaussian (3-21G),
	and four couple of valence orbitals ϕ 1 2s , ϕ 2 2s , ϕ 2 1 p x , ϕ 2 2 p x , dots, where the ϕ 1 k are
	made of two Gaussian functions, (3-21G), and the ϕ 2

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2 Physical descriptors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3 Kernel Ridge Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4 Dimensionality reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.1 3.5 Generative models: Restricted Boltzmann Machine . . . . . . . . . . . . . . 52

Table 1 : Experimental and reference theoretical enthalpies of adsorption for

 1 the CH 4 and CO 2 molecules in zeolites obtained using different functionals. All values are expressed in kcal/mol.

	System	PBE PBE+D2 SCAN vdW-DF2 SCAN+rVV10	Exp.
	CH 4 @HChab -1.13	-6.08	-2.60	-5.64	-4.77	-4.06 [148]
	CO 2 @SiChab -1.18	-5.48	-2.95	-6.93	-7.11	-5

Table 2 : Enthalpies of adsorption estimates of CH 4 @HChab and CO 2 @SiChab computed

 2

	Production -MD simulations

directly from straightforward MD (Ref. column) as well as using the MLPT method, whereby all functionals have been used as production as well as target methods. Deviations of MLPT results from the reference values (Ref. column) are given in parenthesis. All values are expressed in kcal/mol. The results with deviations above 1.0 kcal/mol are in bold.

13.63 (-8.86) -6.71 (-1.94)

						-5.03 (-0.26)	-4.13 (0.64)	-4.74 (0.03)
		PBE	-1.18	-1.15 (0.02)	-1.47 (-0.29) -1.46 (-0.29) -1.24 (-0.06)	-4.68 (-3.50)
		PBE+D2	-5.48 -5.61 (-0.14)	-5.44 (0.04)	-5.40 (0.08)	-5.60 (-0.12)	-7.56 (-2.08)
	CO 2 @SiChab	SCAN	-2.95	-1.42 ( 1.54)	-2.74 (0.21)	-2.96 (-0.01)	-4.30 (-1.34)	-3.72 (-0.77)
	Exp. -5.02 [126]	vdW-DF2	-6.93 -7.47 (-0.54) -7.57 (-0.64)	-8.11 (-1.17)	-6.94 (-0.00)	-8.90 (-1.96)
		SCAN+rVV10 -7.11	-3.97 ( 3.14)	-5.52 ( 1.59)	-6.01 ( 1.10)	-7.17 (-0.06)	-7.09 (0.02)

Table 3 :

 3 Deviations of the MLPT estimates of the target internal energies of individual systems from reference values. All values are in kcal/mol and all values larger than 1.0 kcal/mol in absolute value are in bold.

	Production -MD simulations

Table 4 : Values of the I w index corresponding to each individual MLPT esti- mate of the ensemble internal energy.

 4 The values in parentheses represent the deviations in kcal/mol from the reference values obtained from straightforward MD. In bold, all deviations large than 1.0 kcal/mol.

	Production -MD simulations

Table 5 : Deviations of the SCAN and SCAN+rVV10 internal energies predicted using MLPT and MLMC from

 5 their respective reference values obtained in straightforward MD simulations (Table1). All values are expressed in kcal/mol.

	System	Prediction	MLPT MLMC
	CH 4 @HChab	SCAN	-4.380	-0.639
		SCAN+rVV10 -5.026	-0.858
	CO 2 @SiChab	SCAN	2.249	-0.564
		SCAN+rVV10	3.395	0.961
	HChab	SCAN	4.142	-0.749
		SCAN+rVV10	3.830	-0.230
	SiChab	SCAN	0.659	-0.345
		SCAN+rVV10	0.246	-0.267
	CH 4	SCAN	0.013	-0.011
		SCAN+rVV10	0.000	-0.014
	CO 2	SCAN	0.054	-0.005
		SCAN+rVV10	0.013	-0.058

Table 6 : Experimental enthalpies of adsorption of CH 4 @HChab and CO

 6

	System	Prediction	Ref.	MLPT	MLMC
	CH 4 @HChab	SCAN	-2.60 -11.13 (-8.54)	-2.48 (0.12)

2 @SiChab and their calculated reference MLPT, and MLMC values. The deviations between the predictions and the reference MD values are given in parentheses. All values are expressed in kcal/mol.

Table 7 : Values of the I w index corresponding to the MLPT predictions of RPA energies from PBE+D2 trajectories

 7

	System	Iw	MLPT	MLMC	Difference
	CH 4 @HChab 0.13 -21398.09 -21398.65	-0.57
	CO 2 @HChab 0.21 -21945.80 -21945.97	-0.17
	CH 4 @SiChab 0.26 -21371.82 -21371.74	-0.08
	CO 2 @SiChab 0.27 -21916.60	-	-
	HChab	0.28 -20214.40 -20215.06	-0.66
	SiChab	0.32 -20188.00 -20188.05	-0.05

for the different systems considered in ref.

[START_REF] Chehaibou | Computing RPA Adsorption Enthalpies by Machine Learning Thermodynamic Perturbation Theory[END_REF]

, together with the MLPT and MLMC total electronic energies, and their difference. Energies are in kcal/mol.

  5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.2.1 Periodic boundary conditions coupled cluster . . . . . . . . . . . . 87 5.2.2 Machine Learning Thermodynamic Perturbation Theory and Monte Carlo resampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

Table 8 : Correlation energy contributions to the adsorption energy (kcal/mol) computed using different levels of theory and basis sets

 8 . N v and N o correspond to the number of virtual and occupied orbitals, respectively.

Table 9 : Enthalpy of adsorption of CO 2 in protonated chabazite (kcal/mol) com- puted using different target and sampling methods.

 9

	Target Method	Sampling Method	Enthalpy (kcal/mol)
	PBE+D2	MD	-9.72 ± 0.27
	MP2	MLPT	-9.50 ± 0.24
	CCSD(T)	MLPT	-8.32 ± 0.28
	CCSD(T)	MLMC	-8.09 ± 0.71
	Experiment [147] adsorption isotherms	-8.41

  6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.2 Methodological approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.2.1 Configuration interaction approach . . . . . . . . . . . . . . . . . . 100 6.2.2 Generative model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.2.3 Direct generation of single and double excitations . . . . . . . . . . 103 6.2.4 Computational complexity . . . . . . . . . . . . . . . . . . . . . . . 106 6.2.5 Computational details . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Figure 23: Effect of the temperature in the generative procedure: energy

  convergence for water in the cc-pVDZ basis with one frozen orbital against the number of iteration of the CIgen method, using temperatures T = 1, 5, and 500. The T = 500 curve corresponds to a fully uniform sampling of the connected excited determinants. The chemical accuracy limit and Full CI values are depicted as horizontal lines.However, the generative step (step 4 of the CIgen algorithm) often samples determinants that clearly cannot contribute to the wavefunction and this makes the overall procedure less efficient. For example, it is well known that a fixed set of determinants can only be coupled through the Hamiltonian to new determinants that differ by single and double substitutions. This observation, which has some subtle implications, is exploited by the Monte Carlo Configuration Interaction (MCCI) algorithm proposed by Greer et al.[START_REF] Greer | Monte Carlo configuration interaction[END_REF]. In the MCCI algorithm, new determinants are sampled by randomly proposing single and double excitations of the determinants already included in the wavefunction set. More specifically, a determinant inside the current wave function is randomly chosen, the decision to perform a single or double substitution is chosen at random, and the substituted orbital(s) is (are) randomly chosen. At first sight this could seem a rather inefficient random procedure, equivalent to a uniform sampling among the space of connected determinants, but in practice the important (highly contributing) determinants have often multiple connections to other important determinants through single and double substitutions. Within this sampling method, among all the single and double excited determinants to a given wave function, those that are more connected through single or double excitation are more likely to be sampled, compared to a uniform sampling among them. Consider a wave function made of two determinants a and b, which is connected to four determinants d 1 , d 2 , d 3 , d 4 . Assume that d 1 , d 2 , d 3 are connected to a through single or double excitation, while d 2 , d 3 , d 4

		76.230
		76.232
	Energy [Ha]	76.238 76.236 76.234	CIgen: T=1 CIgen: T=5 CIgen: T=500 ( Random) Chemical Accuracy FCI
		76.240
		76.242
		0	5	10 Iteration number 15	20	25

Table 11 : Total energies (Ha) for the different systems/basis sets: The

 11 -12 , we obtain a total energy of -109.27216 Ha, corresponding to 98.120% of the correlation energy instead of 99.979 % and with only about 156 thousand determinants included in the wavefunction instead of almost 19 millions. CIgen-SD energies are compared to reference FCI and CCSD(T) values. The (N,K) column indicate the size of the active space, where N is the number of correlated electrons and K is the number of active molecular orbitals, and the N conv column indicates the number of iterations to achieve convergence (i.e when the energy in two successive iterations differs by less than 10 -5 Hartree).

	System	Basis	(N,K) E CIgen-SD E CCSD(T )	E FCI	N conv
	C 2	6-31G	(12,18)	-75.64416	-75.64415	-75.64418	12
	N 2	6-31G	(14,18) -109.10824 -109.10635 -109.10842	12
	H 2 O	6-31G	(10,13)	-76.12220	-76.12182	-76.12237	10
	C 2	cc-pVDZ (8,26)	-75.72982	-75.72781	-75.72984	9
	N 2	cc-pVDZ (10,26) -109.27827 -109.27829 -109.27834	13
	H 2 O	cc-pVDZ (8,23)	-76.24192	-76.24131	-76.24195	17

In particular, if there are no interaction between electrons, and if the |ϕ i ⟩ are eigenvectors of the single electron problem, with sorted eigenvalues ϵ i , then the Slater determinant |12 • • • N⟩ is the groundstate to the N electron problem (up to a normalizing constant 1√ N!), with energy E = N i=1 ϵ i .

It is possible to alternatively enforce the wavefunction normalization and keep the freedom of nonorthonormal basis. In that case a general eigenvalue problem Fϕ = ESϕ is found where S ij = ⟨ϕ i |ϕ j ⟩ is the overlap matrix[START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF].

The product of the two non-interacting FCI wave functions will correctly yield the FCI wave function of the total system.

The divergent Coulomb interaction at coinciding positions of particles in the wave function must be compensated by a non-differentiability of the wave function, that will in turn produce an opposite divergent kinetic energy term[START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF].

Note that since the reduced density matrix has positive or zeros eignevalues[START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF], its diagonalization is equivalent to a Singular Value Decomposition[START_REF] Martin | The Extraordinary SVD[END_REF], and keeping only the highest eigenvalues is a dimensionality reduction technique (it is a Principal Component Analysis, see Section 3.4).

In this chapter and the last one, we do not discuss recent linear scaling algorithms, see e.g. ref.[START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF] for DFT, and ref.[START_REF] Péter | Optimization of the Linear-Scaling Local Natural Orbital CCSD(T) Method: Improved Algorithm and Benchmark Applications[END_REF] for CCSD(T).

or alternatively by n(r) = ⟨ n(r)⟩ = ⟨Ψ| i δ(r -r i )|Ψ⟩

If |Ψ 1 ⟩ = |Ψ 2 ⟩, then ( V2 ext -V1 ext ) |Ψ 1 ⟩ = ( T + Vee + V2 ext ) -( T + Vee + V1 ext ) |Ψ 1 ⟩ = (E 2 -E 1 )|Ψ 1 ⟩ so that both potential only differ by a constant in contradiction with the assumption.

This factorization, know as Schmidt decomposition, is a direct application of the singular value decomposition[START_REF] Wouters | A Practical Guide to Density Matrix Embedding Theory in Quantum Chemistry[END_REF], that will be described in Section 3.4

we can always center our dataset, by subtracting the mean of the x i and y i

Consider a model where g(x) = ax 2 , considering the map ϕ(x) = x 2 , we have a linear model in the space defined by ϕ, g(x) = a • ϕ(x).

C = UΣV T (UΣV T ) T = UΣV T VΣ T U T = UΣΣ T U T = UDU T .

Actually X k has the same dimensions of X but is of rank k, we can still pad it with zeros.

Note that the Boltmann machine model is recovered by taking a normal probability for q a (h a ), whose only non-zero cumulant is the second one.

A S S E S S I N G T H E

The content of this chapter is nearly identical to the following published work: [84] "Assessing the Accuracy of Machine Learning Thermodynamic Perturbation Theory: Density Functional Theory and Beyond", by Basile Herzog, Maurício Chagas da Silva, Bastien Casier, Michael Badawi, Fabien Pascale, Tomas Bucko, Sébastien Lebègue and Dario Rocca, published in the Journal of Chemical Theory and Computation 18.3 (2022), pp. 1382-1394.

introduction

Correlated quantum chemical methods could provide an alternative to density functional theory (DFT) in (periodic) materials simulations possibly reaching the threshold of chemical accuracy (1 kcal/mol) with respect to experimental data. Traditional quantum chemical methods, such as Møller-Plesset perturbation theory to second order (MP2) [137] or coupled-cluster theory [START_REF] Bartlett | Coupled-cluster theory in quantum chemistry[END_REF], have been recently implemented for condensed phase materials applications [START_REF] Booth | Towards an exact description of electronic wavefunctions in real solids[END_REF][START_REF] Del Ben | Second-Order Møller-Plesset Perturbation Theory in the Condensed Phase: An Efficient and Massively Parallel Gaussian and Plane Waves Approach[END_REF][START_REF] Dixit | Communication: A novel implementation to compute MP2 correlation energies without basis set superposition errors and complete basis set extrapolation[END_REF][START_REF] Marsman | Second-order Møller-Plesset perturbation theory applied to extended systems. I. Within the projectoraugmented-wave formalism using a plane wave basis set[END_REF][START_REF] Pisani | Periodic local MP2 method for the study of electronic correlation in crystals: Theory and preliminary applications[END_REF]] . An alternative particularly suitable for condensed matter applications is represented by the random phase approximation (RPA) [START_REF] Bohm | A Collective Description of Electron Interactions: III. Coulomb Interactions in a Degenerate Electron Gas[END_REF][START_REF] Dobson | Successful Test of a Seamless van der Waals Density Functional[END_REF][START_REF] Furche | Molecular tests of the random phase approximation to the exchange-correlation energy functional[END_REF][START_REF] Gell-Mann | Correlation Energy of an Electron Gas at High Density[END_REF][START_REF] Harl | Cohesive energy curves for noble gas solids calculated by adiabatic connection fluctuation-dissipation theory[END_REF][START_REF] Langreth | The exchange-correlation energy of a metallic surface[END_REF][START_REF] Lu | Ab initio Calculation of van der Waals Bonded Molecular Crystals[END_REF] and its variants which include higher order corrections [START_REF] Bates | Communication: Random phase approximation renormalized many-body perturbation theory[END_REF][START_REF] Colonna | Correlation energy within exact-exchange adiabatic connection fluctuation-dissipation theory: Systematic development and simple approximations[END_REF][START_REF] Dixit | Improving the accuracy of ground-state correlation energies within a plane-wave basis set: The electron-hole exchange kernel[END_REF][START_REF] Dixit | Improving the Efficiency of Beyond-RPA Methods within the Dielectric Matrix Formulation: Algorithms and Applications to the A24 and S22 Test Sets[END_REF][START_REF] Grüneis | Making the random phase approximation to electronic correlation accurate[END_REF][START_REF] Hellgren | Beyond the random phase approximation with a local exchange vertex[END_REF][START_REF] Olsen | Extending the random-phase approximation for electronic correlation energies: The renormalized adiabatic local density approximation[END_REF]. Due to the significantly high computational cost of these approaches, their use is limited, especially in (finite-temperature) molecular dynamics (MD) simulations.