Despite much efforts towards the use of renewable energies, fossil fuels still dominate the global energy production leading to substantial greenhouse-gas emissions to the atmosphere. The use of highly porous, chemically tunable solid adsorbents such as metalorganic frameworks (MOFs) has attracted much interest in the past 20 years for their potential application in gas separation and carbon capture technologies. In this work, we propose to computationally design MOFs whose high affinity for guest molecules can be modified under temperature treatment. With a suitable choice of ligands and metal centers, MOFs exhibiting a thermally-induced spin crossover (SCO) and a concomitant change in the adsorption properties can be developed to possibly yield more energy-efficient gas capture-and-release. From a computational point of view, the challenge in simulating the SCO phenomenon using ab initio electronic structure methods is represented by the evaluation of the thermodynamics of SCO. To address this issue, we develop a DFT+Ubased scheme to achieve an accurate description of spin-state energy differences in Fe(II) complexes. A critical assessment of the Hubbard U approach in the description of spinstate energetics is presented and insights from this analysis are used to propose a practical and efficient way to overcome them. This approach, which is tested and validated against experiments and coupled cluster-corrected CASPT2 results, is then used to propose the proof-of-concept in silico of temperature-induced SCO-MOF for efficient gas release.

Résumé en français

Malgré de nombreux efforts en faveur de l'utilisation des énergies renouvelables, les combustibles fossiles continuent de dominer la production énergétique mondiale, ce qui entraîne des émissions importantes de gaz à effet de serre dans l'atmosphère. L'utilisation d'adsorbants solides très poreux et chimiquement accordables, tels que les cadres métalloorganiques (MOF), a suscité beaucoup d'intérêt au cours des 20 dernières années en raison de leur application potentielle dans les technologies de séparation des gaz et de capture du carbone. Dans ce travail, nous proposons de concevoir par calcul des MOFs dont la haute affinité pour les molécules invitées peut être modifiée par traitement thermique. Avec un choix approprié de ligands et de centres métalliques, les MOF présentant un croisement de spin induit thermiquement (SCO) et un changement concomitant des propriétés d'adsorption peuvent être développés pour éventuellement produire une capture et une libération de gaz plus efficaces sur le plan énergétique. D'un point de vue informatique, le défi de la simulation du phénomène de SCO en utilisant des méthodes de structure électronique ab initio est représenté par l'évaluation de la thermodynamique du SCO.

Pour résoudre ce problème, nous développons un schéma basé sur la DFT+U pour obtenir une description précise des différences d'énergie à l'état de spin dans les complexes de Fe(II). Une évaluation critique de l'approche Hubbard U dans la description de l'énergie à l'état de spin est présentée et les résultats de cette analyse sont utilisés pour proposer une manière pratique et efficace de les surmonter. Cette approche, qui est testée et validée par rapport aux expériences et aux résultats de CASPT2/CC, est ensuite utilisée pour proposer la preuve de concept in silico du SCO-MOF induit par la température pour une libération efficace du gaz.

Chapter 1

Introduction 1.1 Motivation

Since the first industrial revolution, human activity has been characterized by a massive production of industrial waste. Their harmful effects on human health and on the environment have become nowadays clear and unquestionable. However, if on one hand, a complete transition in the direction of non-polluting energies and industrial production is desirable, on the other hand, the associated technologies cannot ensure a quick energy revolution and a stage where green and polluting energies coexist together must be foreseen. In this context, gas adsorption and separation can play a crucial role in limiting toxic emissions as well as global warming.

Metal-organic frameworks (MOFs) are an emerging class of nanoporous materials with unique structural and electronic properties that make them promising candidates to be applied in gas adsorption and separation. The use of computational methods to design new MOFs with selected features turned out to be an efficient strategy to facilitate the screening process and to highlight the microscopic nature of the gas adsorption mechanisms. Recent studies [1,2,3,4,5] have highlighted the importance of a full quantum mechanical description of the interaction between gas and MOF especially in frameworks characterized by the presence of open-metal sites (OMS). For these MOFs the electronic configuration of the OMS determines the binding mechanism with the guest molecules.

These findings have paved the way for the development of new strategies where a change in the MOF electronic structure upon gas adsorption can modify the binding mechanism itself and consequently potentially achieve unprecedented performances in gas capture and release. In 2017, Jeffrey R. Long and collaborators [START_REF] Reed | A spin transition mechanism for cooperative adsorption in metal-organic frameworks[END_REF] have demonstrated how Fe open metal sites can undergo a gas pressure-induced cooperative spin-transition upon gas ad-Introduction sorption with a consequent change in the metal-gas binding mechanism and a resulting exceptional working capacity.

Motivated by this recent experimental work and by previous computational studies [1] in our group demonstrating the importance of the electronic configuration for the binding mechanism, the scope of my PhD thesis was to provide a proof-of-concept, in silico, for the use of spin-crossover (SCO) MOFs where the binding mechanism with the guest molecules can be modified through a temperature swing as a consequence of the change in spin state of the system. SCO materials are transition-metal compounds whose spin state can change as a result of an external perturbation such as temperature, light irradiation, pressure or external magnetic fields; the induced spin-transition is called SCO transition.

In this thesis, I investigate the possibility of inducing a temperature-driven SCO transition on the OMS upon gas adsorption. In this way a change in the occupation of the σ * antibonding orbitals between the OMS and the molecule upon SCO will lead to a change of the binding strength and, possibly, to the desorption of the gas. The practical implementation as well as the technological relevance of this approach depends on the specific value of the SCO transition temperature. The idea is that once the gas is adsorbed, a SCO occurring at temperatures slightly larger than the adsorption temperature will potentially result in capture-and-release processes with reduced energy penalties and thus reduced cost. It follows that the methodologies employed to investigate SCO MOFs must be able to correctly capture the thermodynamic of SCO in order to predict as precisely as possible the SCO transition temperature.

Achieving an accurate ab initio description of SCO transition in transition-metal complexes represents a great challenge for modern quantum chemistry and solid-state first principle methods. Chemical accuracy can be achieved only for small complexes using accurate quantum chemistry methods (full Configuration Interaction (full-CI), complete active space multiconfigurational method followed by second-order perturbation theory (CASSCF/CASPT2) or coupled-cluster (CC)), but these are extremely computational demanding and cannot be applied to study MOFs which contain hundreds of atoms in the unit cell. Thus, during my PhD I have developed, tested and validated a new strategy based on the density-functional theory (DFT) and the Hubbard U-correction method to efficiently compute the spin-state energetics of Fe(II) molecular and crystalline complexes with good accuracy compared with experimental values or with higher level theory (depending on the availability) at the cost of standard DFT calculation.

1.2 Thesis outline

Thesis outline

In chapter 2 the SCO phenomenon is presented. After a brief historical introduction to contextualize the subject, the relevant scientific aspects are reviewed: the underlying physics of SCO is approached by means of ligand field theory and thermodynamics. The chapter is concluded by presenting the state-of-the-art for what concerns the computational investigation and simulations of spin-transitions as well as the problematics concerning the evaluation of spin state energetics in transition metal complexes.

Chapter 3 is an overlook of MOFs and their application for gas adsorption purposes.

In particular, I discuss the influence of the OMS electronic structure on the adsorption properties. The thermodynamics of gas adsorption and some basic concepts of molecular orbitals (MOs) theory are presented in order to prepare the ground for the understanding of the results discussed in the following chapters.

In chapter 4 the computational methods employed to produce all the original results presented in this thesis are reviewed.

Each subsequent chapter is a study published (chapter 5), submitted for publication (chapter 6), or ready for publication (chapter 7).

Chapter 5 is a detailed analysis of the performances of DFT+U method in computing spin energetics in octahedrally coordinated Fe(II) molecular complexes. Comparing the spin gap computed with DFT and DFT+U with the accurate CASPT2/CC method, I show that, within the DFT+U framework, low-spin states are systematically penalized over high-spin states. The Hubbard energy functional introduces an energy penalty for lowspin states whose extent and origin is discussed in detail. I show that the linear response DFT+U method, despite its dramatic failure in computing the spin state energetics, leads to a systematic improvement of the electronic charge density reducing the spurious curvature, due to the self-interaction error (SIE), of the DFT energy as a function of the total number of electrons. In this chapter MOs arguments are used to rationalize the observed bias and to explain the metal-to-ligand charge transfer found upon U correction.

Motivated by the previous study, in chapter 6 I adopt a new approach consisting of a Hubbard U density-corrected DFT where the PBE functional is evaluated on the Hubbard U density, using a linear-response U. By doing so, the energy bias introduced by the Hubbard energy term is removed while keeping an improved electronic density. The adiabatic energy differences computed using this method for six molecular complexes yield a very good agreement with the reference CASPT2/CC values resulting in a main absolute error smaller than the best performers among the standard DFT functionals (TPSSh and Introduction M06-L). I then apply it to study several crystalline compounds for which experimental data are available and again an excellent agreement is found.

In chapter 7 the DFT+U density-corrected method presented and validated in chapter 6 is used to design MOFs employing a SCO transition to efficiently desorb the adsobed gas. This is done by computing the adiabatic energy differences between different spin states in presence of adsorbate molecule, and the associated change in binding energy.

A screening study is performed on M(CN) 4 fragment (M=Fe(II), Fe(III), Ni(II), Co(II))

with adsorbate CO, CO 2 , N 2 , H 2 and H 2 O molecules. Two different mechanisms behind the change in binding energy upon spin transition are observed. The first is linked with the (de)population of the metal d z 2 orbital upon SCO transition and the second, never analyzed in detail before, results in enhanced electrostatic interaction between the adsorbate and the metal center when the metal d x 2 -y 2 orbital passes from being empty to be singly occupied. I conclude my thesis by reporting calculations on the Hofmann-clathrate Fe(bpac)Fe(CN)4 and I show that this is a promising candidate for SCO-assisted gas capture/release applications, thus providing the proof-of-concept for our idea.

Chapter 2

Spin-crossover phenomena

Introduction

In 1931 L. Cambi and L. Szegö reported for the first time unusual magnetic properties of several dithiocarbamates compounds based on Fe(III) ions [START_REF] Cambi | Über die magnetische susceptibilität der komplexen verbindungen[END_REF]. A change from paramagnetic to the diamagnetic state was observed as a result of a temperature variation. A few years later, Pauling and his collaborators found similar magnetic behavior in ferrihemoglobin [START_REF] Coryell | The magnetic properties and structure of ferrihemoglobin (methemoglobin) and some of its compounds[END_REF]. The underlying physics of these observations was not immediately grasped and it took around 25 years for the scientific community to understand it. The experimental observations are rationalized on the basis of the ligand field-theory and the Boltzmann statistics [START_REF] Ballhausen | Some comments on the anomalous magnetic behavior of certain ni(II) complexes[END_REF][START_REF] Stoufer | UNUSUAL MAG-NETIC PROPERTIES OF SOME SIX-CO ÖRDINATE COBALT(II) COM-PLEXES1-ELECTRONIC ISOMERS[END_REF]: when the energy gap that separates two spin states at 0K is small, the thermal energy k b T can induce a change in the electronic population of these states and a modification of the magnetic state of the system. The relative energy difference between two spin states depends primarily on the electronic configuration of the transition-metal ion of the system and on its coordination environment (e.g. the nature of the axial ligands).

Such a dependency was later investigated in 1961 by Busch et al. [START_REF] Figgins | Complexes of iron(II), cobalt(II) and nickel(II) with biacetyl-bis-methylimine, 2-pyridinal-methylimine and 2, 6-pyridindial-bismethylimine[END_REF], who also reported for the first time a spin transition in Co(II) compounds. In 1964, Baker and Bobonich observed for the first time a spin transition in Fe(II) complexes, [Fe(phen) 2 (NCX) 2 ] (X = S, Se; phen= 1,10-phenanthroline) and [Fe(bpy) 2 (NCS) 2 ] (bpy=2,2'-bipyridine) [START_REF] Baker | Magnetic properties of some high-spin complexes of iron(II)[END_REF]. Soon after, in 1967, König and Madeja extensively investigated these complexes by employing magnetic measurements, and electromagnetic and Mössbauer spectroscopy [START_REF] Koenig | 5t2-1a1 equilibriums in some iron(II)-bis(1, 10phenanthroline) complexes[END_REF] and confirmed that the nature of the spin-crossover (SCO) phenomenon is linked with the ability of certain transition metal complexes to undergo a change in their spin configuration as a result of external stimuli (temperature, pressure or light).

Starting from the 70's the number of studies on SCO, especially on Fe(II) compounds, Spin-crossover phenomena has rapidly increased [START_REF] Gütlich | Spin crossover in transition metal compounds I[END_REF], and the technological applications of these bistable materials are found today in molecular spintronic, molecular electronics, sensors, and actuators [START_REF] Kumar | Emerging trends in spin crossover (SCO) based functional materials and devices[END_REF][START_REF] Molnár | Molecular spin crossover materials: Review of the lattice dynamical properties[END_REF][START_REF] Resines-Urien | A switchable ironbased coordination polymer toward reversible acetonitrile electro-optical readout[END_REF][START_REF] Resines-Urien | Covalent post-synthetic modification of switchable iron-based coordination polymers by volatile organic compounds: a versatile strategy for selective sensor development[END_REF].

Ligand-field theory and spin transition

In transition-metal complexes, the relative equilibrium between different spin states depends on the transition-metal element itself as well as on the coordination geometry around it and the nature of the metal-to-ligand interactions. The ligand-field theory is used to rationalize SCO. A complete and detailed presentation of this theory can be found in ref [START_REF] Hauser | Ligand field theoretical considerations[END_REF]. I restrict here my discussion to octahedrally coordinated transition metal compounds and I consider the metal ion to be Fe(II). This choice is motivated by the fact that: 1) most of the results presented in this thesis concern octahedrally coordinated Fe(II) compounds; 2) it allows me to easily introduce all the basic concepts needed for a satisfactory understanding of the spin-crossover phenomena.

In figure 2.1 the energy diagram for the 3d electronic levels is shown for an isolated atom 2.2 Ligand-field theory and spin transition 11 and for an ion in an octahedrally coordinated configuration. Because of the electrostatic interaction between the metal and the ligands, the 5-fold degenerate 3d levels are split into two energy levels e g and t 2g , which are 2-fold and 3-fold degenerate, respectively. The energy gap between the e g and the t 2g levels is called ligand-field splitting, ∆, and it depends on the metal ion, the metal-to-ligand distances, and the nature of the metal-ligand interaction. The way in which the electrons fill the e g and t 2g levels depend on the competition between ∆ and the electronic repulsion, Π, which results when two electrons sit on the same energy level. Hence it follows that when the electrostatic interaction between the iron and the ligands is weak (i.e. ∆ < Π) the electrons fill the energy levels following the Hund's rule. In this case, the system is in the paramagnetic high-spin (HS) state, 5 T 2g (t 4 2g e 2 g ). On the other hand, when the interaction with the ligand is strong due to the presence of a strong ligand field (i.e. ∆ > Π), the energy cost of filling the e g levels is too large and the t 2g levels are doubly occupied. The diamagnetic electronic configuration resulting from the condition ∆ > Π is low-spin state (LS), 1 A 1g (t 6 2g e 0 g ). If the energy difference between HS and LS state is of the order of the thermal energy, k b T , then an SCO transition may occur upon heating. The competition between Π and ∆ is schematically represented in the Tanabe-Sugano's diagram reported in figure 2.2. Here, the energy of the different iron 3d spin configurations, expressed through their molecular-symmetry labels, are plotted as a function of the ligand-field splitting ∆. This diagram shows how for ∆ < Π the 5 T 2g state (HS) is energetically more stable compared to other symmetries.

For increasing ligand fields the 1 A 1g (LS) is gradually stabilized to the critical point ∆ c where a crossover between the HS and LS energetics occurs and the magnetic state of the system changes. Such a behavior is at the origin of the term: spin crossover.

The SCO transition occurs along with a change in the volume of the complex. Because the e g orbitals have an anti-bonding character between the iron and the ligands, their occupation results in lowering the metal-ligand bond order with concurrent elongation of the metal-ligand distance, r. This change has a significant effect on the vibrational properties and, therefore, the potential energy surface E p in the two spin states assumes different shapes. This is due to the different force constant of the vibrations f through which E p is expressed within the harmonic approximation as E p,i (r) = 1 2 (f i r 2 i ) (where i = LS, HS). In figure 2.3 a schematic representation of E p as a function of the metal-to-ligand distance r is shown. Inside the two parabolic E p profiles, horizontal lines represent the energy levels associated with the symmetric harmonic modes of the complex. The quantity ∆E H-L is the internal energy difference at zero temperature that, within the Born-Oppenheimer approximation or adiabatic approximation (see section 4.2), is given by the sum of the (2.1)

The schematic diagram of figure 2.3 represents a situation in which the ground-state of the system is LS at zero temperature. As the temperature increases, the thermal energy brings the system to occupy higher vibrational states. Thus, HS is gradually stabilized due to the greater increase in its entropic term compared to that of LS. This is because HS has larger spin degeneracy and higher vibrational density of states. The crossing point of the two curves, where the geometries associated with LS and HS states have the same metal-to-ligand distance r, corresponds to the ∆ c point in the Tanabe-Sugano's diagram which represents an energetically unstable stage where the SCO complex undergoes a LS→HS spin-transition. LS↔HS transition can also be induced by light irradiation or by the action of pressure. In the first case, the population of the two spin states can be reversed through a process of excitation-relaxation induced by irradiating the sample. In the second case, a pressure (chemical or mechanical) applied to the system results in a stronger ligand field for the metal ion with a consequent destabilization of HS state. Both light-and pressure-induced spin transitions do not have a central role in this thesis and they will not be further discussed, see ref. [START_REF] Chastanet | Light-induced spin crossover-solution and solid-state processes[END_REF][START_REF] Gaspar | Pressure effect investigations on spin-crossover coordination compounds[END_REF] for more details.

Jahn-Teller distorsion

The change in the metal-to-ligand distance r upon spin transition is due to the change in the occupation of the anti-bonding e g orbitals. The volume expansion observed in octahedrally coordinated Fe(II) complexes when passing from LS to HS state is accompanied by a geometrical distortion that reduces the total symmetry. This effect is called Jahn-Teller distortion and it is named after the two scientists who, in 1937, postulated a theorem stating that "orbital electronic degeneracy and stability of the nuclear configuration are incompatible unless all the atoms of a molecule lie on a straight line" [22]. This means that when orbital degeneracy occurs in octahedrally coordinated Fe(II) complexes, the system lowers its energy by distorting the coordination environment of the ion to lift the orbital degeneracy. This is the case for the Fe(II) HS configuration, where the degeneracy arises because of the double occupation of one of the three t 2g states (see figure 2 shown above (the compression is here exaggerated to illustrate qualitatively the effect of the distortion).

Thermodynamics of spin-crossover

The macroscopic manifestation of SCO phenomena is the result of the collective behavior of a group of metal SCO centers interacting with an external environment. To correctly describe the nature of the spin transition at this scale, we need to approach the problem from a thermodynamical point of view. We consider an ensemble of N noninteracting SCO centers, that can be either in HS or LS state, in contact with a thermal bath T and kept at a constant pressure P. We define the fraction of HS centers as γ HS = N HS /N , where N HS is the number of SCO centers in HS state. The LS→HS transition is represented by a change of γ HS from 0 to 1. The thermodynamics of the system is described by the

Thermodynamics of spin-crossover

Gibbs free energy G = H -T S, with H and S enthalpy and entropy, respectively. The spin state of the system is thus determined by the relative equilibrium between the HS and LS phases, which is expressed through the Gibbs free energy difference:

∆G = G HS -G LS = ∆H -T ∆S, (2.2)
where ∆H = H HS -H LS and ∆S = S HS -S LS represent enthalpy and entropy variations.

We define the transition temperature T 1/2 as the temperature at which the two phases are at equilibrium, i.e. when ∆G = 0. The transition temperature can be obtained as:

T 1/2 = ∆H ∆S . (2.
3)

The Gibbs free energy difference in eq. 2.2 can be expressed in terms of the HS fraction γ HS through the equilibrium constant K eq as:

∆G = -RT ln K eq = -RT ln γ HS 1 -γ HS , (2.4)
where R is the gas constant. When T < T 1/2 the enthalpic term is dominant in eq. 2.2, ∆G is negative and the LS state is favoured (γ HS < 1/2). On the contrary for T > T 1/2 , G HS < G LS because T ∆S > ∆H and the HS state is favoured (γ HS > 1/2). The temperature-induced SCO transition is an entropy-driven phase conversion.

The enthalpy and entropy differences can be computed using electronic structure methods in terms of temperature-independent (or electronic) and temperature-dependent contributions as:

∆H(T ) = ∆E HL + ∆E vib (T ) (2.5)
and ∆S(T ) = ∆S el + ∆S vib (T ) + ∆S rot (T ) + ∆S trans (T ).

(2.6)

The electronic contribution to ∆H is the adiabatic energy difference. ∆S el is the electronic entropy variation due to the different spin degeneracy of LS and HS states and it can be expressed in terms of the total spin momenta S j (j = HS, LS) as:

∆S el = R ln 2S HS + 1 2S LS + 1 . (2.7)
In the case of LS-HS transition in Fe(II) complexes between S LS = 0 and S HS = 2 we have ∆S el = R ln 5 = 13.38Jmol -1 K -1 . The temperature-dependent terms in ∆S(T ) are the vibrational entropy variation ∆S vib (T ), the rotational entropy difference ∆S rot (T ), and the translation entropy variation ∆S trans (T ). The two latter terms are absent in solid state. The expression for the vibrational energy E vib and the vibrational entropy S vib can Spin-crossover phenomena be derived within the statistical thermodynamic theory from the partition function of an harmonic oscillator. For a given normal modes vibrational spectrum {ν i }, the vibrational contributions are:

E vib = R ν i hν i k B 1 2 + 1 e hν i /k B T -1 (2.8)
and

S vib = R ν i hν i k B T 1 e (hν i /k B T ) -1 -ln 1 -e -hν i /k B T .
(2.9)

Thus, combining eq. 2.2-2.8 it is possible to compute γ HS (T ), as well as the transition temperature T 1/2 with ab initio approaches.

It is worth making a final remark about the SCO model just presented: the assumption of non-interacting SCO centers allows for an accurate evaluation of the transition temperature but it is not expected to reproduce the correct shape of the transition curve γ HS (T ) observed experimentally [START_REF] Cirera | Theoretical prediction of spin-crossover temperatures in ligand-driven light-induced spin change systems[END_REF]. This is because the real cooperative behavior of SCO centers observed in bulk materials is not considered. Such an effect is responsible, for instance, for the first-order spin transition measured experimentally and accompanied by large hysteresis [START_REF] Brooker | Spin crossover with thermal hysteresis: practicalities and lessons learnt[END_REF].

Calculation of ∆E H-L using electronic structure methods

The ab initio description of SCO phenomena represents an unquestionable challenge for ab initio electronic structure methods. Even in the assumption of non-interacting SCO centers, the evaluation of the quantities that enter eq. 2.2 require an accurate description of electronic correlations, dispersion interactions, relativistic effects and vibrational properties [START_REF] Kepp | Consistent descriptions of metal-ligand bonds and spin-crossover in inorganic chemistry[END_REF]. The vibrational contribution to the enthalpy and entropy differences, ∆H vib and ∆S vib , can be computed with good accuracy within the harmonic approximation [START_REF] Cirera | Theoretical prediction of spin-crossover temperatures in ligand-driven light-induced spin change systems[END_REF][START_REF] Cirera | Computational modeling of transition temperatures in spincrossover systems[END_REF][START_REF] Kepp | Theoretical study of spin crossover in 30 iron complexes[END_REF]. On the other hand the calculation of the remaining term, the adiabatic electronic energy difference ∆E H-L , represents a true challenge and it is the main obstacle to calculate the SCO transition temperature T 1/2 from first principles.

The difficulty of a correct evaluation of ∆E H-L stems from the lack of error cancellation when computing energy differences between spin states exhibiting different types and amount of electronic correlations. Electron correlations are due to the complex interelectron interactions and different electronic structure methods treat them differently. In quantum chemistry, the correlation energy E corr is defined as the difference between the total nonrelativistic energy E and the total self-consistent Hartree-Fock (HF) energy E HF 2.5 Calculation of ∆E H-L using electronic structure methods 17 [START_REF] Shavitt | The method of configuration interaction[END_REF]:

E corr = E -E HF , (2.10)
where the HF energy takes into account only the Coulomb and exchange interactions between electrons within a main-field single Slater-determinant picture. For small systems, the exact solution of the Schrödinger equation, and hence the evaluation of the total energy E (for a given basis set), can be obtained with full configuration interaction (CI) approach [START_REF] Bytautas | The range of electron correlation between localized molecular orbitals. a full configuration interaction analysis for the NCCN molecule †[END_REF][START_REF] Olsen | Surprising cases of divergent behavior in mo/ller-plesset perturbation theory[END_REF]. The computational cost of full-CI method grows exponentially with the size of the system [START_REF] Olsen | Passing the one-billion limit in full configuration-interaction (FCI) calculations[END_REF][START_REF] Rossi | A full-configuration benchmark for the n2 molecule[END_REF] and its use to study Fe(II) SCO complexes is computationally infeasible. Other (approximate) computational methods can be used to capture the essential features of the electronic interactions needed for a good description of the relevant physical quantities under investigation.

The electronic correlations, absent in the HF wavefunction by definition, can be di- The amount of static correlation is thus proportional to the number and the strength of the occupied σ bonding MOs which result in an electronic charge transfer from the ligand to the metal and a consequent increase in the multiconfigurational character of the electronic ground-state wave function [START_REF] Domingo | Spin crossover in fe(II) complexes: An ab initio study of ligand σ-donation[END_REF][START_REF] Pierloot | The caspt2 method in inorganic electronic spectroscopy: from ionic transition metal to covalent actinide complexes[END_REF]. The LS state has two doubly occupied σ bonding orbitals and two empty σ * orbitals. In HS, the two anti-bonding σ * orbitals are singly occupied resulting in lower metal-ligands bond order and lower static correlation contribution. Also dynamical correlation enters differently in the HS and the LS states due to the different number of parallel and antiparallel electrons. In addition, static correlation also depends on the metal-ligand strength, and a single Slater-determinant picture may become less and less accurate in describing Fe(II) SCO complexes moving along the Spin-crossover phenomena spectrochemical series.

vided
Dynamical correlation depends on the electronic repulsion that is generally higher for antiparallel electrons sitting on the same energy level (antiparallel intraorbital pairs) and lower for antiparallel electrons on different orbitals (antiparallel interorbital pairs); the lowest repulsion is found for parallel interorbital pairs instead. In LS there are three antiparallel intraorbital pairs and three antiparallel interorbital pairs while, in HS, the total number of antiparallel electron pairs is five with only one antiparallel intraorbital pair.

Thus the dynamical correlation is expected to be higher in LS compared to HS.

For a method to be able to correctly calculate ∆E H-L in Fe(II) complexes it has to capture correctly all these correlations, both in HS and LS and for different ligand field strengths, otherwise the different treatment in the two spin states does not result in error cancellation and biased total energies are found. In HF method the exchange term keeps parallel interorbital electrons far apart resulting in an underestimation of their Coulomb interaction compared to antiparallel electrons. The result is a strong overstabilization of HS over LS. On the other hand, DFT within its standard local density and generalized gradient approximations tend to overstabilize LS over HS [START_REF] Ganzenmüller | Comparison of density functionals for differences between the high-(t2g5) and low-(a1g1) spin states of iron(II) compounds. IV. results for the ferrous complexes [fe(l)('NHS4')[END_REF]. By adding a fraction of exact exchange, by means of hybrid functionals one can reduce such overstabilization [START_REF] Pinter | Conceptual insights into DFT spin-state energetics of octahedral transition-metal complexes through a density difference analysis[END_REF].

The main drawback in using hybrid functionals is that the amount of exact exchange needed to correctly compute ∆E H-L is, in general, system-dependent. Among the post-HF methods, coupled-cluster singly and double theory with the perturbative inclusion of triple excitations, i.e. CCSD(T), has been successfully implemented to calculate ∆E H-L

for systems where the electronic wavefunction has a small multiconfigurational character, i.e. in weak field ligand complexes [START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF]. For strong field ligand complexes, CASPT2

with a suitable choice of the active space can be implemented to capture both dynamical and static correlation [38]. In a recent study [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF], it has been shown that CASPT2

tends to underestimate ∆E H-L because of a lack of electrons core 3s3d treatment. These 3s3p correlation effect on the computed adiabatic energy difference can be computed with CCSD(T) and then added to the CASPT2 result. This combination of CCSD(T) and CASPT2 is called CASPT2/CC method [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined CASPT2/CC approach[END_REF] and it is the method implemented in this thesis to compute the reference values used to benchmark DFT functionals. For adiabatic energy difference calculations, some author also suggested DMC method [START_REF] Fumanal | Diffusion monte carlo perspective on the spin-state energetics of [fe(NCH)6]2+[END_REF][START_REF] Song | Benchmarks and reliable DFT results for spin gaps of small ligand fe(II) complexes[END_REF].

Which method should be taken as a reference is a matter of debate in the scientific community and a generally recognized "golden standard" does not exist [START_REF] Römer | A review of density functional models for the description of fe(II) spin-crossover complexes[END_REF]. In fact, the three abovementioned "accurate" quantum chemistry methodologies, i.e. CCSD(T), CASPT2/CC and DMC, when compared together on the same systems may show, as a consequence of the different approximations, a significant deviation of the computed adiabatic energy differences ∆E H-L (see, for example, table 6.1 in chapter 6). For some Fe(II)

complexes ∆E H-L has been extracted from experimental values of enthalpy differences ∆H once the environmental experimental effects are removed [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF][START_REF] Radoń | Benchmarking quantum chemistry methods for spin-state energetics of iron complexes against quantitative experimental data[END_REF][START_REF] Mata | Benchmarking quantum chemical methods: Are we heading in the right direction?[END_REF][START_REF] Daku | Assessment of density functionals for the high-spin/low-spin energy difference in the low-spin iron(II) tris(2, 2 -bipyridine) complex[END_REF]. We do make use also of these extracted values for our analysis in chapter 6 but an internal reference is still needed to investigate the performances of electronic structure methods when weak field ligand systems are considered (e.g.

[Fe(H 2 O) 6 ] +2 , [Fe(NCH) 6 ] +2 , [Fe(NH 3 ) 6 ] +2 ) and
for which experimental values do not exist.

Chapter 3

Metal-organic frameworks for efficient gas separation

Introduction

Metal-organic frameworks (MOFs) are a class of nanoporous crystalline materials made of metal ions centers coordinated by organic linkers to form 3D-structures [START_REF] Li | Design and synthesis of an exceptionally stable and highly porous metal-organic framework[END_REF] (See figure 3.1). The large number of possible combinations of organic linkers and metal sites results in a high tunability of physical and chemical properties [START_REF] Kuppler | Potential applications of metalorganic frameworks[END_REF]. This feature has resulted in more than 70000 synthesized MOFs published in CSD [START_REF] Moghadam | Development of a cambridge structural database subset: A collection of metal-organic frameworks for past, present, and future[END_REF] with a wide range of applications like catalysis [START_REF] Lee | Metal-organic framework materials as catalysts[END_REF][START_REF] Huang | Multifunctional metal-organic framework catalysts: synergistic catalysis and tandem reactions[END_REF][START_REF] Herbst | MOF catalysts in biomass upgrading towards value-added fine chemicals[END_REF][START_REF] Bagheri | A MoO3-metal-organic framework composite as a simultaneous photocatalyst and catalyst in the PODS process of light oil[END_REF], gas storage and separation [START_REF] Li | Recent advances in gas storage and separation using metal-organic frameworks[END_REF][START_REF] Li | Selective gas adsorption and separation in metal-organic frameworks[END_REF] drug delivery [START_REF] Horcajada | Metal-organic frameworks in biomedicine[END_REF], sensing [START_REF] Bagheri | A MoO3-metal-organic framework composite as a simultaneous photocatalyst and catalyst in the PODS process of light oil[END_REF][START_REF] Gamonal | Divergent adsorption-dependent luminescence of amino-functionalized lanthanide metal-organic frameworks for highly sensitive NO2 sensors[END_REF] , among others [START_REF] Yuan | Stable metal-organic frameworks: Design, synthesis, and applications[END_REF]. The term 'MOF' has been used in the literature since 1995 [START_REF] Yaghi | Hydrothermal synthesis of a metal-organic framework containing large rectangular channels[END_REF] to refer to a broad class of porous coordination polymers (PCPs) and porous coordination networks (PCNs) but it is only in 2012 that IUPAC suggested the following definition:" MOF is a Coordination Polymer (or alternatively Coordination Network) with open framework containing potential voids" [START_REF] Batten | Coordination polymers, metal-organic frameworks and the need for terminology guidelines[END_REF].

In this thesis, we are interested in two aspects of MOFs: the first is the application of MOFs in gas adsorption and separation; the second, which is motivated by the first, is the influence of the electronic structure on the adsorption properties.

MOFs for gas capture and storage

Efficient gas adsorption and separation can be possibly achieved in MOFs by exploiting their extremely large internal surface area and the presence of active metal sites with high affinity for gas molecules. The internal surface area of MOFs, generally evaluated applying [62].

the Brunauer, Emmett, and Teller (BET) theory to N 2 adsorption isotherms (measured at 77 K) [START_REF] Walton | Applicability of the BET method for determining surface areas of microporous metal-organic frameworks[END_REF], is of the order of thousands of meters square per gram. Hence, the gas can easily flow within the framework and interact, eventually, with the metal sites. The interaction between the gas molecules and the MOF can be largely enhanced by the presence of coordinatively unsaturated metal centers, called open-metal sites (OMSs) [START_REF] Kökçam-Demir | Coordinatively unsaturated metal sites (open metal sites) in metal-organic frameworks: design and applications[END_REF] (see figure 3.2). In this case, the gas adsorption process can happen via the formation of a coordination bond between the adsorbed molecule and the active site, with a consequent enhanced affinity for different sorbate species compared to MOF without OMSs [START_REF] Yoon | Selective nitrogen capture by porous hybrid materials containing accessible transition metal ion sites[END_REF]. Most of the technological efforts towards the development of MOFs suitable for gas storage/separation have been made to enhance the adsorption properties of greenhouse gases (CO 2 ), harmful gases (CO, NH 3 , NO 2 and SO 2 ) and energy-related gases (H 2 and CH 4 ) [START_REF] Li | Recent advances in gas storage and separation using metal-organic frameworks[END_REF].

Gas adsorption in these materials can occur via two possible binding mechanisms between the gas (adsorbate) and the MOF (adsorbent). The adsorption process can happen via van der Waals interactions (physical adsorption) or through the formation of chemical bondings between the adsorbate and the adsorbent (chemisorption). The distinction between these two interactions, however, is not always easy to establish. When a gas is adsorbed onto an unloaded MOF at fixed temperature (T), then the amount of adsorbed gas per unit mass of the solid (Q) is a function only of the applied gas pressure (P).

The evaluation of Q as a function of P is called adsorption isotherm and it is often used to investigate the adsorption properties of a system. Most adsorption isotherms can be described by the Langmuir-type function, i.e: where Q m represents the saturation capacity, i.e. the maximum value of Q obtained when the adsorbent internal surface is completely covered with a monolayer of the adsorbate.

Q(P ) = Q m bP 1 + bP (3.1)
The b parameter is given by the b = k a /k b , with k a and k b desorption and adsorption rate constants. Q m and b are constant for a given system and two regimes of Q(P) can be identified:

• Low pressure regime (bP 1). The eq. 3.1 assumes the linear form Q = KP (with K constant) and the adsorption isotherm describes a low adsorbate concentration regime, in which each of the adsorbed molecules "sees" the adsorbent as if there were no other adsorbates participating in the adsorption process. This part of the Q(P) diagram is also called Henry region.

• High pressure regime (bP 1). The adsorption sites in the adsorbent are gradually saturated and the molecules have fewer available binding sites as the pressure increases. Q approaches asymptotically the saturation capacity Q m .

For hydrogen and methane capture, MOFs are engineered to obtain high uptake values Q at room temperature and low pressure, thus avoiding the use of expensive cryogenic or high-pressure techniques to store the gas [START_REF] Rowsell | Strategies for hydrogen storage in metal-organic frameworks[END_REF][START_REF] Mason | Evaluating metal-organic frameworks for natural gas storage[END_REF]. In applications such as carbon dioxide post-combustion capture in thermal power plants, the target CO 2 adsorbate represents only the 14-16% of the whole flue gas that is mainly composed of N 2 (73-77 %). In this case, high uptake values have to be achieved along with selective adsorption of CO 2 over N 2 [START_REF] Hu | CO2 capture in metal-organic framework adsorbents: An engineering perspective[END_REF]. Selectivity properties are also desired in CO purification processes where the Metal-organic frameworks for efficient gas separation carbon monoxide molecule, produced in a variety of chemical processes [START_REF] Dutta | Developments in CO separation[END_REF], has to be separated from other chemical products like H 2 , CO 2 , N 2 and hydrocarbons [START_REF] Cheng | A review of PEM hydrogen fuel cell contamination: Impacts, mechanisms, and mitigation[END_REF].

I give now two examples of how the same MOF structure, in this case MOF-74, can be modified in order to tune its adsorption and selectivity properties.

• After the molecule is adsorbed, it is desirable to efficiently remove the adsorbate and regenerate the adsorbent so that it can be reused several times. An efficient adsorption/desorption cycle is of crucial importance, for example, in post-combustion carbon capture and storage (CCS) processes in thermal power plants, where the CO 2 capture process is systematically followed by the regeneration of the adsorbent and the storage of the carbon dioxide [START_REF] Bhown | Analysis and status of post-combustion carbon dioxide capture technologies[END_REF].

3.2 MOFs for gas capture and storage

25

In a temperature swing adsorption/desorption process, a gas (mixture) is first put in contact with the adsorbent at low pressure and low temperature. The resulting gas uptake depends on the pressure and temperature as explained before. The subsequent desorption is triggered through an increase in temperature that brings the system to a different adsorption isotherm characterized by a lower gas uptake and higher gas pressure. The difference between the adsorbate uptake at low temperature (adsorption) and high temperature (desorption) corresponds to actual working capacity (see figure 3.4).

For an adsorption/desorption strategy to be technologically relevant, the working ca- The difference between the gas uptake at (T low ,P ads ) and at (T high ,P des ) defines the working capacity. A smaller working capacity is obtained if the system is heated at T med <T high . Figure modified from ref. [START_REF] Reed | A spin transition mechanism for cooperative adsorption in metal-organic frameworks[END_REF].

pacity has to be maximized while keeping the temperature swing as small as possible.

The temperature variation is in fact directly linked with the cost for the implementation of the specific adsorption/desorption scheme. In the context of CCS in thermal power plants this cost is expressed by the energy penalty, i.e. the fraction of fuel that must be dedicated to CCS activities for a given quantity of input fuel used in the thermal power plant for energy production [START_REF] House | The energy penalty of post-combustion CO2 capture & storage and its implications for retrofitting the u.s. installed base[END_REF]. In 'phase-change' adsorbents the gas pressure can induce a structural modification that can greatly enhance the efficiency of the adsorption/desorption process achieving a large working capacity with a small energy penalty [START_REF] Mcdonald | Cooperative insertion of CO2 in diamine-appended metal-organic frameworks[END_REF]. In this situation, the Lagmuir-type isotherms do not correctly describe This work demonstrated how adsorption strategies based on the modification of the MOF electronic configuration upon gas adsorption can yield to the development of new technological solutions for efficient gas capture and storage. In the next section, I will present the main mechanisms through which gas uptake can be modified by a change in the electronic configuration of the open-metal site.

Q(P)

Role of the electronic configuration at the OMS on the gas adsorption mechanism

The interaction between gas and MOFs can happen through van der Waals forces, shortrange and long-range Coulomb mechanisms, exchange forces, and orbital rehybridization [4]. The interplay of these effects makes it difficult to rigorously separate each contribution to the binding energy. However, the analysis of the MOs diagram resulting from the hybridization between the metal, the molecule and the ligand orbitals is particularly relevant to establish a connection between the electronic configuration of the metal and the binding energy of the gas molecule. Therefore, I will only discuss this aspect here.

The covalent bond formation between the metal and the molecule can be understood by using MO diagrams between couples of valence atomic orbitals. The interaction of two orbitals A and B with non-zero overlap S results in the formation of two new MOs, one bonding AB and one anti-bonding AB * orbital. The first is lower in energy and the second Metal-organic frameworks for efficient gas separation is higher in energy than the two original A and B orbitals. The stabilization (destabilization) of the AB (AB * ) is proportional to S 2 /δ [START_REF] Jean | Molecular Orbitals of Transition Metal Complexes[END_REF], where δ is the energy difference between A and B. A strong interaction between orbitals thus occurs when the two interacting levels are close in energy and they have a big overlap. Two main interactions characterize the transition-metal complexes, the σ and the π interactions. In σ-bonding, there is axial overlap, S σ , while in π-bonding there is lateral (or sideways) overlap, S π . In general, S σ >S π and σ interactions are stronger compared to the π ones. When an empty (or singly occupied) metal e g orbital interacts with a doubly occupied orbital of the ligand (which can be the adsorbate gas molecule), the resulting σ interaction is maximal if the bonding σ orbital is doubly occupied while the anti-bonding σ * is empty. The occupation of σ bonding orbitals produces an electronic charge transfer from the molecule to the metal because of the hybridization. In the case of π interaction, there are two possible scenarios depending on the electronic occupation of the considered molecular orbitals as illustrated in figure 3.6:

• The metal (M) orbital is empty (or singly occupied) and the ligand (L) orbital is doubly occupied. This case is similar to the σ-bonding case and the occupation of the bonding π-bonding MO strengthens the metal-molecule σ-bond and electrons are transferred from the molecule to the metal (see figure 3.6, left).

• The M orbital is occupied and the L orbital is empty (typically the case of π * MOs of the ligand such as CO); the two orbitals are close in energy and they have a nonzero overlap S π . A π-backbonding interaction between the metal and the molecule can occur. This interaction strengthens the metal-molecule bond and it results in an electronic charge transfer from the M to the L (thus called π-backbonding or π-backdonation).

The MOs diagrams are a powerful tool to understand the trend of binding energies of adsorbates in MOFs. I present below two particularly significant cases, taken from the literature, where the MOs analysis has been used to rationalize the trends in binding energy.

• σ interaction. In a computational study from 2014, the authors of ref. • π-backbonding. Lee et al. [2] found that the presence of large π-backbonding interaction between N 2 and the V OMS in V-MOF-74 is at the origin of the observed strong N 2 /CH 4 selectivity of this MOF [START_REF] Wu | High-capacity methane storage in metal-organic frameworks m2(dhtp): The important role of open metal sites[END_REF][START_REF] Valenzano | Computational and experimental studies on the adsorption of CO, n2, and CO2 on mg-MOF-74[END_REF] (see figure 3.8). The LUMO of N 2 has π-like symmetry. When the overlap S π with the π-like metal orbital is non-zero then a π-backbonding interaction occurs. On the contrary, methane does not have an appropriate π-like empty orbital and S π is zero with a consequent weaker interaction with the V metal site. When the vanadium open-metal site is replaced with an HS Fe Metal-organic frameworks for efficient gas separation Figures reproduced from ref. [2].

By explaining these two specific cases in more detail, my intent is to set the ground for a better understanding of the results discussed in chapter 7 where a gas capture/release 3.3 Role of the electronic configuration at the OMS on the gas adsorption mechanism 31 SCO based strategy is presented. For this task we consider the Hofmann-type coordination polymer Fe(bpac)Fe(CN) 4 ( bpac=bis(4-pyridyl)acetylene)[79] (see figure 3.9) that represents ideal candidate for our purposes for the following reasons:

1. Presence of open-metal sites. On the plane, the structure of this MOF presents a series of Fe cyano-bridged metal sites with an alternate octahedral and squareplanar coordination environment. The square-planar coordinated Fe OMS is the target adsorption site.

2. Large pores for gas adsorption and diffusion. The 2D metal-cyanide sheets are connected by long bpac ligand along the z-direction (see figure 3.9).

3. Temperature-driven SCO transition near room temperature and with small hysteresis. A temperature-induced spin transition of the closed Fe(II) site has been experimentally observed with transition temperatures ranging between 225 K and 325 K according to the nature of the OMS and on the experimental conditions. The reported SCO curves are characterized by small hysteresis.

Once these properties are accomplished, the SCO based adsorption/desorption strategy is predicted to be able to achieve a large working capacity with a small energy penalty. The former is because of the expected step-shaped isotherms due to the spin-induced change in the gas-metal binding mechanism. The small energy penalty derives from the narrow temperature swing around room temperature needed to induce the spin-transition and that would require a limited amount of thermal energy to be triggered.

The working hypothesis of this work is that by achieving a temperature-induced SCO on the metal-site where the molecule is adsorbed, the energy penalty associated with the regeneration process is reduced owing to the decrease in binding energy thus potentially leading to more cost-effective technologies. Chapter 4

Computational Methods

Introduction

For a given Hamiltonian Ĥ containing all the relevant information about the system under investigation, the ground-state wavefunction Ψ(x) and its ground-state total energy E can be obtained solving the time-independent nonrelativistic Schrödinger equation:

Ĥ(x)Ψ(x) = EΨ(x), (4.1)
where x is a general position variable. For our purposes, we need to evaluate E from first principles (ab initio) and in this chapter we review some methodologies that allow us to do this in an approximate way.

Eq. 4.1 can be written equivalently as:

E = Ψ * (x) ĤΨ(x)dx Ψ * (x)Ψ(x)dx (4.2)
and the following inequality holds

E ≤ Ψ * (x) ĤΨ(x)dx Ψ * (x)Ψ(x)dx . (4.3)
Accordingly, if the total energy can be uniquely expressed as a functional of some suitable function X as E[X], then it is possible to obtain a good approximation of the ground-state total energy E by looking for the X that gives the lowest values of E[X]. This approach is called variational method and it is an extremely powerful tool to obtain approximate ground-state solutions of the time-independent Schrödinger equation. Despite that the variational approach is widely used in electronic structure methods, this does not represent the only possible approach to calculate the total energy of many-electron systems. Other approaches, based for example on Green-function formalism [START_REF] Galitskii | The energy spectrum of a non-ideal fermi gas[END_REF], exist but they will not be discussed in the next sessions of this thesis.

Born-Oppenheimer approximation

In most ab initio approaches, the underlying assumption is that the electronic degrees of freedom can be treated separately from the ions. Such approximation is called the Born-Oppenheimer approximation [START_REF] Born | Zur quantentheorie der molekeln[END_REF] and it is a consequence of the big mass difference between the electrons and the nuclei. Being much lighter, the electrons move in the solid faster than the ions and they relax rapidly to their ground-state for a given instantaneous position of the ions. It follows that the total Hilbert space of the system can be separated into the electronic Hilbert space and the ionic Hilbert space, consequently the total system wavefunction Ψ(R, r) can be written as:

Ψ(R, r) = Φ(R)ψ R (r), (4.4)
where Φ and ψ R are the ionic and electonic wavefunctions respectively. R = {R I } are the nuclear positions of the N ions ions and r = {r i } are the N electron positions. Note that the electronic wavefunction ψ R depends parametrically on the ion positions R and, moreover, it depends on the electronic spin degree of freedom that, for the sake of clarity, is not made explicit in the notation. As a consequence of eq. 4.4 the ionic and electronic wavefunctions can be obtained separately from eq. 4.1. The ionic wavefunction Φ(R) is the solution of the Schrödinger equation:

- I h2 2M I ∂ 2 ∂R 2 I + E(R) Φ(R) = Φ(R) (4.5)
where M I is the mass of the ion I and E(R) is the potential energy surface (PES) that represents the eigenvalue of the electronic Schrödinger equation solved for the ionic positions A first, rough, variational solution of eq. 4.6 can be obtained with HF method in which 4.3 Hartree-Fock method 35 the electrons are treated in a single-particle picture and they interact with each other through an effective potential defined by the configuration of the surrounding particles.

R, i.e.   - i h2 2m ∂ 2 ∂r 2 i + e 2 2 i =j 1 |r i -r j | - iI Z I e 2 |r i -R I | + e 2 2 I =J Z I Z J R I -R J   ψ R (r) = E(R)ψ R (r).

Hartree-Fock method

In HF method [82], the N-electrons wavefunction ψ R (x) is written as a single Slaterdeterminant:

ψ HF ({x i }) = 1 √ N ! φ 1 (x 1 ) φ 2 (x 1 ) • • • φ N (x 1 ) φ 1 (x 2 ) φ 2 (x 2 ) • • • φ N (x 2 ) . . . . . . . . . . . . φ 1 (x N ) φ 2 (x N ) • • • φ N (x N ) . (4.7)
The {φ i {x j }} are single-particle orbitals and x j are generalized coordinates that include spatial as well as spin degrees of freedom. The HF wavefunction ψ HF still depends on the ionic positions R but we drop the subscript for clarity. To show how the time-independent electronic Schrödinger equation reads for a wavefunction in the form of eq. 4.7 we write the electronic Hamiltonian operator as the sum of single-particle operators ĥ1 and pair interaction operators ĥ2 :

Ĥe = i ĥ1 (x i ) + 1 2 i =j ĥ2 (x i , x j ), (4.8)
where for fixed spin:

ĥ1 (r i ) = - h2 2m ∂ 2 ∂r 2 i - I Z I e 2 |r i -R I | (4.9) and ĥ2 (r i , r j ) = e 2 |r i -r j | . (4.10)
Notice that in the electronic Hamiltonian Ĥe the last left hand term in eq. 4.6 is not included. Within the Born-Oppenheimer approximation, this term contributes to the total energy with a constant value E ion due to the Coulomb interactions between the ions.

The Hartree-fock energy E HF can be written as:

E HF = ψ HF | Ĥe |ψ HF + E ion . (4.11)
The single electron operator contribution to the total energy is just the sum over the total number of electrons of the expectation value of ĥ1 calculated on the single-particle wavefunctions |φ i :

ψ HF | i ĥ1 (x i ) |ψ HF = i φ i | ĥ1 (x i ) |φ i . (4.12)
Computational Methods

The two-particle energy term is:

ψ HF | 1 2 i =j ĥ2 (x i , x j ) |ψ HF = 1 2 i,j φ i φ j | ĥ2 (x i , x j ) |φ i φ j -φ j φ i | ĥ2 (x i , x j ) |φ i φ j , (4.13)
where the first right hand term is the Coulomb interaction between pairs of electrons and the second right hand term is the exchange term. It is important to notice that for i = j the two terms cancel out. This is going to be an important difference between HF and density-functional theory. In the HF method, each electron feels the presence of the other electrons through an effective potential and each |φ i can be determined upon knowing all the others N -1 single-electron wavefunctions |φ j through Ĥe . HF is a mean-field theory for electrons and the HF wavefunction |ψ HF has to be computed self-consistently, i.e. using the variational principle in eq. 4.3 by minimizing the total energy with respect to some suitable set of parameters. These parameters are the coefficients that enter the expression of the single-particle wavefunctions {φ i }.

Basis set approximation

The possibility of expanding the single-particle wavefunctions {φ i } as a finite linear combination of known basis functions is at the core of any ab initio method (that use classical bits). A finite basis expansion is an approximation to the complete basis set (infinite) needed to express a quantum-mechanical wavefunction. The choice of the basis set depends on the properties investigated and the required accuracy; the bigger the number of basis functions the more computation expensive is the calculation. For a given basis set {χ i }, each single-particle orbital (often called also one-electron MO) φ i (r) can be expressed as:

φ i (r) = j b ij χ j (r). (4.14)
The basis function {χ i } can be expressed in different ways (plane waves, localized functions, augmented plane waves,...) [START_REF] Martin | Electronic Structure[END_REF]. In the following, the expressions of {χ i } constructed through three different methods are shown as an example:

• Slater-type orbital (STO):

χ i (r, θ, θ ) = N ST O r n-1 e -α i r Y l,m (θ, θ ). ( 4

.15)

• Gaussian-type orbital (GTO):

χ i (r, θ, θ ) = N GT O r 2n-2-l e -α i r 2 Y l,m (θ, θ ). (4.16)

Full Configuration Interaction wavefunction

• Plane-wave basis:

χ i (r) = N P W e -G i •r . (4.17)
where N ST O , N GT O , and N P W are normalized constants. STO and GTO types are expressed in polar coordinates (r, θ, θ ), Y l,m is a spherical harmonic function and α i is the orbital constant. In plane-wave basis set, G i is the reciprocal lattice vector.

STO and GTO basis functions are generally (but not always) used in calculations where periodic boundary conditions (PBC) are not imposed. On the contrary, when periodic systems are considered, the electronic wavefunctions in a band theory picture are more efficiently described as a superposition of plane-waves (free electrons).

Full Configuration Interaction wavefunction

The HF theory can capture almost 99% of the total electronic energy of a system, nevertheless, the remaining 1% contribution due to the many-body nature of the electronic wavefunction is often crucial when real systems are investigated. As already discussed in section 2.5, in quantum chemistry the correlation energy is defined as the difference between the total exact total energy E and the HF energy. To fully recover the correlation effects, it is necessary to move from the HF single-particle theory to a proper many-body description of the electronic wavefunction. This can be achieved through the full-CI expansion of the wavefunction.

Let's consider the N -electrons HF single Slater-determinant wavefunction |ψ HF built with M single-particle wavefunctions {φ i } (M > N ). There are N occupied states and M -N unoccupied or virtual states. To facilitate the discussion, |ψ HF is now expressed writing explicitely the occupied single-particle orbitals as:

|ψ HF = |φ a φ b φ c . . . φ N . (4.18)
A single excitation in which an electron is promoted from an occupied |φ a to a virtual |φ r MO can be written as:

|ψ r a = |φ r φ b φ c . . . φ N . (4.19)
In the same way, a double excitation in which two electrons are excited from their initial |φ a and |φ b single-particle orbitals to the |φ r and |φ s virtual orbitals is:

|ψ rs ab = |φ r φ s φ c . . . φ N (4.20)
Computational Methods and similarly, higher excitation orders can be defined. At each excitation we can associate an excited Slater-determinant of the form |ψ r,s,... a,b,... . The total number of possible excited determinants depends on the number of electrons N and the number of virtual orbitals M-N. We can define the CI operator ĈI such that its effect on a given Hartree-Fock wavefunction |ψ HF is that it generates all the I th excited Slater-determinants. The full-CI wavefunction is then expressed as a linear combination of all the possible excited determinants: The full-CI wavefunction is reduced in size imposing that |Ψ CI is an eigenstate of the angular momentum and spin operators. After this symmetry operation, the remaining terms in 4.21 are called configuration state functions (CSFs) and the total CI ground-state energy E CI {c r,s,... a,b,... } is variationally computed optimizing the CI coefficients {c r,s,... a,b,... }. The full-CI approach can recover all the correlation energy within the selected basis set.

|Ψ CI = 1 + N I=1 ĈI |ψ HF = |ψ HF + N I=1 1 (I!) 2

Coupled cluster method

In CC theory, the electronic wavefunction |Ψ CC is expressed as an exponential ansatz:

|Ψ CC = e T |ψ HF (4.22)
where T is the cluster operator which can be written as the sum of cluster operators of different excitation levels:

T = T1 + T2 + • • • + TN . (4.23)
When all the terms in 4.23 are considered, the cluster operator e T is equivalent to the full-CI operator 1 + N I=1 ĈI of eq. 4.21, i.e.:

e T = 1 + T + T 2 2! + T 3 3! + . . . = 1 + Ĉ1 + Ĉ2 + • • • + ĈN (4.24)
and the relationships with the CI operators are:

Ĉ1 = T1 , (4.25) Ĉ2 = T2 + 1 2 T 2 1 , (4.26)
Ĉ3 = T3 + T1 T2 + 1 6 T 3 1 , (4.27)
4.7 Complete active space self-consistent field theory 39 and so on. This means that single excitations are described by T1 . Double excitations can be obtained in two ways: with a simultaneous two electrons excitation T2 or with two independent single excitations T 2 1 . The triple excitations can be generated as three simultaneous excitations T3 , three independents excitations T 3 1 or as an independent excitation and a double excitation T1 T2 . It is important to note that for a complete description of the ĈI operators all the cluster operator terms, Ti , have to be considered and, for example, the double independent excitations T 2 1 can capture only some of the double exitations Ĉ2 . Instead of considering the full cluster operator it is common practice to keep only few orders, for example truncating the expansion at the double excitation operator T2 we obtain the CCSD theory (where SD stands for single and double excitations). The main advantage in using the exponential anzats derives from the fact that even at CCSD level of theory also some triple, quadruple and higher excitations are treated. This is clear if we write the CCSD cluster operator e T1 + T2 using eq. 4.24 , i.e.:

e T1 + T2 = 1 + T1 + T2 + T 2 1 2 + T 2 2 2 + T1 T2 + . . . . (4.28)
With the inclusion of many excited Slater-determinants, the coupled cluster method, especially in its CCSD(T) formulation (where the triple excitations are included perturbatively), is able to capture most of the electrons dynamical correlation and it is often referred to the "gold standard" of quantum chemistry [START_REF] Watts | Coupled-cluster methods with noniterative triple excitations for restricted open-shell hartree-fock and other general single determinant reference functions. energies and analytical gradients[END_REF][START_REF] Knowles | Coupled cluster theory for high spin, open shell reference wave functions[END_REF].

When systems with significant degeneracy in the ground-state electronic configuration are treated with CCSD(T) theory, some important contributions to the correlation energy may not be capture (static correlation). In these cases a single Slater-determinant employed for the reference wavefunction cannot describe these electronic correlations and a combination of multiple Slater determinants should be used instead.

Complete active space self-consistent field theory

Systems that are not well described in a single Slater-determinant picture can be better represented using a linear combination of different Slater-determinants. This is the idea behind the multiconfigurational self-consistent field (MCSCF) theory in which the electronic wavefunction |Ψ M CSCF is initially written as a linear combination of CSFs, i.e:

|Ψ M CSCF = i c i |Ψ CSF i . (4.29)
|Ψ M CSCF is then parametrized through the CSFs CI coefficients {c i } plus the singleparticle MO coefficients of eq. 4.14, i.e. the set

{b i }. The total MCSCF energy E M CSCF Computational Methods is obtained by minimizing E M CSCF [{c i }, {b i }].
The main drawback of MCSCF is the choice of the CSFs to include in |Ψ M CSCF . When selecting the set of CSFs to include in the calculation, one should consider the specific chemical situation under investigation as well as the fact that the calculation becomes computationally prohibitive as the number of CSFs increases. A strategy to facilitate the choice of the CSFs to include in an MCSCF calculation has been developed by Roos et al. [START_REF] Roos | A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach[END_REF]. The method avoids the manual selection of the relevant CSFs which are, instead, generated with a full-CI expansion within a selected subset of MOs. The single-particle orbitals {φ i } used to construct the CSFs are divided into inactive (doubly occupied), active (occupation determined through the CI expansion) and virtual orbitals (unoccupied). The active orbitals can be empty, singly, or doubly occupied and are selected by hand according to the specific chemical situation. The MCSCF method where the CSFs are selected through this procedure, i.e. by embedding full-CI in the HF mean field, is called complete active space self-consistent field (CASSCF).

For the calculation of excited electronic states, CASSCF can be applied individually on a specific electronic state or, alternatively, by considering N SA electronic states within a given symmetry and performing a state-average CASSCF (SA-CASSCF). In the latter, a single set of MOs is optimized for the N SA states while the CI coefficients are optimized on each state independently. The minimization process is applied on the weighted sum of each state energy. More reliable results are found by weighing equally each of the N SA CASSCF state {|Ψ i CASSCF }. When a single electronic state is considered, SA-CASSCF simply describes a single-state (SS) CASSCF method. CASSCF method can efficiently capture static correlation. On the other hand, dynamical correlations effects need a large number of CSFs to be correctly treated and they converge much more slowly with the size of the active space compared to the static correlations. Thus a more efficient way to include dynamical correlations is by performing complete active space second-order perturbation (CASPT2) theory in which dynamical correlation effects are perturbatively added on the converged CASSCF wavefunction.

Complete active space second-order perturbation theory

CASPT2 method [START_REF] Andersson | Second-order perturbation theory with a complete active space self-consistent field reference function[END_REF] can be used on top of a CASSCF wave function |Ψ CASSCF to obtain a new wavefunction |Ψ CASP T 2 , with relative energy E CASP T 2 , where the dynamical effects are recovered perturbatively in a way similar to that used in the Møller and Plesset 4.8 Complete active space second-order perturbation theory method [START_REF] Møller | Note on an approximation treatment for many-electron systems[END_REF].

In perturbation theory, the Hamiltonian operator Ĥ is written as the sum of the unperturbed Hamiltonian Ĥ0 plus a small perturbation V :

Ĥ |Ψ = ( Ĥ0 + λ V ) |Ψ = E |Ψ . (4.30)
Where λ is a dimensionless parameter ranging continuously from 0 (no perturbation) to 1 (full perturbation). Since the perturbation is small, the wavefunction |Ψ and the energy E can be expanded as

|Ψ = |Ψ (0) + λ |Ψ (1) + λ 2 |Ψ (2) + . . . (4.31)
and

E = E (0) + λE (1) + λ 2 E (2) + . . . . (4.32)
The Schrödinger equation 4.30 combined with eq. 4.31-4.32 leads to the following expressions for the first three energy orders

E (0) = Ψ (0) | Ĥ0 |Ψ (0) (4.33)
E (1) = Ψ (0) | V |Ψ (0) (4.34)
E (2) = Ψ (0) | V |Ψ (1) (4.35)
Møller and Plesset showed that when the unperturbed Hamiltonian Ĥ0 is considered to be the sum of one-electron Fock operators, the first order correction to the energy, E (1) , is already included in the Hartree-Fock energy and so, only the second order correction E (2)

has to be calculated in a second order perturbation theory.

In state-specific, or single-state, (SS) CASPT2, the second order perturbation theory is applied considering as reference wavefunction a single SA-CASSCF wavefunction |Ψ i CASSCF describing one electronic state, i, or a SS-CASSCF wavefunction. The relative CASPT2 wavefunction |Ψ i CASP T 2 is written as

|Ψ i CASP T 2 = |Ψ i CASSCF + |Ψ i (1)
(4.36)

where |Ψ i (1) is the first order contribution to |Ψ i CASP T 2 . The second order energy contribution reads:

ii = Ψ i CASSCF | V |Ψ i (1) . (4.37)
Once the second order correction energy ii is computed, the SS-CASPT2 energy for the electronic state i is:

E i SS-CASP T 2 = E i CASSCF + ii (4.38)
Computational Methods An extension of CASPT2 has been suggested by Finley and collaborators [START_REF] Finley | The multi-state CASPT2 method[END_REF], where several SA-CASSCF functions are considered to represent a multi-dimensional reference space in multi-state (MS) CASPT2. In MS-CASPT2 the zeroth-order wavefunction |Ψ (0)

(also called perturbation modified CAS (PMCAS) reference function) is expressed as a linear combination of SA-CASSCF reference wavefunctions. The effective Hamiltonian matrix elements are

H ef f ij = E i CASSCF δ ij + 1 2 ( ij + ji ) (4.39)
where the diagonal terms elements H ef f ii are the SS-CASPT2 energies for the state i, while the off-diagonal elements couple the electronic states i and j. The total MS-CASPT2 wavefunction for the electronic state i is

|Ψ i M S-CASP T 2 = a C ia |Ψ a CASSCF + |Ψ i (1)
(4.40)

and the coefficients C ia are obtained diagonalizing the effective Hamiltonian in eq. 4.39.

MS-CASPT2 (as well as its extended formulation XMS-CASPT2 [START_REF] Granovsky | Extended multi-configuration quasi-degenerate perturbation theory: The new approach to multi-state multi-reference perturbation theory[END_REF]) better captures dynamical correlation compared with SS-CASPT2 when different SA-CASSCF states are mixed together and a single SA-CASSCF wavefunction does not represent a good reference.

Density-functional theory

In DFT the total electronic ground-state energy, E, of an electronic system is written as a functional of the total electron density ρ, i.e. ∼ E[ρ]. The total electronic charge density ρ(r) associated with a generic N-electrons ground-state wavefunction ψ R ({r i }) can be written as:

ρ(r) = N |ψ R (r 1 , r 2 , . . . , r N )| 2 dr 1 . . . dr N . (4.41)
The electronic-charge density is non negative, continuous and ρ(r)dr = integer. The practical advantage in expressing the ground state energy as a functional of the charge density (that depends only on one spatial variable) is the significant reduction of the problem degrees of freedom with respect to expressing the ground-state energy using the full many-body wavefunction (it depends on the spatial variables of all the electrons).

Formally, the basis of the DFT is expressed by the Hohenber-Kohn theorems [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] which establish a one-to-one relationship between the ground-state electronic charge density ρ(r)

and the external potential V ext (r) acting on the electrons, and prove the existence of a universal functional F [ρ(r)], independent of V ext (r), such that the energy functional

E[ρ(r)] = F [ρ(r)] + V ext (r)ρ(r)dr (4.42)
has its minimum at the correct ground-state energy. Nondegeneracy is required to ensure the uniqueness of the ground-state charge density ρ(r) associeted with the ground-state wavefunction ψ R (r).

The universal functional F [ρ(r)] can be decomposed as the sum of the kinetic functional T [ρ(r)], the "classical" Coulomb interaction and the exchange-correlation term E xc [ρ(r)],

i.e:

F [ρ(r)] = T [ρ(r)] + 1 2 ρ(r)ρ(r ) |r -r | drdr + E xc [ρ(r)]. ( 4
ĤKS ψ KS i (r) = - h2 ∇ 2 2m + V KS (r) ψ KS i (r) = i ψ KS i (r) (4.44)
and the ground-state electronic charge density can be written as:

ρ(r) = i f i |ψ KS i (r)| 2 (4.45)
where f i are the orbital occupation numbers. Within the Kohn-Sham formulation of DFT, it is straightforward to calculate the kinetic energy for the non-interacting electrons as:

T s [ρ(r)] = - i f i (ψ KS i (r)) * h2 ∇ 2 2m ψ KS i (r)dr. (4.46)
The only unknown term left in the expression of F [ρ(r)] is the exchange-correlation energy functional E xc [ρ(r)], which enters eq. 4.44 through the exchange-correlation potential v xc (r) in the effective KS potential, i.e.: where ε hom xc (ρ(r)) is the exchange-correlation energy density per particle. The exchange energy for the homogeneous electron gas has an analytic expression [START_REF] Dirac | Note on exchange phenomena in the thomas atom[END_REF] while the correlation energy has been calculated with great accuracy with Monte Carlo methods [START_REF] Ceperley | Ground state of the electron gas by a stochastic method[END_REF].

V KS (r) = V ext (r) + e 2 ρ(r ) |r -r | dr + v xc (r) (4.47) v xc (r) = δE xc [ρ(r)] δρ(r) . ( 4
The LDA approximation thus holds especially in systems where the charge density varies smoothly in the space. It tends to predict shorter equilibrium bond distances and higher bonding energies when compared to the experiments.

• Generalized Gradient Approximation (GGA). To improve over the LDA, the exchange-correlation functional can be derived as a function not only of the local electronic charge density (as in LDA) but also of the local density gradient ∇ρ(r) (gradient approximation). In GGA the exchange-correlation functional reads:

E GGA xc [ρ] = ε GGA xc (ρ(r), ∇ρ(r))ρ(r)dr. (4.50)
Accordingly to the expression of ε GGA xc several GGA functionals have been proposed, e.g. BP86 [START_REF] Perdew | Density-functional approximation for the correlation energy of the inhomogeneous electron gas[END_REF], BLYP [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the colle-salvetti correlationenergy formula into a functional of the electron density[END_REF], PW91 [START_REF] Perdew | Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation[END_REF] and PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] among others. GGA functionals provide a substantial improvement over LDA for, almost, the same computational cost.

• Meta-Gradient Approximation. In meta-gradient approximation functionals (meta-GGA) the exchange-correlation functional depends on the kinetic energy density τ (r) and/or the Laplacian of the density ∇ 2 ρ(r) in addition to the local density ρ(r) and the density gradient ∇ρ(r), i.e.:

E meta-GGA xc [ρ] = ε meta-GGA xc ρ(r), ∇ρ(r), ∇ 2 ρ(r), τ (r) ρ(r)dr (4.51)
where the kinetic energy density is:

τ (r) = 1 2 i f i |∇ψ KS i | 2 . (4.52)
Meta-GGA functionals, like TPSS [START_REF] Tao | Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids[END_REF] and M06-L [START_REF] Zhao | A new local density functional for main-group thermochemistry, transition metal bonding, thermochemical kinetics, and noncovalent interactions[END_REF], are non-local functionals of the density in contrast with the local and semilocal density dependency of LDA and 4.9 Density-functional theory

45

GGA functionals [START_REF] Perdew | Jacob's ladder of density functional approximations for the exchangecorrelation energy[END_REF][START_REF] Perdew | Accurate density functional with correct formal properties: A step beyond the generalized gradient approximation[END_REF]. Nevertheless, meta-GGA functionals depend semilocally on the KS orbitals and thus they are not much more computational demanding compared to GGA functionals.

The systematic improvements of GGA over LDA are not found when GGA and meta-GGA are compared. Bond lengths of molecules are in some cases worse in meta-GGA compared to GGA [START_REF] Adamo | The meta-GGA functional: Thermochemistry with a kinetic energy density dependent exchange-correlation functional[END_REF] but in other cases, for example, lattice constants for solids [START_REF] Kurth | Molecular and solid-state tests of density functional approximations: LSD, GGAs, and meta-GGAs[END_REF] or metal surface energies [START_REF] Yan | Density-functional versus wave-function methods: Toward a benchmark for the jellium surface energy[END_REF], meta-GGA functionals perform better.

• Hybrid functionals. Hybrid functionals result from the combination of one or more of the abovementioned exchange-correlation functionals (LDA,GGA, and meta-GGA) with some portion of the HF exchange [START_REF] Becke | A new mixing of hartree-fock and local density-functional theories[END_REF]. HF exchange E HF x (see eq. 4.13) is a non-local functional of the occupied orbitals and when it is mixed with others exchange-correlation functionals in DFT can significantly improve the description of several electronic properties [START_REF] Corà | Principles and applications of density functional theory in inorganic chemistry ii[END_REF][START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. One of the most widely used DFT hybrid functional [START_REF] Kepp | Consistent descriptions of metal-ligand bonds and spin-crossover in inorganic chemistry[END_REF] is B3LYP [START_REF] Lee | Development of the colle-salvetti correlationenergy formula into a functional of the electron density[END_REF][START_REF] Hertwig | On the parameterization of the local correlation functional. what is becke-3-LYP?[END_REF] that can be written as:

E B3LY P xc = (1 -a)E LDA xc + aE HF x + bE B88 x + cE LY P c + (1 -c)E V W N c (4.53)
where the B88 [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF] GGA exchange functional, the LYP [START_REF] Lee | Development of the colle-salvetti correlationenergy formula into a functional of the electron density[END_REF] GGA and the VWN [START_REF] Vosko | Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis[END_REF] LDA correlation functionals are used. The three parameters a = 0.20, b = 0.72, and c = 0.81 are obtained empirically. Other popular hybrid functionals are PBE0 [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF],

HSE [START_REF] Heyd | Hybrid functionals based on a screened coulomb potential[END_REF] ,M06 [START_REF] Zhao | The m06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four m06class functionals and 12 other functionals[END_REF] and TPSSh [START_REF] Staroverov | Comparative assessment of a new nonempirical density functional: Molecules and hydrogen-bonded complexes[END_REF]. Hybrid functionals are more computationally expensive compared to LDA, GGA, and meta-GGA functionals.

Because approximate KS DFT can capture part of the dynamical correlation that is completely missed in the HF theory, it globally improves over HF. Most of the inaccuracy of approximate KS DFT can be traced back to the so-called self-interaction error (SIE) [START_REF] Perdew | Self-interaction correction to density-functional approximations for many-electron systems[END_REF]. SIE arises because of the Coulomb term (for a singly occupied orbital) in the Coulomb energy expression does not cancel out with the corresponding term in approximate exchange-correlation functionals, as it happens in the HF theory. This results in an electronic charge density that feels an excessive Coulomb repulsion that pushes electrons apart, thus resulting in delocalization. This manifests itself in the incorrect behavior of the approximate functionals KS DFT total energy between integer number of electrons.

As it has been shown by Perdew et al. [START_REF] Perdew | Density-functional theory for fractional particle number: derivative discontinuities of the energy[END_REF], the trend of the exact electronic energy of a system calculated between integer number of electrons follows a straight line by virtue of the discrete nature of electrons. The delocalization error lowers the energy for fractional number of electrons resulting in a convex behavior of the energy (see figure 4.1). However, also HF energy deviates from the correct energy trend. In this case the lack of correlation results in an overlocalization of the electron charge density that occurs in a concave behavior of the energy.

For certain systems is of central importance to describe the localized character of the charge density, such as in strongly-correlated systems. For these, SIE may heavily affect the computed electron charge density and total energy and thus it can compromise the description of the ground-state character of the systems (e.g. Mott insulators) A possible improvement in this direction is represented by the DFT+U extension of DFT.

DFT+U approach

In strongly-correlated systems the correlations are mainly due to a specific electronic subshell of the system (for example d or f localized electrons) and the DFT energy functional can be shaped on the electronic Hubbard Hamiltonian which in its single-band formulation reads:

H Hub = t <i,j>,σ (c † i,σ c j,σ + h.c.) + U i n i,↑ n i,↓ (4.54)
where c † i,σ , c j,σ and n i,σ are electron creation, annihilation and number operators , respectively, < i, j > denotes nearest-neighbour atom sites, and σ =↑, ↓ the electron spin.

The Hubbard Hamiltonian describes electrons jumping from one site to the other with a probability t and subject to a Coulomb on-site repulsion U when they occupy the same atomic site.

In DFT+U the total energy functional assumes the form:

E DF T +U [ρ(r)] = E DF T [ρ(r)] + E Hub [{n Iσ mm }] -E dc [{n Iσ }] (4.55)
where E DF T is the energy functional of an arbitrary approximate DFT functional (generally LDA or GGA), E hub is the Hubbard term used to incorporate the Hubbard Hamiltonian model in eq. 4.54 into the DFT framework, and E dc is the double-counting term which avoids the inclusion of the description of the selected sub-shell twice because it is already treated in DFT via main field. To determine the electronic levels treated at the DFT+U level, a localized basis set, {φ I m }, describing the relevant localized electronic orbitals {m} of the considered atomic site, I, has to be defined (and the choice is not unique nor irrelevant). The elements n Iσ mm are the occupations of such a basis defined as the projections of the KS orbitals {ψ KS i,σ } on {φ I m }, i.e.:

n Iσ m,m = i f σ i ψ KS i,σ |φ I m φ I m |ψ KS i,σ (4.56)
where f σ i are the Fermi-Dirac occupations of the KS orbitals. The form of the Hubbard energy functional E U = E Hub -E dc has evolved during the years, from the pioneering works of Anisimov and coworkers [START_REF] Anisimov | Band theory and mott insulators: HubbardUinstead of StonerI[END_REF][START_REF] Anisimov | Density-functional theory and NiO photoemission spectra[END_REF][START_REF] Solovyev | Corrected atomic limit in the local-density approximation and the electronic structure ofdimpurities in rb[END_REF] to the simplified formulation proposed by Dudarev et al. [START_REF] Dudarev | Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+u study[END_REF] and based on the rotationally invariant DFT+U reformulation of Liechtenstein and Anisimov [START_REF] Liechtenstein | Density-functional theory and strong interactions: Orbital ordering in mott-hubbard insulators[END_REF]. Here I present the simplified version [START_REF] Dudarev | Electron-energy-loss spectra and the structural stability of nickel oxide: an LSDA+u study[END_REF] because 1) it represents the formulation adopted for all the DFT+U calculations of this thesis and 2) it allows one to implement a linear-response scheme to calculate the parameter U ab initio (vide infra) [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in theLDA+u method[END_REF].

In the simplified rotationally-invariat formulation, the Hubbard energy functional reads:

E U [{n Iσ mm }] = E hub [{n Iσ mm }] -E hub [{n Iσ mm }] = = I U I 2 (n I ) 2 - σ T r[(n Iσ ) 2 ] - I U I 2 n I (n I -1) = = I,σ U I 2 T r[n Iσ (1 -n Iσ )] (4.57)
where n Iσ is the occupation matrix whose elements are the occupations 4.56, n I = m,σ n Iσ m and {n Iσ m } are the diagonal elements of the occupation matrix, i.e. n Iσ m = n Iσ mm . The Hubbard term U I for the atomic site I entering eq. 4.57 represents an effective Hubbard potential U ef f = U I -J I , with J exchange interaction. We assume hereafter Computational Methods U I = U ef f . The KS potential V KS DF T +U in the DFT+U formulation takes the form:

V KS DF T +U |ψ KS i,σ = V KS |ψ KS i,σ + I,m U I 1 2 -n Iσ m |φ I m φ I m |ψ KS i,σ (4.58)
where V KS is the KS potential of the DFT approximation adopted to build the DFT+U functional. From eq. 4.58, the effect of the DFT+U potential (second term on the right side) on the KS orbitals becomes clear: when the occupation elements n Iσ m are greater then 1/2 the potential is attractive, while for n Iσ m < 1/2 it is repulsive. During the self-consistent cycle the electronic levels defined by the localized basis set {φ I m } are "emptied" or "filled" according to the relative occupation matrix elements {n Iσ m }. When self-consistency is reached, the total DFT+U energy is the sum of the DFT energy calculated on the converged DFT+U electronic density plus the Hubbard U energy E U [{n Iσ mm }}] which destabilizes the total energy proportionally to the non-zero occupation matrix elements, {n Iσ m }, and the Hubbard parameter U I . So far, the Hubbard parameter U entering DFT+U is treated as an external parameter and the specific value has to be determined for any practical purpose. It is common practice to extract U by fitting the experimental values of certain investigated properties with the same computed with DFT+U [START_REF] Wang | Electronic structure: Wide-band, narrowband, and strongly correlated systems-oxidation energies of transition metal oxides within the gga+ u framework[END_REF][START_REF] Lebègue | Modeling spin-crossover compounds by periodic dft+ u approach[END_REF][START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF]. Alternatively, U can be calculated in an ab initio fashion noticing that the Hubbard U energy in eq. 4.57 represents the contribution to the total energy curvature due to fractional values of the occupation number n of the localized sub-shell, i.e. U = δ 2 E DF T /δn 2 [START_REF] Gunnarsson | Density-functional calculation of the parameters in the anderson model: Application to mn in CdTe[END_REF][START_REF] Anisimov | Band theory and mott insulators: HubbardUinstead of StonerI[END_REF]. In this way the Hubbard correction is introduced to recover the linearity of the DFT energy between integer number of electrons mitigating the SIE effects on the electronic charge density. It is important to notice that the curvature shown in figure 4.1 arises because of the behavior of the total energy with respect to fractional values of the total electron number (global curvature), while U is the curvature because of fractional values of the selected localized sub-shell (local curvature).

Nevertheless, the Hubbard U correction has been shown to alleviate both local and global curvature [START_REF] Zhao | Global and local curvature in density functional theory[END_REF] (see also chapter 5 of this thesis). To show how to calculate U as the DFT local energy curvature I follow the formulation proposed by Cococcioni and de Gironcoli [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in theLDA+u method[END_REF] inspired by the work of Pickett and collaborators [START_REF] Pickett | Reformulation of theLDA+u method for a local-orbital basis[END_REF]. The on-site Hubbard U I for the atomic site I is written as:

U I = ∂ 2 E DF T ∂(n I ) 2 - ∂ 2 E DF T ∂(n I 0 ) 2 (4.59)
where the second right-hand term represents the noninteracting contribution to the local curvature arising from the rehybridization of the electronic states induced by a change in the occupation number n I 0 in the noninteracting system. Computing numerically eq. 4.59 by varying n I and collecting the relative variation of the DFT energy is not very practical.

A more efficient procedure consists in replacing the change in n I with a perturbation α I

to the total DFT KS potential. This perturbation shifts the KS orbitals having a nonzero overlap with the selected localizes basis set {φ I m } thus resulting in a change in the occupation numbers n I . The perturbed KS potential reads:

V KS α = V KS + α P I (4.60)
where P I projects onto the localized basis set {φ I m }:

P I = m |φ I m φ I m | (4.61)
The DFT minimization problem with the potential of eq. 4.60 can be seen as a constrained minimization problem in the form

E DF T [α I ] = min ρ {E DF T [ρ] + α I n I } (4.62)
and for a given value of the perturbation α I , the relative occupation n I can be obtained.

In order to switch the independent variable from α I to n I , and thus calculate

U I = ∂ 2 E DF T /∂(n I ) 2
, we use a Legendre trasformation to write:

E DF T [n I ] = min α I {E DF T [α I ] -α I n I }. (4.63)
Noticing that dE DF T [α I ]/dn I = 0 we can write:

dE DF T [n I ] dn I = -α I (4.64) d 2 E DF T [n I ] d(n I ) 2 = - dα I dn I (4.65)
Practically, the quantity that is varied "by hand" is the perturbation α I . By "collecting" the change in occupation n I , one can calculate the response functions as χ II = dn I /dα I and χ 0 II = dn I 0 /dα I . The latter gives the abovediscussed noninteracting curvature, ∂ 2 E DF T /∂(n I 0 ) 2 , and the occupations n I 0 are computed without allowing the electrons reorganize after the potential V KS α is applied on the unperturbed charge density. Finally the Hubbard U I can be written as:

U = (χ -1 0 -χ -1 ) II . ( 4

Main paper

The past decade has seen a substantial effort from the computational chemistry community in predicting the spin-state energetics of molecular complexes containing transition-metal atoms [START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF][START_REF] Domingo | Spin crossover in fe(ii) complexes: An ab initio study of ligand σ-donation[END_REF][START_REF] Radoń | Benchmarking quantum chemistry methods for spin-state energetics of iron complexes against quantitative experimental data[END_REF][START_REF] Cirera | Theoretical prediction of spin-crossover temperatures in ligand-driven light-induced spin change systems[END_REF]134,[START_REF] Swart | Accurate spin-state energies for iron complexes[END_REF][START_REF] Droghetti | Assessment of density functional theory for iron(ii) molecules across the spin-crossover transition[END_REF][START_REF] Cirera | Benchmarking density functional methods for calculation of state energies of first row spin-crossover molecules[END_REF]. Achieving such a description with chemical accuracy appears today as a significant challenge for ab initio methods with strong implications in future assistance and guidance of experimental efforts devoted to the (i) development and optimization of catalytic reactivity and the (ii) design of novel spin crossover (SCO) complexes. The vast majority of SCO materials, which are of great interest for applications such as molecular spintronics, molecular electronics, sensors and actuators [START_REF] Kumar | Emerging trends in spin crossover (sco) based functional materials and devices[END_REF][START_REF] Molnár | Molecular spin crossover materials: Review of the lattice dynamical properties[END_REF][START_REF] Resines-Urien | A switchable iron-based coordination polymer toward reversible acetonitrile electro-optical readout[END_REF][START_REF] Resines-Urien | Covalent post-synthetic modification of switchable iron-based coordination polymers by volatile organic compounds: a versatile strategy for selective sensor development[END_REF], are octahedrally-coordinated Fe(II)-based molecular complexes. For these complexes, the thermodynamics of the spin transition (i.e. the transition temperature T 1/2 ) between the low-spin (LS) and high-spin (HS) states is dominated by the adiabatic energy difference, ∆E H-L =E HS -E LS , i.e. the energy difference between the two spin states computed at their corresponding geometry. Within density functional theory, large deviations in the value of ∆E H-L are found among different families of exchange and correlation functionals [START_REF] Swart | Validation of exchange-correlation functionals for spin states of iron complexes[END_REF]134,[START_REF] Swart | Accurate spin-state energies for iron complexes[END_REF][START_REF] Droghetti | Assessment of density functional theory for iron(ii) molecules across the spin-crossover transition[END_REF][START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF][START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF]. Semilocal functionals, such as GGA for example, overstabilize LS [START_REF] Swart | Validation of exchange-correlation functionals for spin states of iron complexes[END_REF]145,[START_REF] Radoń | Revisiting the role of exact exchange in dft spin-state energetics of transition metal complexes[END_REF][START_REF] Mortensen | Spin propensities of octahedral complexes from density functional theory[END_REF][START_REF] Ioannidis | Towards quantifying the role of exact exchange in predictions of transition metal complex properties[END_REF][START_REF] Kulik | Perspective: Treating electron over-delocalization with the dft+u method[END_REF], while Hartree-Fock overstabilizes HS [START_REF] Rehier | Reparameterization of hybrid functionals based on energy differences of states of different multiplicit[END_REF]. Thus, depending on the system of choice, accurate ∆E H-L may be obtained by tuning the amount of exact exchange in global hybrids [START_REF] Rehier | Reparameterization of hybrid functionals based on energy differences of states of different multiplicit[END_REF][START_REF] Prokopiou | Spin-state energetics of fe complexes from an optimally tuned range-separated hybrid functional[END_REF], or by adoping density-corrected approaches [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF].

Since its original formulation [START_REF] Anisimov | Band theory and mott insulators: Hubbard u instead of stoner i[END_REF][START_REF] Anisimov | Densityfunctional theory and nio photoemission spectra[END_REF], the capability of Hubbard U -corrected DFT [START_REF] Czyżyk | Local-density functional and on-site correlations: The electronic structure of la 2 cuo 4 and lacuo 3[END_REF][START_REF] Solovyev | Corrected atomic limit in the localdensity approximation and the electronic structure of d impurities in rb[END_REF] in solving the deficiencies of standard semilocal functionals, mostly owing to self-interaction error [START_REF] Perdew | Self-interaction correction to density-functional approximations for many-electron systems[END_REF][START_REF] Perdew | Density-functional theory for fractional particle number: derivative discontinuities of the energy[END_REF][START_REF] Perdew | Physical content of the exact kohn-sham orbital energies: Band gaps and derivative discontinuities[END_REF][START_REF] Zhang | A challenge for density functionals: Self-interaction error increases for systems with a noninteger number of electrons[END_REF], has allowed it to provide a quantitatively correct description of varying properties of transition metal complexes [START_REF] Cramer | Density functional theory for transition metals and transition metal chemistry[END_REF][START_REF] Himmetoglu | Hubbard-corrected dft energy functionals: The lda+u description of correlated systems[END_REF][START_REF] Poloni | Probing the electric field-induced doping mechanism in yba 2 cu 3 o 7 using computed cu kedge x-ray absorption spectra[END_REF][START_REF] Tolba | The dft+u: Approaches, accuracy, and applications[END_REF].

Main paper

53

Its performance in predicting spin-state energetics has also been studied, to some extent. DFT+U with a self consistent U averaged among different spin states was shown to significantly improve upon DFT [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF][START_REF] Kulik | Systematic study of first-row transition-metal diatomic molecules: A self-consistent dft+u approach[END_REF]. In Ref. [START_REF] Hsu | Spin-state crossover and hyperfine interactions of ferric iron in mgsio 3 perovskite[END_REF], by using a U computed at different spin-state configurations the authors correctly predicted the ground state and the experimentally-observed pressure-induced SCO for (Mg,Fe)(Si,Fe)O 3 perovskites. Several studies have reported quantitative agreement with experiment by adopting a value of U of 4 eV [START_REF] Scherlis | Simulation of heme using dft+u: A step toward accurate spin-state energetics[END_REF][START_REF] Zhang | Intra-and interatomic spin interactions by the density functional theory plus u approach: A critical assessment[END_REF][START_REF] Ali | Electronic structure, spin-states, and spin-crossover reaction of heme-related fe-porphyrins: A theoretical perspective[END_REF]; other studies, instead, show the necessity to parametrize U against experiment for a quantitative description of the energetics [START_REF] Lebègue | Modeling spin-crossover compounds by periodic dft+u approach[END_REF][START_REF] Panchmatia | Halide ligated iron porphines: A dft+u and ub3lyp study[END_REF][START_REF] Maldonado | Two-step spin-switchable tetranuclear fe(ii) molecular solid: Ab initio theory and predictions[END_REF][START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF]. Vela and coworkers recently showed that the linear response U -approach yields an adiabatic energy difference for a well-known SCO molecular crystal, Fe(phen) 2 (NCS) 2 , deviating almost 1.8

eV from the experimentally-extracted reference value [START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF].

Here, we report a comprehensive analysis of the performance of the DFT+U method to predict adiabatic energy differences of 6-fold-coordinated Fe(II) complexes and provide an explanation of the apparent contrasting findings mentioned above, thus clarifying the limitations of this approach. We compute the ∆E H-L between the 1 A 1g LS and the 5 T 2g HS states of six Fe(II) octahedral complexes with varying ligand-field strengths by employing LDA, GGA and their Hubbard U -corrected version. Specifically, we use LDA, LDA+U sc , PBE and PBE+U sc , where the linear response-Hubbard U [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF], named U sc , is computed for each spin state. We then report ∆E H-L computed with a PBE+U approach where the same U is used for both LS and HS, in the same spirit of previous studies [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF][START_REF] Scherlis | Simulation of heme using dft+u: A step toward accurate spin-state energetics[END_REF][START_REF] Lebègue | Modeling spin-crossover compounds by periodic dft+u approach[END_REF][START_REF] Kulik | Systematic study of first-row transition-metal diatomic molecules: A self-consistent dft+u approach[END_REF][START_REF] Panchmatia | Halide ligated iron porphines: A dft+u and ub3lyp study[END_REF][START_REF] Zhang | Intra-and interatomic spin interactions by the density functional theory plus u approach: A critical assessment[END_REF][START_REF] Ali | Electronic structure, spin-states, and spin-crossover reaction of heme-related fe-porphyrins: A theoretical perspective[END_REF][START_REF] Maldonado | Two-step spin-switchable tetranuclear fe(ii) molecular solid: Ab initio theory and predictions[END_REF][START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF]. As a reference method for ∆E H-L we compute coupled cluster-corrected CASPT2 energies, i.e. CASPT2/CC, a method recently proposed by Pierloot and coworkers [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF][START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF]. Such an approach allows for the improvement of the description of electronic correlations in semicore 3s3p electrons by adopting CCSD(T) thus removing the well-known bias of CASPT2 towards high-spin states [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF][START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF].

We charge transfer whose extent and sign depends on both the spin state and the ligand-field strength. The term hybridization is here employed to describe mixing of metal and ligand atomic orbitals.

Recall that the Hubbard U term enters a DFT calculation as an extra potential acting self-consistently on the Kohn-Sham wave functions |Ψ σ i of orbital index i and spin σ:

V DFT+U = V DFT |Ψ σ i + m U 1 2 -n σ m |φ m φ m |Ψ σ i (5.1)
where the projectors {|φ m } are here atomic-like functions (the five 3d atomic functions) and the occupations, {n σ m }, are the eigenvalues of the 5×5 occupation matrix

n σ mm = i f σ i Ψ σ i |φ m φ m |Ψ σ i
, with f σ i being Fermi-Dirac occupations. Equation 5.1 shows that the Hubbard potential is repulsive for n <1/2, attractive for n >1/2 and zero for n =1/2. The DFT+U energy functional reads:

E DFT+U [ρ(r), {n σ m }] = E DFT [ρ(r)] + m,σ U 2 [n σ m (1 -n σ m )] (5.2)
For each molecule at either spin state the following strategy is adopted: (1) U is first computed with linear response [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF] by employing a DFT (without U ) optimized geometry;

then (2) this U is used to optimize again the geometry and (3) another U is calculated on the new structure. This process is iterated until convergence is reached (see SI for more details). This self-consistent and structurally-consistent U is named U sc [START_REF] Hsu | First-principles study for low-spin lacoo 3 with a structurally consistent hubbard u[END_REF]. This procedure is performed using both LDA and PBE by employing Quantum Espresso [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF]. The reference coupled cluster-corrected CASPT2 adiabatic energy differences, ∆E CASPT2/CC , are computed as follows, by adopting the procedure of Pierloot and coworkers [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF]:

∆E CASPT2/CC = ∆E CASPT2 + ∆E 3s3p,CCSD(T) (5.3) ∆E 3s3p,CCSD(T) = ∆E +3s3p -∆E no-3s3p (5.4)
Extended multi-state (XMS) CASPT2 calculations were performed using Bagel [START_REF] Shiozaki | Bagel: Brilliantly advanced general electronic-structure library modular and open-source software project for quantum simulations of materials[END_REF]177] on the TPSSH [178, 179]-optimized geometries for both HS and LS by using aug-cc-pVQZ-DK basis sets and an active space of 10 electrons on 12 orbitals, (10e,12o). This includes the 3d electrons of Fe(II), the two ligand-e g molecular orbitals plus the Fe 4d double-shell [134,[START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF], and their corresponding electrons (Figure 5.6 and 5.7). The CASPT2 calculations used for the reference set were performed without any ionisation potential-electron affinity (IPEA) shift to the zeroth-order Hamiltonian. Because this does affect energy differences of different spin states [START_REF] Ghigo | A modified definition of the zerothorder hamiltonian in multiconfigurational perturbation theory (caspt2)[END_REF][START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF], we report also the ∆E CASPT2 computed with an IPEA=0.25 a.u. in Table 5.4. For [Fe(H 2 O) 6 ] +2 we were unable to converge a (10e,12o)

active space where the two ligand-e g orbitals remained in the active space for LS. The Fe 3s orbital consistently rotated into the active space replacing one of the ligand e g orbitals.

Thus, for water the ∆E CASPT2 is taken from (10e,12o) calculations by Gagliardi and coworkers [START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF]. See computational methods in the SI for more details. To add the metal semicore 3s3p correlations, ∆E 3s3p,CCSD(T) , to the CASPT2 energy difference, CCSD(T) was performed using orca [START_REF] Neese | The orca program system[END_REF][START_REF] Neese | Software update: The orca program system, version 4.0[END_REF] on the same TPSSH-optimized geometries by including and freezing the 3s3p electrons [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF] (equation 5.4). The aug-cc-pwCVTZ-DK and cc-pVDZ basis sets were used for Fe and the ligand atoms, respectively. The whole set of CASPT2, 3s3p contribution and final CASPT2/CC adiabatic energy differences is reported in Table 5.2. CASPT2 values obtained with an active space of (6e,10o) and the Biased spin-state energetics of Fe(II) molecular complexes within density functional theory and the linear-response Hubbard-U correction ANO-RCC basis sets are also reported for comparison in Table 5.3.

We note that accurate adiabatic energy differences from diffusion Monte Carlo (DMC)

-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0 0.5 1 S=0 E dev (eV) Fe (3-q)+ 0 0.5 1 S=0 S=2 Fe (3-q)+ H 2 O NH 3 NCH PH 3 CO CNH Figure 5
.2: Deviation from linearity of the total energy with respect to the total charge of the system, from Fe 3+ to Fe 2+ , for the six molecular complexes, computed using PBE.

A lower deviation is found when the Hubbard-U approach is employed. Thick and thin lines refer to DFT and DFT+U sc results. are also available in the literature for Fe(NCH) 2+ 6 [START_REF] Fumanal | Diffusion monte carlo perspective on the spin-state energetics of [fe(nch)6]2+[END_REF]. By adopting different choices of the trial wavefunction and including a multireference treatment, the authors report a best estimate value that varies between -0.55 and -0.95 eV, depending on the geometry. Our CASPT2/CC calculations yield -0.24 eV. For Fe(NH 3 ) 2+ 6 and Fe(H 2 O) 2+ 6 , Ref. [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] reports DMC energy differences of -1.23 eV and -1.78 eV, respectively. We find -0.72 eV and -1.99 eV. The non-negligible difference between the two approaches reflects the limitations associated to each method, mainly the choice of the active space and that of the trial wavefunction. A significantly larger discrepancy, possibly due to the multiconfigurational character of LS, is found for Fe(CO) 2+ 6 , whose ∆E H-L is reported to be +0.59 eV in Ref. [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] using a single-determinant trial wavefunction. We compute a value of +1.97 eV using CASPT2/CC.

The DFT and DFT+U results for ∆E H-L are shown in Figure 5.1 (Tables 5.6 and 5.7) for increasing ligand field strengths together with our CASPT2/CC calculations. LDA and PBE systematically stabilize LS states possibly due to the well-known delocalization error [START_REF] Radoń | Revisiting the role of exact exchange in dft spin-state energetics of transition metal complexes[END_REF][START_REF] Swart | Validation of exchange-correlation functionals for spin states of iron complexes[END_REF][START_REF] Ganzenmüller | Comparison of density functionals for differences between the high-( 5 t 2g ) and low-( 1 a 1g ) spin states of iron(ii) compounds. iv. results for the ferrous complexes [fe(l)(nhs 4 )][END_REF][START_REF] Mortensen | Spin propensities of octahedral complexes from density functional theory[END_REF][START_REF] Ioannidis | Towards quantifying the role of exact exchange in predictions of transition metal complex properties[END_REF][START_REF] Kulik | Perspective: Treating electron over-delocalization with the dft+u method[END_REF]. LDA shows an even larger stabilization of LS possibly due to greater electronic delocalization as compared to PBE [START_REF] Cohen | Challenges for density functional theory[END_REF]. Conversely, DFT+U sc overstabilizes HS with an increasing deviation from reference CASPT2/CC results for stronger-field ligands. Remarkably, for weak-field ligands (i.e. for ∆E CASPT2/CC < 0), among PBE, PBE+U sc , LDA and LDA+U sc , the best agreement with the reference calculations is found for LDA+U sc , (see in Table 5.1), with a mean absolute error (MAE) of 0.644 eV, in agreement with findings by Zhang and coworkers [START_REF] Zhang | Intra-and interatomic spin interactions by the density functional theory plus u approach: A critical assessment[END_REF]. As we that PBE+U sc performs worse than PBE despite that the U sc approach leads to a better electronic density with reduced self-interaction error [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF][START_REF] Zhao | Global and local curvature in density functional theory[END_REF]. This is shown in Figure 5.2 where the deviation from linearity of the total energy with respect to fractional occupancy, which arises from the self-interaction error [START_REF] Perdew | Density-functional theory for fractional particle number: derivative discontinuities of the energy[END_REF][START_REF] Sham | Density-functional theory of the energy gap[END_REF][START_REF] Perdew | Physical content of the exact kohn-sham orbital energies: Band gaps and derivative discontinuities[END_REF][START_REF] Mori-Sánchez | The derivative discontinuity of the exchangecorrelation functional[END_REF], is reported for the six molecular complexes, for PBE and PBE+U sc . Kronik and Baer demonstrated that the average deviation from linearity is given by the energy difference between the LUMO of N-electron system and the HOMO of the N+1-electron system [START_REF] Stein | Curvature and frontier orbital energies in density functional theory[END_REF]. The deviation as a function of the fractional charge, q, is computed as the difference between the energy E(q) and the linear interpolation between the q=0 LUMO and the q=1 HOMO, which correspond to the N and N+1-electron system, here Fe 3+ and Fe 2+ , respectively. E(q) is computed by interpolating E(q=0) and E(q=1) using a cubic function which depends through q upon the total energy difference ∆E q =E(q=0)-E(q=1), and the q=1 HOMO

(Fe 2+ ), HOM O q=1
, and q=0 LUMO (Fe 3+ ), LU M O q=0

, eigenvalues [START_REF] Johnson | Extreme density-driven delocalization error for a model solvated-electron system[END_REF][START_REF] Janak | Proof that ∂e ∂n i = ε in density-functional theory[END_REF]:

E dev (q) =[( LUMO q=0 -∆E q )(1 -q)+ + (∆E q -HOMO q=1
)q]q(1 -q).

(5.5)

Biased spin-state energetics of Fe(II) molecular complexes within density functional theory and the linear-response Hubbard-U correction This function approximates the energy versus fractional occupancy curve without the need for explicit fractional-charge calculations [START_REF] Johnson | Extreme density-driven delocalization error for a model solvated-electron system[END_REF]. Figure 5.2 shows that, overall, DFT+U sc decreases the SIE, as demonstrated by the systematic decrease of E dev , upon U correction, for both LS and HS. Thus, while the density overall improves U , the DFT+U energetics do not: this is because the functional form of the Hubbard energy within the DFT+U approach, E U , leads to a bias toward the high-spin state due to the larger penalization of LS than HS [START_REF] Janet | Density functional theory for modelling large molecular adsorbate-surface interactions: a mini-review and worked example[END_REF][START_REF] Chen | Density functional plus dynamical mean-field theory of the spin-crossover molecule fe(phen) 2 (ncs) 2[END_REF].

m,σ [n σ m (1 -n σ m )]
, is plotted separately from the occupations of the e g and the t 2g states: (1) the summation term increases for increasing ligand-field strengths as a result of the larger M/L hybridization [START_REF] Kulik | Perspective: Treating electron over-delocalization with the dft+u method[END_REF]; (2) this term is larger for LS than HS, as already pointed out [START_REF] Janet | Density functional theory for modelling large molecular adsorbate-surface interactions: a mini-review and worked example[END_REF]; (3) the difference in this term between LS and HS increases for strong-field ligands (see Figure 5.3); (4) the difference between the U sc of LS and HS also increases as a function of the ligand field strength. Hybridization is here quantified by performing the Hirshfeld polulation analysis in order to decompose the active molecular orbitals from ). Thus, fractional {n σ m } arise from the projection of these occupied ligand e g onto the atomic basis because these states do exhibit some metal-like character. For HS states, the e * g are empty only for one spin channel, yielding to a total summation term that is lower than in LS, for which both spin channels are summed up in the summation.

Additionally, the larger bond lengths in the HS geometries result from occupied antibonding states that yield a smaller hybridization. This, as it is well known, contributes to a further decrease in the fractional occupation and is supported by the Hirshfeld population analysis (Tables 5.8-5.12).

Similar arguments hold to explain the larger fractional occupations of the t 2g states in LS than in HS. The hybridization between the metal non-bonding orbitals and the ligand states via π-bonding gives rise to three sets of molecular orbitals [START_REF] Jean | Molecular Orbitals of Transition Metal Complexes[END_REF] LDA+U sc gives the same qualitative trend across the molecular series as PBE+U sc but with a larger deviation from the LDA result (Figure 5.1). The ∆E U,H-L increases more in LDA than in PBE across the molecular series due to the larger values of U sc (Table 5.6) and larger fractional values for occupations {n σ m } computed in LS using LDA as compared to PBE (Figure 5.8). The better performance of LDA+U compared to PBE+U, claimed previously for weak ligands [START_REF] Zhang | Intra-and interatomic spin interactions by the density functional theory plus u approach: A critical assessment[END_REF] may now be understood: the bias introduced by E U rigidly downshifts in energy the LDA-and PBE-computed ∆E H-L . In the case of [Fe(H 2 O) 6 ] 2+ , the overstabilization of LS by LDA (approx. 2 eV) and the bias introduced by LDA+U cancel out yielding a ∆E H-L in remarkably good agreement with CASPT2/CC data, with a deviation smaller than 0.1 eV.

The larger values of U sc computed for LS compared to HS may derive from the overestimation of U sc for LS as already discussed in the literature [START_REF] Scherlis | Simulation of heme using dft+u: A step toward accurate spin-state energetics[END_REF][START_REF] Huang | The screened pseudo-charge repulsive potential in perturbed orbitals for band calculations by dft+u[END_REF][START_REF] Linscott | Role of spin in the calculation of hubbard u and hund's j parameters from first principles[END_REF]. Such an overestimation is here reflected by the overcorrection of the density produced by U sc in LS and is shown by the concave E dev in Figure 5.2.

From the above observations, since the penalizing Hubbard term, E U , is proportional to the m,σ [n σ m (1 -n σ m )] and U , the larger bias towards HS predicted for strong-field lig-

Main paper

61

ands can be mitigated by adopting a value of U smaller than U sc . To illustrate this, we report the adiabatic energy differences computed by adopting the same U value for both LS and HS for the whole series, as performed in previous studies on different compounds [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF][START_REF] Scherlis | Simulation of heme using dft+u: A step toward accurate spin-state energetics[END_REF][START_REF] Kulik | Systematic study of first-row transition-metal diatomic molecules: A self-consistent dft+u approach[END_REF]. Geometries have been optimized in each case. The computed ∆E H-L are shown in Figure 5.4 for PBE and values of U =[2,4,6,8] eV (and Table 5.6). Additionally, we report ∆E H-L computed using the average U value between the linear-response U computed for HS and LS, for each case. This procedure allows us to assess the validity of an approach where no bias is imposed by the use of different U -values while keeping an ab initio-derived parametrization of the Hubbard term. The difference ∆E U,H-L reported in the lower panel of Figure 5.4 shows, as expected, a lower bias towards HS when small values of U are used. Thus, for weak-field ligands U =4 eV yields the best agreement with CASPT2/CC, with a MAE of 0.300 eV (Table 5.1), in agreement with previous studies adopting the same value of U also for weak-field molecules [START_REF] Scherlis | Simulation of heme using dft+u: A step toward accurate spin-state energetics[END_REF][START_REF] Zhang | Intra-and interatomic spin interactions by the density functional theory plus u approach: A critical assessment[END_REF][START_REF] Ali | Electronic structure, spin-states, and spin-crossover reaction of heme-related fe-porphyrins: A theoretical perspective[END_REF]. The strongligand field complexes exhibit the best agreement with reference data for U =2 eV with a MAE of 0.354 eV (Table 5.1). Overall, the best performer is U =2 eV with a MAE of 0.365 eV. These results are consistent with U values recently obtained by mapping DFT+U ∆E H-L onto (experimentally-extracted) adiabatic energy differences and ranging between 2 and 3 eV for molecular crystals with intermediate ligand field-strengths [START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF].

The average-U calculations yield energies differences largely deviating from the reference and consistent with U values between 7 and 8 eV, thus showing the same qualitative result as the linear response-U for each spin state. Because the analysis of the density (shown in Figure 5.2) reveals a reduction of the deviation from the exact-density behavior upon Ucorrection [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF][START_REF] Zhao | Global and local curvature in density functional theory[END_REF], we employ again molecular orbital theory to rationalize the changes in electron density betweed DFT and DFT+U for the six molecular complexes. Besides, a better description of electronic correlations has been shown to affect the metal/ligand charge transfer [START_REF] Fumanal | Diffusion monte carlo perspective on the spin-state energetics of [fe(nch)6]2+[END_REF]. By applying U , a ligand↔metal charge transfer is observed whose sign and extent depends on both the electronic occupation and the ligand-field strength.

To illustrate this, the charge transfer computed using Löwdin charges for both LS and HS are reported in Figure 5. (they are shifted up in energy) whose extent depends on the projections. The d-character increases for the unoccupied e * g and decreases for the occupied e g molecular orbitals thus resulting in a more ionic-like character of the metal-ligand σ-bond and in a metal→ligand charge transfer (upper panel of Figure 5.5). This can be seen as an upshift of the metal states with respect to the ligand's, resulting in a larger relative energy difference between these (i.e. δE, blue arrow in Figure 5.5), and, thus, in a smaller M/L hybridization [START_REF] Jean | Molecular Orbitals of Transition Metal Complexes[END_REF].

Conversely, in LS the Hubbard potential is attractive for KS states that project onto the For HS, the whole set of occupations {n σ m } > 1/2 and the Hubbard potential is attractive for the whole majority spin manifold resulting in a L→M charge transfer throughout the molecular series. For the minority spin channel, only one n m >1/2 and overall the charge transfer is M→L.

In conclusion, we show and explain the performance of DFT+U sc approach in describing spin-state energetics. For weak-field molecules, such as [Fe(H 2 O) 6 ] 2+ , the overstabilization of LS by standard PBE or LDA functionals may cancel out with the the bias introduced by the Hubbard correction thus yielding to a good agreement with the reference values. For stronger ligand-fields, the larger bias towards HS state yields to deviations from the reference CASPT2/CC values by up to several eV. As a semiempirical approach, lower values of U may be adopted to mitigate this bias, thus explaining why U =4 eV, often employed in the literature, may perform better. Based on the analysis provided here, we anticipate a reduced energy-bias for spin-state energy differences involving a reduced change in the summation term upon spin transition. This is the case for weak-field ligand-molecules, as shown here, and when the lowest number of molecular orbitals is involved in the change of occupation upon spin transition, for example two molecular orbitals [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF][START_REF] Kulik | Systematic study of first-row transition-metal diatomic molecules: A self-consistent dft+u approach[END_REF] rather than four as in the present case. Finally, we rationalize the density change upon application of U in terms of change in metal/ligand hybridization. This analysis offers a general tool to interprete the U -induced charge redistribution and allows for the explaination of previous results [START_REF] Gani | Where does the density localize? convergent behavior for global hybrids, range separation, and dft+u[END_REF].

Acknowledgement

Supplementary information

Computational details 5.4.1 DFT

The DFT+U calculations are performed by using the Quantum espresso [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF] package. We use the GBRV ultrasoft pseudopotentials [START_REF] Garrity | Pseudopotentials for high-throughput dft calculations[END_REF][START_REF] Bennett | A systematic determination of hubbard u using the gbrv ultrasoft pseudopotential set[END_REF] and wavefunction and charge density cutoffs of 80 Ry and 800 Ry, respectively. Geometrical optimization is performed for PBE, PBE+U, LDA and LDA+U until the interatomic forces are smaller than 10 -4 Ry/a.u. For LDA and PBE a single optimization cycle is applied while for LDA+U and PBE+U multiple optimization cycles are performed in order to obtain a structurallyconsistent Hubbard-U parameter [START_REF] Hsu | First-principles study for low-spin lacoo 3 with a structurally consistent hubbard u[END_REF]. This is computed via an iterative procedure: starting from the U =0 optimized geometry, linear response U sc is computed and a new DFT+U sc structure is then calculated. The DFT+U sc geometry usually differs from the starting DFT structure and, therefore, a new calculation of linear response U sc leads to a different value of the Hubbard parameter that is subsequently used for the next optimization cycle. This procedure is iterated until the new DFT+U sc atomic structure does not differ from the previous one (i.e. forces change less than the tolerance given above). The whole set of geometries is listed at the end of this document. We employ atomic d-like orbitals for the projection of the Kohn-Sham states. The self-consistent correction scheme of the electrostatic potential developed by Martyna and Tuckerman [START_REF] Martyna | A reciprocal space based method for treating long range interactions in ab initio and force-field-based calculations in clusters[END_REF] was employed in all calculations to remove periodic boundary conditions effects.

CASPT2

In order to calculate the adiabatic energy difference ∆E H-L,CASPT2/CC for the six iron(II) complexes, we followed the approach proposed by Pierloot and coworkers [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF][START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF]. CASSCF/CASPT2 energy differences are computed using bagel [START_REF] Shiozaki | Bagel: Brilliantly advanced general electronic-structure library modular and open-source software project for quantum simulations of materials[END_REF] by adopting an 5.4 Supplementary information 65 active space of 10 electrons in 12 orbitals (10e,12o), i.e., including the 3d and (empty) 4d shells of Fe and the two occupied ligand-e g orbitals. CASPT2 calculations are peformed using the extended-multistate approach (XMS) with a real vertical shift of 0.2 a.u. and a Cholesky decomposition of the two-electron integrals. The three lowest states were state avaraged and treated with XMS-CASPT2 for HS while only one states was computed for LS.

Using the DFT optimized geometry from this work with no symmetry constraint, we were unable to converge a (10e,12o) active space for [Fe(H 2 O) 6 ] 2+ where the two ligand e g orbitals remained in the active space after orbital optimization for the LS state. The Fe 3s orbital consistently rotated into the active space replacing one of the e g orbitals. Attempts with other software packages were made and similar results were obtained. Therefore, we expect this is due to water being a weak field ligand and the ligand e g orbitals are sufficiently low in energy that they are doubly occupied. This implies that there is no driving force for their presence in the active space (the CASSCF energy is not lowered if these orbitals are in the active space versus the inactive space). This is supported by the fact that previously published calculations on this molecule [START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF], although on a different geometry, using the (10e,12o) active space give the same for ∆E H-L,CASPT2 as a (6e,10o) active space (the Fe 3d and 4d orbitals are active), specifically, 2.147 eV and 2.155 eV, respectively [START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF]. This implies that for this particular molecule using the smaller active space does not impact the HS-LS energy splitting significantly. This is not the case for the strong field ligands. The literature [START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF] value will be used so all reference data has the same active space.

The XMS-CASPT2 ∆E H-L results reported in the manuscript and used as reference are those computed using the Dunning aug-cc-pVQZ-DK basis set. These are reported in Table 5.3. As a comparison, adiabatic energy differences were also computed using the cc-pVTZ-DK and the ANO-RCC-TZVP basis sets, and the multistate approach for CASPT2 and results are reported in Table 5.3. These calculations are performed without any correction applied to the zeroth-order Hamiltonian used for CASPT2, i.e. the so-called ionisation potential-electron affinity (IPEA) shift. Because the adiabatic energy differences for Fe(II) complexes are known to be affected by the inclusion of the IPEA shift [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF],

, we also report values computed using =0.25. Table 5.4 reports MS-CASPT2 energy differences computed with =0 and =0.25 a.u., for comparison, using ANO-RCC-TZVP basis sets and by employing molcas [199]. The IPEA shift always increases the ∆E H-L , as expected [START_REF] Ghigo | A modified definition of the zerothorder hamiltonian in multiconfigurational perturbation theory (caspt2)[END_REF], thus stabilizing LS with respect to HS. The deviation ranges between 96 meV for NH 3 and 173 meV for NCH, consistent with recent results by Pierloot and theory and the linear-response Hubbard-U correction coworkers [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF]. The calculations in molcas also employ a 0.2 a.u. imaginary shift. The second order DKH Hamiltonian is used as well. The integrals are computed using Cholesky decomposition in combination with local exchange screening.

In order to quantify hybridization (covalency), Hirshfeld population analysis is used to decompose the active molecular orbitals into their atomic compositions, as implemented in the MultiWFN package [START_REF] Lu | Multiwfn: A multifunctional wavefunction analyzer[END_REF]. The atomic densities are computed using the built-in densities in MultiWFN and the normalized contributions to the molecular orbital are reported in Tables 5.8-5.12.

The metal semicore (3s3p) correlation effects are computed with coupled cluster CCSD(T) and the two energy differences (the CASPT2 and the 3s3p contribution computed using CC) are added to yield the final ∆E H-L,CASPT2/CC (see below). For open-shell systems, CCSD(T) calculations are performed on unrestricted Hartree-Fock wave functions (UHF).

We note that no significant spin contamination is found (see Table 5.5) for the quintet spin states of the six systems under investigation. All CASPT2 and CCSD(T) calculations were performed on the same TPSSH-optimized geometries and scalar relativistic effects were included using a second-order Douglas-Kroll-Hess (DKH) Hamiltonian while spinorbit coupling was not considered [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF]. Geometrical optimization was performed with nwchem [START_REF] Valiev | Nwchem: A comprehensive and scalable open-source solution for large scale molecular simulations[END_REF] using a cc-pVTZ basis set along with a total energy target accuracy of the numerical integration of the XC potential of 10 -8 . The default convergence criteria are enforced. These geometries are reported at the end of this document.

CCSD(T)

Single-point CCSD(T) calculations were performed for all the complexes and spin states employing an aug-CC-pwCVTZ-DK basis set for the iron center and a CC-pVDZ for the ligand atoms. The orca [START_REF] Neese | The orca program system[END_REF] package was used. As demonstrated recently by Pierloot and coworkers [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF] these basis sets allow for a correct description of the ∆E 3s3p . We worked in the frozen core approximation (FC) to correlate only the outer-core electrons while freezing the core electrons of each atom: 1s2s2p(3s3p) for Fe; 1s2s2p for P; 1s for C,N and O.

The 3s3p-correlation contribution ∆E 3s3p,CCSD(T) is calculated as the energy difference between calculations with the 3s3p Fe electrons explicitely included (∆E +3s3p ) and calculations with the 3s3p electrons kept frozen (∆E no-3s3p ). ∆E 3s3p,CCSD(T) is then obtained as: 5.14: Eigenvalues of the occupation matrix {n σ j } obtained with PBE+U sc for each spin state (HS, LS) and spin channel (σ=maj(ority),min(ority)). The difference between the sum of the eigenvalues calculated with PBE+U sc and PBE is also reported. The five occupation matrix eigenvalues refer to the occupation of the five d levels of the iron center.

∆E 3s3p,CCSD(T) = ∆E +3s3p -∆E no-3s3p (
Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density eV). We take advantage of the computational efficiency of this approach and compute the adiabatic energy differences of five molecular crystals using PBE[U] with periodic boundary conditions. The results show, again, an excellent agreement (MAE=0.07 eV)

with experimentally-extracted values and a superior performance compared with the best performers TPSSh (MAE=0.08 eV) and M06-L (MAE=0.31 eV) computed on molecular fragments.

Main paper

The accurate description of spin-state energetics of transition metal complexes represents a great challenge for electronic structure ab initio methods [START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF][START_REF] Domingo | Spin crossover in fe(ii) complexes: An ab initio study of ligand σ-donation[END_REF][START_REF] Radoń | Benchmarking quantum chemistry methods for spin-state energetics of iron complexes against quantitative experimental data[END_REF][START_REF] Cirera | Theoretical prediction of spin-crossover temperatures in ligand-driven light-induced spin change systems[END_REF]134,[START_REF] Swart | Accurate spin-state energies for iron complexes[END_REF][START_REF] Droghetti | Assessment of density functional theory for iron(ii) molecules across the spin-crossover transition[END_REF][START_REF] Cirera | Benchmarking density functional methods for calculation of state energies of first row spin-crossover molecules[END_REF]]. Yet, the accurate prediction of spin-state energy differences are of critical importance for the understanding of spin crossover phenomena relevant for example for spintronics, molecular elecronics and sensors [START_REF] Kumar | Emerging trends in spin crossover (sco) based functional materials and devices[END_REF][START_REF] Molnár | Molecular spin crossover materials: Review of the lattice dynamical properties[END_REF][START_REF] Resines-Urien | A switchable iron-based coordination polymer toward reversible acetonitrile electro-optical readout[END_REF][START_REF] Resines-Urien | Covalent post-synthetic modification of switchable iron-based coordination polymers by volatile organic compounds: a versatile strategy for selective sensor development[END_REF] and for the catalytic reactivity of biological systems [START_REF] Shaik | Exchange-enhanced reactivity in bond activation by metal-oxo enzymes and synthetic reagents[END_REF]. This challenge stems from the lack of error cancellation when computing energy differences, using approximate electronic structure methods, between spin states exhibiting different types and amounts of electronic correlations (dynamic and non-dynamical) [START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF]. Because Hartree-Fock (HF) only treats exchange correlations, for example, it tends to stabilize high-spin states over low spin states due to the absence of dynamical correlation that would stabilize doubly occupied orbitals [START_REF] Swart | Accurate spin-state energies for iron complexes[END_REF][START_REF] Rehier | Reparameterization of hybrid functionals based on energy differences of states of different multiplicit[END_REF]. On the contrary, local and semilocal functionals within DFT tend to overstabilze low spin states [START_REF] Swart | Validation of exchange-correlation functionals for spin states of iron complexes[END_REF]145,[START_REF] Mortensen | Spin propensities of octahedral complexes from density functional theory[END_REF][START_REF] Ioannidis | Towards quantifying the role of exact exchange in predictions of transition metal complex properties[END_REF][START_REF] Ganzenmüller | Comparison of density functionals for differences between the high-( 5 t 2g ) and low-( 1 a 1g ) spin states of iron(ii) compounds. iv. results for the ferrous complexes [fe(l)(nhs 4 )][END_REF][START_REF] Salomon | Assertion and validation of the performance of the b3lyp* functional for the first transition metal row and the g2 test set[END_REF] and by adding a fraction of exact exchange one can, in most cases [START_REF] Radoń | Revisiting the role of exact exchange in dft spin-state energetics of transition metal complexes[END_REF], reduce such overstabilization [START_REF] Pinter | Conceptual insights into DFT spin-state energetics of octahedral transition-metal complexes through a density difference analysis[END_REF][START_REF] Prokopiou | Spin-state energetics of fe complexes from an optimally tuned range-separated hybrid functional[END_REF]. Thus, global hybrids can provide a reasonable decription of spin-state energetics depending on the system of choice and the amount of exact exchange [START_REF] Rehier | Reparameterization of hybrid functionals based on energy differences of states of different multiplicit[END_REF][START_REF] Prokopiou | Spin-state energetics of fe complexes from an optimally tuned range-separated hybrid functional[END_REF].

Song at al. showed that a HF density-corrected DFT approach can yield adiabatic energy differences in good agreement with diffusion Monte Carlo (DMC) calculations [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF]. The DFT+U approach has also been investigated in this respect in a few studies and we refer the reader to the introduction of Ref. [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF] for a recent summary on the topic. The present authors have shown and discussed the bias towards high spin states imposed by the Hubbard term in the total energy and how it can be mitigated by adopting values of U smaller than the computed self-consistent value, U sc [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. Despite the energetics being pathologically wrong for strong-field ligands, the electronic density exhibits a systematic improvement with respect to local and semilocal functionals for both low spin and high spin states and for all systems when adopting a self-consistent U [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF][START_REF] Zhao | Global and local curvature in density functional theory[END_REF][START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. In this work we merge the above ideas and adopt a new approach consisting of a Hubbard U densitycorrected DFT where the PBE functional is evaluated on the Hubbard U density, using a linear-response U [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF] computed self-consistently [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF]. We show that this method allows one to obtain adiabatic energy differences for a series of six Fe(II) molecular complexes in excellent agreement with the chosen reference set. The molecular complexes include varying ligand field strengths from the weak H 2 O ligand, whose reference ∆E H-L is -1.83 eV, to the strong CNH one, whose ∆E H-L is 2.87 eV (see Tab. 6.1), thus allowing for a better assessment of the validity of this approach. Although the choice of the reference method is still matter of debate, we choose the coupled-cluster corrected CASPT2 approach proposed by Pierloot and coworkers [START_REF] Pierloot | Spin state energetics in first-row transition metal complexes: Contribution of (3s3p) correlation and its description by second-order perturbation theory[END_REF][START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF]. This approach reduces the overstabilization Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density of high spin state by treating the 3s and 3p semicore electrons using CCSD(T) and can be used in principle in systems with non-negligible multiconfigurational character such as the strong-ligand field molecules studied here. Its accuracy has been recently further validated by Radoń by comparing with ∆E H-L values extracted from experiments [START_REF] Radoń | Benchmarking quantum chemistry methods for spin-state energetics of iron complexes against quantitative experimental data[END_REF]. We show that for Fe(II) complexes exhibiting a weak ligand strength, our result compare very well also with CCSD(T) [START_REF] Domingo | Spin crossover in fe(ii) complexes: An ab initio study of ligand σ-donation[END_REF] and recent DLPNO-CCSD(T) results [START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF]. Larger deviations are found with respect to DMC results (vide infra).

6.2 Main paper 85 -3 -2 -1 0 1 2 3 4 H 2 O NH 3 NCH PH 3 CO CNH ΔE H-L (eV) PBE PBE+U sc PBE+U sc +V sc PBE[U sc +V sc ] PBE[U sc ] PBE[HF] CASPT2/CC -3 -2 -1 0 1 2 3 4 H 2 O NH 3 NCH PH 3 CO CNH ΔE H-L (eV) PBE PW91 BLYP TPSS SCAN M06L CASPT2/CC -3 -2 -1 0 1 2 3 4 H 2 O NH 3 NCH PH 3 CO CNH ΔE H-L (eV) B3LYP PBE0 M06 M06-2X TPSSh TPSS0 CASPT2/CC
We then apply this approach to compute the spin crossover energies of seven compounds, either crystalline or molecular, for which the adiabatic energy differences have been extracted from experiments and we find again very good agreement. In light of this accuracy, this approach can be adopted to study molecular crystals very efficiently with any DFT code including a DFT+U implementation thus avoiding the use of hybrid functionals.

We recall that the DFT+U total energy can be written as:

E DFT+U [ρ(r)] = E DFT [ρ(r)] + E U [n] (6.1)
In the above formula the term E DFT represents the unperturbed DFT energy functional and the E U is the Hubbard term containing the Hubbard correction for the electronic repulsion within a given subshell and a double-counting term that removes the interactions that are already counted within the DFT term via mean-field. For a clear review we refer the reader to Ref. [START_REF] Himmetoglu | Hubbard-corrected dft energy functionals: The lda+u description of correlated systems[END_REF]. The E U depends on the density through the occupation numbers n computed from the projection of the occupied Kohn-Sham eigenfunctions onto a localized basis set. For projection numbers close to 1/2 the summation term that enters E U and that mutliplies U is the largest. In our recent work we showed that the DFT+U energy yields a systematic bias towards high spin due to the E U term being systematically larger for low spin states thus resulting in a destabilization of the latter with respect to the former. This bias increases as a function of the ligand field strength: for stronger field ligands the more covalent bonding between Fe and the ligand yields more fractional occupations thus resulting in larger penalizing summation terms [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. While this penalizing term is necessary to recover the localization of electrons and stabilize the insulating phase in Mott physics, here it results in a systematic unphysical overstabilization of high spin which further increases for molecular complexes with larger covalent character, such as the CO and CNH strong field ligands.

Cococcioni and coworkers implemented an extended Hubbard model in DFT through the inclusion an inter-site effective interaction V within the Hubbard energy. Such a generalized scheme, named DFT+U+V [START_REF] Cococcioni | Extended dft+u+v vmethod with on-site and inter-site electronic interactions[END_REF], aims at an improved treatment of electronic -0.09 [START_REF] Daku | Accurate spinstate energetics of transition metal complexes. 1. ccsd(t), caspt2, and dft study of [m(nch) 6 ] 2 (m = fe, co)[END_REF] DLPNO-CCSD(T1) -1.44 [START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF] -0.49 [START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF] -0.38 [START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF] DMC -1.78 [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] -1.23 [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] -1.17 [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] 0.59 [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] -2.60 [START_REF] Droghetti | Assessment of density functional theory for iron(ii) molecules across the spin-crossover transition[END_REF] -1.55 [START_REF] Droghetti | Assessment of density functional theory for iron(ii) molecules across the spin-crossover transition[END_REF] -1.37 [START_REF] Droghetti | Assessment of density functional theory for iron(ii) molecules across the spin-crossover transition[END_REF] -0.85,-0.95 [START_REF] Fumanal | Diffusion monte carlo perspective on the spin-state energetics of [fe(nch)6]2+[END_REF] Table 6.1: ∆E H-L (in eV) computed using different DFT approaches (upper left table).

These are reported in order of descreasing total MAE computed with respect to the CASPT2/CC reference (see text). The MAE computed separately for weak-, strong-ligand and for the whole set (total) are also reported (upper right table). The ∆E H-L computed using varying wavefunction methods and taken from the literature are also reported (lower table). a refers to the CASPT2/CC values computed in this work.

Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density correlations. The new Hubbard potential includes two terms of opposite sign: the first on-site term is attractive for Kohn-Sham states that exhibit a localized character (the standard on-site U term) whereas the second inter-site term stabilizes hybridized states. Thus, a competition between these two opposing behaviors should allow for a more balanced description of electronic correlations and thus improved structural and electronic properties [START_REF] Cococcioni | Extended dft+u+v vmethod with on-site and inter-site electronic interactions[END_REF][START_REF] Ricca | Self-consistent DFT + u + v study of oxygen vacancies in srtio 3[END_REF]. For the six Fe(II) molecular complexes computed here, i.e. [Fe(H

2 O) 6 ] +2 , [Fe(NH 3 ) 6 ] +2 , [Fe(NCH) 6 ] +2 , [Fe(PH 3 ) 6 ] +2 , [Fe(CNH) 6 ]
+2 and [Fe(CO) 6 ] +2 , the geometries optimized using the TPSSh functional are taken from Ref. [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF] and used for all calculations, i.e. DFT, CCSD(T) and CASPT2. All DFT calculations, except for the DFT+U, were performed using ORCA [START_REF] Neese | The orca program system[END_REF][START_REF] Neese | Software update: The orca program system, version 4.0[END_REF]. The DFT+U and DFT+U+V calculations were performed using Quantum ESPRESSO [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF] by adopting a linear-response approach [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF] for the self-consistent calculation of U [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF], i.e. U sc , and U +V [START_REF] Timrov | Hubbard parameters from densityfunctional perturbation theory[END_REF], i.e.

U sc +V sc , respectively. We stress that in what follows, DFT+U or DFT+U+V always refer to self-consistent calculations, unless otherwise specified (e.g. in the results discussed later in Figure 6.4). See SI for more details. Unlike our recent work where a few geometrical optimizations were performed upon calculation of the linear-response U to yield a structurally consistent U [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF], here U and U +V are computed on the TPSSH geometries for LS and HS separately. Because in what follows we report errors computed as deviations from the reference values, we intend to avoid including effects arising from different geometries. The effect of the employed geometry on the ∆E H-L =E HS -E LS has also been investigated and will be discussed below. The projections for the Hubbard term are performed using orthonormalized atomic wavefunctions. This yields ∆E H-L systematically smaller than those computed with non-orthogonal non-normalized atomic projectors (as those we employed previously [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]), as shown in Tab. 6.1 and Tab. 6.4.

Extended multi-state CASPT2 calculations were performed employing BAGEL [START_REF] Shiozaki | Bagel: Brilliantly advanced general electronic-structure library modular and open-source software project for quantum simulations of materials[END_REF]177] using an active space of 10 electrons in 12 orbitals, (10e,12o). This includes the 3d electrons of Fe(II), the two ligand-e g molecular orbitals plus the Fe 4d double-shell [134,[START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF], and their corresponding electrons. Density fitting was used for all calculations by employing the fitting basis set cc-pV5Z-JKFIT and no symmetry constraints were imposed. The CASPT2 calculations used for the reference set were performed without any ionisation potential-electron affinity (IPEA) shift to the zeroth-order Hamiltonian. Because of the well established slow convergence of the CASPT2 energy with respect to basis set size, we perform the extrapolation of the spin-state energies to the complete basis set (CBS) limit. This is done separately for CASSCF and the CASPT2 energies, by adopting the three-point extrapolation method described in Refs. [START_REF] Feller | Application of systematic sequences of wave functions to the water dimer[END_REF][START_REF] Feller | The use of systematic sequences of wave functions for estimating the complete basis set, full configuration interaction limit in water[END_REF][START_REF] Helgaker | Basis-set convergence of correlated calculations on water[END_REF]. The cc-pVTZ-DK, cc-pVQZ-DK, and cc-pV5Z-DK basis sets were used for this and the corresponding CASPT2 ∆E H-L are reported in Tab. 6.3.

For [Fe(H 2 O) 6 ] +2 we were unable to converge a (10e,12o) active space where the two ligand-e g orbitals remained in the active space for LS. The Fe 3s orbital consistently rotated into the active space replacing one of the ligand e g orbitals. This implies that for this particular molecule using the smaller active space should not impact the HS-LS energy splitting significantly, as reported by Gagliardi and co-workers [START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF] who reported CASPT2 values of -2.14 eV with (6e,10o) and -2.15 eV with (10e,12o). Differences of the order of 0.1 eV are reported in Ref. [134]. Thus, for water the extrapolation to CBS is performed using a (6e,10o) active space. See computational methods in the SI for more details.

The Fe semicore 3s3p correlation energy is computed using CCSD(T) by including and freezing the 3s3p electrons [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF], using ORCA. This correction is then added to the CASPT2 energy difference to yield the CASPT2/CC energy difference. The aug-cc-pwCVTZ-DK and cc-pVDZ basis sets were used for Fe and the ligand atoms, respectively (see details in SI). Extrapolation to the CBS limit is not required here as demonstrated by Pierloot and coworkers [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF].

The adiabatic energy differences, ∆E H-L , computed using several choices of DFT func- throughout the molecular series [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF] due to the penalizing Hubbard term being larger for LS and for strong-field ligands [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. A minor improvement of DFT+U+V as compared to DFT+U is found, possibly due to the values of V being too low to correct for the bias towards HS (see Tab. 6.5). For [Fe(PH 3 ) 6 ] +2 , we were unable to converge the DFT+U+V calculations for HS and thus the corresponding ∆E H-L is omitted.

Despite yielding erroneous spin-state energetics for the molecular series reported here, the DFT+U with a linear-response U approach systematically improves the electronic density, regardless of the spin state, with a reduction of the energy bowing as a function of fractional occupations which is a manifestation of self-interaction error [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF][START_REF] Zhao | Global and local curvature in density functional theory[END_REF][START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]]. Song at al. [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] discuss the case of spin gaps in Fe(II) octahedrally-coordinated complexes in terms of calculations affected by large errors in the density: the error that arises from the approximation of the exchange-correlation functional is comparable or smaller than the error introduced by the use of the approximate density [START_REF] Sim | Quantifying density errors in DFT[END_REF]. In this respect, the density-corrected DFT approach, discussed in detail in Refs. [START_REF] Wasserman | The importance of being inconsistent[END_REF][START_REF] Vuckovic | Density functional analysis: The theory of density-corrected DFT[END_REF], consists in employing approximate density functionals on a density different than the self-consistent one and possibly closer to the exact one. This approach implemented using the Hartree-Fock density has been shown to improve over the self-consistent DFT results the description of many properties such as reaction barriers [START_REF] Janesko | Hartree-fock orbitals significantly improve the reaction barrier heights predicted by semilocal density functionals[END_REF][START_REF] Verma | Increasing the applicability of dft i: Nonvariational correlation corrections from hartree-fock DFT for predicting transition states[END_REF], weak intermolecular forces [START_REF] Gordon | Theory for the forces between closed-shell atoms and molecules[END_REF], bond energies [START_REF] Kim | Halogen and chalcogen binding dominated by density-driven errors[END_REF] and the binding properties of anions [START_REF] Kim | Communication: Avoiding unbound anions in density functional calculations[END_REF][START_REF] Kim | Ions in solution: Density corrected density functional theory (dc-dft)[END_REF]. The authors of Ref. [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] showed systematically improved results for spin-state splittings of fours Fe(II) molecular complexes computed using the DFT[HF] approach. Our working hypothesis is that the Hubbard U -corrected density should yield more accurate results compared to a HF density, since the latter only includes exchange correlations while neglecting dynamic and non-dynamical correlations. We employ a density-corrected DFT by adopting a standard semilocal functional, such as PBE, evaluated on the Hubbard U density. By doing so, we remove the energy bias introduced by the E U term discussed above while keeping an improved electronic density. Practically, we perform self-consistent DFT+U calculations by computing U using linear-response theory and then remove the E U term from the total energy to compute ∆E H-L . This is not, strictly speaking, a non-self consistent, density-corrected DFT calculation as the kinetic term is computed using the DFT+U orbitals. However, our assumption is that the kinetic energy computed using DFT orbitals matches closely the DFT+U case so that this approach can be seen as a non-self consistent density-corrected DFT method. This assumption is motivated by a recent study employ- Hartree-Fock one are negligibly different when computed on the same HF density [START_REF] Nam | Measuring density-driven errors using kohn-sham inversion[END_REF].

The results of the Hubbard U sc -corrected density employed using a PBE functional are shown in Figure 6.1 and are named PBE[U] henceforth. In the same figure we also show the PBE[HF] results, i.e. using the Hartree-Fock density. We stress that the PBE+U and the PBE values are slightly different compared to those reported in our previous work [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF] because of the different geometries employed and the atomic basis used for the projections.

The ), however, our conclusion on the accuracy of DFT[HF] is somehow different owing to the difference in the corresponding reference values. In particular, the DMC values in Ref. [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF] are systematically lower compared to CASPT2/CC values and the largest deviation is found for the CO and NCH. See Tab. 6.1 for the whole list of ∆E H-L computed either here or in previous studies using wavefunction methods.

We note the reasonably good agreement between our CASPT2/CC reference values and published CCSD(T) [START_REF] Domingo | Spin crossover in fe(ii) complexes: An ab initio study of ligand σ-donation[END_REF][START_REF] Daku | Accurate spinstate energetics of transition metal complexes. 1. ccsd(t), caspt2, and dft study of [m(nch) 6 ] 2 (m = fe, co)[END_REF] and DLPNO-CCSD(T) [START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF] values for weak-field molecules (see Tab. 6.1).

The performance of varying DFT functionals for the calculation of adiabatic energy differences has been reported in the literature by several authors [134, 227,

. Thus, we do compute the ∆E H-L using a few DFT functionals in order to establish a comparison for the performance of DFT[U], however, we refer the reader to these articles for a more detailed discussion. GGA functionals overstabilize the LS state, although BLYP does so to a lesser extent compared to PBE and PW91. Excellent results have been reported in the past [START_REF] Swart | Accurate spin-state energies for iron complexes[END_REF][START_REF] Pierloot | Relative energy of the high-( 5 t 2g ) and low-( 1 a 1g ) spin states of the ferrous complexes [fe(l)(nhs 4 )]: Caspt2 versus density functional theory[END_REF] using the optimized OPTX exchange proposed by Handy and Cohen [START_REF] Handy | Left-right correlation energy[END_REF]. By adopting global hybrids with increasing admixtures of exact exchange HS is systematically stabilized. PBE0 and B3LYP with a 25% and 20% admixture of exact exchange added respectively to PBE and BLYP functionals [START_REF] Mortensen | Spin propensities of octahedral complexes from density functional theory[END_REF] show an overcorrection and overall overstabilize HS. Among the meta-GGAs, M06-L performs well in comparison to other meta-GGA functionals as already observed in previous studies [START_REF] Ioannidis | Towards quantifying the role of exact exchange in predictions of transition metal complex properties[END_REF][START_REF] Ioannidis | Ligand-field-dependent behavior of meta-GGA exchange in transition-metal complex spin-state ordering[END_REF][START_REF] Cirera | Assessment of the SCAN functional for spin-state energies in spin-crossover systems[END_REF][START_REF] Flöser | Detailed pair natural orbital-based coupled cluster studies of spin crossover energetics[END_REF].

Among the studied functionals, the smallest error is found for the meta-hybrid TPSSh (15% of exact exchange) in agreement with several recent studies [START_REF] Cirera | Theoretical prediction of spin-crossover temperatures in ligand-driven light-induced spin change systems[END_REF][START_REF] Cirera | Theoretical modeling of spin crossover in metal-organic frameworks: [fe(pz) 2 pt(cn) 4 ] as a case study[END_REF][START_REF] Cirera | Benchmarking density functional methods for calculation of state energies of first row spin-crossover molecules[END_REF]. Climb-Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density ing up the DFT Jacob's ladder other functionals like double-hybrid [START_REF] Wilbraham | Communication: Evaluating nonempirical double hybrid functionals for spin-state energetics in transition-metal complexes[END_REF][START_REF] Alipour | Appraising spin-state energetics in transition metal complexes using double-hybrid models: Accountability of sos0-pbescan0-2(a) as a promising paradigm[END_REF] and rangeseparated hybrid functionals [START_REF] Prokopiou | Spin-state energetics of fe complexes from an optimally tuned range-separated hybrid functional[END_REF] Specifically, regardless of the functional used to compute ∆E H-L , the PBE+U geometry always yields the largest decrease in ∆E H-L . This is consistent with the fact that the PBE+U geometry computed using a structurally consistent approach deviates the most from the TPSSh optimized geometry as shown in Tabs. 6.8 and 6.9. Because the effect of U is larger for LS [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF] due to the larger values of E U as compared to HS, any other functional would destabilize LS more than HS. For stronger-field molecules this effect is more pronounced as confirmed by the larger increase in metal-ligand bond distances in the PBE+U LS geometry with respect to PBE (Tabs. 6.8 and 6.9).

To understand why the Hubbard U -corrected PBE density yields significantly improved results compared with the PBE density and the HF density, we compare densities from PBE, PBE+U, PBE+U+V (again computed with U sc and V sc ), and Hartree-Fock density with that obtained from the relaxed CASPT2 spin-density matrix in BAGEL [START_REF] Shiozaki | Hyperfine coupling constants from internally contracted multireference perturbation theory[END_REF]. The relaxed spin-density matrix is obtained by adding orbital and configurational relaxation contributions due to dynamical correlation to the unrelaxed density matrix using the CASPT2

Lagrangian [START_REF] Shiozaki | Hyperfine coupling constants from internally contracted multireference perturbation theory[END_REF]. We study the LS case of [Fe(NH 3 one indicating charge depletion and accumulation, respectively, for PBE compared with CASPT2. We note that CASPT2 calculations with the inclusion of bonding metal-ligand e g states and the 4d double shell have been shown to account for non-dynamical correlation [START_REF] Pierloot | The caspt2 method in inorganic electronic spectroscopy: from ionic transition metal to covalent actinide complexes[END_REF]134,[START_REF] Pierloot | Computational organometallic chemistry[END_REF][START_REF] Pierloot | Transition metals compounds: Outstanding challenges for multiconfigurational methods[END_REF]. DFT does not account for non-dynamical correlation, however the self-interaction error arising from the implementation of approximate density-functionals yields an overdelocalization of the charge density along the chemical bonds (and less charge near the atom) and a more diffuse character of the electron cloud around the nuclei, as shown in Figure 6.2, that can actually mimic these effects [START_REF] Polo | Some thoughts about the stability and reliability of commonly used exchange?correlation functionals ? coverage of dynamic and nondynamic correlation effects[END_REF][START_REF] Gräfenstein | The impact of the self-interaction error on the density functional theory description of dissociating radical cations: Ionic and covalent dissociation limits[END_REF][START_REF] Harvey | On the accuracy of density functional theory in transition metal chemistry[END_REF] sometimes called left-right and radial correlations, respectively. As shown in Figure 6.2, and as discussed in the literature, these effects are exaggerated in PBE. The PBE+U sc density is qualitatively similar to the PBE density but with a smaller deviation from the reference one within the Fe region. We chose not plot the PBE+U+V density because it yields a plot visually identical to the PBE+U one. The Hartree-Fock density exhibits the opposite behavior near the Fe, i.e. charge density accumulates and depletes with the e g and t 2g orbitals, respectively. This is consistent with the lack of explicit non-dynamical correlations and absence of self-interaction error. Thus, the effect of the Hubbard U term on the density is qualitatively similar to the case found when increasing the exact exchange admixture in global hybrid functionals [START_REF] Pinter | Conceptual insights into DFT spin-state energetics of octahedral transition-metal complexes through a density difference analysis[END_REF].

In what follows we attempt to quantify the error on the density by considering, again, only the region around the Fe. For each method, we extract the electronic charge density distribution around the iron centre, ρ F e x (r), by employing the Bader scheme [START_REF] Henkelman | A fast and robust algorithm for bader decomposition of charge density[END_REF][START_REF] Sanville | Improved grid-based algorithm for bader charge allocation[END_REF][START_REF] Tang | A grid-based bader analysis algorithm without lattice bias[END_REF][START_REF] Yu | Accurate and efficient algorithm for bader charge integration[END_REF]. We evaluate the error on the density, ∆ρ x , as a deviation from the reference, ρ F e CASPT2 (r), within the Bader region (see SI for more details), as follows: This behavior is consistent with the DFT+U approach correcting the density more for strong-field ligand molecules, as shown and discussed previously [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. The trend along the four molecules correlates with trends in non-dynamical correlations. In agreement with previous studies [START_REF] Domingo | Spin crossover in fe(ii) complexes: An ab initio study of ligand σ-donation[END_REF]134,[START_REF] Pierloot | The caspt2 method in inorganic electronic spectroscopy: from ionic transition metal to covalent actinide complexes[END_REF], we find that moving along the spectrochemical series nondynamical correlation becomes more important. The configuration interaction weight of the dominant electronic configuration computed from the CASSCF calculation decreases from 94% to 89% going from NH 3 to CNH (see lower panel of Figure 6.3). This analysis is in line with our results showing HF to perform better for molecular complexes with weakfield ligands and lower non-dynamical correlation. One would thus expect that HF density would overstabilize HS compared to LS more for strong field molecules, which is indeed the case here (see Figure 6.1 and Tab. 6.1). PBE+U (with U sc ) systematically improves the electronic density for both weak and strong-field ligand molecules thus yielding an improved description of the spin-state energetics throughout the spectrochemical series. This is further shown in Figure 6.4 where we report the ∆E H-L computed using the PBE functional evaluated on the PBE+U density, for increasing values of U . We only show the results for four complexes for clarity. Higher values of Hubbard U stabilize HS more compared to LS, as expected, and the deviation of ∆E H-L from the reference value (shown as an horizontal line on the right y-axis) systematically decreases as U increases.

∆ρ x = |ρ F e,Bader
The effect of U -corrected density on the spin energetics is qualitatively similar to the effect observed when adopting densities computed with increasing amounts of exact exchange [START_REF] Radoń | Revisiting the role of exact exchange in DFT spin-state energetics of transition metal complexes[END_REF][START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF][START_REF] Pinter | Conceptual insights into DFT spin-state energetics of octahedral transition-metal complexes through a density difference analysis[END_REF]. It must be noted, however, that the change in ∆E H-L reported here is significantly larger than those computed with a density-corrected approach using hybrid functionals [START_REF] Radoń | Revisiting the role of exact exchange in DFT spin-state energetics of transition metal complexes[END_REF][START_REF] Pinter | Conceptual insights into DFT spin-state energetics of octahedral transition-metal complexes through a density difference analysis[END_REF].

To further test the validity of PBE[U], we compute ∆E H-L for a set of seven Fe(II) compounds for which the HS-LS energy differences have been extracted from experimental data. The first five compounds, Fe(phen) 2 (NCS) 2 [START_REF] Gallois | Structural changes associated with the spin transition in bis(isothiocyanato)bis(1,10-phenanthroline)iron: a single-crystal x-ray investigation[END_REF] (phen=1,10-phenanthroline), Fe(abpt) 2 (NCS) 2 , and Fe(abpt) 2 (NCSe) 2 from Ref. [START_REF] Moliner | Spin-crossover in the [fe(abpt)2(ncx)2] (x=s, se) system: structural, magnetic, calorimetric and photomagnetic studies[END_REF] with abpt=4-amino-3,5-bis(pyridin- (tacn= 1,4,7-triazacyclononane) and [Fe(bpy) 3 ] +2 (bpy=2,2'-bipyridine). The ∆E H-L of [Fe(tacn) 2 ] +2 has been extracted by Radoń [START_REF] Radoń | Benchmarking quantum chemistry methods for spin-state energetics of iron complexes against quantitative experimental data[END_REF] using an approach similar to Ref. [START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF].

Main paper

97

The spin gap of [Fe(bpy) 3 ] +2 has been extracted by Casida et al. from the light-induced population of the high-spin state [START_REF] Daku | Assessment of density functionals for the high-spin/low-spin energy difference in the low-spin iron(II) tris(2, 2 -bipyridine) complex[END_REF]. For all these complexes we first adopt a molecular model to compute the ∆E H-L . For the five molecular crystals, this is done by carving a structure from the fully optimized geometry using periodic boundary conditions, similarly to the procedure adopted in Ref. [START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF]. The geometrical optimization is performed using Quantum Espresso using the PBE functional together with the semiempirical Grimme's D3 correction [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-d) for the 94 elements h-pu[END_REF] combined with the Becke-Johnson (BJ) damping scheme [START_REF] Smith | Revised damping parameters for the d3 dispersion correction to density functional theory[END_REF]. The Hubbard U sc is then computed on the optimized geometry, using periodic boundary condi- with periodic boundary conditions represents the best performers with a MAE of 0.07 eV, i.e. slightly smaller compared to the same calculation performed on molecular fragments.

This results confirms the good performance of PBE[U] established above using ab initio data as reference, and it shows its potential for the effecient calculation of adiabatic energy differences in crystalline complexes.

Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density Kroll-Hess (DKH) Hamiltonian while spin-orbit coupling was not considered [START_REF] Phung | Toward highly accurate spin state energetics in first-row transition metal complexes: A combined caspt2/cc approach[END_REF].

By adopting the TPSSh geometries without any symmetry constraint, we were unable to converge a (10e,12o) active space for [Fe(H 2 O) 6 ] 2+ in LS where the two ligand e g orbitals remained in the active space after orbital optimization. The Fe 3s orbital consistently rotated into the active space replacing one of the e g orbitals. This implies that for this particular molecule using the smaller active space should not impact the HS-LS energy splitting significantly, as discussed in more detail in the SI of Ref. [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. For this reason, the CASPT2 calculations used for the CASPT2/CC reference for water are performed using the (6e,10o) active space.

The extrapolation to the complete basis set (CBS) limit was performed by fitting three energy values computed using the cc-pVTZ-DK, cc-pVQZ-DK, and cc-pV5Z-DK basis sets.

The CASSCF and the CASPT2 energies were fitted separately by adopting the procedure by Feller [START_REF] Feller | Application of systematic sequences of wave functions to the water dimer[END_REF][START_REF] Feller | The use of systematic sequences of wave functions for estimating the complete basis set, full configuration interaction limit in water[END_REF] and Helgaker et al. [START_REF] Helgaker | Basis-set convergence of correlated calculations on water[END_REF], respectively. The CASSCF energy is extrapolated using the a+bexp(-cX) function while for the second-order energy correction we use a+bX -3 . The ∆E H-L computed using CASPT2 and the cc-pVTZ-DK, cc-pVQZ-DK, and cc-pV5Z-DK basis sets, respectively, are reported in Tab. 6.3. All the calculations are performed without any empirical correction applied to the zeroth-order Hamiltonian used for CASPT2, i.e. the so-called ionisation potential-electron affinity (IPEA) shift. In the SI of our previous study we have reported a comparison between calculations performed without and with an IPEA shift of 0.25 a.u. The ∆E H-L always decrease upon inclusion of the IPEA shift and the largest deviation was found to be 173 meV for NCH [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. for CASSCF and CASPT2) to the limit of complete basis set.

6.4 Supplementary information 101

CCSD(T)

The 3s3p correlation contribution computed using CCSD(T) to yield the final CASPT2/CC results are taken from our previous work [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. We therefore refer the reader to the SI of that work for a detailed discussion.

Projectors in DFT+U

Within the DFT+U approach the energy is written as:

E DFT+U [ρ(r), {n σ m }] = E DFT [ρ(r)] + m,σ U 2 [n σ m (1 -n σ m )] (6.3)
where {|φ m } are the chosen atomic-like functions and the occupations, {n σ m }, are the eigenvalues of the 5×5 occupation matrix n σ mm = i f σ i Ψ σ i |φ m φ m |Ψ σ i , with f σ i being Fermi-Dirac occupations. In Quantum ESPRESSO it is possible to choose among atomic orbitals obtained by solving the Schrödinger equation for the isolated atom, as they are or, alternatively, by applying an orthogonalisation procedure based on Löwdin decomposition in order to obtain a localised basis set with zero overlap with the neighbour atomic basis (ortho-atomic projectors). We find a significant difference in the values of U and, consequently, in the computed ∆E H-L when the atomic and ortho-atomic projects are used (see Tab. 6.4). The ∆E H-L reported in the manuscript are computed using the so called ortho-atomic projectors while in our previous work [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF] we employed the atomic ones.

6.4.5 Linear-response U and V DFT+U sc and DFT+U sc +V sc calculations are performed using linear-response to evaluate the Hubbard parameters U and V [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF][START_REF] Campo | Extended dft+u+v method with on-site and inter-site electronic interactions[END_REF][START_REF] Timrov | Hubbard parameters from densityfunctional perturbation theory[END_REF]. The Hubbard V IJ parameter between atomic sites I and J is computed from the interacting and non-interacting response functions χ IJ and χ 0 IJ as:

V IJ = (χ -1 0 -χ -1 ) IJ , (6.4)
and the on-site U parameter of DFT+U is the special case V II . Within the DFT+U sc calculations we evaluate the U on the iron with an home-made script by computing the response functions χ = ∂n/∂α and χ 0 = ∂n 0 /∂α from a linear regression between the relative change of the screened (self-consistent) and bare occupations, i.e n and n 0 , respectively, and the applied perturbation α, as explained in Ref. [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF]. This is done iteratively. Because the ground state changes upon inclusion of the Hubbard correction, the U computed on the DFT ground state may be different than the one computed on the DFT+U one [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF]. Thus a self-consistent linear-response U calculation is implemented Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density where a first linear-response calculation is performed on the U=0 ground-state, and then the computed U is used to obtained a new DFT+U ground-state density upon which a further linear-response U calculation is run. This is repeated until the U used to obtain the unperturbed charge density is close to the output linear-response U . We impose a convergence threashold of 0.05 eV between input and output U which is usually achieved in 2-3 steps.

The DFT+U+V calculations were performed by adopting the approach reported in

Ref. [START_REF] Timrov | Hubbard parameters from densityfunctional perturbation theory[END_REF] based on the density-functional perturbation theory to computed the on-site V II (i.e. U ) and the inte-site V IJ between the iron site, I, and the six first neighbouring atoms, J, as implemented in the hp.x utility of Quantum ESPRESSO. Because of the Jahn-Teller distortion at the high spin state, the computed inter-site V IJ adopt three different values corresponding to the three M-L lengths. The computed Hubbard parameters used for the DFT+U sc and DFT+U sc +V sc calculations are reported in Tab. 6.4 and 6.5.

Effect of the geometry on the ∆E H-L

Tabs. 6.6 and 6.7 report the adiabatic energy differences, ∆E H-L , computed for [Fe(NH 3 ) 6 ] +2

and [Fe(CO) 6 ] +2 using 11 different functionals. For each functional we use five geometries (PBE+U, PBE, TPSSh, B3LYP and PBE0). For HS, the Jahn-Teller distorsion yields non-identical metal-ligand (M-L) bond lengths. Tab. 6.8 reports the (averaged) M-L distances for (HS) LS at different geometries. For this specific study, the PBE[U] calculations are performed using atomic projectors because the ortho-atomic ones yiels a slower convergence in the case of HS. For each geometry, we compute the corresponding U sc (see Tab. 6.9). The PBE+U geometries are computed using a structurally consistent U and are taken from our previous work [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF].

6.4 Supplementary information

Introduction

Metal-organic frameworks (MOF) are nanoporous materials made of metal ion centers coordinated by organic linkers to for 3D-structures [START_REF] Li | Design and synthesis of an exceptionally stable and highly porous metal-organic framework[END_REF][START_REF] Wang | Metal-organic frameworks as a tunable platform for designing functional molecular materials[END_REF]. The wide range of possible combinations of linkers and metal ions has allowed for the design of MOFs tailor-made for different technological applications [START_REF] Yuan | Stable metal-organic frameworks: Design, synthesis, and applications[END_REF] like catalysis [START_REF] Lee | Metal-organic framework materials as catalysts[END_REF][START_REF] Huang | Multifunctional metal-organic framework catalysts: synergistic catalysis and tandem reactions[END_REF][START_REF] Herbst | MOF catalysts in biomass upgrading towards value-added fine chemicals[END_REF][START_REF] Bagheri | A MoO3-metal-organic framework composite as a simultaneous photocatalyst and catalyst in the PODS process of light oil[END_REF], drug delivery [START_REF] Horcajada | Metal-organic frameworks in biomedicine[END_REF], sensing [START_REF] Bagheri | A MoO3-metal-organic framework composite as a simultaneous photocatalyst and catalyst in the PODS process of light oil[END_REF][START_REF] Gamonal | Divergent adsorption-dependent luminescence of amino-functionalized lanthanide metal-organic frameworks for highly sensitive NO2 sensors[END_REF] and gas storage and separation [START_REF] Li | Recent advances in gas storage and separation using metal-organic frameworks[END_REF][START_REF] Li | Selective gas adsorption and separation in metal-organic frameworks[END_REF]. In these materials, efficient and selective gas adsorption can be achieved due to the large internal surface area and the high affinity of the metal centers with certain gas molecules [START_REF] D'alessandro | Carbon dioxide capture: Prospects for new materials[END_REF]. Specifically, the presence of coordinatively unsaturated metal centers [START_REF] Hwang | Amine grafting on coordinatively unsaturated metal centers of mofs: Consequences for catalysis and metal encapsulation[END_REF], also called open-metal sites (OMS), has resulted in exceptional gas uptake and separation performances [START_REF] Sumida | Hydrogen storage and carbon dioxide capture in an iron-based sodalite-type metal-organic framework (fe-btt) discovered via high-throughput methods[END_REF][START_REF] Herm | Hydrocarbon separations in metal-organic frameworks[END_REF] owing to the formation of selectively strong [START_REF] Lee | Design of a BIBLIOGRAPHY 157 metal-organic framework with enhanced back bonding for separation of n2 and ch4[END_REF]1] coordination bonds between the MOF and the guest molecules [START_REF] Chui | A chemically functionalizable nanoporous material [cu3(TMA)2(h2o)3]n[END_REF][START_REF] Kökçam-Demir | Coordinatively unsaturated metal sites (open metal sites) in metal-organic frameworks: design and applications[END_REF]. Recently, Long and coworkers reported the synthesis of Fe-based OMS-MOFs [START_REF] Reed | Reversible CO scavenging via adsorbate-dependent spin state transitions in an iron(II)-triazolate metal-organic framework[END_REF] exhibiting a cooperative carbon monoxide adsorption mechanism via spin-transition, resulting in unprecedented large working capacities [START_REF] Reed | A spin transition mechanism for cooperative adsorption in metal-organic frameworks[END_REF]. In their work, a cooperative spin-state transition occurs at the Fe centers from high spin (HS) to low spin (LS) as a function of the CO partial pressure, as a result of the strong ligandfield induced by CO upon adsorption. Once the HS→LS is triggered at all Fe centers, the uptake is shown to increase abruptly as a consequence of the significantly larger CO binding energy expected at LS compared to HS [START_REF] Sumida | Hydrogen storage and carbon dioxide capture in an iron-based sodalite-type metal-organic framework (fe-btt) discovered via high-throughput methods[END_REF]. One clear limitation of the reported mechanism is its application to molecules exhibiting a quite strong ligand field so that a large ligand-field splitting, upon binding, can stabilize the LS phase. CO is located at the far right of the spectrochemical series and it exhibits the largest ligand-field strengths, together with CN -. Thus, for most of the gas molecules relevant for industrial process this mechanism cannot be expected to hold.

In this work, we propose design principles for the development of an alternative strategy that exploits a spin-state transition induced by temperature, rather than the gas pressure, to induce a change in the gas adsorption properties. As first demonstrated by some of us using first-principles [1], then confirmed on other MOFs [START_REF] Mohajeri | Spin crossover as an efficient strategy for controllable gas molecule capturing on open metal sites in ni-btc and cu-btc[END_REF] also using DFT calculations, and shown experimentally for the case of CO later [START_REF] Reed | A spin transition mechanism for cooperative adsorption in metal-organic frameworks[END_REF], large differences in the gas adsorption energy may be found between electronic configurations at the transition-metal atom involving a change in the occupation of the antibonding d z 2 orbital between the gas and MOF. Thus, efficient adsorption and desorption cycles could be achieved by using temperature to change the spin state of the metal center and release the gas upon weakening the coordination bond. In this work, the proof-of-concept for this mechanism is provided Engineering Metal-Organic Frameworks for Spin Crossover-assisted Gas Release: a First-Principles Study MOF Fe(bpac)Fe(CN) 4 are performed using Quantum ESPRESSO [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF] with GBRV ultrasoft pseudopotentials [START_REF] Garrity | Pseudopotentials for high-throughput dft calculations[END_REF][START_REF] Bennett | A systematic determination of hubbard u using the gbrv ultrasoft pseudopotential set[END_REF] and wavefunction and charge density cutoffs of 60

Ry and 600 Ry, respectively. A 2×2×1 Monkhorst-Pack grid is used. The calculations are performed by optimizing the geometry of the bare MOF (without adsorbed gas) at each spin state, starting from the structure already published [79]. For this we employ PBE+U with a Hubbard U computed self-consistently [START_REF] Kulik | Density functional theory in transition-metal chemistry: A self-consistent hubbard u approach[END_REF] using linear-response [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF],

called here U sc . These U sc are then employed for all the subsequent calculations including those with gas molecules. The U sc for the octahedrally-coordinated Fe atom is 8.9 eV and 8.8 eV for LS and HS, respectively; for the Fe atom at the square planar environment U sc is 9.7 eV and 9.1 eV for LS and HS. Ortho-atomic projectors are employed. The binding geometries of CO, CO 2 , N 2 , and H 2 are then computed by optimizing the position of the molecules while keeping the coordinates of the MOF fixed at the bare-MOF case. For these calculations, similar to the molecular case, PBE+U sc is used in conjunction with the Grimme-D3 method [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-d) for the 94 elements h-pu[END_REF] used together with the BJ damping scheme [START_REF] Smith | Revised damping parameters for the d3 dispersion correction to density functional theory[END_REF]. As reported in previous works [4,[START_REF] Rodríguez-Velamazán | Long-range magnetic order in the porous metal-organic framework ni(pyrazine)[pt(CN)4][END_REF], the Hubbard correction is needed in order to correctly describe the transition-metal 3d electrons in MOFs. The different spin states are obtained imposing an initial magnetization and allowing the system to converge to the local energy minimum which corresponds to the desired electronic configuration. The calculations of the adiabatic energy differences are performed by employing the Hubbard U -corrected density scheme. We recently demonstrated that this method allows to compute adiabatic energy differences of Fe(II) molecular and crystalline complexes with excellent agreement compared with coupled cluster-corrected CASPT2 and with experimentally-extracted values [START_REF] Mariano | Improved spin-state energy differences of fe(ii) molecular and crystalline complexes via the hubbard u-corrected density[END_REF]. It consists in performing total energy calculations for each spin state using PBE+U sc including the Grimme-D3+BJ method such that:

E PBE+U+D3 [ρ] = E PBE+D3 [ρ] + E U [n] (7.1)
where ρ is the PBE+U electronic density and E U is the Hubbard energy which depends on the density, ρ, through the projections of the Kohn-Sham orbitals onto a selected basis set, n [START_REF] Cococcioni | Linear response approach to the calculation of the effective interaction parameters in the lda+u method[END_REF]. The adiabatic energy difference between HS and LS, for example, is computed as follows:

∆E H-L [ρ] = E HS PBE+D3 [ρ] -E LS PBE+D3 [ρ], (7.2)
i.e. we keep only the PBE+D3 term of the DFT+U total energy computed on the U density, ρ. By removing the E U HS -E U LS from the adiabatic energy difference, we remove the bias of DFT+U towards HS states, which arises from to the penalizing Hubbard term being larger for LS [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. The binding energy, E bind , is computed by substructing from

Results

111

the total energy of the bound MOF (or fragment)-gas complex, the total energy of the constituents [START_REF] Poloni | Co 2 capture by metal-organic frameworks with van der waals density functionals[END_REF]. To quantify the efficiency of the gas release upon spin crossover we compute the change in binding energy associated with such transition. This change, upon a LS→HS transition, is simply computed as:

∆E bind,H-L = E bind,HS -E bind,LS (7.3)

Results

Fragment calculations

The calculations of the ∆E H-L and the E bind together with the gas-metal bond distance and the electrostatic charge of the metal are is reported in Tab. LS being the ground state. Ideally, in order for a spin crossover-induced gas-release process to be energy-efficient the total energy difference between the two spin states should not be

Results

113

too large. More specifically, the ∆E bind,H-L should be maximized while the ∆E H-L should be the lowest possible to reduce the energy penalty, but consistent with the conditions of the capture process. For example, if we assume that adsorption is performed at room temperature, the binding energy should be the largest at this temperature, and efficient release could be achieved by heating. However, the higher the transition temperature, the larger is the energy penalty associated with the process. exception is CO adsorbed on Co(II). In this case, a significant hybridization between the π-like t 2g 3d metal orbitals and the CO LUMO is observed when the system is in S=1/2.

Such a π-backbonding interaction is weaker when the Co(II) complex is in S=3/2 and the small change in binding in this case is possibly due to the cancellation of the two effects, the π-backbonding and electrostatic Co(II)-CO interaction.

We want to stress that for Fe(II) and Fe(III) the ∆E bind,H-L refer to energy differences between the highest possible spin state and the intermediate spin state. Thus, because the intermediate spin states are already quite high in energy (see Tab.

MOF calculations

The calculations on the fragments show promising results when the first mechanism is computed for the Fe(II), Fe(III) and Ni(II) fragments. In this case the change in binding energy of the adsorbed molecules is significant and the associated ∆E bind,H-L are below 1 eV, except for CO. We recall that for Ni(II), HS corresponds to a stronger interaction with gas molecules because the d z 2 is singly occupied while in LS it is doubly occupied. Thus, in this case gas release is achieved upon cooling down the system. Because the change in binding is lower for Ni(II) (see Figure 7.1) and because thermal effects cannot contribute to further decrease uptake like in standard temperature swing processes, we apply these findings to the Hoffman-like MOF Fe(bpac)M(CN) 4 where M(II) is replaced by Fe(II).

Although interesting, Fe(III) cannot be studied in this MOF.

We compute the ∆E H-L and ∆E bind,H-L between S=1 and S=0 for Fe(bpac)Fe(CN) 4

. The HS S=2 cannot be stabilized in our calculations without imposing a constraint on the spin, suggesting a configuration unlikely reachable upon heating. The octahedrally coordinated Fe is considered in S=2 and S=0, respectively, when the open metal is in S=1 and S=0. This means that at low temperature both Fe atoms are in LS and we assume that the spin transition of the octahedrally coordinated Fe occurs at lower temperature compared to the open-metal site Fe. We have checked nonetheless the effect of the changing the spin state of the octahedral Fe on the computed binding energy and we find an almost negligible difference. The results are shown in Tab. 7.3. We note that the adiabatic energy differences and the change in binding energy somehow differ from the molecular case. This difference results to some extent from the use of different computational methods (basis and HS (S=1) for several gas molecules for Fe(bpac)M(CN) 4 . The ∆E H-L (eV) is also reported.

Discussion and conclusion

To precisely determine the thermodynamics of the spin crossover both the enthalpy and the entropy changes need to be evaluated. Thus, in order to extract the transition temperature, besides the adiabatic energy differences reported above (which are the electronic contribution to the enthalpy change), the calculations of the zero-point energy, the vibrational enthalpy and entropy, and the electronic entropy are also required to compute the free energy change, ∆G=∆H-T∆S [START_REF] Cirera | Theoretical prediction of spin-crossover temperatures in ligand-driven light-induced spin change systems[END_REF]. Despite this fact, some considerations can be made to assess the feasibility of the process starting from the calculation of the (electronic) enthalpy change. Enthalpy and entropy act oppositely: the entropic term is larger for HS due to the larger vibration and electronic contributions (and ZPE term), and, as the temperature increases, this term stabilizes HS more compared to LS, which is the ground state and has larger enthalpy.

We recall that the entropic term in the ∆G is typically well described by standard 7.4 Discussion and conclusion 119 DFT [START_REF] Kepp | Theoretical study of spin crossover in 30 iron complexes[END_REF] while the electronic term is the most challenging one. Motivated by our previous studies showing an improved performance of the Hubbard U -density corrected approach compared to other good performers such as M06-L and TPSSh, we employ this method to accurately compute the ∆E H-L for a series of molecular fragments that can potentially be used in the skeleton of a well-known MOF family. The results obtained from this stage are leveraged to design a Hoffman-type MOF with a Fe(II) in square-planar coordination. In order to comment upon the feasibility of the proposed mechanism it is important to stress that the vibrational contribution to the enthalpy change amounts to a few tens of meV [START_REF] Kepp | Theoretical study of spin crossover in 30 iron complexes[END_REF][START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF]. Thus, if we assume that the computed ∆E H-L reported above well approximate the total ∆H of the spin-state transition, a comparison can be made with experimental values of ∆H reported for spin-crossover compounds. Values of ∆H ranging between 0.1 and 0.3 eV have been reported for several Fe(II)-based molecular [START_REF] Kepp | Theoretical study of spin crossover in 30 iron complexes[END_REF] and crystalline compounds [START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF] resulting in spin-crossover temperature near room temperature. The computed ∆E H-L for Fe(bpac)Fe(CN) 4 upon adsorption of H 2 and CO 2 being 0.298 eV and 0.178 eV, respectively, thus suggest a promising performance of this material and provide the proof-of-concept for the proposed mechanism of spin-crossover induced gas release.

We hope the design principles provided here will assist future computational studies and we are confident that our results for the MOF will motivate future synthetic efforts towards the development of energy-efficient spin crossover-assisted separation technologies.

we find that the MAE for PBE[HF] is significantly larger compared to PBE[U], especially for strong-ligand complexes. As suggested by other authors [START_REF] Radoń | Revisiting the role of exact exchange in DFT spin-state energetics of transition metal complexes[END_REF][START_REF] Pinter | Conceptual insights into DFT spin-state energetics of octahedral transition-metal complexes through a density difference analysis[END_REF][START_REF] Polo | Some thoughts about the stability and reliability of commonly used exchange?correlation functionals ? coverage of dynamic and nondynamic correlation effects[END_REF] and confirmed by our density analysis of LS states, DFT is able to mimic through the SIE some staticcorrelation effects that are, on the contrary, completely missed in HF theory. Nevertheless, the DFT ability to capture some static-correlation effects on the electronic charge density is fortuitous and, in general, exaggerated. We have shown that in PBE+U the electronic density is less delocalized in the spatial region around the metal center and it resembles more the relaxed CASPT2 density, compared to PBE. When PBE+U is compared to HF, the latter method seems to give better electronic densities for weak-ligand complexes and this is consistent with the results of the energetics. However, our analysis is restricted to closed-shell electronic configurations and further investigation on the HS electronic charge distribution is needed to corroborate and confirm the present findings.

In chapter 7 the PBE[U] method has been applied in the calculation of spin energetics of OMS upon gas adsorption in the SCO Hoffman-type MOF. and its adsorption properties upon spin transition. Relatively large changes in binding energies are obtained for all the molecules investigated while values of adiabatic energy differences compatible with a room-temperature spin transition are obtained only for H 2 and CO 2 . In this study we did not calculate the entropic effects that are needed to compute the transition SCO temperature. An accurate analysis of the vibration properties of these systems is needed to calculate the final transition-temperature. A collaboration with

Prof. Lin at the Ohio State University is also ongoing to compute the working capacity using grand canonical Monte Carlo, using our calculations as inputs. This will allow our group to compute the working capacity and to quantitatively assess the efficiency of the proposed mechanism.

In conclusion, my work shows that an SCO-assisted gas adsorption/release strategy is feasible for a suitable combination of OMS and guest molecules. We proposed the Hoffman-type Fe(bpac)Fe(CN) 4 MOF whose synthesis represents an active research project of our experimental collaborators José Sanchez Costa and José Alberto Rodríguez Velamazán from IMDEA Nanociencia (Madrid, Spain) and Institut Laue-Langevin (Grenoble, France).

The experimental realization of this MOF is very challenging because of the presence of square-planar coordinated iron OMS. Yet, we are confident that the proof-of-concept provided in this work will be validated experimentally by our collaborators. Our PBE[U] approach has been shown to represent a reliable method to study Fe(II) spin energetics in crystalline and molecular systems. Nevertheless, the analysis of PBE[U] is restricted to Fe(II) complexes and further investigations are needed to extend our methodology to other transition metal compounds.

Figure 2 . 1 :

 21 Figure 2.1: Schematic representation of the 3d electronic degenerates levels for an isolated ion (left) and for the partially degenerate 3d configuration resulting from an octahedral ligand field (right). ∆E spherical is the energy loss due to the symmetry change and ∆ is the ligand field splitting.

Figure 2 . 2 :

 22 Figure 2.2: Simplified Tanabe-Sugano diagram. The energies of the iron spin states in an octahedrally ligand-field environment (y-axes) are plotted as a function of the ligand field splitting ∆.

Figure 2 . 3 :

 23 Figure 2.3: Schematic representation of the potential energy profile E p for LS (blue) and HS (red) as a function of the metal-to-ligand distance r. r LS and r HS are the equilibrium values of r in LS and HS. The zero-point energy (ZPE), the internal energy difference at zero temperature ∆E H-L , and the adiabatic energy difference ∆E H-L are shown. A schematic representation of the Fe(II) 3d electronic configuration in LS and HS is also shown.

Figure 2 . 4 :

 24 Figure 2.4: Schematic representation of the Jahn-Teller distortion in the case of Fe(II) HS. The degeneracy in the symmetric octahedral coordination environment (left) is lifted lowering the symmetry through a compression (right). The change in the structure is

  into dynamical and static (non-dynamical). The former is due to the instantaneous Coulombic interactions and can be captured by adding multiple excited Slaterdeterminants to the HF single-reference wavefunction. This can be done by employing methods such as configuration interaction (CI), coupled-cluster theory (CC), or manybody perturbation theory (MBPT). On the other hand, static correlations arise from the energetically degenerate character of the ground-state electronic wave function and methods like multiconfigurational complete active space self-consistent field (CASSCF) or diffusion Monte Carlo (DMC) are in principle able to capture their contribution to the total energy. Fe(II) complexes are generally divided into weak field ligand and strong field ligand complexes according to the sign of the adiabatic energy difference ∆E H-L . Particularly for strong field ligand complexes, multireference methods are required to correctly evaluate ∆E H-L . Here, most of the static correlation is due to the metal-to-ligand interactions.

Figure 3 . 1 :

 31 Figure 3.1: Schematic representation of the synthesis of MOF-5 [61] from zinc metal centers and terephthalic acid. Yellow spheres inside the MOF structure are drawn to give a qualitative idea of the void space within the framework. Figure reproduced from ref.

Figure 3 . 2 :

 32 Figure 3.2: Representation of the honeycomblike pore structure of MOF-74 upon CO 2 adsorption on the open-metal sites (left). Magnification of the interaction between the gas and the open-metal site is also shown (right). Figure reproduced from ref. [4].

Figure 3 . 3 : 1 )

 331 Figure 3.3: 1) CO 2 /N 2 adsorption isotherms for different functionalized Mg 2 (dondc) MOFs: a) unfunctionalized MOF; b) en; c) mmen d) ppz. Figure reproduced from ref. [71]. 2) Carbon monoxide isotherms at 298 K • for M 2 (dondc) with M=Ni, Co, Fe, Mg, Mn, Zn. Figure reproduced from ref. [72].

Figure 3 . 4 :

 34 Figure 3.4: Schematic representation of adsorption/desorption process with Langmuirtype isotherms. The adsorption occurs at P ads and T low , the system is heated up reaching gradually the desorption point characterized by a temperature T high and a pressure P des .

3. 3 27 Figure 3 . 5 :

 32735 Figure 3.5: Schematic representation of the CO cooperative adsorption process induced by spin-transition. a) Non-cooperative Langmuir-type adsorption isotherms (see also figure 3.4) and b) step-shaped isotherms. In the latter, a swing in temperature between T low and T med leads to a working capacity larger compared to non-cooperative processes using T low and T high . c) Change in CO uptake due to a collective pressure-induced spin-transition from high-spin to low-spin of the open-metal sites upon gas adsorption. Figure reproduced from ref.[START_REF] Reed | A spin transition mechanism for cooperative adsorption in metal-organic frameworks[END_REF]

[ 1 ] 29 Figure 3 . 6 :

 12936 Figure 3.6: Molecular diagram for a π (left) and a π-backbonding interaction (right) between the metal (M) 3d and the molecule/ligand (L) p electrons. The relative size of the metal and ligand atomic orbitals reflects the contribution to the MO.

Figure 3 . 7 :

 37 Figure 3.7: A) Relationship between the computed electrostatic contribution (Q/r) and CO 2 binding energy in several M-BTT MOFs. B) Schematic view of the change in the V 3d electronic configuration due to the interaction with CO 2 and plot of the CO 2 -metal antibonding orbital resulting from the hybridization between the CO 2 lone pair and the metal d z 2 . Figures reproduced from ref. [1].

Figure 3 . 8 :

 38 Figure 3.8: Ligand-field splitting and electron configuration for vanadium interacting πbackbonding with the π-like (in figure p-like) orbital of N 2 . Figure reproduced from ref. [4]. B) First coordination sphere of Fe in Fe-MOF-74 upon N 2 and CH 4 adsorption.

Figure 3 . 9 :

 39 Figure 3.9: Crystal structure of Fe(bpac)Fe(CN) 4 with CO 2 adsorbed at the apical sites of Fe OMS . Fe closed-metal site (purple), Fe open-metal site (yellow), N (blue), C (grey), hydrogen atoms are omitted for clarity. The polyhedras represent the first coordination of the closed Fe metal sites (purple) and the Fe OMS sites (yellow).

(4. 6 )Z 2 I

 62 I represents the mass of the ion I while m and -e are the mass and the charge of one electron. Eq. 4.6 is the electronic Schrödinger equation within the Born-Oppenheimer approximation. It is the result of the total wavefunction factorization in eq. 4.4 and of the assumption that the action of ionic kinetic termon the electronic wavefunction ψ R is negligible. The time-dependent Schrödinger equation 4.1 is now greatly simplified but still, in most of the situations, the number of degrees of freedom of the problem is too big to have any analytical solutions and further approximations are needed.

  ,... c r,s,... a,b,... |ψ r,s,... a,b,...

. 48 )

 48 Different approximates expressions for the exchange-correlation functional have been proposed over the years. Below we give some examples. • Local Density Approximation (LDA). It is the simplest approximation to the exchange-correlation functional where the exchange-correlation energy E LDA xc [ρ] is evaluated at each point of the space r assuming an homogeneous gas approximation for the electronic density, i.e.: E LDA xc [ρ] = ε hom xc (ρ(r))ρ(r)dr (4.49)

Figure 4 . 1 :

 41 Figure 4.1: Energy as a function of the number of electrons for exact straight line behavior, approximate convex local and semilocal DFT energy functionals and overlocalized concave HF energy. In this, hypothetic, case the three methods give the same electron affinity (A) and ionization potential (I). Figure reproduced from ref. [117].

  show a systematic overstabilization of HS states by DFT+U sc with deviations from the reference CASPT2/CC values reaching up to 4.5 eV for strong-field ligands. Such an overstabilization arises from the penalizing Hubbard energy term in the total energy which is systematically larger in the LS state compared to HS. This effect has two origins: both the summation term that multiplies U in the Hubbard energy term (vide infra) and U sc are always larger in LS. We here explain in detail the origin of the first one. We show that this failure of DFT+U sc can be mitigated by empirically choosing U values that reduce this bias. The density change upon application of U is here rationalized in terms of the change in metal/ligand hybridization effectively resulting in a metal-ligand Biased spin-state energetics of Fe(II) molecular complexes within density functional theory and the linear-response Hubbard-U correction

Figure 5 . 1 :

 51 Figure 5.1: Adiabatic energy differences, ∆E H-L , computed with different approaches and the difference between the E U term of LS and HS computed using PBE+U sc (left panels). The CASPT2/CC reference values are -1.99 eV, -0.72 eV, -0.24 eV, 2.31 eV, 1.97 eV and 2.82 eV, from weak to strong-field ligands. Molecular complexes studied in this work (right figure).

5. 2

 2 Main paper 57 demonstrate below, this agreement is fortuitous and results from a cancellation of errors: the LS overstabilization by LDA is compensated for by the HS overstabilitation by the U -corrected LDA. For strong-field ligands (i.e. for ∆E CASPT2/CC >0) the best agreement is achieved with PBE (MAE=1.115 eV) and the deviation of DFT+U sc from the reference CASPT2/CC increases up to 4.5 eV (see Figure 5.1 and

Figure 5 . 3 :Figure 5 . 1 .

 5351 Figure 5.3: Left panel: molecular-orbital diagram built from both σ and π-bonding together with the corresponding Kohn-Sham orbitals for Fe(CNH) 2+ 6 . Right panel: the m,σ [n σ m (1 -n σ m )] term is plotted for the whole molecular series for increasing strength of the ligand-field and separately for each spin states, σ, and orbital symmetry. The U values (eV) are also reported.

5. 2

 2 Main paper 59 the CASSCF calculation into their corresponding atomic compositions (Tables 5.8-5.12). These trends yield a negative difference of the U energy term between HS and LS, i.e. ∆E U,H-L =E U HS -E U LS , that largely increases as a function of the ligand-field strength, as shown in the lower panels of Figure 5.1. This explains the increasing bias towards HS states for ligands on the right side of the spectrochemical series. As a consequence, the DFT+U-computed adiabatic energy differences are always negative and their dispersion across the series no longer reflects the trend expected for varying ligand field strengths as predicted by CASPT2/CC. While PBE and LDA do follow the CASPT2/CC behavior, the Hubbard-U corrected calculations show no correlation with the reference set. In what follows we explain precisely the origin of the systematically larger summation term for LS than HS, i.e. why the occupations {n σ m } associated to both the t 2g -and e g -like states are systematically closer to 1/2. Since the {n σ m } occupation numbers results from the projection of the occupied Kohn-Sham states, only, onto the atomic basis {|φ m }, fractional values of the {n σ m } occurr when (1) the d orbitals of Fe hybridize with the ligand and (2) the resulting M/L hybridized molecular orbitals are only partially occupied. When all of these molecular orbitals are occupied, even though there is strong hybridization, the {n σ m } are close to 1 (and zero when empty) thus yielding a low m,σ [n σ m (1 -n σ m )]. As illustrated in the left panel of Figure 5.3, the two sets of molecular orbitals arising from the hybridization of the metal and ligand via the σ-bond are the bonding e g and the antibonding e * g . For LS states, the bonding (ligand) e g are occupied while the e * g states, with larger metal-like character, are empty (see also projected density of states in Figures 5.11-5.22

Figure 5 . 4 :

 54 Figure 5.4: Adiabatic energy differences, ∆E H-L , (upper panel) and energy difference of the Hubbard term, ∆E U,H-L , computed at different values of U , by using PBE.

5 .Figure 5 . 5 :

 555 Figure 5.5: Molecular orbital diagram illustrating the change in M/L hybridization upon U correction and difference-plot of the e * g and t * 2g molecular orbitals between U =0 and U sc for LS (upper panel) (red and green correspond to positive and negative isovalues). L→M charge transfer computed from the Löwdin charges for the six molecular complexes using PBE+U sc (lower panel).

Fe 3 63 Figure 5 . 5

 36355 Figure 5.5 and left panel of Figure 5.3). Since the corresponding t *2g are unoccupied, this results in a ligand→metal charge transfer. Thus, the e * g and t * 2g sets respectively increase and decrease their d-character upon U correction as illustrated in Figure5.5 (Figures 5.9

and 5 .

 5 [START_REF] Stoufer | UNUSUAL MAG-NETIC PROPERTIES OF SOME SIX-CO ÖRDINATE COBALT(II) COM-PLEXES1-ELECTRONIC ISOMERS[END_REF]. As the ligand-field strength increases, both effects increase in magnitude with the L→M charge transfer becoming larger than the M→L, due to bigger change, upon application of U , in hybridization of the t 2g set compared to the e g one (Figures 5.11-5.22). This is demonstrated by the change in occupation numbers {n σ m } upon U , shown in Tables 5.13 and 5.14. The M↔L charge transfer computed for LS thus changes its sign along the molecular series (Figure 5.5).

  This work was funded by the French National Agency for Research (ANR-15-CE06-0003-01). Calculations were performed using resources granted by GENCI under the CINES grant number A0020907211. The froggy platform of the CIMENT infrastructure was also theory and the linear-response Hubbard-U correction employed for the calculations. B.V. was supported by the Department of Energy, Basic Energy Sciences (DE-SC0019463). Some of the computations supporting this project were performed on High Performance Computing systems at the University of South Dakota, funded by NSF Award 1626516. Also resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231, were used for this project.

5 :

 5 Deviation of the UHF expectation value < S 2 > from the ideal value S(S +1) = 6 for the quintet states (S = 2).

Figure 5 . 6 : 5 . 4 information 69 Figure 5 . 7 :Figure 5 . 8 :

 5654695758 Figure 5.6: CASPT2 optimized orbitals for LS (left) and HS (right) state configuration of [Fe(NH 3 ) 6 ] +2 .

5. 4 information 71 Figure 5 . 9 :

 47159 Figure 5.9: Electron density difference between PBE+U sc and PBE for [Fe(NCH) 6 ] 2+ . From left to right: HS majority spin, HS minority spin, and LS (maj/min spin). Red and green surfaces represent negative (charge depletion) and positive (charge accumulation) values, respectively, with an isovalue of 0.001 e/Bohr 3 .

Figure 5 . 10 :

 510 Figure 5.10: Electron-density difference between PBE+U sc and PBE for [Fe(CNH) 6 ] 2+ . From left to right: HS majority spin, HS minority spin, and LS (maj/min spin). Red and green surfaces represent negative and positive values, respectively, with an isovalue of 0.001 e/Bohr 3 .

Figure 5 . 13 :

 513 Figure 5.13: Projected density of states of [Fe(NH 3 ) 6 ] 2+ in LS (left) and HS (right) computed using PBE. Solid black line indicates the HOMO level.

Figure 5 . 18 :

 518 Figure 5.18: Projected density of states of [Fe(PH 3 ) 6 ] 2+ in LS (left) and HS (right) computed using PBE+U sc . Solid black line indicates the HOMO level.

-

Figure 6 . 1 :

 61 Figure 6.1: Adiabatic energy differences, ∆E H-L , computed using varying DFT approaches, together with the reference CASPT2/CC energies computed in this work. The values are also reported in Tab. 6.1 for clarity.

Figure 6 . 2 :

 62 Figure 6.2: Density difference plot, δρ x (r), for [Fe(CO) 6 ] +2 between x=PBE, PBE+U sc , and Hartree-Fock and the relaxed CASPT2 density; green and blue correspond to positive and negatives values, respectively. The plot shows values between -0.005 e/bohr 3 and 0.005 e/bohr 3 .

6. 2

 2 Main paper 91 ing Kohn-Sham inversion schemes to show that the Kohn-Sham kinetic energy and the

Figure 6 . 3 :

 63 Figure 6.3: Left figure: charge density difference, ρ F e,Baderx

) 6 ]

 6 2+ , [Fe(NCH) 6 ] 2+ , [Fe(CO) 6 ] 2+ and [Fe(CNH) 6 ] 2+ . In Figure6.2 we plot the difference δρ x (r) = ρ x (r) -ρ CASPT2 (r) between the electronic density obtained with x=[PBE, PBE+U, HF], and the CASPT2 relaxed density, for [Fe(CO) 6 ] 2+ . The same qualitative result is obtained for NCH and the corresponding plots are reported in Figure6.5. A limitation of this analysis is that large density differences are found in the spatial region near the ligand for PBE and PBE+U, while negligible ones are found for the HF density, consistent with a CASSCF active Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density space mostly involving states associated with the Fe, and only marginally associated with the ligand, i.e. the two e g ligand states. Thus, the CASPT2-relaxed density resembles closely the Hartree-Fock one near the ligand, which is the reference method used to get the CASSCF wavefunction. Due to this limitation, in what follows we limit our considerations to the spatial region near the iron. When PBE density is used, the δρ PBE (r) is negative within the spatial region associated with the e g orbitals and positive within for the t 2g

2 )

 2 The charge density difference, ρ F e,Bader x (r)-ρ F e,Bader CAS (r), is plotted in Figure6.3 for a weak and a strong-field ligand case, i.e. NCH and CO. For each molecule the opposite behavior of HF and PBE density is clearly visible, together with a reduced deviation from the reference density of PBE+U sc , in agreement with Figure6.2 and the above considerations.The ∆ρ x computed for the four molecules is plotted in the right panel of Figure6.3. This metric gives a constant error throughout the molecular series for PBE. The error associated with the HF density is smaller compared to PBE, and significantly smaller for weak-field molecules. This result is consistent with the reasonable prediction of ∆E H-L found when employing PBE[HF] for weak-field ligands and with the larger MAE of PBE for both weakand strong-field ligands (see Tab. 6.1). The ∆ρ x increase for molecules with increasing ligand-field strengths when x=HF and the opposite is found for x=[PBE+U, PBE+U+V].

Figure 6 . 4 :

 64 Figure 6.4: Adiabatic energy differences, ∆E H-L , computed with PBE[U] as function of Hubbard U . On the right y-axis the CASPT2/CC reference values are shown by horizontal vertical dashes.

2 -

 2 yl)-1,2,4-triazole, Fe[HB(pz) 3 ] 2 [247] (pz=pyrazolyl), and FeL 2 [BF 4 ] 2 [248] (L=2,6-di(pyrazol-1-yl)pyridine), are molecular crystals for which Vela et al. [126] have extracted the experimental ∆E H-L by removing the (computed) vibrational contribution from the the measured total enthalpy change. The other two are molecular complexes, [Fe(tacn) 2 ] +2

  7.1 for the fragment model M(CN) 4 [M=Fe(II), Fe(III), Co(II), Ni(II)] for CO 2 , CO, N 2 , H 2 , and H 2 O. For Fe(II) and Fe(III) we compute three spin states while for Co(II) and Ni(II) only two spin states are possible. In order to facilitate the understanding of the results, in what follows we divide them in two groups. Each group reflects a specific mechanism that we have identified here as responsible for the change in binding energy upon spin-state transition.The first group involves situations where the change in binding energy mainly results from the change in occupation of the metal-gas molecular orbital d z 2 . Consistent with the analysis reported in Ref.[1], the binding energy always significantly decrease (increase) upon occupation (depopulation) of the antibonding d z 2 , upon spin transition. The second mechanism, largely overlooked until now, involves a change in the occupation of the d x 2 -y 2 molecular orbital, which has an antibonding character between the metal and the in-plane CN -ligands. The latter will be discussed in detail below.Change in the occupation number of the d z 2 . Because the d z 2 has an antibonding character along z, the spin crossover-induced change in the occupation of this molecular orbital yields a change in the bond order of the fragment+gas complex[1]. This case applies to the change between S=1 and S=0 for Fe(II); between S=3/2 and S=1/2 for Fe(III); and between S=1 and S=0 for Ni(II). See Tab. 7.2 for the electronic configurations of the different molecules. The change in binding energy for the five molecules upon the specified spin-state transition is shown in Figure7.

Figure 7 . 1 :Figure 7 . 2 :Figure 7 . 3 :

 717273 Figure 7.1: Computed values of ∆E bind,H-L and ∆E H-L for a spin-state transition involving a change in the occupation of the antibonding metal-gas molecular orbital, d z 2 , (left panels)and the antibonding metal-ligand molecular orbital, d x 2 -y 2 (right panels). The left panels report values for spin-state changes between S=1 and S=0 for Fe(II), between S=3/2 and S=1/2 for Fe(III), and S=1 and S=0 for Ni(II). The right panels report values for spin-state changes between S=2 and S=1 for Fe(II), between 5/2 and 3/2 for Fe(III), and between 3/2 and 1/2 for Co(II). Note that the ∆E bind,H-L are reported with opposite sign for Ni(II) for clarity. For this complex, the binding energies increase upon spin crossover (see text).

7 . 1 )

 71 compared to the ground state, especially for Fe(II), one should expect the transition to high spin to occur at temperatures too high to be relevant for industrial processes. For Co(II) the ∆E H-L refer to the only two spin states accessible but the associated ∆E H-L are also quite large. These considerations, together with the relatively lower change in binding energy computed for this mechanism compared to the previous, allow us to conclude that the most promising scenario is offered by the first set of calculations i.e. a spin transition involving a change in the occupation of the d z 2 .

  Firstly, we have extensively analyzed the effect of spin transition on the computed binding energies of five molecules, i.e. CO, CO 2 , H 2 , H 2 O and N 2 , considering a small fragment M(CN) 4 (M=Fe(II), Fe(III), Co(II), Ni(II)) using the TPSSh DFT functional. This analysis has been used to perform a screening study on the capabilities of this MOF building block for application in SCOassisted gas adsorption/release strategies. The main effect on the binding mechanism is because of the change in the d z 2 metal orbitals which results in a lowering in the bond order between the metal and the molecule and thus leading to large changes in the computed binding energies. However, this study has also highlighted the effect of the change in the population of the d x 2 -y 2 metal orbital on the computed binding energy. This orbital has antibonding character between the metal and the in-plane CN -ligands and, when it is occupied, it results in an electron charge transfer from the metal to the ligand. The resulting partial charge of the metal is reduced with a consequent increase in binding energy with the adsorbate because of the enhanced metal-molecule electrostatic interaction. Finally, we have used the PBE[U] method to investigate the Hoffman-type MOF Fe(bpac)Fe(CN) 4

  Functionalization. Different values of CO 2 uptake and CO 2 /N 2 selectivity have Selective CO adsorption has been achieved using open metal M 2 (dobdc) MOF (M=Mg, Mn, Fe, Co, Ni, Zn; dobdc 4-= 2,5-dioxido-1,4-benzenedicarboxylate), which has resulted in excellent carbon monoxide uptakes and very high CO/H 2 and CO/N 2 selectivities [72] (see figure 3.3, 2).

	been obtained by functionalizing a Mg 2 (dondc) open-metal site MOF (H 4 dondc=1,5-
	dioxido-2,6-naphthalenedicarboxylic acid) with different diamine groups, i.e. en=ethylenediamine,
	mmen=N,N' or ppz=piperazine [71] (see figure 3.3, 1).

• Metal substitution.

  .[START_REF] Römer | A review of density functional models for the description of fe(II) spin-crossover complexes[END_REF] The unknown expressions for T [ρ(r)] and E xc [ρ(r)] do not allow any practical implementation of the variational principle to E[ρ(r)] and further approximations are needed. The

problem of determining the kinetic term T [ρ(r)] can be simplified by mapping the original many-body problem into an equivalent non-interacting particle model in which the ground-state charge density corresponds to that of the actual many-body system. This method is called Kohn-Sham method and it is named after W. Kohn and L. J. Sham who developed this theory in 1965

[START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]

. The single-particle Kohn-Sham wavefunctions {ψ KS i (r)} can be obtained diagonalizing the Kohn-Sham Hamiltonian:

  .66) Biased spin-state energetics of Fe(II) molecular complexes within density functional theory and the linear-response Hubbard-U correction with accurate coupled cluster-corrected CASPT2 results. We show that DFT+U fails in correctly capturing the ground state for strong field-ligands yielding ∆E H-L that are almost constant throughout the molecular series. This bias towards high spin together with the metal/ligand charge transfer upon U correction are here quantified and explained using molecular orbital diagrams involving both σ-and π-bonding interactions. With increasing ligand-field strengths this bias also increases owing to the stronger molecular character of the metal/ligand Kohn-Sham orbitals thus resulting in large deviations from the reference larger than 4 eV. Smaller values of U can be employed to mitigate this effect and recover the right energetics.

Table 5 .

 5

	1). We note

Table 5 .

 5

1: Mean absolute error (MAE) in eV for weak-field, strong-field ligand molecules, and the total value.

  which exhibit π bonding, non-bonding and antibonding character for increasing energy, and are here referred to as t 2g , t 2g -n and t Biased spin-state energetics of Fe(II) molecular complexes within density functional theory and the linear-response Hubbard-U correction set as expected from the smaller hybridization resulting from a π-bond versus a σ-bond (Tables 5.8-5.12 and Figures 5.11-5.22). The correlation between ligand-field strength

	H-L (eV)	0 2 4	2eV 4eV 6eV 8eV Average CASPT2/CC
	ΔE
		-2
		-4 0
	U,H-L (eV)	-2
	ΔE	-4
			H 2 O	NH 3	NCH	PH 3	CO	CNH

* 2g , respectively (Figure 5.3 and Figures 5.11-5.22). Overall, the m,σ [n σ m (1 -n σ m )] associated to the t 2g set is significantly smaller than for the e g

Table 5 .

 5

	5.6)

2: Couple cluster-computed 3s3p contribution, ∆E 3s3p,CCSD(T) , CASPT2computed adibatic energy differences ( a is taken from Ref.

[START_REF] Wilbraham | Multiconfiguration pair-density functional theory predicts spin-state ordering in iron complexes with the same accuracy as complete active space second-order perturbation theory at a significantly reduced computational cost[END_REF]

), ∆E H-L,CASPT2 , taken from XMS-CASPT2 with aug-cc-pVQZ-DK basis set, and the final CASPT2/CC adiabatic energy differences obtained by summing up the first two terms. All values are in eV.

Table 5 .

 5

Table 5 .

 5 [START_REF] Cambi | Über die magnetische susceptibilität der komplexen verbindungen[END_REF]. The e g and t 2g refers to occupations associated to the σ and π molecular orbitals, respectively.

	Biased spin-state energetics of Fe(II) molecular complexes within density functional
			theory and the linear-response Hubbard-U correction
			∆E H-L (eV)
	Functional	Hubbard U (eV) H 2 O NH 3 NCH PH 3	CO CNH
		0	-1.11 0.08 1.16 2.67 3.65 4.12
		2	-1.50 -0.57 0.25 1.59 2.08 3.04
	PBE	4 6	-1.91 -1.16 -0.62 0.03 1.27 1.85 -2.28 -1.70 -1.39 -0.69 0.02 0.62
		8	-2.61 -2.17 -2.01 -1.74 -1.18 -0.59
		U sc	-2.59 -2.12 -2.11 -2.34 -2.13 -1.82
		U average	-2.47 -1.91 -1.83 -1.41 -1.10 -0.70
	LDA	0 U sc	-0.17 1.16 2.60 4.14 5.42 5.89 -2.06 -1.37 -1.45 -1.63 -1.21 -0.81
	CASPT2/CC		-1.99 -0.72 -0.24 2.31 1.97 2.82

Table 5 .

 5

			U sc (eV)
	Functional Spin State H 2 O NH 3 NCH PH 3	CO	CNH
	PBE	HS LS	6.76 6.33 6.70 6.12 6.51 7.44 7.37 7.87 8.57 9.22	6.68 9.45
	LDA	HS LS	5.81 6.57 6.72 6.54 6.83 7.37 7.53 8.68 9.56 10.29 10.48 6.91

6: Adiabatic energy differences ∆E H-L (eV) for LDA, LDA+U, PBE and PBE+U.

The values of U sc are computed within the linear response approach as decribed in the manuscript and are listed in Table

5

.7. The reference CASPT2/CC values from Table

5

.2 are also reported for comparison.

Table 5 .

 5

7: Computed values of U sc (eV) used in this work.

  Figure 5.11: Projected density of states of [Fe(H 2 O) 6 ] 2+ in LS (left) and HS (right) computed using PBE. Solid black line indicates the HOMO level.Figure 5.12: Projected density of states of [Fe(H 2 O) 6 ] 2+ in LS (left) and HS (right) computed using PBE+U sc . Solid black line indicates the HOMO level.

			Biased spin-state energetics of Fe(II) molecular complexes within density functional
								theory and the linear-response Hubbard-U correction
							Tot d-xz d-xy d-yz d-z2 d-x2y2								Tot d-xz d-xy d-yz d-z2 d-x2y2
	DOS								DOS
	-20	-18	-16	-14	-12	-10	-8	-6	-22	-20	-18	-16	-14	-12	-10	-8
				E(eV)									E(eV)
							Tot d-xz d-xy d-yz d-z2 d-x2y2								Tot d-xz d-xy d-yz d-z2 d-x2y2
	DOS								DOS
	-20	-18	-16	-14	-12	-10	-8	-6	-22	-20	-18	-16	-14	-12	-10	-8
				E(eV)									E(eV)
							Tot								Tot
							d xz d xy d yz d z 2 d x 2 -y 2								d xz d xy d yz d z 2 d x 2 -y 2
	DOS								DOS
	-20	-18	-16	-14	-12	-10	-8	-6	-20	-18	-16	-14	-12	-10	-8	-6
				E(eV)									E(eV)

3

.

The optimized geometries can be found in https://pubs.acs.org/doi/10.1021/acs.jctc.0c00628.

  Figure 5.14: Projected density of states of [Fe(NH 3 ) 6 ] 2+ in LS (left) and HS (right) computed using PBE+U sc . Solid black line indicates the HOMO level.Figure 5.15: Projected density of states of [Fe(NCH) 6 ] 2+ in LS (left) and HS (right) computed using PBE. Solid black line indicates the HOMO level.Figure 5.16: Projected density of states of [Fe(NCH) 6 ] 2+ in LS (left) and HS (right) computed using PBE+U sc . Solid black line indicates the HOMO level. Biased spin-state energetics of Fe(II) molecular complexes within density functional theory and the linear-response Hubbard-U correction Figure 5.17: Projected density of states of [Fe(PH 3 ) 6 ] 2+ in LS (left) and HS (right) computed using PBE. Solid black line indicates the HOMO level.

				5.4 Supplementary information															73
													Tot	Tot												Tot	Tot
												d xz d xy d yz d z 2 d x 2 -y 2	d xz d xy d yz d z 2 d x 2 -y 2												d xz d xy d yz d z 2 d x 2 -y 2	d xz d xy d yz d z 2 d x 2 -y 2
	DOS	DOS													DOS	DOS
	-20	-20	-18	-18	-16	-16	-14	-14	-12	-12	-10	-10	-8	-8	-6 -20	-20	-18	-18	-16	-16	-14	-14	-12	-12	-10	-10	-8	-8
								E(eV)	E(eV)													E(eV)	E(eV)
													Tot	Tot												Tot	Tot
												d xz d xy d yz d z 2 d x 2 -y 2	d xz d xy d yz d z 2 d x 2 -y 2												d xz d xy d yz d z 2 d x 2 -y 2	d xz d xy d yz d z 2 d x 2 -y 2
	DOS	DOS													DOS	DOS
	-20	-20	-18	-18 -16		-16 -14		-14 -12		-12 -10		-8	-10	-8 -20		-20 -18		-18 -16	-16 -14	-14 -12	-12	-10	-10	-8	-8	-6
								E(eV)	E(eV)													E(eV)	E(eV)
														Tot														Tot
														d xz d xy d yz d z 2 d x 2 -y 2														d xz d xy d yz d z 2 d x 2 -y 2
		DOS														DOS
		-20		-18		-16		-14		-12			-10	-8		-20		-18	-16	-14	-12		-10		-8	-6
									E(eV)														E(eV)

  Figure 5.19: Projected density of states of [Fe(CO) 6 ] 2+ in LS (left) and HS (right) computed using PBE. Solid black line indicates the HOMO level.Figure 5.20: Projected density of states of [Fe(CO) 6 ] 2+ in LS (left) and HS (right) computed using PBE+U sc . Solid black line indicates the HOMO level.Figure 5.21: Projected density of states of [Fe(CNH) 6 ] 2+ in LS (left) and HS (right) computed using PBE. Solid black line indicates the HOMO level.Figure 5.22: Projected density of states of [Fe(CNH) 6 ] 2+ in LS (left) and HS (right) computed using PBE+U sc . Solid black line indicates the HOMO level. Biased spin-state energetics of Fe(II) molecular complexes within density functional theory and the linear-response Hubbard-U correction

			5.4 Supplementary information										75
											Tot									Tot
											d xz d xy d yz d z 2 d x 2 -y 2									d xz d xy d yz d z 2 d x 2 -y 2
		DOS										DOS
		-22		-20		-18		-16		-14	-12		-22		-20	-18		-16	-14	-12	-10
							E(eV)				NH 3						E(eV)
											HS
				Atom		e g		e g		3d t 2g	3d t 2g	3d t 2g	3d e g		3d e g
				Fe	10.41% 10.56% 95.87% 95.85% 95.91% 92.73% 92.91%
				N	6.28%	18.83% 0.42%	0.70%	0.48%	1.49%		0.50%
				N	6.51%	Tot 18.10% 0.46%	0.68%	0.45%	1.43%		0.53%		Tot
				N N	6.61% 6.65%	d xz d xy d yz d z 2 d x 2 -y 2 18.82% 0.51% 18.30% 0.41%	0.68% 0.72%	0.45% 0.48%	1.46% 1.45%		0.54% 0.52%		d xz d xy d yz d z 2 d x 2 -y 2
		DOS		N N	24.21% 0.14% 24.18% 0.14%	0.69% 0.68%	0.20% DOS 0.20%	0.64% 0.65%	0.07% 0.07%		1.87% 1.87%
											LS
				Atom		e g		e g		3d t 2g	3d t 2g	3d t 2g	3d e g		3d e g
		-20	-18	-16 Fe	-14 15.17% 15.19% 92.55% 92.35% 92.31% 88.79% 88.80% -12 -10 -8 -20 -18 -16 -14 -12	-10	-8
				N	E(eV) 5.67%	17.97% 0.82%	0.84%	1.17%	2.31%		E(eV) 0.77%
				N	23.30% 0.26%	1.16%	1.25%	0.41%	0.10%		2.98%
				N	23.30% 0.27%	1.16%	1.26%	0.40%	0.10%		2.98%
				N	5.45%	18.16% 0.79%	0.86%	1.18%	2.34%		0.74%
				N	6.43%	17.18% 0.81%	0.75%	1.27%	2.21%		0.87%
				N	6.67%	Tot 16.97% 0.84% Tot	0.72%	1.27%	2.19%		0.89%	Tot	Tot
										d xz d xy d yz d z 2 d x 2 -y 2	d xz d xy d yz d z 2 d x 2 -y 2								d xz d xy d yz d z 2 d x 2 -y 2	d xz d xy d yz d z 2 d x 2 -y 2
	DOS	DOS									DOS	DOS
		22 -20	-20 -18	-16	-18	-14	-16	-12	-14	-10	-12 -8	-22	-20	-20	-18	-18 -16	-16 -14	-14 -12	-12 -10	-10 -8
					E(eV)	E(eV)									E(eV)	E(eV)

Table 5 .

 5 8: Hirshfeld decomposition of the active molecular orbitals into their atomic compositions for HS and LS NH 3 . H atoms are neglected for clarity.

	5.4 Supplementary information

Table 5 .

 5 9: Hirshfeld decomposition of the active molecular orbitals into their atomic compositions for HS and LS NCH. H atoms are neglected for clarity.

	Biased spin-state energetics of Fe(II) molecular complexes within density functional
				theory and the linear-response Hubbard-U correction
					PH 3
					HS
	Atom	e g	e g	3d t 2g	3d t 2g	3d t 2g	3d e g	3d e g
	Fe	15.07% 15.33% 94.92% 94.82% 94.82% 89.45% 89.70%
	P	17.72% 6.21%	1.02%	0.65%	0.65%	0.81%	2.08%
	P	0.60%	23.55% 0.30%	1.07%	1.08%	2.79%	0.15%
	P	0.60%	23.57% 0.30%	1.07%	1.08%	2.79%	0.15%
	P	18.05% 6.18%	1.02%	0.66%	0.65%	0.80%	2.10%
	P	17.90% 6.35%	1.02%	0.66%	0.66%	0.82%	2.09%
	P	17.56% 6.35%	1.01%	0.64%	0.66%	0.82%	2.07%
					LS
	Atom	e g	e g	3d t 2g	3d t 2g	3d t 2g	3d e g	3d e g
	Fe	35.15% 35.23% 81.40% 81.16% 81.17% 73.25% 73.24%
	P	14.13% 5.28%	2.59%	2.87%	2.35%	2.12%	5.67%
	P	1.06%	17.57% 2.38%	1.96%	3.50%	7.15%	0.41%
	P	0.15%	1.07%	17.56% 2.38%	1.97%	3.49%	7.15%
	P	2.10%	13.10% 4.78%	2.75%	2.97%	2.08%	1.92%
	P	2.09%	12.81% 5.07%	2.74%	3.04%	2.01%	2.01%
	P	2.07%	13.78% 5.62%	2.58%	2.94%	2.30%	2.22%

Table 5 .

 5 10: Hirshfeld decomposition of the active molecular orbitals into their atomic compositions for HS and LS PH 3 . H atoms are neglected for clarity.

	5.4 Supplementary information

Table 5 .

 5

	11: Hirshfeld decomposition of the active molecular orbitals into their atomic
	compositions for HS and LS CO.

Table 5 .

 5 12: Hirshfeld decomposition of the active molecular orbitals into their atomic compositions for HS and LS CNH. H atoms are neglected for clarity.

	Molecule	Spin State	d z 2	d x 2 -y 2	d xy	d xz	d yz
		HS(maj)	0.996 0.996 0.986 0.994 0.994
	[Fe(H 2 O) 6 ] +2	HS(min)	0.131 0.142 0.044 0.048 0.964
		LS(maj/min) 0.222 0.222 0.971 0.978 0.978
		HS(maj)	0.996 0.996 0.983 0.989 0.990
	[Fe(NH 3 ) 6 ] +2	HS(min)	0.178 0.197 0.011 0.011 0.956
		LS(maj/min) 0.294 0.294 0.955 0.957 0.957
		HS(maj)	0.997 0.997 0.971 0.983 0.983
	[Fe(NCH) 6 ] +2	HS(min)	0.181 0.196 0.033 0.033 0.891
		LS(maj/min) 0.381 0.381 0.903 0.903 0.903
		HS(maj)	0.988 0.990 0.961 0.974 0.975
	[Fe(PH 3 ) 6 ] +2	HS(min)	0.280 0.347 0.011 0.011 0.887
		LS(maj/min) 0.547 0.547 0.853 0.857 0.858
		HS(maj)	0.995 0.996 0.959 0.977 0.977
	[Fe(CO) 6 ] +2	HS(min)	0.259 0.260 0.010 0.011 0.860
		LS(maj/min) 0.540 0.540 0.851 0.851 0.851
		HS(maj)	0.993 0.995 0.952 0.973 0.973
	[Fe(CNH) 6 ] +2	HS(min)	0.283 0.286 0.017 0.017 0.829
		LS(maj/min) 0.550 0.550 0.845 0.845 0.845

Table 5 .

 5 [START_REF] Koenig | 5t2-1a1 equilibriums in some iron(II)-bis(1, 10phenanthroline) complexes[END_REF]: Eigenvalues of the occupation matrix {n σ j } obtained with PBE for each spin state (HS, LS) and spin channel (σ=maj(ority),min(ority)). The five occupation matrix eigenvalues are referred to the occupation of the five d levels of the iron center.

	Biased spin-state energetics of Fe(II) molecular complexes within density functional
			theory and the linear-response Hubbard-U correction
	Complex	Spin State	d z 2	d x 2 -y 2	d xy	d xz	d yz	j n sc j -j n 0 j
		HS(maj)	0.998 0.998 0.991 0.996 0.996	0.013
	[Fe(H 2 O) 6 ] +2	HS(min)	0.081 0.089 0.019 0.021 0.979	-0.140
		LS(maj/min) 0.155 0.155 0.985 0.988 0.988	-0.100
		HS(maj)	0.997 0.998 0.990 0.994 0.994	0.019
	[Fe(NH 3 ) 6 ] +2	HS(min)	0.102 0.117 0.006 0.006 0.975	-0.147
		LS(maj/min) 0.212 0.212 0.977 0.978 0.978	-0.100
		HS(maj)	0.997 0.998 0.984 0.991 0.991	0.030
	[Fe(NCH) 6 ] +2	HS(min)	0.111 0.120 0.013 0.013 0.956	-0.121
		LS(maj/min) 0.286 0.286 0.951 0.951 0.951	-0.046
		HS(maj)	0.993 0.994 0.980 0.987 0.988	0.054
	[Fe(PH 3 ) 6 ] +2 PH 3	HS(min)	0.172 0.193 0.005 0.005 0.947	-0.214
		LS(maj/min) 0.450 0.450 0.915 0.917 0.917	-0.013
		HS(maj)	0.996 0.997 0.980 0.987 0.987	0.043
	[Fe(CO) 6 ] +2	HS(min)	0.130 0.173 0.004 0.005 0.944	-0.144
		LS(maj/min) 0.458 0.458 0.909 0.909 0.909	0.010
		HS(maj)	0.995 0.996 0.975 0.986 0.986	0.052
	[Fe(CNH) 6 ] +2	HS(min)	0.170 0.174 0.007 0.007 0.925	-0.149
		LS(maj/min) 0.469 0.469 0.904 0.904 0.904	0.015
	Table

  PBE[U] results are in excellent agreement with the reference CASPT2/CC set and provide a systematically improved description of ∆E H-L compared with PBE[HF]. Specifically, PBE[HF] yields a reasonable prediction for weak-field ligands but it performs poorly for strong field ligands. Our computed values of PBE[HF] energies are similar to those reported in Ref. [144] on the same molecular complexes (i.e. [Fe(H 2 O) 6 ] +2 , [Fe(NCH) 6 ] +2 , [Fe(NH 3 ) 6 ] +2 and [Fe(CO) 6 ] +2

  have been tested. Kronik et al.[START_REF] Prokopiou | Spin-state energetics of fe complexes from an optimally tuned range-separated hybrid functional[END_REF] employed optimally tuned range-separated hybrid functionals to study [Fe(H 2 O) 6 ] +2 , [Fe(NCH) 6 ] +2 , [Fe(NH 3 ) 6 ] +2 and [Fe(bpy) 3 ] +2 and found good agreement with the available CCSD(T) and CASPT2 reference values. In Refs.[START_REF] Wilbraham | Communication: Evaluating nonempirical double hybrid functionals for spin-state energetics in transition-metal complexes[END_REF] and[START_REF] Alipour | Appraising spin-state energetics in transition metal complexes using double-hybrid models: Accountability of sos0-pbescan0-2(a) as a promising paradigm[END_REF] the PBE0-DH-based double-hybrid and SOS0-PBESCAN0-2(a) double-hybrid were employed, respectively, to study [Fe(H 2 O) 6 ] +2 , [Fe(NCH) 6 ] +2 , [Fe(NH 3 ) 6 ] +2 and [Fe(CO) 6 ] +2 . The authors reported a good agreement when comparing with the DMC reference values of Ref.[START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF]. If these values are compared with our CASPT2/CC results, a systematic underestimation of ∆E H-L is observed.In Tab. 6.1 we report the ∆E H-L of each method tested here and the MAE computed with respect to CASPT2/CC. Among the employed approaches, PBE[U] is the best performer with the lowest total MAE of 0.44 eV. It follows TPSSh with a MAE of 0.49 eV and M06-L with MAE of 0.51. PBE[U] with values computed using atomic projectors represents the forth best performer. For the weak-field molcules (i.e. H 2 O, NH 3 and NCH), TPSS0, PBE0 and B3LYP yield the best agreement with the reference with MAEs of 0.01 eV, 0.15 eV and 0.16 eV, respectively. These are followed by M06-L (0.28 eV)and PBE[U] (0.30 eV). We note that PBE[HF] yields a MAE similar to TPPSh for theseThe role of the geometry on the computed ∆E H-L was also investigated. For each of the 11 functionals tested above we compute the ∆E H-L using the geometries optimized with PBE+U, PBE, TPSSh, B3LYP and PBE0. We do this for [Fe(NH 3 ) 6 ] 2+ and [Fe(CO) 6 ] 2+ .The PBE+U geometry is optimized using a procedure that iteratively computes U and then relaxes the geometry with this U until convergence is achieved. Overall, a nonnegligible effect of the geometry on the ∆E H-L is found (see Tabs. 6.6 and 6.7). The ∆E H-L change at most by 0.08 eV for [Fe(NH 3 ) 6 ] 2+ , and 0.25 eV for [Fe(CO) 6 ] 2+ , if the PBE+U geometry is excluded. When the PBE+U geometry is considered, the largest deviation is 0.18 eV and 0.78 eV for the weak-and strong-field ligand molecules, respectively.

	three molecules (0.57 eV). For strong-field ligands (i.e. PH 3 , CO and CNH), the best
	performers are TPSSh (MAE=0.45 eV) and PBE[U] (MAE=0.57 eV) followed by M06-L
	and PBE[U] with atomic projectors, while PBE[HF] is the second worst performer after
	M06-2X (MAE=1.87 eV).

Table 6 .

 6 Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density 2: Adiabatic energy difference (eV), mean signed error (MSE) and mean absolute error (MAE) (in eV) computed with different DFT methods. The reference values are extracted from experimental (exp.) data. a : Ref. [126]; b : Ref. [46]; c : Ref. [44]. The reference value used to calculate MSE and MAE for [Fe(bpy) 3 ] 2+ is 0.589 eV.

		∆E H-L / Periodic		∆E H-L / Gas phase
	Complex
		exp.	PBE[U sc ]	exp.	PBE[U sc ] TPSSh M06-L PBE[HF]
	Fe(phen) 2 (NCS) 2 0.155 a	-0.065	0.093 a	-0.117	0.372 a -0.151 a	-0.887
	Fe(abpt) 2 (NCS) 2	0.137 a	0.086	0.156 a	-0.032	0.433 a 0.121 a	-1.006
	Fe(abpt) 2 (NCSe) 2 0.150 a	0.159	0.184 a	-0.009	0.491 a 0.115 a	-0.950
	Fe[HB(pz) 3 ] 2	0.223 a	0.179	0.363 a	0.251	0.722 a 0.428 a	-0.757
	FeL 2 [BF 4 ] 2	0.198 a	0.196	0.208 a	0.191	0.574 a 0.150 a	-0.671
	[Fe(tacn) 2 ] 2+			0.165 c	0.166	0.443	0.171	-0.727
	[Fe(bpy) 3 ] 2+			0.434-0.744 b	0.466	0.858	0.513	-0.626
	MSE		-0.062		-0.120	0.305	-0.059	-1.055
	MAE		0.065		0.121	0.305	0.079	1.055

  tons. More details are reported in the SI. For [Fe(tacn) 2 ] +2 and [Fe(bpy) 3 ] +2 the structure is optimized using TPSSh with ORCA. The gas phase calculations of ∆E H-L computed using PBE[U ], PBE[HF], TPSSh, and M06-L are reported in Tab. 6.2 together with the experimentally-extracted reference value. For the five crystals, the TPSSh and M06-L results are taken from Ref.[START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF]. For the rest of the calculations (i.e. PBE[HF] and PBE[U] on the seven molecules and TPSSh and M06-L on the last two) we add the D3 correction (similar to Vela et al.[START_REF] Vela | Thermal spin crossover in fe(ii) and fe(iii). accurate spin state energetics at the solid state[END_REF]) with the BJ damping scheme, except for M06-L for which this is not implemented. We note the use of four significant digits in Tab. 6.2, compared to three in Tab. 6.1: the choice in Tab. 6.1 was made for consistency with the approximation reported in the values taken from the literature. We choose however to add a significant figure in Tab. 6.2 because the reported values are closer. PBE[U] and M06-L are the best performers with a MAE of 0.12 eV and 0.08 eV, respectively. They both slightly underestimate the ∆E H-L resulting in negative values of the mean signed error (MSE). TPSSh systematically overestimates the adiabatic energy differences with a MAE and MSE of 0.31 eV. PBE[HF] yields the largest error with a MAE of 1.06 eV. Consistent with the study of the six molecular complexes reported above, PBE[HF] systematically underestimate the ∆E H-L for these intermediate-/strong-field molecules resulting in the wrong prediction of the ground state for the whole set under study. The five molecular crystals were also studied using a full periodic approach using PBE[U] within the D3+BJ approximation for the dispersion forces. Compared to the gas phase calculations, the only difference is the molecular versus periodic model because the energy functional and the U computed on the LS and HS periodic geometries (vide supra) are the same. PBE[U]

Table 6 .

 6 3: CASPT2 values of ∆E H-L in eV for computed using the three basis sets. The last column shows the final energies obtained upon extrapolation (performed separately

			∆E H-L / CASPT2
	Complex	Basis Set cc-pVTZ-DK/DKH cc-pVQZ-DK/DKH cc-pV5Z-DK/DKH CBS
	H 2 O	-2.161	-2.064	-2.021	-1.990
	NH 3	-1.053	-0.942	-0.896	-0.846
	NCH	-0.518	-0.382	-0.323	-0.265
	PH 3	1.642	1.982	2.154	2.308
	CO	1.505	1.650	1.717	1.783
	CNH	2.355	2.524	2.595	2.661

Table 6 .

 6

	103

4: U sc and ∆E H-L (in eV) computed using ortho-atomic and atomic projectors. In parentheses the values computed using PAW pseudopotentials to extract the all-electrons charge density (see text).

Table 6 .

 6 5: U sc and V sc (in eV) used for the DFT+U+V calculations computed using the density-functional perturbation theory as reported in Ref.[START_REF] Timrov | Hubbard parameters from densityfunctional perturbation theory[END_REF] using the hp.x routine.

	In parentheses the values computed using the PAW pseudopotentials to extract the all-
	electron charge density.

Table 6 .

 6

				∆E H-L
	Geometry → PBE+U PBE TPSSh B3LYP PBE0
	Functional ↓
	PBE[U]	2.41	3.19	3.07	2.94	3.10
	PBE	3.06	3.45	3.63	3.39	3.44
	TPSSh	2.01	2.23	2.25	2.23	2.24
	B3LYP	1.16	1.23	1.25	1.27	1.24
	PBE0	1.10	1.34	1.34	1.32	1.35
	TPSS	2.79	3.05	3.07	3.04	3.06
	TPSS0	0.97	1.07	1.10	1.11	1.10
	M06	0.88	0.91	1.03	0.99	0.97
	M06-L	1.68	1.78	1.78	1.83	1.81
	M06L-2X	-1.10	-1.64	-1.54	-1.37	-1.55
	PBE[HF]	0.27	0.66	0.58	0.52	0.60

[START_REF] Reed | A spin transition mechanism for cooperative adsorption in metal-organic frameworks[END_REF]

: ∆E H-L (eV) of Fe(NH 3 ) 6 ] 2+ computed with different functionals using five geometries (PBE+U, PBE, TPSSh, B3LYP, PBE0).

Table 6 .

 6

	6.4 Supplementary information

[START_REF] Cambi | Über die magnetische susceptibilität der komplexen verbindungen[END_REF]

: ∆E H-L (eV) of Fe(CO) 6 ] 2+ computed with different functionals using five geometries (PBE+U, PBE, TPSSh, B3LYP, PBE0).

Table 6 .

 6 9: U sc computed at each geometry.

	Complex	a	b	c	α	β	γ
	Fe(phen) 2 (NCS) 2	LS 12.73 10.14 17.09 90.00 90.00 90.00 HS 13.14 10.02 17.14 90.00 90.00 90.00
	Fe(abpt) 2 (NCS) 2	LS 8.52 10.05 15.96 89.66 91.86 89.72 HS 8.60 10.41 15.81 89.81 91.12 89.76
	Fe(abpt) 2 (NCSe) 2	LS 8.59 10.06 16.06 89.71 91.17 89.76 HS 8.66 10.37 15.90 89.86 90.50 89.80
	Fe[HB(pz) 3 ] 2	LS 9.84 17.23 12.78 90.73 95.75 91.10 HS 9.71 17.57 13.17 89.94 97.31 90.06
	FeL 2 [BF 4 ] 2	LS 8.61 HS 8.55	8.64 18.33 89.98 95.39 90.00 8.61 18.63 90.00 94.93 90.00

Table 6 .

 6

10: Lattice parameters computed using PBE+D3+BJ (see text) for the molecular crystals.

  1; the associated ∆E H-L in presence of gas molecule is reported in the same figure in the lower panel. We use the term ∆E H-L to refer

	Engineering Metal-Organic Frameworks for Spin Crossover-assisted Gas Release: a
								First-Principles Study
	Metal center	gas	S	∆E Charge M-L M-mol Ebind
			0	0	1.62	1.91	1.87	2.159
		CO	1	2.00	1.84	1.96	2.40	0.660
			2	2.69	2.02	2.22	2.27	0.879
			0	0.06	1.80	1.92	2.06	0.960
		CO2	1	0	2.01	1.96	2.42	0.487
	Fe(II)		2 0	0.66 0	2.22 1.75	2.19 1.93	2.24 1.91	0.698 1.343
		N2	1	0.87	1.87	1.96	2.44	0.410
			2	1.70	2.08	2.19	2.34	0.563
			0	0	1.49	1.92	1.67	1.131
		H2	1	0.77	1.70	1.96	2.22	0.247
			2	1.77	1.94	2.18	2.21	0.314
			0	0	1.71	1.91	2.03	1.741
		H2O	1	0.39	1.85	1.96	2.26	1.043
			2	0.90	2.06	2.21	2.14	1.361
			1/2	0	1.82	1.97	1.93	2.270
		CO	3/2 1.21	2.11	2.01	2.38	1.330
			5/2 1.86	2.36	2.19	2.31	1.550
			1/2 0.16	2.11	1.98	1.98	1.599
		CO2	3/2	0	2.27	2.01	2.14	1.344
	Fe(III)		5/2 0.45 1/2 0	2.50 1.98	2.20 1.98	2.06 1.95	1.584 1.507
		N2	3/2 0.40	2.15	2.01	2.30	0.981
			5/2 1.10	2.42	2.19	2.22	1.165
			1/2	0	1.73	1.97	1.71	1.093
		H2	3/2 0.41	2.01	2.01	2.19	0.552
			5/2 1.27	2.27	2.17	2.12	0.700
			1/2	0	1.88	2.01	2.10	2.554
		H2O	3/2 0.21	2.13	2.02	2.12	2.115
			5/2 0.67	2.38	2.21	2.05	2.429
		CO	1/2 3/2 1.30 0	1.69 1.96	1.91 2.16	2.26 2.17	0.755 0.778
	Co(II)	CO2 N2	1/2 3/2 1.15 0 1/2 0 3/2 1.24	1.91 2.17 1.75 2.01	1.90 2.14 1.91 2.14	2.41 2.27 2.34 2.30	0.484 0.581 0.449 0.498
		H2	1/2 3/2 1.34 0	1.75 2.01	1.91 2.13	2.10 2.14	0.291 0.292
		H2O	1/2 3/2 0.86 0	1.71 1.99	1.90 2.15	2.28 2.12	1.004 1.244
		CO	0 1	0.03 0	1.19 1.84	1.87 2.09	2.91 2.13	0.318 0.953
	Ni(II)	CO2 N2	0 1 0 1	0 0.64 0 0.48	1.82 2.10 1.66 1.94	1.87 2.08 1.87 2.09	2.84 2.23 2.97 2.19	0.307 0.609 0.210 0.591
		H2	0 1	0 0.66	1.47 1.78	1.87 2.08	2.81 2.00	0.104 0.398
		H2O	0 1	0 0.09	1.61 1.90	1.87 2.09	2.61 2.14	0.656 1.235

more generally to the relative HS versus LS. Specifically, in the case of Fe(II) and Fe(III) we consider intermediate spin versus low spin (see Tab. 7.1). Using the convention of positive binding energies for attractive interactions, positive values of ∆E bind,H-L indicate a lowering of the bond order upon LS→HS transition. Positive values of ∆E H-L are associated with

Table 7 .

 7

	1: Binding energies (eV), metal-gas bond distance ( Å), electrostatic charge (elec-
	trons), and the adiabatic energy difference ∆E H-L (eV) (computed with respect to the

ground state set at 0 energy), for the four molecular complexes at different spin states, S.

Table 7 .

 7 xz ) 2 (d yz ) 2 (d xy ) 2 (d z 2 ) 0 (d x 2 -y 2 ) 0 1 (d xz ) 2 (d yz ) 2 (d xy ) 1 (d z 2 ) 1 (d x 2 -y 2 ) 0 2 (d xz ) 2 (d yz ) 1 (d xy ) 1 (d z 2 ) 1 (d x 2 -y 2 ) 1 xz ) 2 (d yz ) 2 (d xy ) 1 (d z 2 ) 0 (d x 2 -y 2 ) 0 3/2 (d xz ) 2 (d yz ) 1 (d xy ) 1 (d z 2 ) 1 (d x 2 -y 2 ) 0 5/2 (d xz ) 1 (d yz ) 1 (d xy ) 1 (d z 2 ) 1 (d x 2 -y 2 ) 1 Co(II) 1/2 (d xz ) 2 (d yz ) 2 (d xy ) 2 (d z 2 ) 1 (d x 2 -y 2 ) 0 3/2 (d xz ) 2 (d yz ) 2 (d xy ) 1 (d z 2 ) 1 (d x 2 -y 2 ) 1 xz ) 2 (d yz ) 2 (d xy ) 2 (d z 2 ) 2 (d x 2 -y 2 ) 0 1 (d xz ) 2 (d yz ) 2 (d xy ) 2 (d z 2 ) 1 (d x 2 -y 2 )1 2: Schematic of the electronic configuration of all the spin states and metal complexes considered in this work. These occupations confirmed by looking at the Kohn-Sham orbitals are consistent with the molecular orbital diagram of an elongated octahedral field.

	7.3 Results

Table 7 . 3

 73 Engineering Metal-Organic Frameworks for Spin Crossover-assisted Gas Release: a First-Principles Study sets, pseudopotentials, exchange and correlation functional) but also to the presence of the MOF environment. The binding mechanism and binding geometry of the molecules in MOF is similar to the molecular case, except for CO 2 and H 2 in HS which have a negative binding energy and move towards the center of the pore. Expect for CO, for which the ∆E H-L is larger than 2 eV, our calculations predict that significant changes in binding energy are associated with ∆E H-L values below 1 eV, compatible with transition temperatures near room temperature.

		LS	HS
				∆E H-L
		E bind M-L	E bind	M-L
	CO 0.787 2.031 0.033 2.762 1.967
	N 2	0.499 2.119 0.181 2.890 0.954
	H 2	0.379 1.974 -0.030 2.850 0.298
	CO 2 0.525 2.402 -0.048 3.459 0.178
	H 2 O 0.714 2.202 0.227 2.946 0.469

: Binding energies (eV) and metal-gas bond distance ( Å) computed for LS (S=0)

Acknowledgements I would like to express my very great appreciation to Roberta Poloni who guided me during these years. Thanks to her expertise, her enthusiasm and her genuine respect for me and my ideas I have learned how to do research. I am particularly grateful for the assistance given by Bess Vlaisavljevich and Marc de Boissieu. I wish to thank Claudio Attaccalite and Gabriele D'Avino for the interesting discussions. Thanks to all the people I met at SIMaP laboratory, especially Angel, Mike, Imad, Ali, Charline, Amandine, Zélie, Damian, Solène, Chaima, Manuel. A special thank to Aseem with whom I shared the whole Ph.D. experience. My special thanks are extended to all my friends and my family who represent the most precious thing life ever gave me. Dulcis in fundo, my sincere thanks to Emilie who supported me with her love in the most important moments of this wonderful adventure.

Copyright 2020 American Chemical Society.

Abstract

The spin-state energetics of six Fe(II) molecular complexes are computed using the linearresponse Hubbard U approach within DFT. The adiabatic energy differences, ∆E H-L , between the high spin (S=2) and the low spin (S=0) states are computed and compared

Abstract

We recently showed that the DFT+U approach with a linear-response U yields adiabatic energy differences biased towards high spin [Mariano et In conclusion, we show that the PBE[U] approach consisting of adopting the PBE functional evaluated on the PBE+U density, with a self-consistent approach for the calculation of U , represents a reliable and computationally efficient method for the calculation of spin gaps of both molecular complexes and molecular crystals. We show that for the six Fe(II) molecular complexes ranging from weak-(H PBE+U+V ∆E H-L calculations are performed by using Quantum ESPRESSO [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF]. We use the GBRV ultrasoft pseudopotentials [START_REF] Garrity | Pseudopotentials for high-throughput dft calculations[END_REF][START_REF] Bennett | A systematic determination of hubbard u using the gbrv ultrasoft pseudopotential set[END_REF] with wavefunction and charge density 6.4 Supplementary information 99 cutoffs of 80 Ry and 800 Ry, respectively. The Makov-Payne correction [START_REF] Makov | Periodic boundary conditions in ab initio calculations[END_REF] is added to the total energy calculations of isolated charged systems. The other DFT single point energy calculations are performed using ORCA [START_REF] Neese | The orca program system[END_REF] with a tight integration grid (grid6)

and aug-cc-pVTZ basis sets. For the density analysis, we use PAW pseudopotentials [START_REF] Blöchl | Projector augmented-wave method[END_REF] from the PSlibrary [START_REF] Corso | Pseudopotentials periodic table: From h to pu[END_REF] to compute the all-electron charge density. In this case we use 130 Ry and 900 Ry for the wavefunction and charge density cutoff respectively. U sc and

V sc are recomputed in this case (see Tabs. set for Fe and cc-pVTZ for the other atoms together with a tight integration grid. For all these calculations we add the D3 correction with the BJ damping scheme, except for M06-L for which this is not implemented. In this case we use D3 without the BJ damping.

The PBE[U] calculations of [Fe(tacn) 2 ] +2 and [Fe(bpy) 3 ] +2 and the five gas phase version of the molecular crystals are instead performed using Quantum ESPRESSO, using large unit cells to avoid interaction between images.

CASPT2

The computational procedure adopted here is similar to Ref. [START_REF] Mariano | Biased spin-state energetics of fe(II) molecular complexes within density-functional theory and the linear-response hubbard u correction[END_REF]. The CASPT2 calculations are peformed using the extended-multistate approach (XMS) with a real vertical shift of 0.2 a.u.and a Cholesky decomposition of the two-electron integrals. The three lowest states were state-avaraged and treated with XMS-CASPT2 for HS while only one state was computed for LS. Scalar relativistic effects were included using a second-order Douglas-Improved Spin-State Energy Differences of Fe(II) molecular and crystalline complexes via the Hubbard U -corrected Density

∆ρ analysis

We extract the PBE, PBE+U sc and PBE+U sc +V sc all-electron charge density using Quantum ESPRESSO pp.x utility program. The generated Gaussian cube files contain voxel volume elements with side length of 0.103925 Bohr. The same exact voxel volumes are used in the HF and CASPT2 cube files generated using Multiwfn [START_REF] Lu | Multiwfn: A multifunctional wavefunction analyzer[END_REF] from the files in Molden format written by ORCA and BAGEL [START_REF] Shiozaki | Bagel: Brilliantly advanced general electronic-structure library modular and open-source software project for quantum simulations of materials[END_REF], respectively. We employ the Bader charge analysis code by Henkelman et al. [START_REF] Henkelman | A fast and robust algorithm for bader decomposition of charge density[END_REF][START_REF] Sanville | Improved grid-based algorithm for bader charge allocation[END_REF][START_REF] Tang | A grid-based bader analysis algorithm without lattice bias[END_REF][START_REF] Yu | Accurate and efficient algorithm for bader charge integration[END_REF] to extract the Fe electron charge density ρ F e (r). by reporting values of adiabatic energy differences computed between different spin states and by studying them in conjunction with changes in binding energies.

We employ density-functional theory calculations to predict the change in binding energy upon spin-state transitions among different spin configurations and the adiabatic energy difference associated with the spin crossover, in presence of adsorbed gas, in order to assess the feasibility of the process. We first perform a screening study by computing a small molecular fragment consisting of a square planar M(CN)

Computation Methods

The DFT calculations of the molecular fragments are performed by employing the TPSSh functional along with aug-cc-pVTZ basis sets as implemented in ORCA [START_REF] Neese | The orca program system[END_REF]. The choice of this functional is motivated by a large number of computational studies showing that this yields accurate values of adiabatic energy differences, i.e. the total energy difference computed between different spin states each at the corresponding geometry. Good agreements are found against higher level theory [264] and experimentally-extracted values, whenever available [START_REF] Cirera | Theoretical prediction of spin-crossover temperatures in ligand-driven light-induced spin change systems[END_REF][START_REF] Cirera | Theoretical modeling of spin crossover in metal-organic frameworks: [fe(pz) 2 pt(cn) 4 ] as a case study[END_REF][START_REF] Cirera | Benchmarking density functional methods for calculation of state energies of first row spin-crossover molecules[END_REF]. The semiempirical Grimme's D3 correction [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-d) for the 94 elements h-pu[END_REF] is employed in conjunction with the Becke-Johnson (BJ) damping scheme [START_REF] Smith | Revised damping parameters for the d3 dispersion correction to density functional theory[END_REF] to account for van der Waals forces which are critically important for the correct description of the adsorption mechanism [START_REF] Poloni | Co 2 capture by metal-organic frameworks with van der waals density functionals[END_REF]4]. For each calculations we place two molecules on each side of the metal atom and perform a full geometrical optimization. The molecules are always bound in end-on configuration except for H 2 molecule for which the side-on geometry is energetically favored, in agreement with previous studies [3,[START_REF] Sun | Effect of spin state on the dihydrogen binding strength to transition metal centers in metal-organic frameworks[END_REF][START_REF] Kim | Enhanced dihydrogen adsorption in symmetry-lowered metal-porphyrin-containing frameworks[END_REF]. For the charge analysis, we employ the Bader charge method [START_REF] Henkelman | A fast and robust algorithm for bader decomposition of charge density[END_REF][START_REF] Sanville | Improved grid-based algorithm for bader charge allocation[END_REF][START_REF] Tang | A grid-based bader analysis algorithm without lattice bias[END_REF][START_REF] Yu | Accurate and efficient algorithm for bader charge integration[END_REF] to extract the charge associated with the metal cation. The calculations of spin energetics and binding energies the Hoffman-like First-Principles Study Chapter 8

Conclusions

In the present thesis SCO MOFs and their gas adsorption properties have been investigated by means of computational electronic structures methods in order to establish a connection between OMS spin state and binding mechanism with the guest adsorbate.

MOFs are nanoporous materials whose chemical and physical properties can be tuned to achieve large gas uptake and high affinity with different gas molecules. However, for a capture/release gas strategy to be technologically attractive the energy cost associated with the gas release (energy penalty) has to be as small as possible. The release process is often performed by heating the material, and the use of such temperature swing may lead to prohibitively high energy penalties. The release process can be greatly enhanced if the microscopic binding mechanism between the adsorbate and the adsorbent can be "switched" through some external stimulus thus possibly leading to large working capacity with minimum energy penalty. Motivated by the studies of J. Long et al. [START_REF] Reed | A spin transition mechanism for cooperative adsorption in metal-organic frameworks[END_REF] and by some previous computational studies in the group [1], where a cooperative spin-transition at the metal centers, upon adsorption, lead to exceptionally large working capacities, we have investigated a novel capture/release strategy where a temperature-driven SCO transition can occur on the OMS upon gas adsorption so that a modification of the binding mechanism can happen.

The computational description of SCO phenomena in molecular and, even more, in crystalline periodic systems (like MOFs) is extremely challenging and part of this thesis has been devoted to the development of a DFT strategy that allows us to calculate with enough precision the adiabatic energy differences (∆E H-L ) of Fe(II) SCO complexes. A reliable DFT approach is in fact desirable in order to describe large systems in their crystalline configuration where highly-accurate quantum chemistry methods cannot be applied. If on the one hand ∆E H-L is not the only quantity needed for a complete description of the energies of HS and LS, i.e. ∆E U = E HS U -E LS U , is at the origin of the observed bias. For self-consistent values of U, ∆E U is up to -6 eV for the strongest ligand-field complex, i.e [Fe(CNH) 6 ] +2 , resulting in an error in the adiabatic energy difference of around 4.6 eV, much larger than in the case of pure PBE. The systematic shift introduced by the Hubbard energy has been rationalized by means of MO analysis: the different coordination environment of the Fe(II) center in HS and LS results in higher fractional occupancies of the localized basis set in LS compared to HS thus leading to a systematic penalization of closed-shell electronic configurations. Despite this, DFT+U within the linear-response approach is found to be able to correct the spurious DFT energy curvature resulting from the SIE effects on the electronic charge density.

The analysis of the DFT+U performances in computing ∆E H-L has served us to develop a strategy in which the semilocal PBE functional is used to calculate the total energy using as input density the one obtained with the self-consistent linear-response DFT+U method. In this way the bias introduced by ∆E U is removed while keeping the improved description of the DFT+U electronic charge density. This method is referred to with the notation PBE[U] and it has been used in chapter 6 to compute the adiabatic energy difference for the same set of six Fe(II) molecular complexes used for the analysis of DFT+U functional plus a set of seven crystalline and molecular Fe(II) compounds for which the reference ∆E H-L has been extracted from experimental measurements of enthalpy difference. We obtain excellent agreement with the reference CASPT2/CC values and with the experimentally extracted adiabatic energy differences. However, the average error obtained in crystalline systems is of the order of ∼ 0.07 eV which would lead to the evaluation of the transition temperature still with a substantial imprecision (around 100 K accordingly to our qualitative analysis). This has to be put in a more general context considering that the chemical accuracy, often considered as a good level of approximation, is around 0.04 eV. So the role of ab initio approaches in SCO simulations is not that of giving precise quantitative results but rather to set qualitative guidelines for experimental investigations. For these calculations, we found that PBE[U] overall outperforms most of the commonly used DFT functionals and, moreover, it allows us to efficiently compute ∆E H-L in periodic boundary condition calculation. In a recent study [START_REF] Song | Benchmarks and reliable dft results for spin gaps of small ligand fe(ii) complexes[END_REF], Burke et al. proposed to use PBE energy functional in a non-self-consistent way where the input density is taken from HF method (PBE[HF]). This approach, like PBE[U], aims at correct the SIE effects on the electronic density thus possibly yielding improved spin energetics with the assumption that these are strongly affected by density-driven error. However,