Tzu - Yao

transfert d'électrons dans des systèmes réalistes présentée par

First of all, I would like to attribute this thesis to my supervisors, Mathieu Salanne and Guillaume Jeanamiret. From the beginning of the recruitment process for this thesis project until the very last moment before the thesis defence, they have constantly provided me with advice to guide me along the winding path to completion. Without their firm and enlightening support, I would never have been able to accomplish such a large amount of work and progress.

Next, I would like to thank the people who have been working with me. Carlo Pierleoni and Michele Ruggeri from the partner lab, Maison de la Simulation, as well as Kyle Reeves from my lab, PHENIX, provided me with the quantum mechanical calculation data that allow me to study the solvation effect at the electrochemical interface between Lithium-ion aqueous solution and carbon electrode surface. Roxanne Berthin and Alessandra Serva from my lab provided me with classical molecular dynamics simulation data that allow me to benchmark my results of anthraquinone solvated in acetonitrile. The scientific contributions of our co-authored articles could not be as important without your collaboration.

PHENIX lab is the first working place that I fully committed to in France. People here are so nice and supportive, no matter if you are permanent or non-permanent personnel, having a glimpse of time and work overlap with me, or, like Roxanne Berthin, Camille Bacon, and Jeanne Decayeaux, being a PhD student in the same year as me. It is a great pleasure to learn technically, to grow professionally, and to do science in such a positive working environment created by all of you.

Finally, I want to thank my family and my friends at Paris as well as those in Taiwan. With your emotional supports, I was able to overcome difficulties of studying abroad. Without you, this thesis may never have been completed or even be started.

List of Figures

without and (dashed) with the PMV correction and (dotted) using the WDA bridge functional as a function of the charge of the fictitious solutes. The value of the PMV correction is also displayed as a function of the charge of the fictitious solutes. e is the elementary charge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.3 Average vertical energy gap computed using the HNC (black solid line) and WDA functional (red solid line) as a function of the atomic charge q. The dotted and dashed lines are the tangent to the HNC curve at q = -1 and q = 1, respectively. . . . . . . . . . . . . . . . . 76 6.4 2 Molecular structures of anthraquinone (AQ) and ionic liquid (IL) grafted with anthraquinone (AQ) [START_REF] Reeves | A first-principles investigation of the structural and electrochemical properties of biredox ionic species in acetonitrile[END_REF] . . . . . . . . . . . . . . . . . . 84 7.3 Two dimensional maps of acetonitrile (MeCN) density and polarization maps in the aromatic plane of AQ computed with ab-initio molecular dynamics (AIMD), classical molecular dynamics (classical MD) and MDFT. Grey, blue, and red circles represent carbon, hydrogen, and oxygen atoms, respectively. Polarization is represented by arrows oriented in the direction of the projection of P and whose length is proportional to ||P || . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 7.4 Solvent density isosurfaces of (A) AQ and (B) 1,4-benzoquinone (BQ) at 10 -2 × n 0 calculated by MDFT. Oxygen, hydrogen, and carbon atoms are displayed in red, cyan and white respectively. . . . . . . . . 88 7.5 Free energy profiles of electron transfer half reactions of (A) AQ and anthraquinone anion (AQ -) calculated by classical MD and MDFT, (B) 1,4-benzoquinone (BQ) and 1,4-benzoquinone anion (BQ -), (C) duroquinone (DQ) and duroquinone anion (DQ -) calculated by MDFT.

As described in the text, the vertical energy gaps were shifted by the same constant value, which was chosen in order to align the minima for the AQ/anthraquinone anion (AQ -) redox couple. Solid lines are fit of MDFT results using Marcus theory. . . . . . . . . . . . . . . . . 90

x LIST OF FIGURES 7.6 Primary molecular structure of (A) 2-R-BQ, (B) 2,3-R-BQ, (C) 2,6-R-BQ, (D) 2,5-R-BQ, (E) 2,3,5-R-BQ, and (F) 2,3,5,6-R-BQ derivatives with R's represent the functional groups. R's in the figure represent the sites to be substituted by functional groups. . . . . . . . . . . . . 91 7.7 Reorganization free energies of functionalized BQ derivatives as a function of substitution positions. The notation of primary molecular structures follows those in Figure 7.6. The black dashed line is the reorganization free energy of BQ. . . . . . . . . . . . . . . . . . . . . 92 7.8 Top: Reorganization free energies of AQ, BQ and BQ derivatives with respect to the inverse cavity radius, 1/r c . The solid line is a linear regression of the data. Bottom: Reorganization free energies of mono-functionalized BQ derivatives with respect to the volume of the functional group estimated as the difference of cavity volume between BQ and the mono-functionalized BQ derivatives. . . . . . . . . . . . . 93 7.9 Ideal (black), external (blue), and excess (red) contributions to the total (green) reorganization free energy of BQ and its monofunctionalized derivatives computed using Equations 3.

Abstract

The power output of electrochemical energy storage devices are closely related to the underlying redox reaction kinetics. Marcus theory provides a theoretical framework to understand the reaction mechanism and predict the rate constant. However, some key properties introduced by Marcus theory cannot be easily studied by experiments and one may want to resort to molecular simulations. Molecular density functional theory (MDFT) is a theory developed based on classical density functional theory (cDFT). MDFT possesses the advantage of numerical efficiency for the study of phenomena that are dominated by solvent effects. This makes MDFT an appropriate tool to study the kinetics of redox reactions involving electron transfer.

Developments and verification of MDFT is done in this thesis for the purpose of studying realistic electrochemical electron transfer reactions. Precise predictions of the electron transfer properties at electrochemical interfaces using MDFT require the application of pressure correction methods. The correctness of the two pressure correction strategies to MDFT proposed in the literature has been examined in the context of electron transfer. It has been shown that only those correction methods that respect the variational principle can be used for studying electron transfer. On the other hand, we shown that MDFT reproduces very well the results of solvation structure and electron transfer properties of quinones solvated in acetonitrile (MeCN) that are calculated with molecular dynamics (MD) simulations using much more sophisticated force fields. We then take advantage of the numerical efficiency of MDFT to calculate the reorganization free energies of numerous 1,4-benzoquinone (BQ) derivatives. We found that the electron transfer of all the BQ derivatives solvated in MeCN can be well described by Marcus theory. Finally, a constant potential electrode model has been implemented in MDFT. The adsorption free energy profiles of lithium ion solvated in water on a graphite electrode is studied based on the development. It has be found that, upon including the water solvent effects, the energy minimum obtained in gas phase calculation disappears. The results of this thesis pave the wave for studying realistic electrochemical electron transfer reactions of a wide range of systems using MDFT. In recent years, more and more attention has been paid to the urgent need of reducing carbon emission in order to mitigate the climate change. Energy is one of the important sectors of activities that should be transformed to attain net zero carbon emission. The recent development of infrastructures to exploit renewable energy, for instance wind power and solar energy, is rapidly growing. However, renewable energy suffers from production intermittence making 100% reliance on them impossible without energy storage.

Electrochemical energy storage devices

Electrochemical energy storage devices are among the most widely used techniques for storing electric energy from renewable sources. There is a range of energy storage solutions based on electrochemical energy storage devices for different scales 1.1. ELECTROCHEMICAL ENERGY STORAGE DEVICES Figure 1.1: Schematic representation of a rechargeable lithium ion battery during discharge process. Reprinted from Ref. [START_REF] Tarascon | Issues and challenges facing rechargeable lithium batteries[END_REF]. of stored energy. These solutions may make use of electrochemical reactions in different ways.

Rechargeable batteries

Rechargeable batteries are the most well known and most developed electrochemical energy storage devices. A rechargeable battery is composed of two electrodes, an anode and a cathode, internally connected by an electrolyte and externally connected to the load. During the discharge process, the chemical energy stored in the battery is converted to electric energy through the oxidation and reduction reactions that take place separately at the surface of anode and cathode, respectively. Both redox reactions involve the change of electrode materials. The cell potential is the measure of the potential difference between the two redox half reactions in an electrochemical cell. For a battery, it is the difference between the potentials of the anode and cathode materials. The amount of chemical energy converted is the integral of the cell potential multiplied by the current during the discharge process. The changes of ionic concentration and ionic charge in the electrolyte are balanced through the movement of ions within it. During the reverse process, namely the charge process, electric energy is used to reverse the above mentioned process and converted to chemical energy stored in the materials.

Lithium ion battery is a well-known example of rechargeable battery. Thanks to its high energy density which is the amount of energy stored in an energy storage device divided by its own volume, lithium ion battery is nowadays the preferred choice for the power source of mobile devices [START_REF] Tarascon | Issues and challenges facing rechargeable lithium batteries[END_REF][START_REF] Dunn | Electrical energy storage for the grid: a battery of choices[END_REF]. A schematic representation of a rechargeable lithium ion battery is shown in Figure 1.1. The two electrodes of lithium ion batteries are composed of materials that allow lithium ion to be inserted in their structure. However, despite the cost reduction at industrial production scale, the production cost of lithium ion batteries is still high. Furthermore, the power output of batteries throughout the discharge is not homogeneous enough to adequately regulate the renewable energy production [START_REF] Dunn | Electrical energy storage for the grid: a battery of choices[END_REF][START_REF] Yang | Electrochemical energy storage for green grid[END_REF]. As a result, it is challenging to meet the need of larger capacity, stationary energy storage solutions with rechargeable lithium ion batteries.

Redox flow batteries

Redox flow battery is another type of electrochemical energy storage device that stores chemical energy in the electrolytes [START_REF] Gentil | Aqueous organic and redox-mediated redox flow batteries: a review[END_REF][START_REF] Leung | Recent developments in organic redox flow batteries: a critical review[END_REF][START_REF] Yuan | Membranes in non-aqueous redox flow battery: a review[END_REF][START_REF] Potash | On the benefits of a symmetric redox flow battery[END_REF]. Two electrolytes either containing oxidant or reductant, respectively called catholyte and anolyte, are stored in two separate tanks. During the discharge process, electrolytes are pumped to the reacting cell which is composed of an ion-selective membrane and two electrodes. The two redox half reactions take place separately at the two electrode/electrolyte interfaces, reduction on the catholyte side and oxidation on the anolyte side, resulting in the change of the oxidation state of the redox active species. Unlike common batteries, the electrodes only serve as current collectors. The ion-selective membrane allows supporting electrolyte ions flowing through it to balance the ionic charge but prohibit the crossover of redox active molecules. During the charge process, the external input electricity is used to reverse the redox reactions and is converted back to chemical energy in the form of redox potential difference between the two redox active species.

Vanadium redox flow battery is a typical and commercialized flow battery. Its catholyte and anolyte are composed of pervanadyl/vanadyl ion pair (VO 2 + and VO 2+ ) and vanadium ions (V 3+ and V 2+ ), respectively. A schematic representation of a redox flow battery based on vanadium is shown in Figure 1.2.

Organic molecules for energy storage

Organic redox flow battery

Commercially available redox flow batteries are usually built with aqueous electrolytes containing metallic elements, such as vanadium. However, the production cost of these systems is also too high for widespread applications. To reduce the cost, it is proposed to instead use organic redox active molecules that can be synthesized from inexpensive commodity chemicals [START_REF] Park | Investigation of charge transfer kinetics at carbon/hydroquinone interfaces for redox-active-electrolyte supercapacitors[END_REF]. Another great advantage of organic molecules is their highly versatile molecular structure. Not only the redox potential [START_REF] Yu | Tuning the performance of aqueous organic redox flow batteries via first-principles calculations[END_REF][START_REF] Huynh | Quinone 1 eand 2 e-/2 h+ reduction potentials: identification and analysis of deviations from systematic scaling relationships[END_REF][START_REF] Er | Computational design of molecules for an all-quinone redox flow battery[END_REF][START_REF] Bachman | Investigation of the redox chemistry of anthraquinone derivatives using density functional theory[END_REF][START_REF] Schwan | Substituent pattern effects on the redox potentials of quinone-based active materials for aqueous redox flow batteries[END_REF] but also the solubility [START_REF] Er | Computational design of molecules for an all-quinone redox flow battery[END_REF][START_REF] Pineda Flores | Bio-inspired electroactive organic molecules for aqueous redox flow batteries. 1. thiophenoquinones[END_REF], possible number of oxidation state changes, chemical stability [START_REF] Wedege | Organic redox species in aqueous flow batteries: redox potentials, chemical stability and solubility[END_REF][START_REF] Tabor | Mapping the frontiers of quinone stability in aqueous media: implications for organic aqueous redox flow batteries[END_REF][17] and reaction kinetics can be finely tailored by grafting functional groups to the main molecular skeleton [18][19][20][21].

Compared to common rechargeable batteries that store energy in solid electrode materials, it is much easier to scale up redox flow batteries that store energy in liquid phase. The amount of stored energy and power output increase as the size of electrolyte tanks and surface area of electrodes increase, respectively. The cycling stability of common rechargeable batteries is another important issue that should be addressed for long-term usage. Due to the nature of energy storage mecha-1.2. ORGANIC MOLECULES FOR ENERGY STORAGE Figure 1.2: Schematic representation of a vanadium redox flow battery during a discharge process. Reprinted with permission from [START_REF] Yang | Electrochemical energy storage for green grid[END_REF]. Copyright 2011 American Chemical Society. nism, rechargeable batteries suffer from numerous degradation mechanisms in solid phase [22]. For example, dendrite formation at lithium metal anode may cause short circuit and catastrophic outcome [23]. On the contrary, the capacity fade mechanisms of aqueous organic redox flow batteries are mostly associated to the decomposition of redox active molecules in liquid electrolyte [17]. Due to the different device working mechanisms, not only low capacity fade rate is achievable for aqueous organic redox flow batteries, but also the system maintenance is relatively easy. The scalability, low production cost, and low maintenance cost thus make them suitable for stationary applications [START_REF] Dunn | Electrical energy storage for the grid: a battery of choices[END_REF].

The energy density of electrochemical energy storage devices is closely related to the cell potential. Despite that we can design organic redox active molecules pair exhibiting high potential difference, the maximum cell potential is yet limited by the electrochemical window of the solvent. For aqueous redox flow batteries, it is around 1.23 V under standard condition due to water electrolysis [24]. Although aqueous electrolytes benefit from higher level of safety and lower cost, non-aqueous redox flow battery using organic solvent can work at a larger cell potential. However, non-aqueous systems often suffer from slow transport properties in the electrolytes and bad chemical stability [25]. Other than cell potential, the solubility and possible number of oxidation state changes of redox species also directly influence the energy density of redox flow battery. These properties as well as the reaction kinetics and chemical stability may be different for every combinations of redox species and solvent [26]. Thus, it is important to understand the relationships between the molecular structure and its chemical properties in order to sort out promising candidates.

Quinones as organic redox active materials

Quinones is a family of organic molecules that can undergo reversible redox reactions [27][28][29][30][31]. It has recently been demonstrated in the literature that they can serve as the redox active materials for redox flow battery applications [32][33][34][START_REF] Son | Quinone and its derivatives for energy harvesting and storage materials[END_REF][START_REF] Ding | Molecular engineering of organic electroactive materials for redox flow batteries[END_REF][START_REF] Leung | A new aqueous all-organic flow battery with high cell voltage in acidic electrolytes[END_REF][START_REF] Han | Organic quinones towards advanced electrochemical energy storage: recent advances and challenges[END_REF]. The main structural feature of quinones are the aromatic rings with a pair of hydrogen atoms substituted by oxygen atoms. Their chemical properties can be modified by the substitution of the rest of the hydrogen atoms with functional groups [9-11, 14-16, 18]. Therefore, it is possible to design suitable quinones molecules for the specific solvent used to build the redox flow batteries.

Let's take simple quinones without functional groups as an example to illustrate these features. 1,2-benzoquinone and 1,4-benzoquinone (BQ) are the simplest quinones and are derived from a benzene ring with two hydrogen atoms substituted by oxygen forming two carbonyl groups. Their molecular structures are displayed in Figure 1.3. The π bonds of these carbonyl groups also participate in the conjugation. The main structural feature of quinones may also comprise 2 or more aromatic rings. Naphtoquinone and anthraquinone (AQ) are two typical quinones comprising respectively 2 and 3 aromatic rings. Their molecular structures are displayed in Figure 1.4. Quinones are able to uptake the 2 electrons to form lone pairs on oxygen. In aprotic solvent environment, e.g. acetonitrile (MeCN), reduced quinones exist as radical anions. In contrast, in protic solvent environment, e.g. water, reduced quinones may or may not, depending on the pH value, combine with one or two protons to form hydroquinone anions or hydroquinones [START_REF] Yu | Tuning the performance of aqueous organic redox flow batteries via first-principles calculations[END_REF][START_REF] Huynh | Quinone 1 eand 2 e-/2 h+ reduction potentials: identification and analysis of deviations from systematic scaling relationships[END_REF]25,28,30]. The Thorough review on the organic redox flow batteries can be found in the references [START_REF] Leung | Recent developments in organic redox flow batteries: a critical review[END_REF]24].

Furthermore, quinones have recently been shown capable of improving the capacity of electric double layer capacitors [START_REF] Mourad | Biredox ionic liquids: electrochemical investigation and impact of ion size on electron transfer[END_REF][START_REF] Mourad | Biredox ionic liquids with solid-like redox density in the liquid state for high-energy supercapacitors[END_REF]. The energy storage mechanism of normal electric double layer capacitors does not involve electrochemical reaction [START_REF] Merlet | Influence of solvation on the structural and capacitive properties of electrical double layer capacitors[END_REF][START_REF] Jeanmairet | Microscopic simulations of electrochemical double-layer capacitors[END_REF]. During the charge process, ions in the liquid electrolyte drift under the applied potential and adsorb on the electrode surface to form electric double layers. The electric energy is stored in the capacitors that are composed of charges of electric double layers on one side and counter charges in the electrodes on the other. Electric double layer capacitors have high power density, the amount of power an energy storage device can deliver divided by its own volume. Nonetheless, they have poor energy density and short discharge duration compared to electrochemical energy storage devices despite the improvements in capacity achieved by employing porous electrodes to increase the active surface area and ionic liquids to allow for a higher cell potential. Further improvements of energy density and discharge duration could expand the range of application of electric double layer capacitors. It is thus proposed to gain capacity from redox reactions by grafting quinones onto ionic liquids [START_REF] Mourad | Biredox ionic liquids: electrochemical investigation and impact of ion size on electron transfer[END_REF][START_REF] Mourad | Biredox ionic liquids with solid-like redox density in the liquid state for high-energy supercapacitors[END_REF][START_REF] Reeves | A first-principles investigation of the structural and electrochemical properties of biredox ionic species in acetonitrile[END_REF][START_REF] Berthin | Solvation of anthraquinone and TEMPO redox-active species in acetonitrile using a polarizable force field[END_REF]. The so-formed devices are called biredox ionic liquid supercapacitors.

Reaction kinetics

The power output of electrochemical energy storage devices are closely related to the kinetics of the underlying redox reactions. During charge and discharge process of a redox flow battery or a biredox ionic liquid supercapacitor, redox species in solution grabs or donates electrons from or to current collector, resulting in the change of their oxidation state. These reactions can be categorized as redox half reactions involving electron transfer at the electrode surface. A half-cell reduction CHAPTER 1. MOTIVATION reaction can be expressed as

Ox + ne -= Red (1.1)
where Ox and Red denote the single redox active species in its oxidant and reductant form, respectively, n denotes the number of electrons transferred, and e -stands for the electron.

1.3 Theoretical studies of electron transfer

1.3.

Marcus theory of electron transfer

Marcus theory is a widely accepted model of electron transfer in solution and for electrochemical reactions [START_REF] Marcus | On the theory of oxidation-reduction reactions involving electron transfer. I[END_REF][START_REF] Marcus | Electrostatic free energy and other properties of states having nonequilibrium polarization. I[END_REF][START_REF] Marcus | On the theory of oxidation-reduction reactions involving electron transfer. II. applications to data on the rates of isotopic exchange reactions[END_REF][START_REF] Marcus | On the theory of oxidation-reduction reactions involving electron transfer. III. applications to data on the rates of organic redox reactions[END_REF][START_REF] Marcus | Exchange reactions and electron transfer reactions including isotopic exchange. Theory of oxidation-reduction reactions involving electron transfer. Part 4.-A statistical-mechanical basis for treating contributions from solvent, ligands, and inert salt[END_REF][START_REF] Marcus | Free energy of nonequilibrium polarization systems. II. homogeneous and electrode systems[END_REF][START_REF] Marcus | ON THE THEORY OF OXIDATION-REDUCTION REACTIONS INVOLV-ING ELECTRON TRANSFER. V. COMPARISON AND PROPERTIES OF ELECTROCHEM-ICAL AND CHEMICAL RATE CONSTANTS 1[END_REF][START_REF] Marcus | On the theory of electron-transfer reactions. VI. unified treatment for homogeneous and electrode reactions[END_REF][53]. It provides a conceptual framework for the interpretation of measured rate constants of redox reactions. The theory can also be applied in an inverse manner to predict the rate constants. However, it is not easy to measure the parameters that appear in the equations of Marcus theory with experiments. The modern formulation of Marcus theory provides a guideline for applying molecular simulations to study electron transfer [54][55][56]. Introductions to Marcus theory and the applications of molecular simulations to electron transfer are given in Chapter 4.

Molecular simulations on electron transfer

Molecular simulations have evolved to tools that complement theories and experimental works for the development of new technologies. They can be used to help us understand complex material systems and predict their properties. Especially, they provide a point of view at molecular scale that may explain properties observed yet hard to investigate by experiments. For instance, Mourad et al. shown that biredox ionic liquid supercapacitors can have doubled capacity compared to the electric double layer capacitors built without grafting quinones onto the ionic liquids [START_REF] Mourad | Biredox ionic liquids with solid-like redox density in the liquid state for high-energy supercapacitors[END_REF]. The underlying mechanism of the dramatic improvement is now under study by our group using molecular dynamics (MD) simulations [START_REF] Berthin | Solvation of anthraquinone and TEMPO redox-active species in acetonitrile using a polarizable force field[END_REF].

Theories are commonly established with simplifications or approximations. To what extent a theory is valid can be examined by molecular simulations. The results calculated by molecular simulations can also provide insights for the improvement of the theory. For example, MD simulations have been used in our group to study the confinement effects on electron transfer in nano-porous carbon electrodes. It has been found that redox properties of iron cations are strongly influenced by the breakdown of solvation shell in the nano-pores [57].

Since the development of the modern formulation of Marcus theory, MD has long been the simulation method of choice for studying electron transfer. It consists in distributing particles inside a simulation box of finite size, specifying their initial momentum as well as the interaction between them and letting the system to evolve as time passes according to the Newton's equations of motion. An introduction to the working principles of MD and the application of MD to study electron transfer 1.3. THEORETICAL STUDIES OF ELECTRON TRANSFER is given in Section 4.4. It should be noted that properties of the system at a specific time of a simulation is not directly the properties that we can measure with experiments. To calculate measurable properties, one should make use of the equations derived from statistical physics. An introduction to concepts of statistical physics, which may help in understanding other sections of the thesis, is given in Chapter 2.

There are several possible scenarios in which the basic model of Marcus theory may break down due to effects at molecular scale [58]. MD is also applied to help understand molecular effects that are ignored by Marcus' macroscopic description and to develop more sophisticated models based on Marcus theory. However, it is a heavy computational task to study electron transfer with MD due to the required simulation time to reduce statistical error. The computational cost of MD simulations makes the studies of large amount of redox active molecules very difficult.

Molecular density functional theory (MDFT)

Recently, our group has followed Marcus theory to develop a numerically more efficient alternative approach based on molecular density functional theory (MDFT) to study electron transfer [59,60]. MDFT is a flavor of classical density functional theory (cDFT) developed to study solvation problems. The solvent is described by its density field ρ(r, Ω) that depends on the position r and orientation Ω of the solvent molecules and a free energy functional of this density F[ρ(r, Ω)]. To see the solvation effects, the homogeneous bulk solvent with no externally applied potential is taken as the thermodynamic reference state. The solute is then introduced as an external perturbation that causes solvent to deviate from the bulk density ρ b . According to the variational principle of cDFT for grand-canonical ensemble, for a given solute, there exists only one equilibrium solvent density field ρ eq (r, Ω) and the free energy functional of this equilibrium density F[ρ eq (r, Ω)] is exactly the minimum of the free energy functional F[ρ eq (r, Ω)] = F min . Furthermore, the free energy functional at minimum F min is the solvation free energy of the solute that perturbs the solvent. Therefore, one can find the solvation free energy of a solute as well as the corresponding equilibrium solvent density field by minimizing the free energy functional with respect to the solvent density [61,62]. A complete derivation of MDFT is given in Chapter 3.

The application of MDFT to electron transfer

Compared to MD in which properties are calculated by accumulating data and taking their average, the minimization procedure mentioned above is much more efficient. The computational cost of studying one system is thus largely reduced. It is then possible to study electron transfer properties of a series of molecules with the approach based on MDFT. In the thesis, the electron transfer properties of a series of quinones solvated in MeCN solvent are studied with the approach based on MDFT. Before doing so, the validity of the application of the approach on quinones is examined by comparing the results of the solvation structure and the electron transfer half reaction of AQ solvated in MeCN calculated with MDFT to those calculated with MD simulations. (The MD simulations are performed by one of our CHAPTER 1. MOTIVATION group member.) The derivation of the MDFT based approach is given in Section 4.7. The validation of the approach and the results of electron transfer properties of quinones are discussed in Chapter 7.

Pressure correction to MDFT

Following the derivation of MDFT in Chapter 3, we can see in Section 3.1.5 that a part of the free energy functional should be approximated in order to make it computer tractable [61]. Later in Section 3.3, we also see that the approximation leads to the overestimation of the solvation free energy [63]. Several attempts were made to fix this issue, either based on simple a posteriori correction [64][65][66] or by introducing additional terms to the approximated functional [63,67,68]. In the course of studying electron transfer, we found that the application of simple a posteriori correction leads to artifacts in the calculated results. In Chapter 6, the correctness of these two types of corrections in the study of electron transfer is assessed by studying two model half reactions in water, Cl → Cl + and Cl → Cl -[69].

Constant potential electrode in MDFT

In order to simulate electrochemical reactions, it is necessary to incorporate electrode model in MDFT. Constant potential electrode is a realistic model that allows the atomic charges on electrode atoms to fluctuate according to the electric field generated by the electrolyte [70][71][START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF][START_REF] Coretti | MetalWalls: simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes[END_REF]. In a recent paper, our group has demonstrated the theoretical background for incorporating constant potential electrode model in MDFT [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF]. Introductions to the constant potential electrode model and the theoretical background for the incorporation in MDFT are given in Chapter 5.

The MDFT code equipped with constant potential electrode model is then applied to calculate the adsorption free energy profile of a lithium ion on graphite electrode. The adsorption free energy profile is the adsorption free energy as a function of distance between the lithium ion and carbon electrode. The adsorption free energy includes the solvation free energy of lithium ion as well as the carbon atoms and the free energy of interactions between lithium ion and the carbon atoms. The results are discussed in Chapter 8 [START_REF] Ruggeri | Multi-scale simulation of the adsorption of lithium ion on graphite surface: From quantum monte carlo to molecular density functional theory[END_REF]. This chapter is dedicated to a short introduction of concepts of statistical physics [START_REF] Pathria | Statistical mechanics[END_REF]. The introduction is intended to help understand the derivation of molecular density functional theory (MDFT), Marcus theory of electron transfer, the working principles of molecular dynamics (MD), and the applications of MD and MDFT to the study of electron transfer.

Chapter 2

Statistical physics

Entropy and temperature

Let's consider a physical system composed of N identical particles confined to a space of volume V having a total energy E. In thermodynamic limit, N → ∞, V → ∞, and the particle density ρ = N/V is fixed at a preassigned value. The total energy E of the system equals to the sum of the kinetic energy of individual particles, the potential energy between particles, and the energy due to externally applied potential. A combination of actual values of N , V , and E defines a macrostate of the system.

For a given macrostate, there is a large number of possible distributions of energies at microscopic level. Consider an ideal system with no interaction between particles and no external applied potential as an example, the total energy E is the sum of the kinetic energy of individual particles. There are many possible ways of distributing the total energy E among all the particles. Each different distribution of energy at microscopic level defines a microstate of the system and the total number of microstates may be denoted by Ω(N, V, E). The fundamental assumption in statistical mechanics is that, at any time, the system in a given macrostate with fixed N , V , and E is equally likely to be in any one of numerous possible microstates that corresponds to the macrostate.

ENTROPY AND TEMPERATURE

Now let's regard the system as composed of two subsystems A 1 and A 2 that can only exchange energy between them. The subsystem A 1 is equally likely to be in any one of Ω 1 (E 1 ) microstates while the subsystem A 2 is equally likely to be in any one of Ω 2 (E 2 ) microstates. Note that the dependencies of

Ω 1 on N 1 , V 1 and Ω 2 on N 2 , V 2 are dropped due to fixed N 1 , V 1 , N 2 and V 2 .
The whole system has a total energy E equal to the sum of the energy of the subsystems E 1 + E 2 and is equally likely to be in any one of

Ω 1 (E 1 )Ω 2 (E 2 ) = Ω 1 (E 1 )Ω 2 (E -E 1 ) = Ω(E, E 1 ) (2.1)
microstates. Now, the question that we want to answer to is that, at which E 1 , the system attains equilibrium? According to the fundamental assumption, macrostates having larger number of corresponding microstates are more probable macrostates. Actually, for a typical system, the most probable macrostate has orders of magnitude more corresponding microstates than all the other macrostates. In another word, only the most probable macrostate has a non-negligible possibility to be found when the system is in equilibrium. Thus, we can search for the value of E 1 which maximizes Ω(E, E 1 ) by taking the total derivative of Ω(E, E 1 ) with respect to E 1

∂Ω 1 (E 1 ) ∂E 1 E 1 = Ē1 Ω 2 ( Ē2 ) + Ω 1 ( Ē1 ) ∂Ω 2 (E 2 ) ∂E 2 E 2 = Ē2 • ∂E 2 ∂E 1 = 0 (2.2)
where Ē1 and Ē2 denote the equilibrium value of E 1 and E 2 , respectively. Since ∂E 2 /∂E 1 = -1, Equation 2.2 can be expressed as

∂ ln Ω 1 (E 1 ) ∂E 1 E 1 = Ē1 = ∂ ln Ω 2 (E 2 ) ∂E 2 E 2 = Ē2
( 2.3) or simply β 1 = β 2 with β defined by

β ≡ ∂ ln Ω(N, V, E) ∂E N,V,E= Ē .
(2.4) Equation 2.3 reveals that β is somehow related to the thermodynamic temperature T . Recall the thermodynamic formula

∂S ∂E N,V = 1 T (2.5)
where S is the entropy of the system. Through the following formula which links the absolute value of the entropy of the system in a macrostate to the total number of microstates accessible to it

S = k B ln Ω(N, V, E) (2.6)
where k B is the Boltzmann constant, one can easily see

β = 1 k B T . (2.7)

Ensemble theory and micro-canonical ensemble

In the framework of classical statistical physics, the total energy of the system can be calculated with classical mechanics. Consider a classical system with fixed number of particles N , volume V , and total energy E. A microstate of the system can be characterized by the full set of instantaneous position and momentum of each constituent particle. It thus requires 3N independent variables r N ≡ {r 1 , r 2 , ..., r N } to specify all the particle positions and additional 3N independent variables p N ≡ {p 1 , p 2 , ..., p N } to specify the momenta. Note that r a and p a where a = 1, 2, ..., N are three dimensional vectors that respectively denote the position and momentum of one specific particle. The 6N -dimensional space spanned by the vector {r N , p N } is called phase space. A point specified by a specific set of coordinates in the phase space is called a phase point.

We can regard an ensemble of systems as a bunch of virtual replica of the same system in the same macrostate yet each of them resides in different microstate. In another word, an ensemble of systems is a collection of all the distinct microstates that point to the same given macrostate. It should be mentioned that an ensemble does not include all the phase points in the phase space but only those pointing to the same macrostate. With the constraint on total energy, only those microstates whose Hamiltonian H(r N , p N ) equals to E is included in the ensemble.

We have mentioned the fundamental assumption in Section 2.1 which states that, for a given macrostate with fixed N , V , and E at any time, the system is equally likely to be in any one of the corresponding microstates. It is equivalent to say that the probability f N (r N , p N ) is uniform within the phase space enclosed in the ensemble. Namely, the probability of finding the system in any one of the phase points included in the ensemble is the same. f N (r N , p N ) is then called the probability density. We call the ensemble of systems with the same N , V , and E micro-canonical ensemble. The integral of probability density f N (r N , p N ) over the whole relevant phase space is restricted to 1 by introducing the partition function Q that serves as the normalization factor,

f N (r N , p N ) = 1 Q . (2.8)
The partition function Q of a micro-canonical ensemble is simply the total number of all possible microstates. The micro-canonical ensemble average of a physical quantity B(r N , p N ) which may be different for systems in different microstates is simply

⟨B⟩ = B(r N , p N )f N (r N , p N )dr N dp N = Q -1 B(r N , p N )dr N dp N .
(2.9)

CANONICAL ENSEMBLE

Canonical ensemble

It is not practical for real world experiments to fix the total energy E of the system due to the fact that it's both not straightforward to measure or control it. A better approach is to fix temperature of the system instead of energy. Such an ensemble of systems with the same N , V , and T yet variable energy is called canonical ensemble.

In thermodynamic limit, the available energy states are continuous due to its large volume. The partition function Q N of a canonical ensemble is defined as .10) where h is the Planck's constant. The N ! in the denominator in Equation 2.10 is included to account for the fact that particles are indistinguishable and h 3N is added to compensate the gain of unit upon integration. The factor h 3 is related to the smallest volume of a phase point limited by Heisenberg's uncertainty principle.

Q N = 1 h 3N N ! exp -βH(r N , p N )dr N dp N . ( 2
The probability density of finding a system in canonical ensemble having energy equal to H(r N , p N ) at thermodynamic equilibrium can be expressed as

f N 0 (r N , p N ) = Q -1 N exp -βH(r N , p N ). (2.11)
The canonical ensemble average of a physical quantity B(r N , p N ) is

⟨B⟩ = B(r N , p N ) exp -βH(r N , p N )dr N dp N exp -βH(r N , p N )dr N dp N = 1 h 3N N ! B(r N , p N )f N 0 (r N , p N )dr N dp N .
(2.12)

For example, the average energy of the system U is .13) Recall the thermodynamic relationships between entropy S and Helmholtz free energy .14) Express U as a partial differential function of only T and F , .15) The relationship between Helmholtz free energy and the canonical partition function can be derived by comparing equations 2.13 and 2.15

U = ⟨H(r N , p N )⟩ = H(r N , p N ) exp -βH(r N , p N )dr N dp N exp -βH(r N , p N )dr N dp N = - ∂ ln Q N ∂β . ( 2
F (= U -T S), S = - ∂F ∂T N,V . ( 2
U = F -T ∂F ∂T N,V = -T 2 ∂ ∂T F T N,V = ∂(F/T ) ∂(1/T ) N,V . ( 2
F = -k B T ln Q N .
(2.16)

Grand canonical ensemble

The grand-canonical ensemble is composed of systems having the same values of volume, V , temperature, T , and chemical potential of particles, µ. A system in grand-canonical ensemble may have any possible number N of particles. The Hamiltonian of the system, H(r N , p N ; N ), as well as its total energy, H(r N , p N ; N )-µN , thus also depend on the number of particles. Therefore, all the microstates having any possible N but the same V , T , and µ should be taken into account in the grand partition function. The grand partition function is defined as

Ξ = ∞ N =0 1 h 3N N ! exp -β(H(r N , p N ; N ) -µN )dr N dp N = Tr cl exp -β(H -µN ) (2.17)
where

Tr cl = ∞ N =0 1 h 3N N ! dr N dp N
is the classical trace and the dependency of H(r N , p N ; N ) is dropped to simplify the notation. The probability density f (r N , p N ; N ) of finding a system in grandcanonical ensemble having N particles and energy H at thermodynamic equilibrium can be expressed as

f (r N , p N ; N ) = Ξ -1 exp -β(H -µN ). (2.18)
Such a definition of the probability density satisfies that

∞ N =0 1 h 3N N ! f (r N , p N ; N )dr N dp N = 1. (2.19)
The grand-canonical ensemble average of a physical quantity B is

⟨B⟩ = ∞ N =0 1 h 3N N ! B(r N , p N ; N )f (r N , p N ; N )dr N dp N = Tr cl (f B) (2.20)
where the dependency of f (r N , p N ; N ) is dropped to simplify the notation. The similar quantity to the Helmholtz free energy for the grand-canonical ensemble is the grand potential Ω. The grand potential Ω is defined as (2.21) and linked to Helmholtz free energy via This chapter is dedicated to the introduction of state-of-the-art molecular density functional theory (MDFT). MDFT is a flavor of classical density functional theory (cDFT) developed for studying solvation problems. Comprehensive lectures of cDFT can be found in the book Theory of simple liquids written by Hansen and McDonald [77].

Ω = -k B T ln Ξ,
Ω = F -µN. ( 2
In Section 3.1.1, we first demonstrate that, in grand-canonical ensemble, it is possible to write down a functional of particle density for spherical mono-atomic ideal gas. The minimum of the functional can be found analytically and is equal to the grand potential. In Section 3.1.2, we incorporate two additional terms that account for the external potentials and particle interactions in the functional. The critical principle of cDFT, variational principle for grand-canonical ensemble, is then derived in Section 3.1.3 [78].

CLASSICAL DENSITY FUNCTIONAL THEORY

In the framework of MDFT, the solvent is described by its density and the solute molecules (or confining walls, electrode atoms) are considered as the external potential that cause the solvent density to deviate from the bulk density. The functional that we want to evaluate, usually called the free energy functional, is the difference between the functional of solution and the one of the homogeneous bulk solvent. Within such formalism, the minimized free energy functional is the solvation free energy of the solute. The exact mathematical form of the free energy functional is discussed in Section 3.1.4. To properly describe the solvation of molecular solutes in molecular solvent, one should take into account their molecular orientation. The free energy functional of molecular density that depends on both position and orientation of solvent molecules is discussed in Section 3.2. The part that accounts for particle interactions, the excess free energy functional, is yet complicated to be evaluated. In Section 3.1.5, a common approximation to simplify the excess free energy functional of mono-atomic solvent is introduced. The approximation is also valid and applied for molecular solvent. Finally, correction methods proposed to fix problems due to the approximation is discussed in Section 3.3. A part of the results of the thesis is on assessing the correctness of the correction methods in the study of electron transfer. The results will be discussed in Chapter 6.

Classical density functional theory

The density functional for ideal gas

Let's start with considering again a classical system in canonical ensemble containing N identical spherical particles, confined in a volume V , having a fixed temperature T . As mentioned in Section 2.1, the energy of a system H(r N , p N ) can be separated into 3 parts,

H(r N , p N ) = K(p N ) + U(r N ) + Φ(r N ) (3.1)
where K(p N ) is the kinetic energy, U(r N ) is the potential energy between particles, and Φ(r N ) is the total external energy that accounts for the energy due to externally applied potential ϕ(r). Note that K(p N ) only depends on the momentum of the particles while both U(r N ) and Φ(r N ) only depend on the position of particles.

If we consider a canonical ensemble of systems containing ideal gas with no externally applied potential, both U and Φ vanish. The kinetic energy of a system is simply the sum of the kinetic energy of each particle

K(p N ) = N i=1 p 2 i 2m (3.2)
where m is the mass of particles. The partition function defined in Equation 2.10 is CHAPTER 3. MOLECULAR DENSITY FUNCTIONAL THEORY reduced to the partition function of ideal gas

Q N,id = 1 h 3N N ! dr N exp -β N i=1 p 2 i 2m dp N = V N N ! 2mπk B T h 2 3N (3.3) = 1 N ! V Λ 3 N (3.4)
where the Gaussian integral is used to obtain the factor √ 2mπk B T in Equation 3.3, and Λ is the thermal de Broglie wavelength. Λ 3 is the smallest volume of a phase point limited by Heisenberg's uncertainty principle. We can regard Equation 3.4 as counting the total number of possible arrangements of N identical particles in a space of volume V divided into V /Λ 3 number of elementary volumes. By inserting Equation 3.4 into Equation 2.16 and applying Stirling's approximation (ln N ! ≈ N ln N -N ), we get the Helmholtz free energy of ideal gas

F id = k B T [N ln (n id Λ 3 ) -N ] (3.5)
where the equilibrium particle density of ideal gas is defined as

n id = N V (3.6)
With the thermodynamic relation that defines the chemical potential,

µ = ∂F ∂N T,V , (3.7)
we can find the chemical potential of ideal gas,

µ id = k B T ln n id Λ 3 . (3.8)
Let's now look at systems containing ideal gas in grand-canonical ensemble. The grand-partition function defined in Equation 2.17 is reduced to the grand-partition function of ideal gas 3.9) where z = e βµ /Λ 3 , defined as activity. We can see from Equation 3.8 that z = n id for ideal gas but, for grand-canonical ensemble,

Ξ id = ∞ N =0 e βµN h 3N N ! dr N exp -β N i=1 p 2 i 2m dp N = ∞ N =0 e βµN V N N !Λ 3N = ∞ N =0 (zV ) N N ! = e zV (
n id = ⟨N ⟩ V .
(3.10)

CLASSICAL DENSITY FUNCTIONAL THEORY

The grand potential of ideal gas is

Ω id = -k B T zV = -k B T n id V. (3.11)
Here, we follow Equation 2.22 and 3.5 to write down a functional of arbitrary particle number density n(r) which is itself a function of the coordinate r in the 3-dimensional space denoted as R 3 ,

Ω id [n(r)] = k B T R 3 n(r)(ln [n(r)Λ 3 ] -1)dr -µ R 3 n(r)dr = F id [n(r)] -µ R 3 n(r)dr (3.12)
where another functional F id [n(r)] is defined. For the equilibrium particle density of ideal gas n id , F id [n id ] is equal to the Helmholtz free energy of ideal gas defined in Equation 3.5 but with n id = ⟨N ⟩/V . The derivative of Ω id [n(r)] is (see Appendix B for the definition of functional derivative) ,

δΩ id [n(r)] δn(r) = k B T ln n(r)Λ 3 -µ. (3.13)
It can then be shown that when n(r) = n id = z = e βµ /Λ 3 , Ω id [n(r)] reaches the minimum. By inserting n id back into Equation 3.12 and using Equation 3.11, we get

Ω id [n(r)] = k B T R 3 n id (ln [n id Λ 3 ] -1)dr -µ R 3 n id dr = k B T n id V (βµ -1) -µn id V = Ω id (3.14)
Thus, the grand potential as well as the equilibrium particle density of the ideal gas can be found by minimizing the functional with respect to the particle density.

The density functional for non-ideal fluid

Consider a canonical ensemble of systems containing interacting particles with no external applied potential. Due to the fact that U(r N ) is independent of K(p N ), the partition function of the ensemble can be expressed as .15) where Q N,id is defined in Equation 3.4 and

Q N = 1 N ! 1 Λ 3N exp -βU(r N )dr N = Q N,id Z N V N . ( 3
Z N = exp -βU(r N )dr N (3.16)
is called the configuration integral. On taking the logarithm of Equation 3.15, it can be shown that the Helmholtz free energy of the ensemble is naturally separated into ideal and excess parts, .17) where the excess part F exc accounts purely to the interactions between particles. Now take into account the externally applied potential ϕ(r) that perturbs the system to deviate from homogeneous bulk density. The external energy is the integral of the product of the external potential ϕ(r) and the microscopic particle density n(r), Φ(r N ) = R 3 ϕ(r)n(r)dr.

F id -k B T ln Z N V N = F id + F exc . ( 3
(3.18)

The microscopic particle density n(r) is defined as

n(r) = N i=1 δ(r -r i ) (3.19)
where δ(rr i ) is the Dirac delta distribution. The equilibrium particle density n eq (r) is the ensemble average of the microscopic particle density, .20) Note that, for a fluid composed of interacting particles with no externally applied potential, the equilibrium particle density should still be homogeneous. Nevertheless, it could not be the case if an external potential is present. The average value of external energy is

n eq (r) = ⟨n(r)⟩ = Tr cl [f n(r)]. ( 3
⟨Φ(r N )⟩ = R 3 ϕ(r)n eq (r)dr. (3.21)
From Equation 3.18, we can see that Φ(r N ) is also independent of K(p N ) and U(r N ). Thus, the Helmholtz free energy F ext associated to the external potential ϕ(r) can be regarded as the third independent contribution to the Helmholtz free energy.

F = F id + F exc + F ext (3.22)
Let's turn our attention to grand-canonical ensemble. The thermodynamic relations for grand-canonical ensemble state that the infinitesimal change in Helmholtz free energy can be related to the infinitesimal change in temperature, volume and number of particles of the system. The external potential ϕ(r) presented above includes the potential that confines the fluid inside the volume V . As a result, the infinitesimal change in Helmholtz free energy can be reformulated as

∂F = -S∂T + R 3
∂ϕ(r)n eq (r)dr + µ∂N.

(3.23)

Here, we define the intrinsic free energy F int that only contains the ideal and excess parts of the Helmholtz free energy with the aid of Equation 3.21 as

F int = F - R 3 ϕ(r)n eq (r)dr (3.24)
The infinitesimal change in intrinsic free energy is thus

∂F int = -S∂T + R 3
ψ(r)∂n eq (r)dr. (3.25) where ψ(r) = µϕ(r) is defined as the intrinsic chemical potential. Inserting Equation 3.23 into the thermodynamic relation in Equation 2.22, we get

∂Ω = -S∂T - R 3
∂ψ(r)n eq (r)dr. (3.26) By taking the functional derivative of Ω and F int with respect to ψ and n, respectively, we can show δΩ δψ(r)

= -n eq (r) (3.27) and δF int δn eq (r)

= ψ(r).

(3.28)

Note that Equation 3.28 is derived with thermodynamic relations and therefore δn eq (r) means the infinitesimal change in the equilibrium state as well as the corresponding equilibrium particle density. From Equation 3.27 and 3.28, we can see that it is natural to express the intrinsic free energy as a functional of n eq (r) and the grand potential as a functional of ψ(r).

For a grand-canonical ensemble, its grand potential can also be separated into ideal, excess and external parts. Accordingly, we can write down another functional of particle density by adding two separated terms which account for particle interactions and external potential to Equation 3.12,

Ω ϕ [n(r)] = k B T R 3 n(r)(ln [n(r)Λ 3 ] -1)dr -µ R 3 n(r)dr + F exc [n(r)] + R 3
ϕ(r)n(r)dr (3.29) where the third term F exc [n(r)] is the excess part that accounts for particle interactions potential. We name the functional Ω ϕ [n(r)] full functional since it takes into account all the three contributions to the Hamiltonian of the systems. Note that n(r) is an arbitrary particle density that is not necessary equal to the equilibrium particle density n eq (r). The exact form of F exc [n(r)] will be discussed in Section 3.1.4. The last term of Equation 3.29 is written in such a form so that it reduces to the ensemble average of external energy ⟨Φ⟩ defined in Equation 3.21. In addition, its functional derivative with respect to n(r) is the external potential ϕ(r). Equation 3.29 can be further simplified to

Ω ϕ [n(r)] = F id [n(r)] + F exc [n(r)] - R 3 ψ(r)n(r)dr = F int [n(r)] - R 3 ψ(r)n(r)dr (3.30)
where F int [n(r)] is the intrinsic part of the full functional In the following section, we will show that the full functional Ω ϕ [n(r)] shown in Equation 3.29 with the exact form of F exc [n(r)] given in Section 3.1.4 has its minimum value that is equal to the grand potential Ω. The density field that gives rise to the minimum value is the equilibrium density n eq (r). Furthermore, for a given external potential ϕ(r), there is only one equilibrium density n eq (r) and one grand potential.

Variational principle for grand-canonical ensemble

Gibbs' inequality

We have introduced the equilibrium probability density f (r N , p N ; N ) in Equation 2.18. Let's now consider also all the other probability densities, all denoted as f ′ (r N , p N ; N ), with the only condition that they are normalized, i.e. Tr cl f ′ = 1. A functional of probability density can be written as an ensemble average

Ω[f ′ ] = Tr cl [f ′ (H -µN + k B T ln f ′ )].
(3.31)

At equilibrium where f ′ = f , the functional Ω[f ′ ] is equal to the grand potential,

Ω[f ′ ] = Ω[f ] = Tr cl [f (H -µN -k B T ln Ξ -H + µN )] = ⟨Ω⟩ = Ω (3.32)
where the last equality is true because the grand potential of the system Ω is a scalar. In general cases where f ′ ̸ = f , because both f ′ and f are normalized,

Tr cl f ′ = Tr cl f (3.33)
We can multiply both side of the equation by Ω,

Tr cl (f ′ Ω) = Tr cl (f Ω). (3.34)
Using Equation 2.18 and 2.21, we can replace Ω on the left hand side of the above equation, giving

Tr cl [f ′ (H -µN + k B T ln f )] = Ω. (3.35)
With Equation 3.31, the left hand side can be further transformed, resulting in .36) According to Equation 3.33, we can add Tr cl f and subtract Tr cl f ′ on the right hand side while conserving the equality. Equation 3.36 then becomes

Ω[f ′ ] = Ω + k B T • Tr cl [f ′ ln f ′ f ]. ( 3
Ω[f ′ ] = Ω + k B T • Tr cl [f f ′ f ln f ′ f ] = Ω + k B T • Tr cl [f f ′ f ln f ′ f - f ′ f + 1 ]. (3.37) Let x = f ′ /f , we can rewrite the part in parenthesis f ′ f ln f ′ f - f ′ f + 1 ≡ x ln x -(x -1) (3.38) By taking derivative, ∂(x ln x -(x -1)) ∂x = ln x (3.39)
it can be shown that x ln x ≥ (x -1) for all x ≥ 0 and that only when x -1 = 0, or say f ′ = f , x ln x = (x -1). The inequality can also be seen by plotting the two functions as shown in the Figure 3.1. As a consequence,

Ω[f ′ ] ≥ Ω (3.40)
where the equality only satisfies when f ′ = f . One-to-one relations

In previous section, we have demonstrated a key lemma of cDFT: the minimum of the functional Ω[f ′ ] is reached for the equilibrium probability density f ′ = f and is equal to the grand potential Ω. Here, we want to show that, for given µ, V , and T , the intrinsic free energy is a unique functional of the equilibrium particle density n eq (r).

The intrinsic free energy defined in Equation 3.24 can be expressed with the aid of Equation 2.22 as

F int = Ω + µ R 3 n eq (r)dr - R 3
ϕ(r)n eq (r)dr.

(3.41)

In Equation 3.41, the number of particles averaged over systems in the grandcanonical ensemble is implicitly defined as

⟨N ⟩ = R 3
n eq (r)dr.

(3.42)

The intrinsic free energy can also be expressed as an ensemble average .43) We first demonstrate that there is an one-to-one relation between external potential ϕ(r) and the equilibrium particle density n eq (r). To do so, we apply the proof by contradiction. Assume that there exists another external potential ϕ ′ (r) ̸ = ϕ(r) that gives rise to the same n eq (r). The Hamiltonian associated to ϕ ′ (r) is H ′ = K+U +Φ ′ and the associated equilibrium probability density is f ′ eq . Using Equation 3.31 and the above derived lemma, we can find the inequality,

F int = Ω + µ⟨N ⟩ -⟨Φ⟩ = ⟨k B T ln f + K + U⟩ = Tr cl [f (k B T ln f + K + U)]. ( 3
Ω = Tr cl [f (H -µN + k B T ln f )] < Tr cl [f ′ eq (H -µN + k B T ln f ′ eq )].
(3.44)

We then introduce H ′ in above equation and get, .45) On replacing the second term in Equation 3.45 by the definition in Equation 3.21 and using the assumption that ϕ ′ (r) ̸ = ϕ(r) gives rise to the same n eq (r), we arrive at

Ω < Tr cl [f ′ eq (H ′ + Φ -Φ ′ -µN + k B T ln f ′ eq )] = Ω[f ′ eq ] + Tr cl [f ′ eq (Φ -Φ ′ )]. ( 3
Ω < Ω[f ′ eq ] + R 3
[ϕ(r)ϕ ′ (r)]n eq (r)dr. (3.46) We can also derive the inequality in an inverse way giving that

Ω[f ′ eq ] = Ω[f ′ eq ] < Ω + R 3 [ϕ ′ (r) -ϕ(r)
]n eq (r)dr. (3.47) By adding the two inequality, we get

Ω[f ′ eq ] + Ω < Ω[f ′ eq ] + Ω (3.48)
which is a contradiction. Hence there is no ϕ ′ (r) ̸ = ϕ(r) that gives rise to the same n eq (r). It follows that there is no f ′ eq ̸ = f that gives rise to the same n eq (r). Namely, both n eq (r) and f are unique functionals of the external potential ϕ(r). This important conclusion implies that any quantity that can be wholly expressed as a functional of f is a unique functional of ϕ(r). From the expression of intrinsic free energy in Equation 3. [START_REF] Reeves | A first-principles investigation of the structural and electrochemical properties of biredox ionic species in acetonitrile[END_REF], we can see that the intrinsic free energy is a unique functional of ϕ(r) since K and U are independent of f . Note that the intrinsic free energy does not explicitly depends on the external potential ϕ(r). Finally, from Equation 3.41, we can see that the grand potential is also a unique functional of ϕ(r). The one-to-one relations between the external potential, the intrinsic free energy, the equilibrium probability density, the equilibrium density field, and the grand potential can be written as

ϕ(r) ↔ F int ↔ f ↔ n eq (r) ↔ Ω. (3.49)
The Variational principle

If we require that the form of the excess part F exc [n(r)] in Equation 3.30 makes F int [n(r)] evaluated at n(r) = n eq (r) being equal to the intrinsic free energy F int , then, according to Equation 3.41,Ω ϕ [n(r)] reduces to the grand potential Ω at n(r) = n eq (r). Now we want to show that, Ω ϕ [n(r)] in such a form follows the theorem derived above.

Let the arbitrary particle density n(r) appeared in Equation 3.29 be associated with a phase space probability density f ′ . The full functional in the from described above can also be written as

Ω ϕ [n(r)] = Tr cl [f ′ (H -µN + k B T ln f ′ )].
(3.50)

CLASSICAL DENSITY FUNCTIONAL THEORY

At equilibrium where f ′ = f and n(r) = n eq (r)

,

Ω ϕ [n eq (r)] = Tr cl [f (H -µN + k B T ln f )] = Ω. (3.51)
According to the one-to-one relation between the equilibrium probability density and the equilibrium particle density, the Gibbs' inequality derived with the probability density holds also for grand potentials as a functional of particle density. It follows that the full functional complies to the theorem derived above. For a given external potential, any value of the full functional in the form of Equation 3.29 evaluated at a non-equilibrium particle density is always greater than the one evaluated at the equilibrium particle density n eq (r). The minimum of the full functional is the grand potential Ω. This conclusion can be expressed as

Ω ϕ [n(r)] ≥ Ω ϕ [n eq (r)] = Ω (3.52)
On differentiating the Equation 3.30 with respect to n(r), we can find that for a given external potential ϕ(r), Ω ϕ [n(r)] is stationary with respect to variations in n(r) around the equilibrium density n eq (r),

δΩ ϕ [n(r)] δn(r) n(r)=neq = δF int [n(r)] δn(r) n(r)=neq -µ + ϕ(r) = 0. (3.53)
The second equality in the above equation is fulfilled by virtue of Equation 3.28. Equations 3.29, 3.52 and 3.53 together allow one to estimate the grand potential and the equilibrium particle density of an inhomogeneous fluid through a variational procedure. The quality of estimation depends on the applied minimization algorithm.

In the following section, the exact form of F exc [n(r)] is discussed.

Free energy functionals

To find the exact expression for F exc [n(r)], we first express the intrinsic chemical potential as the functional derivative of F int [n(r)] and separate it into ideal and excess part.

βψ(r) = β δF id [n(r)] δn(r) + β δF exc [n(r)] δn(r) = ln [n(r)Λ 3 ] -c (1) (r; n(r)) (3.54)
In the second equality of Equation 3.54, the functional derivative of F id [n(r)] defined in Equation 3.12 and the definition of the one body direct correlation function,

c (1) (r; n(r)) = -β δF exc [n(r)] δn(r) , (3.55)
are used. The two body direct correlation function is defined as the functional derivative of c (1) (r; n(r)),

c (2) (r, r ′ ; n(r)) = -β δ 2 F exc [n(r)] δn(r)δn(r ′ ) . (3.56) CHAPTER 3. MOLECULAR DENSITY FUNCTIONAL THEORY
Similarly, higher order direct correlation functions are defined as the functional derivative of the direct correlation function that is one order lower. Rearrange Equation 3.54 and use the definition of activity, one can express the equilibrium particle density as,

n eq (r) = z exp [-βϕ(r) + c (1) (r; n eq (r))] (3.57)
Comparing to the equilibrium particle density of ideal gas n id = z, we can see that the effects of external potential and particle interactions on the density field are wholly contained in two separated terms, exp [-βϕ(r)] and c (1) (r; n eq (r)).

In an attempt to find a useful expression for F exc [n(r)], we integrate Equation 3.55 with respect to density and choose a linear integration path between the state that we are interested in having equilibrium density n eq (r) and a reference equilibrium state having equilibrium density n 0 (r). The linear integration path can be expressed as n χ (r) = n 0 (r) + χ∆n eq (r) (3.58) where ∆n eq (r) = n eq (r)n 0 (r) and χ is the integration parameter varies from 0 to 1 with χ = 0, 1 designating the reference state and the targeted state, respectively. The integration gives the following expression for

F exc [n eq (r)] = F exc , F exc -F exc [n 0 (r)] = -k B T 1 0 R 3 ∂n χ (r) ∂χ c (1) (r; n χ (r))dr dχ = -k B T 1 0 R 3
∆n eq (r)c (1) (r; n χ (r))dr dχ. (3.59) where F exc [n 0 (r)] is the excess Helmholtz free energy of the reference state. We can evaluate the one body direct correlation function c (1) (r; n χ (r)) in a similar way by integrating the two body direct correlation function in Equation 3.56,

c (1) (r; n χ (r)) -c (1) (r; n 0 (r)) = χ 0 R 3 ∆n eq (r)c (2) (r, r ′ ; n χ ′ (r))dr ′ dχ ′ . (3.60)
Finally, by inserting Equation 3.60 back into Equation 3.59, we get

F exc = F exc [n 0 (r)] -k B T R 3
∆n eq (r)c (1) (r; n 0 (r))dr

-k B T 1 0 (1 -χ) R 3 R 3 ∆n eq (r)c (2) (r, r ′ ; n χ (r)) ∆n eq (r ′ )drdr ′ dχ (3.61)
with the fact that

1 0 χ 0 f (χ ′ )dχ ′ dχ = 1 0 (1 -χ)f (χ)dχ. (3.62)
It should be emphasized that no matter which integration path is chosen, the result remains the same since F exc is a unique functional of equilibrium particle density n eq (r).

We can now make the functional in Equation 3.61 depend on an arbitrary density n(r),

F exc [n(r)] = F exc [n 0 (r)] -k B T R 3 ∆n(r)c (1) (r; n 0 (r))dr -k B T 1 0 (1 -χ) R 3 R 3 ∆n(r)c (2) (r, r ′ ; n χ (r)) ∆n(r ′ )drdr ′ dχ (3.63)
with ∆n(r) = n(r)n 0 (r). When n(r) = n eq (r), F exc [n(r)] reduces to the excess part of Helmholtz free energy F exc . The full functional Ω ϕ [n(r)] that contains this functional F exc [n(r)] complies to the variational principle derived above. The full functional can be explicitly written as follow with the aid of equations 3.29 and 3.63,

Ω ϕ [n(r)] = k B T R 3 n(r)(ln [n(r)Λ 3 ] -1)dr - R 3 ψ(r)n(r)dr + F exc [n 0 (r)] -k B T R 3 ∆n(r)c (1) (r; n 0 (r))dr -k B T R 3 R 3 ∆n(r)C (2) (r, r ′ )∆n(r ′ )drdr ′ , (3.64) with C (2) (r, r ′ ) = 1 0 (1 -χ)c (2) (r, r ′ ; n χ (r))dχ.
(3.65)

If we take the homogeneous bulk fluid having the same chemical potential as the targeted fluid and the equilibrium particle density n b as the reference state, the grand potential of the homogeneous reference fluid is

Ω b = k B T R 3 n b [ln (n b Λ 3 ) -1]dr -µ R 3 n b dr + F exc [n b ]. (3.66)
On replacing n eq (r) by n b in Equation 3.54, we can show that

µ = k B T ln (n b Λ 3 ) -k B T c (1) (r; n b ). (3.67)
Then using the above equation, the full functional in Equation 3.64 with homogeneous bulk fluid as the reference state can be written in a more compact form,

Ω ϕ [n(r)] = Ω b + k B T R 3 n(r) ln n(r) n b -∆n(r) dr + k B T R 3 ϕ(r)n(r)dr -k B T R 3 R 3 ∆n(r)C (2) (r, r ′ )∆n(r ′ )drdr ′ (3.68) = Ω b + F id [n(r)] + F ext [n(r)] + F exc [n(r)] (3.69)
where ∆n(r) = n(r)n b here. The last three terms of Equation 3.69 are defined respectively as the ideal free energy functional, external free energy functional and excess free energy functional with their formulas given in the following equations

F id [n(r)] = k B T R 3 n(r) ln n(r) n b -∆n(r) dr (3.70) F ext [n(r)] = R 3 ϕ(r)n(r)dr (3.71) F exc [n(r)] = -k B T R 3 R 3 ∆n(r)C (2) (r, r ′ )∆n(r ′ )drdr ′ (3.72)
The free energy functional is defined as

F[n(r)] = Ω ϕ [n(r)] -Ω b = F id [n(r)] + F ext [n(r)] + F exc [n(r)].
(3.73)

Homogeneous reference fluid approximation

The ideal and external free energy functional defined respectively in equations 3.70 and 3.71 can be calculated analytically. However, the excess free energy functional defined in Equation 3.72 is hard to be evaluated because it involves the function C (2) (r, r ′ ). In order to calculate C (2) (r, r ′ ), it requires the knowledge of the two body direct correlation functions c (2) (r, r ′ ; n χ (r)) at all the density fields n χ (r) along the integration path (between the density of the homogeneous reference fluid and the density of the targeted state). To simplify the calculation of this function, we approximate the two body direct correlation functions at all the density fields by that of the homogeneous reference fluid, i.e., c (2) (r, r ′ ; n χ (r)) ≈ c (2) (r, r ′ ; n 0 (r)).

(3.74)

Using Equation 3.65, we get the approximated C (2) (r, r ′ ),

C (2) (r, r ′ ) ≈ 1 2 c (2) (r, r ′ ; n 0 (r)). (3.75)
The homogeneous reference fluid approximation for the two body direct correlation functions is actually equivalent to the hyper netted chain (HNC) closure relation for the key equation of liquid theory, the Orstein-Zernike relation. The two body direct correlation function in Equation 3.56 can be defined equivalently via the Ornstein-Zernike relation,

h (2) (r, r ′ ) = c (2) (r, r ′ ) + c (2) (r, r ′′ )n eq (r ′′ )h (2) (r ′′ , r ′ )dr ′′ (3.76)
where the explicit notation of the dependence on n(r) is dropped to simplify the expression and h (2) (r, r ′ ) is the two body total correlation function. Equation 3.76 separate the two body total correlation function into a direct correlation described by c (2) (r, r ′ ) and an indirect correlation. However, Equation 3.76 only relates two unknown properties h (2) (r, r ′ ) and c (2) (r, r ′ ) and to solve it, one should propose a 3.2. MOLECULAR FREE ENERGY FUNCTIONALS closure relation. h (2) (r, r ′ ) and c (2) (r, r ′ ) are also related via the following equation,

g (2) (r, r ′ ) = exp (-βv (2) (r, r ′ ) + h (2) (r, r ′ ) -c (2) (r, r ′ ) + b (2) (r, r ′ )) (3.77)
where g (2) (r, r ′ ) = h (2) (r, r ′ ) + 1 is pair distribution function, v (2) (r, r ′ ) is pair potential, and b (2) (r, r ′ ) is formally defined here as bridge function. Nonetheless, b (2) (r, r ′ ) is usually unknown. The HNC closure relation simply neglects the bridge function, setting b (2) (r, r ′ ) = 0.

We can see the HNC closure relation in another way. Starting from Equation 3.63, we can again express the two body direct correlation function c (2) (r, r ′ ; n χ (r)) by higher order direct correlation functions just like what we've done in Equation 3.60. On applying the definitions of all the higher order direct correlation functions and integrate them with respect to the density, one can write an exact expression for the function c (2) (r, r ′ ; n χ (r)) as an infinite series of higher order direct correlation functions of the reference fluid.

c (2) (r, r ′ ; n χ (r)) = c (2) (r, r ′ ; n 0 (r)) + χ R 3 ∆n(r 3 ) δc (2) (r, r ′ ; n 0 (r)) δn(r 3 ) dr 3 + ∞ n=2 χ n n! R 3 • • • R 3 δ n c (2) (r, r ′ ; n 0 (r)) δn(r 3 ) • • • δn(r n+2 ) ∆n(r 3 ) • • • ∆n(r n+2 )dr 3 • • • dr n+2 = c (2) (r, r ′ ; n 0 (r)) + χ R 3 ∆n(r 3 )c (3) (r, r ′ , r 3 ; n 0 (r))dr 3 + ∞ n=2 χ n n! R 3 • • • R 3 c (n) (r, r ′ , r 3 , • • • , r n+2 ; n 0 (r))∆n(r 3 ) • • • ∆n(r n+2 )dr 3 • • • dr n+2 (3.78)
The HNC approximation consists in restricting the Ornstein-Zernike equation to second order correlations and consider only the two body direct correlation function of the reference fluid in the above equation. We then return to the approximation of Equation 3.74.

The formal functional differentiation of Equation 3.73 with homogeneous reference fluid approximation leads to

n(r) = n b exp [-βϕ(r) + γ(r)] (3.79) with γ(r) = R 3 c (2) (r, r ′ ; n 0 (r))∆n(r ′ )dr ′ (3.80)
describing the indirect correlation. The integral Equation 3.79 is also the HNC approximation expressed in another way.

Molecular free energy functionals

Until now, we have derived the variational principle for spherical particles, proposed a free energy functional of spherical particle density F[n(r)] that can be applied in the variational procedure and introduced the homogeneous reference fluid approximation to simplify the functional. In the study of solvation problems, we can describe solvent molecules by its density field and regard the solute molecules as the external potential ϕ(r) that couples with the density and perturbs the solvent. In such a formalism, the minimized free energy functional is then the difference of grand potential between the equilibrium solution and the homogeneous bulk solvent. Therefore, by definition, it is the solvation free energy of the solute. Other sources of external potential such as confining walls or externally applied electric field can also be introduced to perturb the solvent. To properly describe the coupling between solvent molecule density and external potential, we should locate the solvent molecules in space with not only a position vector r but also an orientation Ω with respect to the fixed laboratory reference frame. The orientation of the molecule Ω is defined by the three Euler angles Ω ≡ (θ, ϕ, ψ) defined in Figure 3.2. Thus we should now perform the variational procedure on the free energy functional F[ρ(r, Ω)] of density that depends on both position and orientation ρ(r, Ω). Note that ρ(r, Ω) and n(r) are related via the following integral,

π θ=0 2π ϕ=0 2π ψ=0 ρ(r, Ω)dψdϕdθ = 8π 2 ρ(r, Ω)dΩ = n(r).
(3.81)

The above derived expression for functional F[n(r)] can be generalized to depend also on the orientation F[ρ(r, Ω)] [START_REF] Ramirez | Density functional theory of solvation in a polar solvent: Extracting the functional from homogeneous solvent simulations[END_REF][START_REF] Ramirez | Density functional theory of solvation and its relation to implicit solvent models[END_REF][START_REF] Ramirez | Direct correlation functions and the density functional theory of polar solvents[END_REF][START_REF] Gendre | Classical density functional theory of solvation in molecular solvents: Angular grid implementation[END_REF][START_REF] Zhao | Molecular density functional theory of solvation: From polar solvents to water[END_REF][START_REF] Levesque | Scalar fundamental measure theory for hard spheres in three dimensions: Application to hydrophobic solvation[END_REF]. The free energy functional can still be separated into 3 parts

F[ρ(r, Ω)] = F id [ρ(r, Ω)] + F ext [ρ(r, Ω)] + F exc [ρ(r, Ω)].
(3.82)

The ideal and external free energy functionals are redefined as follow

F id [ρ(r, Ω)] = k B T R 3 8π 2 ρ(r, Ω) ln ρ(r, Ω) ρ b -∆ρ(r, Ω) dΩdr (3.83) F ext [ρ(r, Ω)] = R 3 8π 2 ϕ(r, Ω)ρ(r, Ω)dΩdr.
(3.84)

CORRECTIONS TO THE APPROXIMATION

Note that the density of reference fluid ρ b is not the same one as n b but the number density per unit solid angle ρ b = n b /8π 2 . Equations 3.83 and 3.84 are the results of simply rewriting equations 3.70 and 3.71 with a density that depends also on the orientation. In the framework of homogeneous reference fluid approximation, the excess free energy functional can be redefined as

F exc [ρ(r, Ω)] = - k B T 2 R 3 R 3 8π 2 8π 2 ∆ρ(r, Ω) c (2) (r -r ′ , Ω, Ω ′ )∆ρ(r ′ , Ω ′ )dΩdΩ ′ drdr ′ . (3.85)
where c (2) (r -r ′ , Ω, Ω ′ ) is the bulk solvent-solvent molecular direct correlation function, which depends on the vector between the two solvent molecules and five Euler angles characterizing their relative orientation.

The expression for excess free energy functional as it stands in Equation 3.85 cannot be directly calculated for the system under study due to two reasons. On one hand, the c (2) (r -r ′ , Ω, Ω ′ ) should first be calculated by other methods for a given temperature and pressure. The implementation of our lab is to take the results of c (2) (r -r ′ , Ω, Ω ′ ) given by previous extensive Monte Carlo combined with Integral Equation bulk calculations as the input [START_REF] Belloni | Efficient full newton-raphson technique for the solution of molecular integral equations -example of the SPC/E water-like system[END_REF][START_REF] Puibasset | Bridge function for the dipolar fluid from simulation[END_REF]. On the other hand, the evaluation of the integral is numerically very expensive because it requires a double integration over the volume and the unit sphere. An efficient solution using generalized spherical harmonics expansions has been proposed to evaluate the integral and implemented in the home-made Fortran code used in our lab. The exact mathematical method that facilitates the calculation will not be discussed here but all the calculations done during the thesis project benefit from the gain of efficiency thanks to the method. Interested readers are invited to consult the paper written by Ding et al. [62].

Corrections to the approximation

The homogeneous reference fluid approximation (or equivalently the HNC closure relation) is reasonable only if the density field of the system under study remains close to the homogeneous reference fluid. This is definitely not the case when we study the solvation problem where solvent density reduces to zero at the region occupied by the solute molecule.

Overestimation of compressibility route pressure

A consequence of the homogeneous reference fluid approximation is the overestimation of the compressibility route pressure P c of the reference fluid [63,64]. The grand potential of a homogeneous fluid of volume V and pressure P is by definition -P V . For an empty system, it is zero. Thus one can deduce the compressibility route pressure of the reference fluid by evaluating the approximated free energy functional given in Equation 3.82 at zero density,

F[ρ(r, Ω) = 0] = Ω ϕ [ρ(r, Ω) = 0] -Ω b = P c V (3.86) = k B T R 3 n b dr - 1 2 k B T R 3 R 3 n 2 b c (2) b (r -r ′ )drdr ′ (3.87) = n b k B T V - 1 2 n 2 b k B T V ĉ(k = 0) (3.88)
where c

(2) b (r -r ′ ) is the two body direct correlation function of homogeneous fluid thus does not depend on the orientation, k is the Fourier space coordinate and

ĉ(k = 0) is the Fourier transform of c (2) b (r -r ′ ) taken at k = 0. The compressibility route pressure is P c = n b k B T (1 - n b 2 ĉ(k = 0)). (3.89)
With the water model we use in thesis and the temperature set at 298.15 K, we find a pressure around 11450 bar [63,64,67]. The overestimation of the compressibility route pressure of the reference fluid leads to the overestimation of solvation free energy.

Correction strategies

Several attempts were made to fix this issue, either based on simple a posteriori pressure corrections [64][65][66][START_REF] Palmer | Towards a universal method for calculating hydration free energies: a 3D reference interaction site model with partial molar volume correction[END_REF][START_REF] Ratkova | An accurate prediction of hydration free energies by combination of molecular integral equations theory with structural descriptors[END_REF][START_REF] Misin | Communication: accurate hydration free energies at a wide range of temperatures from 3D-RISM[END_REF] or by introducing higher order functional terms [63,67,68,[START_REF] Levesque | Scalar fundamental measure theory for hard spheres in three dimensions: Application to hydrophobic solvation[END_REF][START_REF] Fu | Fast prediction of hydration free energies for SAMPL4 blind test from a classical density functional theory[END_REF][START_REF] Oettel | Integral equations for simple fluids in a general reference functional approach[END_REF][START_REF] Jeanmairet | Molecular density functional theory of water describing hydrophobicity at short and long length scales[END_REF]. The a posteriori pressure correction methods lie in correct the calculated solvation free energy by subtracting a certain amount of free energy from it after the variational procedure. Sergiievskyi et al. proposed to take the partial molar volume (PMV) ∆V as the volume of solute [64]. The solvation free energy calculated with homogeneous reference fluid approximation is reevaluated by subtracting the product of overestimated pressure and the PMV of the solute. They also shown that the addition of an extra n b k B T ∆V correction term can further improve the accuracy of prediction of solvation free energy. However, the physical justification of such an additional term is unclear.

The application of HNC approximation also deteriorates the results of solvation structure. The truncation of the excess free energy functional at the second order causes the elimination of three body as well as higher order correlation terms. The a posteriori pressure correction in nature does not modify the functional and thus the equilibrium solvent density. To recover the badly predicted solvation structure, one should instead try to propose a supplementary functional, often called bridge functional F B [ρ(r, Ω)]. The nomenclature, bridge functional, comes from the bridge function described in Section 3.1.5.

Jeanmairet et al. proposed a third order bridge functional to reinforce the tetrahedral order of solvent structure that is especially important for water due to hydrogen bonding [63]. The bridge functional is constrained to recover the desired compressibility route pressure of the reference fluid. Borgis et al. examined the applicability of two other forms of bridge functional either based on the hard-sphere fluid or the weighted density approximation (WDA) for water. Efforts were put in functionals 3.4. PRACTICAL ASPECTS OF MDFT depending only on the particle number density n(r) which is reasonable for studying hydrophobic solute without partial charge. Parameters of the functionals are chosen to recover the pressure, compressibility, and the liquid-gas surface tension [67,68].

Another deficiency of simple a posteriori pressure corrections is the violation of variational principle. The estimated solvation free energy which is equal to the minimized free energy functional is modified without changing the equilibrium density field. For the studies of solution in equilibrium, the modification would not cause problem. Nevertheless, when using MDFT to study electron transfer, it is necessary to generate solvent densities with biased external potential and evaluate the solvation free energy in these densities. The work of this thesis shown that the violation of the variational principle leads to the inconsistency of the definition of the physical property defined in the framework of the electron transfer theory and the artifacts of the computed property. A discussion of the results is in Chapter 6 [69].

Practical aspects of MDFT

In this section, we discuss about some practical aspects of MDFT that may help understand the implementation of the theory and the fundamental differences between MDFT and other molecular simulation methods such as molecular dynamics (MD).

Solute and solvent models

To evaluate the external potential ϕ(r, Ω) that enters the free energy functional, one should specify the interaction potentials between solute and solvent molecules. Note that although we call ϕ(r, Ω) the external potential, it is actually the external potential energy due to the solute solvent interactions and it has the unit of energy. MDFT is a theory based on classical mechanics and therefore, like other classical molecular simulation methods, the particle interaction potentials are modeled by force fields.

Force field is the analytical expression for interaction potential. Force fields can be classified into two categories, all atom and coarse-grained force fields. All atom force fields assign one interaction site to each atom including hydrogen. United atom force fields is built by the simple coarse-graining method that collects each atom and its bonded hydrogen atoms to one interaction site. Other coarse-grained force fields gather several atoms into one interaction site and allows one to simulate large systems with moderate computational cost. Particle interaction potentials can also be categorized into two types, bonded interactions and non-bonded interactions. The bonded interactions usually include interaction potentials that describe the stretching, bending, and twisting of the chemical bonds. On the other hand, the non-bonded interactions describe the interactions between all the atoms. The non-bonded interactions usually include two types of potentials. One is the electrostatic interaction and the other describes the short-range repulsion and dispersion interaction.

In the framework of MDFT, there is no internal degree of freedom. As a result, we

ϕ LJ (d) = 4ϵ LJ σ LJ d 12 - σ LJ d 6 (3.90)
where d is the inter-particle distance, ϵ LJ is the depth of the potential well, and σ LJ is the distance at which the potential energy is zero. The rigid molecule models for water and acetonitrile (MeCN) solvent used in the thesis are the extended simple point charge (SPC/E) model and the united atom MeCN model developed by Edwards et al. [93], respectively. SPC/E water model has a O-H bond length of 1 Å and a H-O-H bond angle of 109.47 degree. The partial charges q and Lennard-Jones parameters of the SPC/E are listed in Table 3.1. The united atom MeCN model is a 3-sites linear molecular model and has a N-C bond length of 1.17 Å and a C-CH 3 bond length of 1.46 Å. The partial charges and Lennard-Jones parameters of the the united atom MeCN model are listed in Table 3.2. Note that the two body direct correlation function c (2) (r -r ′ , Ω, Ω ′ ) is unique for each solvent model at the specific temperature, pressure, and bulk density. The solute models used in the studies of this thesis will be discussed in the corresponding chapters.

Periodic boundary condition

Despite the significant improvement of computational power since the invention of molecular simulation, it is still not easy to simulate large systems. For a common system size, the portion of the region at the simulation box boundary is nontrivial and, thus, the boundary conditions do have an influence on the results. To avoid the artifacts due to the finite system size, it is common to apply the periodic boundary 3.4. PRACTICAL ASPECTS OF MDFT condition (PBC) [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]. The simulation box is treated as the primitive cell of a lattice. In each dimension with PBC, the replicas of the simulation box are infinitely duplicated and attached to each other. Any one of particles in the simulation box now interacts not only with the particles in the primitive cell but also all other particles in all other periodic cells including the replicas of itself. Note that, since MDFT is originally developed for studying bulk solutions, both Lennard-Jones and electrostatic interactions are calculated with 3 dimensional periodic boundary condition (3D-PBC). The previous calculations that give c (2) (r -r ′ , Ω, Ω ′ ) are also performed with 3D-PBC Nevertheless, due to infinite duplicates of the simulation box, the amount of interactions to be considered becomes infinite too. In practice, a cutoff distance for short-range interactions is given as an input parameter to limit the amount of interactions to be considered [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]. On the other hand, electrostatic interaction potential decreases with squared inter-particle distance and cannot be truncated in a similar way as for short-range interactions. The numerical methods to deal with the electrostatic interaction with 3D-PBC will be discussed in the next section.

In the framework of MDFT, we truncate the Lennard-Jones interaction potential at a cutoff distance d c (always smaller than the half of box length) and correct the error of calculated free energy due to the truncation. The cutoff distance is usually fixed to 10 Å. The solvation free energy is subtracted by an amount of free energy calculated with the following expression after functional minimization [START_REF] Siperstein | Long range corrections for computer simulations of adsorption[END_REF].

16 3 πn b ϵ LJ σ 6 LJ d 3 c .
Note that this term is different from the one usually used in molecular simulations [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF],

8 3 πn b σ 3 LJ ϵ LJ 1 3 σ LJ d c 9 - σ LJ d c 3 .
We don't have a factor of 1/2 due to double counting. The term to the power of nine is negligible compared to the term to the power of three and is therefore ignored.

Numerical evaluation of external potential energy

Due to the long range characteristic of the electrostatic interactions, it is not possible to truncate the electrostatic potential. The most widely used method of calculating the electrostatic potential energy is the Ewald summation [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]. It is used in the constant potential electrode model which will be introduced in Chapter 5. In the framework of MDFT, we use another better adapted method to calculate the electrostatic interaction energy with 3D-PBC.

To calculate the electrostatic part of the external potential energy, we first calculate the electrostatic potential generated by the solute molecule, V c (r), through the Poisson's equation.

-∇ 2 V c (r) = ρ c (r) ϵ 0 (3.91)
where ρ c (r) is the charge density of the solute molecule and ϵ 0 is the vacuum permittivity. The charge density of a solute molecule is defined as

ρ c (r) = Nc i=1 q i δ(r -r i ) (3.92)
where N c is the number of atoms of the solute molecule, and r i and q i are respectively the position and the partial charge of the i th atom of the solute molecule. We can express Poisson's equation in a much simpler form in Fourier space,

k 2 V c (k) = ρ c (k) ϵ 0 (3.93)
where k is the reciprocal lattice vector. Note that, in order to numerically evaluate the free energy functional, we have to use a grid to divide the space into elementary cubes and then evaluate relevant quantities at each point on the grid. As a result, ρ c (r) is interpolated to adjacent grid points and divided by the volume of the elementary cubes. The Lennard-Jones potential generated by the solute molecule is also interpolated to adjacent grid points.

The charge density of a single solvent molecule located at the origin with an orientation Ω is

ρ s (r, Ω) = Ns i=1 q i δ(r -r iΩ ) (3.94)
where N s is the number of atoms of the solvent molecule, r iΩ is the position of the i th site of this single solvent molecule, and q i now is the partial charge of this solvent site. The electrostatic part of the external potential energy is therefore the convolution of the electrostatic potential generated by the solute molecule and the charge density of the solvent molecule located at the origin,

ϕ elec (r, Ω) = V c (r ′ ) ρ s (r -r ′ , Ω)dr ′ . (3.95)
This can be more easily calculated in Fourier space and finally ϕ elec (r, Ω) is recovered by inverse Fast Fourier transform. The subroutines in the FFTW3 library are employed to perform the Fast Fourier transform and inverse Fast Fourier transform with optimized numerical efficiency [START_REF] Frigo | The design and implementation of FFTW3[END_REF].

Computational cost

The major difference between MDFT and the popular molecular simulation methods such as MD is the way of approaching the ensemble average of a property. As will be introduced in Section 4.4, to calculate the ensemble average of a property with MD, one should accumulate data of instantaneous position and momentum of all the particles. The longer the accumulation time is, the smaller the statistical error could be. However, although many physical properties can be calculated with MD, it consumes a great amount of computational resource. On the other hand, with MDFT, we compute directly the equilibrium solvent density and the solvation 3.4. PRACTICAL ASPECTS OF MDFT free energy through a variational procedure which reduces the computational cost by several orders of magnitude with respect to atomistic MD.

Another important difference between the two methods is the object that carries the information. For MD, it is the positions and velocities of particles and, for MDFT, it is the density field. Thus the consumption of the computational resource scales differently. For MD, the cost scales with respect to the number of particles and the simulation time length. The constraint of running large or long MD simulations is the calculation time spent on the central processing units (CPUs). On the other hand, for MDFT, the cost scales with respect to the number of grid points and orientations that divide the space. The principle constraint for MDFT is the available random access memory (RAM). A common choice of the elementary volume of each grid node is 0. Electron transfer is a process by which an electron moves from one atom or molecule to another. It is a key concept in redox chemistry. Redox reactions involving electron transfer can be found in many chemical systems ranging from complex biological processes, such as the electron transport chain in respiration [START_REF] Søballe | Microbial ubiquinones: multiple roles in respiration, gene regulation and oxidative stress management[END_REF][START_REF] Mitchell | Possible molecular mechanisms of the protonmotive function of cytochrome systems[END_REF] and photosynthesis [START_REF] Kurreck | Model reactions for photosynthesis-photoinduced charge and energy transfer between covalently linked porphyrin and quinone units[END_REF]100], to simple electron self-exchange between ions in solution, such as Fe 2+ /Fe 3+ self-exchange in water. In a redox flow battery or a biredox ionic liquid supercapacitor, the redox active molecules in solution undergo redox half reactions with the electrodes and transfer electrons from or to the electrodes during the charging and discharging process. Such redox half reactions can also be described by the theory of electron transfer.

Marcus theory is a widely accepted theory of electron transfer that can be used to predict the reaction rate constant k of a redox reaction involving electron transfer in solution [START_REF] Marcus | On the theory of oxidation-reduction reactions involving electron transfer. I[END_REF][START_REF] Marcus | Electrostatic free energy and other properties of states having nonequilibrium polarization. I[END_REF][START_REF] Marcus | On the theory of oxidation-reduction reactions involving electron transfer. II. applications to data on the rates of isotopic exchange reactions[END_REF][START_REF] Marcus | On the theory of oxidation-reduction reactions involving electron transfer. III. applications to data on the rates of organic redox reactions[END_REF][START_REF] Marcus | Exchange reactions and electron transfer reactions including isotopic exchange. Theory of oxidation-reduction reactions involving electron transfer. Part 4.-A statistical-mechanical basis for treating contributions from solvent, ligands, and inert salt[END_REF][START_REF] Marcus | Free energy of nonequilibrium polarization systems. II. homogeneous and electrode systems[END_REF][START_REF] Marcus | ON THE THEORY OF OXIDATION-REDUCTION REACTIONS INVOLV-ING ELECTRON TRANSFER. V. COMPARISON AND PROPERTIES OF ELECTROCHEM-ICAL AND CHEMICAL RATE CONSTANTS 1[END_REF][START_REF] Marcus | On the theory of electron-transfer reactions. VI. unified treatment for homogeneous and electrode reactions[END_REF][53]. It was developed by Rudolph Marcus since 1956, originally for describing the electron transfer between ions in solution, and later extended to electrochemical systems. The essence of Marcus theory is the separation in time scale of the electronic transition and the reorganization of solvent environment which together constitute an electron transfer reaction in solution. The two-state model is then used to describe electron transfer. It consists in tracking the free energy of the reacting species in two different oxidation states plus the solvent environment with two different free energy profiles (also called free energy curves or free energy landscapes) as a function of a global reaction coordinate. An introduction to fundamental considerations of electron transfer is given in Section 4.1.

Next, Marcus assumed that solvent behaves like a dielectric continuum and reacts linearly to the change of electronic state of the redox species. With this linear response assumption, he derived a simple expression for the activation free energy that enters the expression for rate constant. In Section 4.2, the linear response approximation of Marcus theory and its consequences are introduced. A great success of the Marcus theory is the prediction of the decrease in reaction rate as the reaction driving force increases after exceeding a certain value. This behavior is presented on the free energy profiles as the inverted region and is later discovered by experimental works. However, it is difficult to construct free energy profiles with experiments and one may resort to molecular simulations. The first step of conducting molecular simulations to study electron transfer is the determination of the proper reaction coordinate. The concept of the reaction coordinate and the related free energy profiles are introduced in Section 4. 3 [58, 101].

Vertical energy gap has been shown to be an ideal reaction coordinate for constructing free energy profiles. It was first proposed by Warshel for using molecular dynamics (MD) simulations with free energy perturbation method [54][55][56]. Since then, MD has long been the simulation method of choice for studying electron transfer. A brief introduction to the concepts of MD and its application to electron transfer is given in Section 4.4. The benefits of choosing vertical energy gap as the reaction coordinate are discussed in Section 4.5. The linear response approximation from microscopic viewpoint based on vertical energy gap is discussed in Section 4.6 [58,101].

MD is also used to investigate the influence of molecular effects that may cause deviation from Marcus theory and to build more sophisticated models based on Marcus theory [58]. However, the computational cost of MD is expensive preventing to carry on systematic studies of numerous systems. Our group recently proposed an alternative approach based on molecular density functional theory (MDFT) and applied the approach to study electron transfer half reactions of chlorine in aqueous solution [60]. The approach has shown to considerably reduce the computational cost as compared to MD. The approach based on MDFT to study electron transfer

Fundamentals of electron transfer

Electron transfer mechanism

At the center of Marcus' description of electron transfer lies the separation in time-scale of the electronic transition and the solvent reorganization. Let's take Fe 2+ /Fe 3+ self-exchange electron transfer reaction in water as an example,

Fe 2+ (aq) + Fe 3+ (aq) = Fe 3+ (aq) + Fe 2+ (aq) (4.1)
where the subscript aq denotes species in aqueous solution. A schematic representation of the reactants, products, and the surrounding solvent configurations of reaction 4.1 is shown in Figure 4.1. Note that, it is one of the solvent configurations in equilibrium with the reactants that is depicted on the left hand side of Figure 4.1.

One of solvent configurations in equilibrium with the products is depicted on the right hand side. According to the Born-Oppenheimer approximation, the electronic transition is much faster than the movement of the nuclei. If the electronic transition from ferrous ion to ferric ion happens at the solvent configuration depicted on the left hand side of Figure 4.1, the solvent configuration will be left unchanged and the energy of the system will increase due to the change in the interactions between the redox active species and the solvent molecules. If the electron transfer takes place in the dark (thermally activated), the energy of the system does not conserve and the electronic transition is forbidden. Marcus realized that, to properly describe electron transfer, not only the internal energy of the reacting species but also the interaction energy between the 4.1. FUNDAMENTALS OF ELECTRON TRANSFER Figure 4.2: Schematic representation of diabatic (solid) and adiabatic (dashed) potential energy surfaces of systems containing reactants plus solvent (R) and systems containing products plus solvent (P). The vertical energy gap defined in Equation 4.18 is also depicted. reacting species and the surrounding solvent molecules should be taken into account in the potential energy surface. Potential energy surface is the potential energy of the system as a function of the coordinates of all the atoms in the system. It is the particle-particle interaction potential energy part U(r N ) of the Hamiltonian H(r N , p N ) in Equation 3.1. Nevertheless, it should be mentioned that the degrees of freedom of the solute molecules are included here in U(r N ). Marcus then employed the two-state model to describe the process of electron transfer. The energy of the system is described by two different potential energy surfaces, one corresponds to the reactants plus solvent and the other corresponds to the products plus solvent. A schematic representation of the two potential energy surfaces is shown as solid lines in Figure 4.2 where all the atomic coordinates are collected and simplified to the abscissa in order to illustrate the potential energy surface on a 2-dimensional surface. In the case of weak electronic coupling between the reactants and products, energy conservation is satisfied only on the crossing surface (1 dimension less than the potential energy surfaces) of the two potential energy surfaces. On this crossing surface, the total energy of the reactants plus the solvent molecules is the same as the total energy of the products plus the solvent. It means that, prior to the electronic transition, the solvent configuration should reorganize to one of the configurations corresponding to the crossing surface via thermal fluctuation. Note that the reorganization does not make the system leave the thermal equilibrium but simply move the system to other microstates having higher potential energy. The reaction barrier is thus formed due to the interactions between the redox active species and the solvent molecules.

Global reaction coordinate

Later, Marcus found that it is possible to find a 1-dimensional global reaction coordinate to depict the reaction path. The ideal choice of the global reaction coordinate will be discussed in Section 4.5. The free energy of the reactants plus solvent as well as the one of the products plus solvent can then be defined along this global reaction coordinate. The two potential energy surfaces reduce to the two free energy profiles as a function of the global reaction coordinate. A schematic representation of the free energy profiles as a function of a global reaction coordinate is shown in Figure 4.3. The crossing surface of the two potential energy surfaces also reduces to a single crossing point of the two free energy profiles which represents the so-called transition state. Although there is no real formation of activated complex at the transition state, the reaction barrier introduced above allows the application of the transition state concept and Eyring equation of rate constant to describe the kinetics of electron transfer reactions. The rate constant can be expressed as

k = C • exp(- ∆G ‡ k B T ) (4.2)
where C is the pre-exponential factor that depends on the electron wave functions of the system in reactants state and product states and ∆G ‡ is the activation free energy.

The above description of the potential energy surfaces and free energy profiles is based on the condition that the electronic coupling between reactants and products 4.2. LINEAR RESPONSE APPROXIMATION is weak. It is often true when the electron transfer happens at a large separation distance between the reactants. Such potential energy surfaces and free energy profiles are called diabatic (or non-adiabatic). The occurrence of the electronic transition is restricted to the transition state located at the crossing point of the two free energy curves. However, in some cases, the electronic coupling between reactants and products is non-negligible. The two diabatic potential energy surfaces of reactants and products can mix up to form two adiabatic surfaces and a gap between them. A schematic representation of the two adiabatic potential energy surfaces is shown as dashed lines in Figure 4.2.

Linear response approximation

Inner-sphere and outer-sphere

After describing the fundamental picture of an electron transfer reaction, Marcus went on deriving an expression for the activation free energy. He divided the surrounding solvent environment into an inner-sphere and an outer-sphere. The inner-sphere includes molecules in the first coordination shell of the redox species. Upon electron transfer, the inner-sphere contributes to the energy change through the changes in the bond length and bond angle between the redox species and the molecules in the coordination shell. Outside the inner-sphere, Marcus assumed that the outer-sphere behaves like an unsaturated dielectric continuum and respond linearly to the change of electronic state of the redox species. Finally, he derived the following series of equations that are valid for both electron transfer between ions in solution, often called homogeneous electron transfer, and electron transfer between electrode and the redox species in electrolyte, often called heterogeneous electron transfer at the weak electronic coupling limit. Interested readers are referred to the original paper written by Marcus for the detailed derivation of the equations [START_REF] Marcus | On the theory of electron-transfer reactions. VI. unified treatment for homogeneous and electrode reactions[END_REF].

Reorganization free energy

The linear response approximation leads to a conclusion that the diabatic free energy profiles of both reactants and products are quadratic functions of the reaction coordinate and have the same curvature. Accordingly, a simple relation between the activation free energy, the reaction free energy ∆G, and the reorganization free energy λ,

∆G ‡ = (λ + ∆G) 2 4λ (4.3)
can be established. Note that ∆G and λ are defined on the diabatic free energy profiles which depends on the choice of the global reaction coordinate and ∆G reduces to zero for a self-exchange reaction. The reorganization free energy is defined as the free energy cost to reorganize the solvent configurations corresponding to the minimum of the free energy profile of the reactants to those corresponding to the minimum of the free energy profile of the products while staying on the diabatic potential energy surface of the reactants. On the free energy profiles, the reorganization free energy can be easily defined as the difference between the free energy at Marcus models the inner-sphere contribution as the sum of all the vibrational energy of the different vibrational modes. The outer-sphere term has the following form for electron transfer between two spherical ions in solution,

λ o = (∆e) 2 4πϵ 0 1 2a 1 + 1 2a 2 - 1 R 1 ϵ op - 1 ϵ s . (4.5)
where ∆e is the charge transferred from one reactant to the other, a 1 and a 2 are the ionic radii of reactants including the inner-sphere, R is the center-to-center separation distance, and ϵ op and ϵ s are the optical and static frequency relative permittivity of the solvent. For an outer-sphere electron transfer reaction, the electronic transition can happen at a large separation distance and the inner-sphere contribution is negligible compared to the outer-sphere contribution. For other cases, the innersphere contribution is non-negligible and the reaction is called inner-sphere electron transfer.

STATISTICAL PHYSICS FOR ELECTRON TRANSFER

Electron transfer at electrochemical interface

For electrochemical reactions, equations 4.3 and 4.5 should be slightly modified. In Equation 4.3, ∆G is replaced by a term ∆eE act where ∆e is now the charge transferred between the reactant and the electrode and E act is the activation overpotential. The activation overpotential is the potential difference between metal and solution and is defined at zero when the rate constants for the forward and reverse reactions were equal. In Equation 4.5, the term 1/2a 2 vanishes and R now denotes twice the distance between the center of the reactant's charge and the electrode surface. For outer-sphere electron transfer in which the solute does not interact with electrode, R is large enough to be assumed infinite and Equation 4.5 becomes

λ = (∆e) 2 8πϵ 0 a 1 1 ϵ op - 1 ϵ s . (4.6)
The pre-exponential factor C in Equation 4.2 is now influenced by the electronic coupling between the redox species and electrode. Comprehensive discussions about the pre-exponential factor of heterogeneous electron transfer can be found in the literature [103-105].

Although it is not possible to calculate the pre-exponential factor by classical simulations, one can assume that it is similar for a series of molecules having similar chemical structures exchanging electrons with the same electrode and solvated in the same solvent. The ratio of the rate constant of these reactions can be estimated by comparing the activation free energy of two electron transfer reactions. This concept is applied in the study of electron transfer half reactions of quinones solvated in acetonitrile (MeCN) solvent discussed in Chapter 7 to provide a prediction of rate constant as a function of molecular structure.

Statistical physics for electron transfer

Constructing free energy profiles

Free energy profiles in the two-state model are good pictorial descriptions of an electron transfer reaction. The activation free energy, the reaction free energy, and the reorganization free energy can all be calculated on the diabatic free energy profiles as shown in Figure 4.3. However, it is not easy to construct free energy profiles with experiments. Therefore, one may resort to molecular simulations for deeper understanding of electron transfer reactions.

The first step of constructing diabatic free energy profiles with molecular simulations is the determination of the global reaction coordinate. Once the global reaction coordinate is chosen, we can define the free energy profiles as a function of it. We start from recalling that from Equation 2.16 and the first equality in Equation 3.15, we can write down the free energy of the canonical ensemble of systems containing reactants plus solvent or products plus solvent,

F M = - 1 β ln j 1 N j ! 1 Λ 3N j j exp -βU M (r N )dr N (4.7)
where N j is the number of atoms of species j, Λ j is the thermal de Broglie wavelength of atomic species j determined with its mass m j , M = R denotes the reactants state, M = P denotes the products state and U M (r N ) is the diabatic potential energy surface of the reactants plus solvent (M = R) or the products plus solvent (M = P). Note that, in Equation 4.7, the degrees of freedom of the solute molecules are included in the system. Therefore, the fact that the solute and solvent are different species should be taken into account by the additional factors in the denominator of Equation 4.7. The reaction free energy change can then be related to the ratio of partition functions of reactants and products states,

∆F = F P -F R = - 1 β ln exp -βU P (r N )dr N exp -βU R (r N )dr N . (4.8)
Calculating free energy profiles from potential energy surfaces is actually a process of classifying and collecting the microstates. The definition of the diabatic free energy profiles F M (x), M = R, P, is

F M (x) = - 1 β ln j 1 N j ! 1 Λ 3N j j δ(X(r N ) -x) exp -βU M (r N )dr N (4.9)
where X(r N ) is the physical quantity that we chose to serve as the reaction coordinate, x is a specific value of X(r N ), and the Dirac delta distribution is used to classify the microstates by the value of the physical quantity X(r N ). The presence of the Dirac delta distribution can also be seen in another way as dividing the ensemble of systems according to the value of X(r N ) into numerous sub-ensembles each having a partition function in the form of

Q N (x) = j 1 N j ! 1 h 3N j δ(X(r N ) -x) exp -βH(r N , p N )dr N dp N . (4.10)
In such a formalism, microstates having X(r N ) = x are collected together and are represented by one point on the free energy profiles. The probability distribution for x can be defined as the ratio of the partition function of the sub-ensemble containing systems having X(r N ) = x to the partition function of the whole ensemble,

p M (x) = δ(X(r N ) -x) exp -βU M (r N )dr N exp -βU M (r N )dr N = ⟨δ(X(r N ) -x)⟩ M . (4.11)
Combining equations 4.7, 4.9 and 4.11 we see that

F M (x) = - 1 β ln p M (x) + F M . (4.
12)

The reaction free energy and the reorganization free energy introduced in Equation 4.3 can now be defined mathematically on the free energy profiles as

λ R = F R (x P ) -F R (x R ) (4.13) λ P = F P (x R ) -F P (x P ) (4.14) ∆G = F P (x P ) -F R (x R ) (4.15)

MOLECULAR DYNAMICS ON ELECTRON TRANSFER

where x R and x P correspond respectively to the minimum of the free energy profile of reactants and products and λ R and λ P are the reorganization free energy of the forward and the reverse reactions, respectively. Note that the reaction free energy depends on the choice of the reaction coordinate and is different from the reaction free energy change defined in Equation 4.8. Finally, the activation free energy for forward, ∆G ‡ R , and reverse, ∆G ‡ P , reactions can be defined as

∆G ‡ R = F R (x ‡ ) -F R (x R ) = F (x ‡ ) -F R (x R ) (4.16
)

∆G ‡ P = F P (x ‡ ) -F P (x P ) = F (x ‡ ) -F P (x P ) (4.17)
where x ‡ is the location of the crossing point and thus

F R (x ‡ ) = F P (x ‡ ) = F (x ‡ ).

The ideal global reaction coordinate

We now have the definition of the free energy profile and the related free energies that depends on the reaction coordinate. The remaining question is what is the ideal reaction coordinate? Warshel first suggested that the vertical energy gap of the diabatic potential energy surfaces ∆E(r N ) defined as

∆E(r N ) = U P (r N ) -U R (r N ). (4.18)
is the ideal global reaction coordinate [54]. It is originally proposed for using MD simulations with free energy perturbation method to construct free energy profiles. Note that the vertical energy gap can be positive or negative and that it equals zero on the crossing surface of the two potential energy surfaces. In another word, the crossing surface is completely identified by the configurations having zero vertical energy gap. By using vertical energy gap of diabatic potential energy surface as the global reaction coordinate, the whole crossing surface can be collected to one single point which is the crossing point of the two free energy profiles. Before discussing further the benefits of using the vertical energy gap as the reaction coordinate, the concept of MD should be introduced.

Molecular dynamics on electron transfer

Working principles of molecular dynamics

The procedure for running a MD simulation is to first distribute the particles inside a simulation box according to a predefined density. The initial momentum of each particle as well as the particle interaction potentials should also be specified before starting the simulation. In the second step, we discretize time as composed of elementary time steps and let the system evolve as time passes according to Newton's equations of motion until the system reaches thermodynamic equilibrium. Finally, we let the system continue to evolve in order to sample different microstates of the thermodynamic equilibrium state and we store the instantaneous position and momentum of all the particles. The evolution of the system at equilibrium over time is called a trajectory and, in practice, we only store data per thousands or per tens of thousands of time steps to reduce the memory consumption.

A MD simulation is very similar to a real world experiment. However, physical quantities of the system that we measure with experiments are related to the ensemble average of them, ⟨B(r N , p N )⟩, rather than the value evaluated at a specific time, B(r N , p N ). For example, the diffusion coefficient of an atomic fluid can be calculated from the ensemble average of squared displacement of atoms. Using Einstein's formula,

D = 1 6 lim t→∞ ∂⟨|r i (t) -r i (0)| 2 ⟩ ∂t (4.19)
where r i (t) is the position of i th particle at time t. In order to improve the estimation of the diffusion coefficient and reduce the statistical error, we run the MD simulations over a long enough period of time. This is also what we do in a real word experiment to reduce error. The above mentioned procedure actually corresponds to letting the system at an initial phase point to travel within the phase space as time passes. Thus, one should also specify the ensemble that the simulation sticks to. Without any additional constraint, a MD simulation is naturally conducted in the micro-canonical ensemble. To perform simulations in the canonical ensemble, one should implement thermostats that exchange energy with the system to ensure constant temperature. A thorough discussion about the implementation of the thermostat can be found in the literature [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]106,107].

Application to electron transfer

The ensemble that MD simulations of electron transfer stick to is often set to the canonical ensemble. Since physical quantities evaluated at a specific time of a MD trajectory do not have fixed values but fluctuate as time passes, we can measure not only the ensemble average but also the probability distribution of the physical quantities. In the case of electron transfer, we want to calculate the probability distribution of the value of the reaction coordinate, p M (x), as a function of x. They are both defined in Equation 4.11.

The histogram method is used to construct the probability distribution. We evaluate X(r N ) at each time step along the trajectory and assign each step, according to its X(r N ), to the proper bin of the histogram. We then get the number of steps in each bin as a function of x and the probability distribution is obtained by normalizing the histogram. In order to obtain a good estimation of the probability distribution, one should let the simulations run for a long enough period of time. As the simulation time increases, the statistical error reduces. Finally, we get two probability distributions, one for the system containing reactants plus solvent and the other for the system containing products plus solvent, from the two MD simulations.

According to Equation 4.12, the peak position of probability distributions corresponds to the minimum of the free energy profiles. From Equation 4.12, we can also get the difference between the two free energy profiles which depends on the reaction free energy change, where subscript R and P attached to the brackets denote that the canonical ensemble averages are calculated over all the time steps along the trajectories of the reactants and products state. Note that a MD simulation trajectory is actually a trajectory of a phase point in the phase space or, in another word, a trajectory on the potential energy surface. Due to the thermal fluctuation of the equilibrium state, we can visit a region of the potential energy surface along a trajectory. For some systems, the region of the two potential energy surfaces accessible to and visited by the trajectories of the two simulations may overlap. In the overlapped region, the potential energy surface of reactants has the same vertical energy gap as the potential energy surface of products so that it is also possible to estimate the relative position of the two free energy profiles in the ordinate using equations 4.21 and 4.22. Note that the quality of the estimation of the relative position depends on the quality of the estimation of the reaction free energy change. The overlap of the two potential energy surfaces leads to the overlap of the tails of the two probability distributions. In such cases, the two free energy profiles obtained by using Equation 4.12 cross each other at a specific value of the reaction coordinate. However, for other systems, the region of the two potential energy surfaces sampled by the two simulations as well as the tails of the two probability distributions may not overlap. It is therefore necessary to apply some enhanced sampling techniques, such as the thermodynamic integration or the free energy perturbation method combined with umbrella sampling [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF], to bridge the two free energy profiles [101]. The systems studied with MD simulations in the scope of the thesis are of the first type.

∆F (x) = F P (x) -F R (x) = - 1 β ln p P (x) p R (x) + ∆F. ( 4

Vertical energy gap

To illustrate the benefits of using vertical energy gap as the reaction coordinate, using Equation 4.9, we rewrite the difference between the free energy profiles of the two states as

∆F (x) = - 1 β ln δ(X(r N ) -x) exp -βU P (r N )dr N δ(X(r N ) -x) exp -βU R (r N )dr N = - 1 β ln ⟨δ(X(r N ) -x) exp -β∆E(r N )⟩ R ⟨δ(X(r N ) -x)⟩ R . (4.23)
On choosing X(r N ) = ∆E(r N ), we get

∆F (x) = - 1 β ln ⟨δ(∆E(r N ) -x) exp -β∆E(r N )⟩ R ⟨δ(∆E(r N ) -x)⟩ R . (4.24) CHAPTER 4. ELECTRON TRANSFER
For a given value ϵ of vertical energy gap, the above equation can be simplified to

F P (ϵ) -F R (ϵ) = - 1 β ln [⟨exp -βϵ⟩ R ] = ϵ.
(4.25) Equation 4.25 is a direct consequence of choosing vertical energy gap as the reaction coordinate and leads to the fact that the relative position of the two free energy profiles in the ordinate is exactly the vertical energy gap. Hence, the crossing point of the the two free energy profiles always locates at the zero vertical energy gap, ϵ = 0. Once we decide to use vertical energy gap as the reaction coordinate, the reaction free energy defined in Equation 4.15 can be rewritten as

∆G ϵ = F P (ϵ P ) -F R (ϵ R ) (4.26)
where ∆G ϵ is the reaction free energy defined with respect to vertical energy gap, and ϵ P and ϵ R are the values of vertical energy gap corresponding to the minimum of the free energy profiles of reactants and products, respectively. We can also rewrite the reorganization free energy defined in equations 4.13 and 4.14. Using Equation 4.26,

λ R,ϵ = F R (ϵ P ) -F R (ϵ R ) = ∆G ϵ -ϵ P (4.27) λ P,ϵ = F P (ϵ R ) -F P (ϵ P ) = -∆G ϵ + ϵ R .
(4.28) Equations 4.26, 4.27 and 4.28 imply that if we only want to calculate the reaction free energy and the reorganization free energy, there is no need to construct the whole free energy profiles. It is especially beneficial for the systems whose two probability distributions do not overlap since there is no more need to use enhanced sampling techniques. Finally, the activation free energies defined in equations 4.16 and 4.17 can be rewritten as

∆G ‡ R,ϵ = F R (0) -F R (ϵ R ) = F (0) -F R (ϵ R ) (4.29) ∆G ‡ P,ϵ = F P (0) -F P (ϵ P ) = F (0) -F P (ϵ P ) (4.30)
where F R (0) = F P (0) = F (0). The application of vertical energy gap also provides a benefit in improving the sampling efficiency [101]. Let's take a forward reaction as an example and we perform a MD simulation of the reactants state. A trajectory generated without using any enhanced sampling technique usually only visits the region of the potential energy surface that correspond to the region around the minimum of the free energy profile. Therefore, we get the region around the minimum of the free energy profile of the reactants state with our simulation. Meanwhile, at each time step of the trajectory, we can switch the interaction potential between the redox active species and the surrounding solvent molecules from the one of reactants state to the one of the products state. We evaluate the energy difference upon the switch and it is exactly the vertical energy gap. So from a trajectory of the reactants state, we also get a handful of microstates of the products states by applying the procedure. We can build a probability distribution of these microstates according to Equation 4.11 with X(r N ) = ∆E(r N ). The part of the free energy profile of the products state near the 4.6. MICROSCOPIC LINEAR RESPONSE APPROXIMATION the minimum of the free energy profile of the reactants state which is normally very rarely sampled with the MD simulation of the products state can thus be constructed with the knowledge of the probability distribution and the vertical energy gap. The same procedure can also be applied to the MD simulation of the products state. The sampling efficiency is thus doubled.

Microscopic linear response approximation

From microscopic point of view, the linear response approximation of Marcus theory is equivalent to assuming the probability distribution for vertical energy gap constructed with the simulations of the reactants (M = R) or the products (M = P) state by a Gaussian

p M (ϵ) = 1 σ M √ 2π exp - (ϵ -ϵ M ) 2 2σ 2 M (4.31)
where σ 2 M and ϵ M are the variance and the mean value of the Gaussian of the reactants (M = R) or the products (M = P) state. The mean value of the Gaussian converges to the values of vertical energy gap corresponding to the minimum of the free energy profiles. According to the features of Gaussians,

ϵ M = ⟨∆E(r N )⟩ M (4.32) σ 2 M = ⟨(∆E(r N ) -⟨∆E(r N )⟩ M ) 2 ⟩ M . (4.33)
The free energy profiles built with the Gaussian probability distribution is obtained using Equation 4.12,

F M (ϵ) = F M + k B T 2σ 2 M (ϵ -ϵ M ) 2 + k B T 2 ln (2πσ 2 M ). ( 4

.34)

Now inserting Equation 4.34 into Equation 4.25 leads to

ϵ = ∆F + k B T 2 (ϵ -ϵ P ) 2 σ 2 P - (ϵ -ϵ R ) 2 σ 2 R + ln σ 2 P σ 2 R . (4.35)
In order to guarantee that the above equality is true, the quadratic term of ϵ should be zero. Consequently,

σ 2 P = σ 2 R ≡ 2k B T λ ′ (4.36)
where the linear response reorganization free energy, λ ′ , is defined. So Equation 4.34 can be rewritten as

F M (ϵ) = F M + 1 4λ ′ (ϵ -ϵ M ) 2 + k B T 2 ln (4k B T πλ ′ ). (4.37)
In general, λ ′ is not identical to the reorganization free energy λ R,ϵ and λ P,ϵ defined respectively in equations 4.27 and 4.28. If and only if the probability distributions constructed with the simulations are strictly Gaussian and have the same variance, these three quantities are the same. In this situation,

λ R,ϵ = λ P,ϵ = λ ′ = λ (4.38)
and both free energy profiles have the curvature 1/4λ and we return to Marcus' linear response approximation. Putting back Equation 4.36 into Equation 4.35, we get

ϵ = ∆F + 1 4 (ϵ -ϵ P ) 2 λ - (ϵ -ϵ R ) 2 λ (4.39)
At the minimum of the two free energy profiles,

ϵ R = ∆F + (ϵ R -ϵ P ) 2 4λ (4.40) ϵ P = ∆F - (ϵ P -ϵ R ) 2 4λ ( 4
∆G ϵ = λ + ϵ P = λ + ⟨∆E(r N )⟩ P = -λ + ⟨∆E(r N )⟩ R . (4.46)
While one should perform two simulation runs for the two states to calculate the reaction free energy and the reorganization free energy from equations 4.42 and 4.44, only one simulation run is required from equations 4.36 and 4.46. The reaction free energy and the reorganization free energy can be both obtained directly from the mean and variance of the vertical energy gap fluctuations of a single but long enough trajectory, either reactants or products state. Finally, using Equation 4.29, 4.30, and evaluating Equation 4.37 at ϵ = 0, ϵ R , and ϵ P , we get

∆G ‡ R,ϵ = ϵ 2 R 4λ = (λ + ∆G ϵ ) 2 4λ (4.47) ∆G ‡ P,ϵ = ϵ 2 P 4λ = (λ -∆G ϵ ) 2 4λ . (4.48)
Although many systems can be well described by the Marcus theory derived from the linear response approximation, there are still some circumstances where the approximation could breakdown. A direct way to investigate the deviation from 4.7. MDFT ON ELECTRON TRANSFER the Marcus theory is to construct free energy profiles with molecular simulations and measure the deviation from the ideal two identical curvature parabolas. There are also more sophisticated models that are developed with the aid of molecular simulations to explain the deviation of the free energy profiles [108][109][110][111]. In the following chapters, it will be shown that all the electron transfer reactions studied in the scope of the thesis can be well described by Marcus

MDFT on electron transfer

MD has long been the simulation method of choice for studying electron transfer. However, despite the improved sampling efficiency provided by the usage of vertical energy gap described in Section 4.5, it remains necessary to run long simulations to build free energy profiles. This makes the study of electron transfer with MD rather costly which prevents to carry on systematic studies of large amount of molecules. Our group recently proposed an alternative approach based on MDFT and demonstrated that it is an efficient alternative to tackle this problem [59,60]. This section is dedicated to the description of the approach based on MDFT.

Reexamination of the statistical physics

To illustrate the framework of using MDFT to study electron transfer, we first reexamine the concepts of statistical physics introduced in Section 4.3. Recall that the construction of free energy profiles with vertical energy gap as the reaction coordinate is actually a process of classifying and collecting microstates according to the value of vertical energy gap associated to them. Recall also that, in Equation 4.7, we include the degrees of freedom of solute molecules in the particle-particle interaction potential energy U(r N ). On using vertical energy gap as the reaction coordinate, the Dirac delta distribution in Equation 4.9 only requires the microstates to have the given value of vertical energy gap but different solute geometries are accepted. For example, in Chapter 7, the MD simulations of anthraquinone (AQ)/anthraquinone anion (AQ -) electron transfer half reaction solvated in MeCN solvent are performed with force fields including bonded interactions. As a result, microstates corresponding to the same point on the free energy profile of AQ may have AQ in different geometries. This is also true for the free energy profile of AQ -.

Let's now consider rigid solutes and take the solute molecule as the center of the laboratory frame. In such condition, the degrees of freedom of the solute reduce to zero, r N no longer includes the coordinates of solute molecules, and the internal energy of the solute molecules is constant. The dimension of the potential energy surfaces of the systems considering rigid solutes are thus less than the dimension of the potential energy surfaces that take the degrees of freedom of the solute into account. The part of solute-solvent interactions can be separated from the allparticle interactions U(r N ) and be denoted by Φ(r N ) as in Equation 3.1. This CHAPTER 4. ELECTRON TRANSFER separation of interactions corresponds to the formulation of MDFT which states that U(r N ) in Equation 3.1 only includes the interactions between solvent molecules and that solute molecules are presented as the external potential. Nonetheless, in Chapter 7, it will be shown that the free energy profiles of AQ/AQ -electron transfer half reaction solvated in MeCN solvent constructed with MDFT and with MD using force fields including bonded interactions are very similar.

The consideration of only rigid solute molecules also has other implications. First, since we don't have access to internal degrees of freedom of solute molecules, we assume that the inner-sphere contribution to the free energy profiles of an electron transfer reaction constructed with MDFT to be zero. Such consideration of electron transfer reactions corresponds to the outer-sphere electron transfer and are completely controlled by the solvent response. The free energy profiles become the solvation free energy as a function of a solvent reaction coordinate which is fully characterized by the solvent configurations. In the framework of MDFT, the solvent reaction coordinate is thus directly related to the solvent density field. Secondly, vertical energy gap can be written as

∆E(r N ) = Φ P (r N ) -Φ R (r N ) (4.49)
where Φ M (r N ) is the external energy in Equation 3.1 due to the solute-solvent interactions with solute in reactant states (M = R) or in products state (M = P). However, in the framework of MDFT, we obtain directly the ensemble averaged value of some physical quantity instead of the instantaneous value. Thus we shall use the ensemble averaged vertical energy gap as the solvent reaction coordinate.

Average vertical energy gap

The vertical energy gap is the energy change upon electronic transition while all the atomic coordinates are frozen. We can follow the same concept to define the average vertical energy gap. The average vertical energy gap, ⟨∆E⟩, is defined as the free energy change upon changing the external potential from reactants state to products state while keeping the solvent density field unchanged. The free energy change can be evaluated by using the free energy functional defined in Equation 3.82. We can see that, for a given solvent, the solvent density field is the only variable in equations 3.83 and 3.85. As a consequence, for a given solvent density field, the ideal and excess free energy functional are constants and the average vertical energy gap can be expressed as

⟨∆E⟩ = ρ(r, Ω)[ϕ P (r, Ω) -ϕ R (r, Ω)]drdΩ (4.50)
where ϕ M (r, Ω) is the external potential generated by the solute molecules in reactant states (M = R) or products state (M = P). The average vertical energy gap defined as such in the framework of MDFT actually equals to the ensemble average of the vertical energy gap defined in Equation 4.18 considering rigid solute molecules (or simply the ensemble average of the vertical energy gap defined in the Equation 4.49) at the thermodynamic limit (N → ∞,

MDFT ON ELECTRON TRANSFER

V → ∞). Note that the difference between canonical ensemble and grand-canonical ensemble vanishes at the thermodynamic limit. As a result, all the equations and conclusions derived in sections 4.5 and 4.6 are also true in the framework of MDFT.

Constructing free energy profiles with MDFT

Equation 4.50 implies that, in order to construct free energy profiles, we have to explore solvent density fields other than those two that we can obtained with the two external potentials ϕ R (r, Ω) and ϕ P (r, Ω). The way to do so is to generate a series of biased external potentials, ϕ η (r, Ω), and find the associated equilibrium solvent density fields, ρ η (r, Ω). For example, it is an easy way to linearly interpolate the two external potentials of reactants and products to generate biased external potentials.

ϕ η (r, Ω) = ϕ R (r, Ω) + η[ϕ P (r, Ω) -ϕ R (r, Ω)] (4.51)
where η is the interpolation coefficient. At η = 0,

ϕ 0 (r, Ω) = ϕ R (r, Ω) ρ 0 (r, Ω) = ρ R (r, Ω)
where ρ R (r, Ω) is the equilibrium solvent density of reactants state. At η = 1,

ϕ 1 (r, Ω) = ϕ P (r, Ω) ρ 1 (r, Ω) = ρ P (r, Ω)
where ρ P (r, Ω) is the equilibrium solvent density of products state. It should be mentioned that no matter how we generate the series of biased external potentials, all the corresponding solvent density fields are not in equilibrium with the solute in reactants or products states. As a consequence, the minimum of a free energy profile constructed with MDFT corresponds to the equilibrium solvation state of the solute. All the other solvent density fields represented by the points on the reaction coordinate are out-of-equilibrium with the solute. The solvation free energy of the solute increases as we move away from the minimum of the free energy profile. The more the solvent density differs from equilibrium, the higher the solvation free energy is.

As mentioned above, in the framework of MDFT, the solvent reaction coordinate is directly related to the solvent density field. According to Equation 4.50, for a given biased external potential, ϕ η (r, Ω), average vertical energy gap is uniquely determined by the associated equilibrium solvent density field. We have shown in Section 3.1.3 that there is an one-to-one relation between the external potential and the equilibrium solvent density field. We can follow the same reasoning in Section 3.1.3 to show the one-to-one relation between the external potential and the average vertical energy gap. Let's assume there exists two external potentials ϕ η (r, Ω) and ϕ η ′ (r, Ω) with η ̸ = η ′ that give rise to the same average vertical energy gap. Using Equation 4.51, an inequality similar to inequality 3.46 can be derived,

Ω η < Ω η ′ + (η -η ′ ) ρ η ′ (r, Ω)[ϕ P (r, Ω) -ϕ R (r, Ω)]drdΩ. (4.52)
Similarly to Equation 3.47, we can derived,

Ω η ′ < Ω η + (η ′ -η) ρ η (r, Ω)[ϕ P (r, Ω) -ϕ R (r, Ω)]drdΩ. (4.53)
Sum up the above two inequalities, we get

0 < (η ′ -η) [ρ η (r, Ω) -ρ η ′ (r, Ω)][ϕ P (r, Ω) -ϕ R (r, Ω)]drdΩ. (4.54)
Using the definition of average vertical energy gap in Equation 4.50 and the assumption that the two external potentials give rise to the same average vertical energy gap, we arrive at a contradiction. Accordingly, for the series of biased external potentials generated by linear interpolation in Equation 4.51, there is an one-to-one relation between the interpolation coefficient η, the biased external potential, the average vertical energy gap, and the quantities in Equation 3.49. The free energy functional can thus be equivalently regarded as a function of average vertical energy gap,

F(⟨∆E⟩ η ) ≡ F[ρ η (r, Ω)]. (4.55)
The above conclusion shows that average vertical energy gap is a particularly convenient reaction coordinate. The free energy profiles can be directly constructed by evaluating the free energy functionals at various different solvent density fields, ρ η (r, Ω). Note that although the linear parametrization of the series of biased external potential is chosen for its simplicity, Equation 4.55 and the one-to-one relation are globally true for all strictly increasing continuous parametrization functions. A proof of the argument is given in the Appendix B of reference [60].

Definitions of numerous free energies in MDFT

As stated earlier, all the equations and conclusions derived in sections 4.5 and 4.6 are also true in the framework of MDFT. It's helpful to rephrase the equations in the language of MDFT. The difference between the free energy profiles of the two states in Equation 4.25 can be written as

F P [⟨∆E⟩ η ] -F R [⟨∆E⟩ η ] = F P [ρ η ] -F R [ρ η ] = ⟨∆E⟩ η (4.56)
where F M is the free energy functional with solute in reactants state (M = R) or with solute in products state (M = P). Note that the internal energy change of the solutes upon electronic transition is ignored in Equation 4.56. The reaction free energy in Equation 4.26 can be written as

∆G = F P (⟨∆E⟩ P ) -F R (⟨∆E⟩ R ) = F P [ρ P ] -F R [ρ R ]. (4.57)
The reorganization free energies in equations 4.27 and 4.28 can be written as

λ R = F R (⟨∆E⟩ P ) -F R (⟨∆E⟩ R ) (4.58) = F R [ρ P ] -F R [ρ R ] (4.59) = ∆G -⟨∆E⟩ P λ P = F P (⟨∆E⟩ R ) -F P (⟨∆E⟩ P ) (4.60) = F P [ρ R ] -F P [ρ P ] (4.61) = -∆G + ⟨∆E⟩ R .

MDFT ON ELECTRON TRANSFER

Finally, the activation free energies in equations 4.29 and 4.30 can be written as

∆G ‡ R = F R (0) -F R (⟨∆E⟩ R ) = F(0) -F R (⟨∆E⟩ R ) (4.62) ∆G ‡ P = F P (0) -F P (⟨∆E⟩ P ) = F(0) -F P (⟨∆E⟩ P ) (4.63)
where F R (0) = F P (0) = F(0). Equations derived under the linear response approximation can also be employed with MDFT to calculate free energies [59]. Equations 4.42,4.44,and 4.46

∆G = λ + ⟨∆E⟩ P = -λ + ⟨∆E⟩ R . (4.66)
The practical process of constructing free energy profiles with MDFT is summarized as follow. We first perform the functional minimization for two equilibrium states. Meanwhile, both external potential, ϕ R (r, Ω) and ϕ P (r, Ω), are calculated. Next, we linearly interpolate the external potential field with several values of η and perform the functional minimization to obtain the associated series of equilibrium solvent density field ρ η (r, Ω). Finally, we evaluate the free energy using Equation 3.82 with both external potential, ϕ P (r, Ω) and ϕ R (r, Ω), and the solvent density field ρ η (r, Ω). We get as such a series of points on the profiles, F R [ρ η (r, Ω)] and F P [ρ η (r, Ω)].

Chapter 5

Constant potential electrode model The power density of an electrochemical energy storage device is closely related to the kinetics of the underlying electron transfer reactions at the electrode/electrolyte interfaces. Marcus theory introduced in Chapter 4 provides a theoretical background to understand and predict the electron transfer reaction kinetics at the electrochemical interfaces. Notably, the simple Equation 4.6 is derived for calculating the reorganization free energy of a heterogeneous electron transfer reaction based on the linear response approximation, the assumption of spherical solute, and the assumption of outer-sphere electron transfer. However, this simple picture not only ignores the complex structure of the electric double layer formed at the interfaces [START_REF] Jeanmairet | Microscopic simulations of electrochemical double-layer capacitors[END_REF][START_REF] Merlet | The electric double layer has a life of its own[END_REF] but also assumes the electrode to be a smooth perfect metal. To go beyond and capture the complex behaviour of electrode/electrolyte interfaces, we can resort to molecular simulations. In Section 5.1, we introduce the commonly used electrode and electrolyte models. Then in Section 5.2, we discuss the influence of choice of periodic boundary condition (PBC) on the evaluation of electrostatic interaction energy. The constraints that impose constant potential in electrodes and electroneutrality are also shown. Finally, the strategy and pratical aspects of incorporating constant potential electrode model in molecular density functional theory (MDFT) are discussed in Section 5.3.

ELECTRODE AND ELECTROLYTE MODELS

Electrode and electrolyte models

Molecular dynamics (MD) simulations are intensely used to study the electrochemical interfaces at atomistic resolution [START_REF] Jeanmairet | Microscopic simulations of electrochemical double-layer capacitors[END_REF]113]. A popular way of modelling the electrostatic interaction potentials generated by electrodes is proposed by Siepmann and Sprik [70]. It consists in pinning the electrode atoms in space and assigning to each electrode atom, labelled j, a Gaussian charge distribution,

ρ j (r) = Q j ( σ 2 Q π ) 3 2 e -σ 2 Q (r-r j ) 2 . (5.1)
where r j and Q j are the position and the integral charge of electrode atom j, respectively, and σ Q is the Gaussian width that can be tuned to account for the metallicity of the electrode [START_REF] Serva | Effect of the metallicity on the capacitance of gold-aqueous sodium chloride interfaces[END_REF]. Given the alternative definition of the electric potential equivalent to Equation 3.91 in Gaussian unit,

V elec (r) = R 3 ρ j (r ′ ) |r -r ′ | d 3 r ′ , (5.2)
an electrode composed of N Q atoms carrying Gaussian charge distributions generates an electric potential in space that can be expressed as

V Q (r) = N Q j=1 R 3 ( σ 2 Q π ) 3 2 Q j |r -r ′ | e -σ 2 Q (r ′ -r j ) 2 d 3 r ′ = N Q j=1 Q j |r -r j | erf(σ Q |r -r j |). (5.3)
In the framework of classical molecular dynamics (classical MD) simulations, the electrolyte molecules are explicitly presented in the system. A simple way to model the electrostatic interaction potentials generated by the electrolyte is to assign a point charge to each atom, labelled k, of the electrolyte molecules. The charge distribution of a point charge is given by

ρ k (r) = q k δ 3 (r -r k ) (5.4)
where r k and q k are the position and the partial charge of electrolyte atom k, respectively, and δ 3 (r-r k ) is the three dimensional Dirac delta distribution centered at point r k . Following Equation 5.2 and replacing ρ j (r) by ρ k (r), we find the electric potential generated by the electrolyte containing N q of atoms,

V q (r) = Nq k=1 q k |r -r k | .
(5.5)

Periodic boundary condition (PBC)

We mentioned in Section 3.4 that it is common to apply PBC to avoid artifacts due to the finite system size. For Lennard-Jones interactions, we can truncate the potential at a finite distance and correct the energy afterwards. However, we can see that in equations 5.3 and 5.5, the electric potential is a function of 1/r. The long range characteristic of the Coulomb interactions prohibit us to truncate them at a finite distance. A common way to calculate the Coulomb interactions efficiently is to use Ewald summation.

Electrostatic interactions in three dimensional PBC

Ewald summation method is originally developed to calculate the electrostatic interaction energy between point charges, U qq , with 3 dimensional periodic boundary condition (3D-PBC) [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]. We denote the length of an orthorhombic simulation box in the x, y, z directions by l x , l y , l z and the lattice vector of a periodic cell by n ≡ (n x l x , n y l y , n z l z ), where n x , n y , n z ∈ Z. With 3D-PBC, U qq can be expressed as

U qq = 1 2 Nq k=1 Nq k ′ =1 n ′ q k q k ′ |r k -r k ′ + n| (5.6)
where the primed sum indicates that for n = (0, 0, 0) the case k = k ′ is excluded and the summation of n runs over all periodic cells. Equation 5.6 can also be expressed in terms of electric potential,

U qq = 1 2 Nq k=1 n ′ q k V q (r k + n) (5.7)
where V q (r k + n) here is the electric potential at the position r k + n generated by all the periodic images of all the electrolyte atoms except for the atom located at r k + n.

Ewald summation method consists in adding a Gaussian charge distribution of the opposite sign around every atoms to completely cancel its charge and then correct the artifacts due to these additional charge distributions. The technique smooths the electric potential generated by point charges and improve the numerical stability. The concept is illustrated in Figure 5.1. With the above operation, the electrostatic potential energy between point charges can be separated into two parts, one accounts for the short-range interaction energy that is calculated in real space and the other

PERIODIC BOUNDARY CONDITION (PBC)

accounts for the long-range interaction energy that is calculated in Fourier space. In order to treat the long range interaction part in Fourier space properly, one should include the self interaction between each atom and the Gaussian charge distribution around it. This contribution should be computed separately and subtracted from the long range part. Equations of Ewald summation method for calculating the electrostatic interaction energy between point charges with 3D-PBC can be found in reference [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF].

The electrostatic interaction energy of a system containing Gaussian distributed charges and point charges is constituted of three contributions [71, 73]

U elec = U qq + U QQ + U Qq (5.8)
where U QQ is the electrostatic interaction energy between Gaussian distributed charges and U Qq is the cross electrostatic interaction energy between Gaussian distributed charges and point charges. With 3D-PBC, for single electrode, U Qq can be expressed as

U Qq = Nq k=1 N Q j=1 n R 3 ( σ 2 Q π ) 3 2 q k Q j |r k -r j -u + n| e -σ 2 Q u 2 d 3 u (5.9)
where we assign u = r ′r j with r ′ being the one in Equation 5.3. Equation 5.9 can also be expressed in terms of the electric potential given in Equation 5.3,

U Qq = Nq k=1 n q k V Q (r k + n).
(5.10)

In the same way, for single electrode with 3D-PBC, U QQ can be expressed as

U QQ = 1 2 N Q j ′ =1 N Q j=1 n R 3 R 3 ( σ ′ Q σ Q π ) 3 Q j ′ Q j |r j ′ -r j + u ′ -u + n| e -σ 2 Q u ′2 e -σ 2 Q u 2 d 3 u ′ d 3 u
(5.11) where we assign u ′ = r ′′r j ′ , or

U QQ = 1 2 N Q j ′ =1 n R 3 ( σ ′ 2 Q π ) 3 2 Q j ′ V Q (r j ′ + r ′′ + n)e -σ ′2 Q u ′2 d 3 u ′ (5.12)
where V Q (r j ′ + r ′′ + n) here is the electric potential at the position r j ′ + r ′′ + n generated by all the periodic images of all the electrode atoms except for the one located at r j ′ + n. Ewald summation method is also extended to calculate these two contributions with 3D-PBC [71]. Note that the interactions between electrolyte atoms and electrode atoms do not have a self interaction correction term.

Electrostatic interactions in two dimensional PBC

An electrochemical cell is often modelled by two electrodes and an electrolyte between them. For example, Figure 5.2 is a snapshot of a MD simulation of an electric double layer capacitor composed of 322 1-Ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide (EMIM-TFSI) ions pairs and 2 graphite electrodes. The MD simulation is conducted by our group to study the electric double layer of the system under different applied electric potentials. With 2 dimensional periodic boundary condition (2D-PBC) in x and y directions, the two electrodes extend infinitely in x and y directions but remain only two slabs in z direction [71,[START_REF] Gingrich | On the ewald summation of gaussian charges for the simulation of metallic surfaces[END_REF]. This picture resembles more an electrochemical cell than the picture with 3D-PBC where the electrodes are duplicated infinitely in z direction. Especially, the electric potentials at the two ends of the simulation box in z direction don't have to be the same. The potential difference between the two ends can be seen as the cell potential. With 2D-PBC, the equations of Ewald summation method are slightly modified and one additional term appears to account for the long-range interaction energy with reciprocal lattice vector being equal to zero. Complete set of equations and details of implementation in MetalWalls, a MD software developed by our group, can be found in reference [START_REF] Coretti | MetalWalls: simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes[END_REF].

Constant potential and electroneutrality constraints

So now, we have a collection of equations that gives us the electrostatic interaction energy of an electrochemical cell. To simulate realistic electrodes, we should impose two additional constraints. One requires the global electroneutrality and the other fixes the electric potential applied on the electrodes at given values. On using the constant potential electrode model, the integral charges of electrode atoms, Q j , are allowed to fluctuate in response to the change of electric potential generated by the electrolyte [START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF]. The constant potential constraint requires that the electric potential in the regions occupied by the positive or negative electrode, Ω ± , should be uniform and fixed at a given value Ψ ± for the corresponding electrode. This condition can be written in the following form V elec (r) = Ψ ± (5.13) where r here is a point in the regions occupied by the positive or negative electrode, Ω ± . The constant electric potential constraint is imposed via 5.14) depending on which electrode, positive or negative, the atom j is located in.

∂U elec ∂Q j = Ψ ± (
In the framework of MD simulation, the electrode charges are calculated in the Born-Oppenheimer approximation [START_REF] Scalfi | Charge fluctuations from molecular simulations in the constant-potential ensemble[END_REF][START_REF] Coretti | MetalWalls: simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes[END_REF]. At each time step of a simulation, the electrode charges are determined through a minimization procedure of the total energy of the system including U elec with respect to the electrode charges. An optimal set of electrode charges is found for each instantaneous configuration of electrolyte. The conjugate gradient method is often used for the minimization procedure.

Constant potential electrode model in MDFT

We have discussed in Chapter 4 that the computational cost of MD simulations is much greater than MDFT, especially for diluted systems. The high computational cost inhibits the easy access through MD simulations to the properties of electron transfer reactions at electrochemical interfaces. There are only few such studies in the literature [57,[START_REF] Reed | Electrochemical charge transfer at a metallic electrode: a simulation study[END_REF][START_REF] Remsing | Frustrated solvation structures can enhance electron transfer rates[END_REF][118]. Our perspective is to make use of the computational efficiency of MDFT to accelerate the studies of this kind. MDFT has been used to study the influence of a model graphite surface without electrostatic interaction to the electron transfer half reactions of chlorine, Cl → Cl + and Cl -→ Cl, in aqueous phase [60]. However, the electrostatic interactions is critical to the electrochemical interfaces [START_REF] Jeanmairet | Microscopic simulations of electrochemical double-layer capacitors[END_REF]113]. It can be taken into account by incorporating the same constant potential electrode model mentioned above into MDFT.

Approaches to the ensemble average

Recall the fundamental difference between MD and MDFT: one approaches ensemble averaged value with MD simulations by accumulating statistics while with MDFT, one directly obtains the ensemble averaged properties by the minimization procedure. In MD, we look at the atomic charge fluctuation at each time step of a trajectory due to the fluctuation of the positions of electrolyte atoms while, in MDFT, we seek for the ensemble averaged electrode charges and the equilibrium solvent density in the presence of polarized electrodes and solute molecules. Therefore, a different strategy should be followed to find the ensemble averaged properties of an electrochemical cell. The strategy is demonstrated and verified for a watergraphene capacitor in reference [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF]. It is summarized by the flowchart shown in Figure 5.3. In the framework of MDFT, electrode is also considered as a source of the external potential, ϕ(r, Ω). Each electrode atom bears a Gaussian charge distribution just as in MD simulations and a Lennard-Jones interaction site. We iteratively optimize the electrode charges and minimize the free energy functional until attaining a global convergence criterion. We search the optimal electrode charges At the very beginning, we calculate the electrode charges without solvent molecules. If the system includes other solute molecules such as a lithium ion, the optimization of the electrode charges is performed with the presence of these solute molecules. Next, we calculate both the electric potentials generated by the electrode atoms, V Q (r), and by the solute molecules, V q (r), using Ewald summation with 2D-PBC. Recall that in MDFT, space is divided by a grid and therefore, we actually calculate the electric potential on each node of the grid. The electrostatic part of the external potential, ϕ elec (r, Ω), is then calculated using Equation 3.95. Note that our implementation calculates ϕ elec (r, Ω) here in real space instead of using Fast Fourier transform. The external potential in the free energy functional is then the sum of the electrostatic part and Lennard-Jones part. Once the functional minimization is done, we can calculate an inhomogeneous charge distribution generated by the inhomogeneous solvent density ρ eq (r, Ω) via the following equation,

ρ s (r) = R 3 8π 2 ρ eq (r ′ , Ω) ρ s (r -r ′ , Ω)dr ′ dΩ (5.15)
where ρ s (r, Ω) is defined in Equation 3.94. Note that ρ s (r) is also calculated here in real space instead of using Fast Fourier transform. Now we can calculate the electric potential generated by ρ s (r) using Ewald summation with 2D-PBC on regarding each grid point as a point charge. Together with the electric potential generated by other solute molecules, it is taken as an input to optimize again the electrode charges. We continue to iterate between electrode charges optimization and functional minimization. As long as the global convergence criterion is met, the resulting electrode charges and equilibrium solvent density are the ensemble averaged properties.

Practical aspects

It should be mentioned that, we have introduced in Section 3.4 the method adapted to MDFT to calculate the electrostatic interaction energy with 3D-PBC.

Although it provides a convenient numerical route to evaluate the electrostatic interactions, it is not constructed for 2D-PBC. As a result, the above calculations of electric potentials are done using Ewald summation method. Fortran subroutines are implemented during the thesis project in the MetalWalls software package to calculate the electric potential generated by the electrode atoms and by the solute molecules on the grid, and the electric potential generated by ρ s (r). A Python interface that facilitates the calls of Fortran subroutines implemented in MetalWalls and home-made MDFT code, simplifies the data manipulation, and improves the numerical efficiency is also built.

We should also pay attention to the fact that the Lennard-Jones interaction energy is calculated with 3D-PBC in MDFT and the previous calculations that give the two body direct correlation function c (2) (r -r ′ , Ω, Ω ′ ) are also performed with 3D-PBC. Fortunately, these two parameters both have a short range characteristics. A simple way to avoid artifacts due to the discrepancy between 2D-PBC and 3D-PBC is to add thick enough vacuum regions outside the two electrodes in z direction. By doing so, the effects on the Lennard-Jones interaction and the excess free energy functional due to the nearby replicas of simulation box in the z direction are reduced to be negligible. This is accomplished by imposing an infinite external potential within the regions.

Part II

Developments and Applications

Chapter 6

Assessment of pressure corrections in the study of electron transfer The content of this chapter is an adapted version of "Assessing the correctness of pressure correction to solvation theories in the study of electron transfer reactions" published in The Journal of Chemical Physics [69] and has been reproduced here with the permission of the copyright holder.

Following the derivation of the state-of-the-art molecular density functional theory (MDFT) in Chapter 3, we arrived at a computer tractable free energy functional by using homogeneous reference fluid approximation in section 3.2. Then in section 3.3, we discussed in detail the negative consequences of employing the homogeneous reference fluid approximation and the two strategies of correction reported in the literature. In section 4.7, the theoretical background for applying MDFT to study electron transfer was introduced.

One of the works of the thesis is to make use of the computational efficiency of MDFT to study electron transfer half reactions of numerous quinones in acetonitrile (MeCN) solvent. For this purpose, it is necessary to explore out-of-equilibrium solvation which is not considered in the developments of correction methods. In this chapter, we evaluate the impact of the two types of corrections on the free energy profiles and reorganization free energies predicted for the two model electron transfer half reactions in water: Cl → Cl + and Cl -→ Cl.

SELECTED CORRECTION METHODS

Selected correction methods

For the a posteriori correction, we selected the widely used partial molar volume (PMV) correction proposed by Sergiievskyi et al. [64]. It consists in evaluating the partial molar volume of the solute and subtracting an amount of free energy equal to the partial molar volume multiplied by the compressibility route pressure from the calculated solvation free energy. The correction term to the solvation free energy takes the following expression:

P c ∆V = n b k B T (1 - n b 2 ĉ(k = 0)) × R 3 8π 2 ρ eq (r, Ω) -ρ b ρ b dΩdr (6.1)
where the definition of the compressibility route pressure P c can be found in Equation 3.89. Because Equation 6.1 does not depend on ρ(r, Ω), the correction term is not included in the free energy functional. Such a correction does not influence the optimization process but modifies the solvation free energy estimated in the end.

The equilibrium density remains the same as the one obtained using the hyper netted chain (HNC) functional.

For the bridge functional correction, we take the parameter-free angular independent weighted density approximation (WDA) bridge functional recently proposed by Borgis et al.,

F B [n(r)] = k B T n b R 3 a n(r) -n b n b 3 dr (6.2)
where n(r) is the weighted density using a Gaussian weighting function,

n(r) = n(r ′ )w(|r -r ′ |)dr ′ , (6.3) w(r) = (2πσ 2 w ) -3/2 exp (-r 2 /2σ 2 w ). (6.4)
with the Gaussian width σ w = 1 Å. In Equation 6.2, a is a parameter uniquely defined by imposing zero compressibility route pressure. On incorporating the WDA bridge functional, Equation 3.88 becomes

F[ρ(r, Ω) = 0] = n b k B T V - 1 2 n 2 b k B T V ĉ(k = 0) -n b k B T V a. (6.5)
Note that on integrating over the whole 3-dimensional space, we get a factor of V . Thus, a is unambiguously determined by the following equation,

a = 1 - 1 2 n b ĉ(k = 0). (6.6)

Computational details

We study two model electron transfer half reactions Cl → Cl + and Cl -→ Cl in water. The solutes are modeled by one Lennard-Jones site with σ LJ = 4.404 Å STUDY OF ELECTRON TRANSFER and ϵ LJ = 0.4190 kJ/mol and atomic charges q corresponding to the three oxidation states Cl, Cl + , and Cl -. We use a Lennard-Jones cutoff of 10 Å with long range corrections. We generate a series of biased external potentials using Equation 4.50 with ϕ P (r, Ω) being the external potential generated by Cl + and ϕ R (r, Ω) being the external potential generated by Cl. We choose the interpolation interval of -1.2 ≤ η ≤ 2 and a step size of 1/30. It should be mentioned that because the Lennard-Jones parameters are the same for all solutes, the series of biased external potentials is equivalent to the series of external potentials generated by the solutes having atomic charge between -1.2 ≤ q ≤ 2. As a result, at η = -1, ϕ η (r, Ω) is exactly the external potential generated by Cl -.

We use a 40 × 40 × 40 Å 3 box with 120 3 spatial grid points and 196 possible orientations per spatial point with 3 dimensional periodic boundary condition (3D-PBC). We use the type C correction proposed by Kastenholz and Hünenberger due to the interaction between the solute and its periodic replica [119,[START_REF] Kastenholz | Computation of methodology-independent ionic solvation free energies from molecular simulations. II. The hydration free energy of the sodium cation[END_REF]. The temperature is fixed at 298.15 K.

Results and discussion

Assessment of partial molar volume correction

The free energy profiles of the two half reactions calculated using the HNC functional with and without the PMV correction are displayed in figures 6.1A and 6.1B. All minima are shifted to zero to ease the comparison with the curves computed with Hartnig and Koper using MD [121]. The agreement between the free energy profiles obtained using MD and MDFT is satisfactory for the neutral solute. However, important deviations are observed in the case of the cation and the anion. This is a known defect of the HNC functional, which is not able to properly render the hydrogen bond network of water that is essential for the solvation of ions. This is due to the truncation of the expansion of the two body direct correlation function in Equation 3.78 at the first term. The two types of corrections studied here cannot resolve this issue because of their lack of angular dependency.

At first glance, the corrections including the PMV correction do not seem to deeply modify the shape of the curve. However, a closer look reveals that the positions of the minima are shifted. This is especially true for the ions where the minima are shifted toward zero. This might appear surprising since the average vertical energy gap defined in Equation 4.50 solely depends on the density ρ η (r, Ω), which is not modified by the PMV correction.

In fact, the shift of the minima is a consequence of the a posteriori nature of PMV correction. This is visible in Figure 6.2 where the free energy profiles of the ions are depicted as a function of the atomic charge q of the fictitious solutes that generate the series of biased external potentials. The minima of the PMV corrected free energy profiles of Cl + and Cl -no longer correspond to q = 1 and q = -1, respectively. Equivalently, it means that the minimum of the free energy profiles no longer correspond to the equilibrium solvation state of Cl + and Cl -. This is a violation of the variational principle, which is not surprising since the PMV correction is not properly integrated in the classical density functional theory (cDFT) formalism. It is simply an a posteriori correction that does not influence the functional minimization. This might be acceptable when the objective is to reproduce some reference solvation free energies [61,64,66,[START_REF] Misin | Communication: accurate hydration free energies at a wide range of temperatures from 3D-RISM[END_REF][START_REF] Sergiievskyi | Fast computation of solvation free energies with molecular density functional theory: thermodynamic-ensemble partial molar volume corrections[END_REF][START_REF] Misin | Predicting solvation free energies using parameterfree solvent models[END_REF], but it is problematic to study electron transfer. Indeed, the two ways to compute the reorganization free energies, either from the mathematical definition in equations 4.59 and 4.61 or from the graphical definition in equations 4.58 and 4.60 do not coincide anymore.

The reorganization free energies are given in Table 6.1. Note that there are two sets of values for the neutral chlorine since the value of reorganization free energies depends on the choice of the second oxidation state. Reorganization free energies of Cl and Cl + computed with the HNC functional are similar, indicating that this electron transfer half reaction is well described by Marcus theory with linear response approximation, while the Cl -→ Cl reaction deviates from Marcus theory. These results are consistent with the results of MD simulations [121]. When the PMV correction is added, both electron transfer half reactions seem to deviate from Marcus theory. The reorganization free energies of ions differ from the ones of the associated neutral state. Moreover, the reorganization free energies computed using equations 4.59 and 4.61 differ by up to 25 kJ/mol from the ones computed with the graphical definition in equations 4.58 and 4.60.

The evolution of the average vertical energy gap as a function of q is displayed in Figure 6.3. We recall that values of q between 0 and 1 correspond to the Cl → Cl + , while the value between -1 and 0 corresponds to Cl -→ Cl. When an electron transfer is following the linear response approximation, the energy gap should vary linearly. This is the case between q = 0 and q = 1, while the linearity is not and WDA functional (red solid line) as a function of the atomic charge q. The dotted and dashed lines are the tangent to the HNC curve at q = -1 and q = 1, respectively. STUDY OF ELECTRON TRANSFER respected between q = 0 and q = -1. This is consistent with the conclusions drawn from examining the reorganization free energies computed using the HNC functional and with the MD results: the Cl → Cl + half reaction follows the prediction based on linear response approximation while the Cl -→ Cl half reaction does not. Since the energy gap is not modified by the PMV correction, the linear behavior observed for the Cl → Cl + reaction is in contradiction with the different values of λ of Table 6.1. This is another proof of the inappropriateness of the PMV correction to study electron transfer.

Assessment of bridge functional correction

The same electron transfer half reactions were studied using the bridge functional correction with the WDA bridge functional of Equation 6.2. The prediction for the free energy profiles as a function of the average vertical energy gap or the solute charge are given in figures 6.1 and 6.2, respectively. Adding the WDA bridge functional has a limited impact on the free energy profile as compared to the HNC predictions. Nonetheless, more importantly, the minimum of the free energy landscape of the cation (anion) is correctly reached for the solute with charge q = 1 (q = -1). This is because the WDA bridge functional is incorporated in the functional minimization, and the variational principle is respected. The reorganization free energies predicted with the WDA bridge functional are slightly reduced with respect to the ones predicted using HNC functional, but the overall conclusion is recovered: Cl → Cl + follows Marcus theory with linear response approximation while Cl -→ Cl does not. This is also supported by the evolution of the energy gap as a function of the atomic charge in Figure 6.3.

We now attempt to understand why the reorganization free energies of the ions is reduced when the WDA bridge is used. We first show that only the electrostatic contribution is relevant in the comparison of two functionals by breaking down the free energy into individual contributions. Then the reduction of the reorganization free energies of the ions upon adding WDA bridge is discussed through the change of solvent polarization.

Free energy breakdown

We introduce the difference of out-of-equilibrium solvation free energies calculated for a given value of η with the two different functionals for an oxidation state M as a function of η,

∆F M (η) = F WDA M [ρ ′ η ] -F PMV M [ρ η ], (6.7)
where ρ ′ η ̸ = ρ η indicates that the incorporation of the WDA bridge functional modifies the minimization process and the equilibrium solvent density fields. In Equation 6.7, M may indicate any one of the three oxidation states, Cl, Cl + , or Cl -. We can separate the free energy into two contributions, where ∆F ext M (η) is the difference of the external free energy functionals defined in Equation 3.71 and ∆F int (η) groups the difference of the ideal free energy functionals and the excess free energy functionals defined in Equation 3.70 and 3.72, respectively. Note that for a given η, the intrinsic part of out-of-equilibrium solvation free energy F int [ρ η ] is independent of the oxidation state M and so is its difference ∆F int (η).

∆F M (η) = ∆F ext M (η) + ∆F int (η), (6.8)
Using equations 4.59, 4.61, and 6.8, the difference of the reorganization free energies obtained with and without the WDA functional for the oxidation state M can be expressed as

∆λ M = λ WDA M -λ HNC M = ∆F M (0) -∆F M (η) = ∆F ext M (0) -∆F ext M (η) + ∆F int (0) -∆F int (η), (6.9)
where η could be 1 or -1 indicating oxidation states Cl + or Cl -. On the other hand, for the oxidation state Cl,

∆λ 0 = λ WDA 0 -λ HNC 0 = ∆F 0 (η) -∆F 0 (0) = ∆F ext 0 (η) -∆F ext 0 (0) + ∆F int (η) -∆F int (0) (6.10)
where η could be 1 or -1 indicating that ∆λ 0 is calculated on considering half reactions with Cl + or Cl -, respectively. The difference in external, intrinsic and total contributions to the free energy are plotted as functions of the atomic charge q (equivalent to η) used to generate the out-of equilibrium solvent densities for the three oxidation states in Figure 6.4. Recalling that the external potential is the sum of Lennard-Jones and electrostatic contributions, ∆λ M in Equation 6.9 can be STUDY OF ELECTRON TRANSFER further split into .11) where ∆F elec M (η) and ∆F LJ (η) are the electrostatic and the Lennard-Jones part of external functional difference, respectively. Equation 6.10 can be split in the same way, too. Note that ∆F LJ (η) is also independent of the oxidation state due to the fact that Cl, Cl + , and Cl -have the same Lennard-Jones parameters.

∆λ M = ∆F elec M (0) -∆F elec M (η) + ∆F LJ (0) -∆F LJ (η) + ∆F int (0) -∆F int (η). ( 6
Since there is no electrostatics involved in the case of the neutral chlorine, its external contribution ∆F ext Cl (η) is equal to the Lennard-Jones contribution ∆F LJ (η) of Equation 6.11. The variation of ∆F ext Cl (η) remains below 2.1 kJ/mol across the whole range of external potentials and therefore we can assume,

∆F LJ (0) -∆F LJ (η) ≈ 0. (6.
12)

The variation of ∆F int (η) has a more pronounced amplitude than ∆F LJ (η) over the range of external potentials, but it remains below 6 kJ/mol which is considerably lower than the variation of the external (and thus electrostatic) contribution for ions. We can also assume

∆F int (0) -∆F int (η) ≈ 0. (6.13)
The above two assumptions lead to

∆λ M ≈ ∆F elec M (0) -∆F elec M (η). (6.14)

Changes of dipolar polarization of solvent

In Equation 6.14, we assume that since the considered species are ions, the main contribution is due to the electrostatic, i.e., ∆F int (η) and ∆F LJ (η) can be neglected. As long as we are considering spherical ions, the electrostatic potential is spherically symmetric, so is the polarization of the solvent. The equilibrium dipolar polarization density of the fluid can be computed from the equilibrium solvent density as P (r, Ω) µ 0 = ρ eq (r, Ω)ΩdΩ. (6.15) where µ 0 is the molecular dipole of the solvent molecule. The spherically averaged radial component of the solvent equilibrium polarization around each ion is displayed in Figure 6.5.

The WDA functional reduces the polarization of the solvent in the vicinity of the charged solutes, as evidenced by the decrease in the first peaks in Figure 6.5. Since the polarization is reduced, the ion is destabilized by the WDA functional when immersed in its equilibrium density. From Figure 6.4, we can find that both ∆F ext Cl -(-1) and ∆F ext Cl + (1) have positive value. On the contrary, the equilibrium density of the oppositely charged ion is less destabilizing. From Figure 6.4, we can find that both ∆F ext Cl - [START_REF] Tarascon | Issues and challenges facing rechargeable lithium batteries[END_REF] and ∆F ext Cl + (-1) have negative value. Finally, the WDA correction being angular independent, we can expect the polarization around the neutral solute to be almost unchanged, ∆F elec Cl (0) ≈ 0 and the results in Figure 6.4 justifies our expectation. With this simple analysis and using Equation 6.14, we recover the decrease in the reorganization free energies observed in Table 6.1.

Chapter conclusions

In this chapter, we studied the appropriateness of two corrections to the HNC functional to study electron transfer in aqueous medium with MDFT. First, we examine a simple a posteriori correction of the pressure, which has been widely used with different expressions in MDFT [61,64,66,[START_REF] Misin | Communication: accurate hydration free energies at a wide range of temperatures from 3D-RISM[END_REF][START_REF] Sergiievskyi | Fast computation of solvation free energies with molecular density functional theory: thermodynamic-ensemble partial molar volume corrections[END_REF]124,125]. Despite its success to predict solvation free energy in good agreement with reference simulations and experiments, this correction should not be used to study electron transfer. The minimum of the free energy profiles do not correspond to the equilibrium solvent configuration when this correction is used because it modifies the free energy without affecting the functional optimization. The reorganization free energies become illdefined: the graphical definition does not coincide with the mathematical expression. Moreover, we have shown that using any of the two definitions of the reorganization free energy leads to a behavior deviating from the Marcus picture for the Cl → Cl + half reaction in disagreement with MD simulations and HNC MDFT calculations. This deviation from Marcus theory is not recovered when examining the evolution of the average vertical energy gap, which is a final evidence for the inconsistency of the pressure correction.

We then turned to another type of correction, a so-called bridge functional correction, which tries to recover some of the contributions due to the terms beyond the quadratic approximation of the HNC functional. We have chosen to use the simple angular independent WDA bridge functional proposed by Borgis et al. that was shown to properly reproduce the solvation free energies of hydrophobic solutes [68]. Since this bridge has a functional form, the equilibrium density is modified and this approach does not suffer from the flaws of the PMV correction. There are no ambiguities in the definitions of the reorganization free energy, and we recover results consistent with MD simulations and HNC MDFT calculations: Cl → Cl + follows Chapter 7

Electron transfer of quinones In Chapter 1, we have mentioned that quinones are organic molecules that can undergo reversible redox reactions. They have been proposed as an alternative choice to metal-based redox active materials in redox flow batteries [32][33][34][START_REF] Son | Quinone and its derivatives for energy harvesting and storage materials[END_REF][START_REF] Ding | Molecular engineering of organic electroactive materials for redox flow batteries[END_REF][START_REF] Leung | A new aqueous all-organic flow battery with high cell voltage in acidic electrolytes[END_REF][START_REF] Han | Organic quinones towards advanced electrochemical energy storage: recent advances and challenges[END_REF]. A main advantage of such organic molecules is the versatility of their molecular structure. By grafting functional groups onto quinones, it is possible to fine tune not only their chemical stability for reversible energy storage applications [START_REF] Pineda Flores | Bio-inspired electroactive organic molecules for aqueous redox flow batteries. 1. thiophenoquinones[END_REF][START_REF] Wedege | Organic redox species in aqueous flow batteries: redox potentials, chemical stability and solubility[END_REF][START_REF] Tabor | Mapping the frontiers of quinone stability in aqueous media: implications for organic aqueous redox flow batteries[END_REF][17] but also the redox potential [START_REF] Yu | Tuning the performance of aqueous organic redox flow batteries via first-principles calculations[END_REF][START_REF] Huynh | Quinone 1 eand 2 e-/2 h+ reduction potentials: identification and analysis of deviations from systematic scaling relationships[END_REF][START_REF] Er | Computational design of molecules for an all-quinone redox flow battery[END_REF][START_REF] Bachman | Investigation of the redox chemistry of anthraquinone derivatives using density functional theory[END_REF][START_REF] Schwan | Substituent pattern effects on the redox potentials of quinone-based active materials for aqueous redox flow batteries[END_REF], the solubility [START_REF] Er | Computational design of molecules for an all-quinone redox flow battery[END_REF][START_REF] Pineda Flores | Bio-inspired electroactive organic molecules for aqueous redox flow batteries. 1. thiophenoquinones[END_REF], and the kinetics [18] of the associated redox reaction. However, reports discussing reaction kinetics of functionalized quinones, that is related to the power density of a device, remain very scarce [19-21, 27, 29, 126].

Following the derivation of Marcus theory in Chapter 4, we have seen that the reaction rate constant in Equation 4.2 is not solely determined by the redox species but is also impacted by solvent effects. Molecular dynamics (MD) has long been the simulation method of choice to verify the legitimacy of Marcus' linear response approximation and to investigate the influence of molecular effects on electron transfer reactions that may lead to deviations from the approximation. A few reports have applied these methods to study quinones [START_REF] Reeves | A first-principles investigation of the structural and electrochemical properties of biredox ionic species in acetonitrile[END_REF]127,128]. VandeVondele et al. generated trajectories of 1,4-benzoquinone (BQ) and duroquinone (DQ) solvated in acetonitrile (MeCN) using classical molecular dynamics (classical MD). The molecular structure of DQ is displayed in Figure 7.1. They then computed the energy cost to reduce the molecule using electronic density functional theory (eDFT) that enters Equation 4.18 before evaluating the exponential part of Equation 4. 2 [127]. Reeves et al. simulated anthraquinone (AQ) and ionic liquid (IL) grafted with AQ solvated in MeCN using ab-initio molecular dynamics (AIMD) [START_REF] Reeves | A first-principles investigation of the structural and electrochemical properties of biredox ionic species in acetonitrile[END_REF]. The molecu- However, because the associated solvent configurations have a low probability, it is a tedious task to properly sample the region close to the transition state with MD simulations which makes the construction of free energy profiles computationally demanding. In Section 4.7, we have introduced an alternative approach based on molecular density functional theory (MDFT) which is based on functional minimization rather than on statistical sampling, and thus considerably reduces the computational cost as compared to MD [60,69]. The use of MDFT enables us to perform a systematic study of quinones functionalized with different functional groups.

In this chapter, we show how the functionalization of BQ impacts the reorganization free energies of electron transfer half-reactions in MeCN using the approach based on MDFT. We start by validating the application of the method on quinones by comparing the predictions of the solvation shell structure and the free energy profiles to reference classical MD simulations using a recently developed polarizable force field [START_REF] Berthin | Solvation of anthraquinone and TEMPO redox-active species in acetonitrile using a polarizable force field[END_REF] in the case of AQ and its reduced anion anthraquinone anion (AQ -) solvated in MeCN. This serves as a benchmark for MDFT calculations using generic force fields. The results show that MDFT provides a solvation structure in good agreement with classical MD simulations as well as a good estimate of the reorganization free energies. The free energy profiles of one electron transfer half-reactions of BQ and DQ and their reduced anions, 1,4-benzoquinone anion (BQ -) and duroquinone anion (DQ -), solvated in MeCN are then computed using MDFT. Finally, we carry a systematic calculation of reorganization free energies for several functionalized BQ derivatives modelled with a generic force field using MDFT.

Computational details

The results of AIMD simulations were taken from the previous work of our group [START_REF] Reeves | A first-principles investigation of the structural and electrochemical properties of biredox ionic species in acetonitrile[END_REF], where the technical details can be found. Reference classical MD simulations were performed by my colleague Roxanne Berthin using the MetalWalls code [START_REF] Coretti | MetalWalls: simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes[END_REF][START_REF] Marin-Laflèche | Met-alWalls: a classical molecular dynamics software dedicated to the simulation of electrochemical systems[END_REF] and the polarizable force field recently developed by her [START_REF] Berthin | Solvation of anthraquinone and TEMPO redox-active species in acetonitrile using a polarizable force field[END_REF]. The simulation cells contain one redox molecule (either AQ or AQ -) and 95 MeCN molecules. Cubic cells were built using the PACKMOL [130] package. A first run was performed in the constant potential constant temperature ensemble at 298 K and 1 bar to reach the experimental density of MeCN. Then the simulations were carried out in the canonical ensemble at 298 K with cubic cells of 20.553 Å 3 for a production time of 50 ns after 1 ns of equilibration. A timestep of 1 fs was used and the trajectory was saved every 500 fs. The equations of motion were integrated using the velocity-Verlet algorithm [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF], and the temperature was kept constant using a chain 107]. The short-range van der Waals interactions were calculated with a cutoff value of 10 Å, while the Ewald summation method was used for electrostatics interactions (involving the partial charges and the induced dipoles).

MDFT calculations are performed with an in-house developed Fortran code. We use a 30 3 Å 3 cubic simulation box with 125 grid nodes in each direction and 66 orientations per grid point. To compute the density and polarization maps of AQ, a cubic cell of 25 3 Å 3 and 125 3 grid nodes was used to have the same resolution than the maps computed with MD. The non-polarizable model for MeCN [START_REF] Edwards | A computer simulation study of the dielectric properties of a model of methyl cyanide: I. the rigid dipole case[END_REF] is given in Section 3.4. In this study, we don't apply pressure correction to MeCN solvent. Preliminary molecular structure of solutes were first generated via OpenBabel [131] from the canonical SMILES notation of the molecules. In a second step, we applied the semi-empirical quantum model, Austin model 1, together with the bond charge correction (AM1-BCC) [132,133] to optimize the molecular geometry and to find the atomic charges. Lennard-Jones parameters were taken from the second generation generalized amber force field (gaff2) [134] with a 10 Å cutoff radius. Both AM1-BCC modelling tool and gaff2 force field database distributed within the AmberTool software package [135] were used. The force field parameters of quinones anions were chosen applying the strategy mentioned in reference [START_REF] Berthin | Solvation of anthraquinone and TEMPO redox-active species in acetonitrile using a polarizable force field[END_REF]: The structure and the force field parameters are kept identical to the oxidized form but an excess charge of -0.8 e is evenly split between all non-hydrogen atoms. Note that this excess charge is reduced to -0.8 e to implicitly account for the electronic polarizability as is commonly done in non-polarizable force fields [136,137]. The non-dimensional dipolar polarization density is computed using Equation 6.15 We first look at the solvation structure of AQ solvated in MeCN. The solvent density and the polarization in the plane of the aromatic rings computed with MDFT and classical MD are compared in Figure 7.3 together with density results obtained using AIMD trajectories. In the case of MD simulations, these quantities were extracted from the trajectories using Travis [138, 139] using a sampling volume of 0.2 × 0.2 × 4.0 Å 3 , while in MDFT the density is a direct output from the calculation (averaging in a slab of width 4.0 Å in the z direction was made for comparison purposes).

While the density maps obtained by MD and MDFT displayed in Figure 7.3 are very similar, it is quite obvious that the AIMD one is undersampled. The AIMD simulation was run for 18 ps with the same box as classical MD, which is a standard simulation time for this method. This illustrates that the high computational cost associated with AIMD prevents it to access the 3D structure of the solvent around the solute, this is all the more true for polarization properties which also require a sampling of orientations.

A closer look into the MD and MDFT density maps reveals the existence of two solvation shells, the second one showing few variation in intensity. The first solvation shell is more structured with the presence of ten basins of high density. They are all located in the center of the regions delimited by two vicinal C-H or C=O bonds. The agreement between MDFT and MD is excellent considering the force fields used for the acetonitrile solvent and the AQ species differ between the two methods.

We now analyze the polarization density provided by the classical MD and MDFT. In Figure 7.3 it is represented using vectors, which directions show the projection of the polarization in the plane of the AQ aromatic rings and which lengths are proportional to the intensity of the polarization. To improve the readability of the figure, orientations are depicted on a grid 2.5 times looser than the one used for calculation. Again, a very good agreement is found between classical MD and MDFT for the prediction of the polarization maps. In the first solvation shell, MeCN dipoles point towards the oxygen of the carbonyl groups and away from the solute around the C-H bonds. This is an expected behaviour: the negatively charged oxygen atoms attract the positively charged CH 3 fragment of the solvent molecule while the positively charged hydrogen atoms attract the negatively charged nitrile fragment. Between C-H bonds and the carbonyl groups, MeCN dipoles lie almost parallel to the cavity surface created by AQ as shown in Figure 7.4A. This cavity surface is defined as the isosurface on which the solvent density is at 10 -2 × n 0 where n 0 is the bulk density. The value of the polarization is considerably reduced in the second solvation shell because the electric field generated by the solute is screened by the first solvation shell. The vectors are antialigned with the polarization observed in the first solvation shell. This comparison proves that MDFT is able to predict faithfully the structural properties of solvation around AQ. This is of great interest because MDFT is far less computationally demanding: it took more than 100 hours on a 44 core skylake computer to run the 50 ns MD trajectory that was used to compute the density and polarization maps displayed in Figure 7.3B while the MDFT calculation to compute the maps displayed in Figure 7.3C took less than 6 minutes on the same computer.

We now turn to the energetic properties for the electron transfer half-reaction. We first examine the redox half reaction between AQ and AQ -solvated in MeCN. Free energy profiles of the two species are computed with classical MD using the procedure described in Section 4.5. The same reaction is studied with MDFT where we generated a series of potential ϕ η (r, Ω) by varying η in Equation 4.51 between -0.3 and 1.3 by increment of 1/30. Using values of η > 1 and η < 0 allows to compute the branches of the free energy profiles in regions corresponding to the left hand side of the negative minimum (reduced state) and right hand side of the positive minimum (oxidized state).

As mentioned above, the estimation of the free energy profiles from MD allows to calculate the free energy difference between the two states -the same argument also applies to MDFT. However, this free energy difference cannot readily be compared to experimentally measured potentials unless it can be calibrated with another chemical reaction [140]. This was made in the case of AIMD by setting up a computational standard hydrogen electrode [30] and more recently a computational Ag/AgCl reference electrode [141]. In the present work, the use of classical methods prevents us from building such absolute scales, so that we use the AQ/AQ - as a reference. In practice, this is made by assigning an arbitrary energy constant to the free electron in Equation 1.1 so that the minima of the free energy profiles align for this couple -the same shift is then used for all the other couples. Yet it is not possible to comment these free energy differences with respect to experimental values since our calculations do not include the dominant intramolecular energy term, which is properly taken into account in quantum chemistry calculations [START_REF] Huynh | Quinone 1 eand 2 e-/2 h+ reduction potentials: identification and analysis of deviations from systematic scaling relationships[END_REF][START_REF] Er | Computational design of molecules for an all-quinone redox flow battery[END_REF]. In this work we therefore focus on the value of the reorganization free energy. Note that applying a shift to the vertical energy gap values has no influence on this quantity. Moreover, because we concomitantly shift the average vertical energy gap defined in Equation 4.50, the crossing point of the two curves is still at ⟨∆E⟩ η = 0 as it should be by the definition in Equation 4.50.

The free energy profiles obtained for the AQ/AQ -couple using MD and MDFT are compared in Figure 7.5A. Overall, the agreement is excellent both for the curvatures and the position of the minima. To check if this half reaction can be described with Marcus theory we computed the reorganization free energy of AQ and AQ - with Equation 4.61. MD gives a value of 72.3 kJ/mol for both states while MDFT predicts the reorganization free energies of 66.8 and 69.3 kJ/mol for AQ and AQ -, respectively. The values of the reorganization free energies differ by less than 5 % between the two states which implies that this system can be described with a very good approximation using Marcus Theory. Thus, we have fitted the MDFT data using two parabolas with the same curvature (and hence a single λ), as depicted in Figure 7.5.

To further test our methodology, we computed the free energy profiles of the electron transfer involved in the BQ and DQ redox couples solvated in MeCN since these systems were already studied using an ab initio-based approach by Vande-Vondele et al.. The free energy profiles are respectively plotted in Figure 7.5B and 7.5C. For these redox couples as well, they are well fitted by pairs of parabolas with identical curvature indicating that those 2 electron transfer reactions are well described by Marcus theory. The reorganization free energies obtained using MDFT are 86.8 kJ/mol for BQ/BQ -and 74.0 kJ/mol for DQ/DQ -. They are in good agreement with the ab initio-based simulations [127], from which values of 78.3 and 69.4 kJ/mol were reported for BQ/BQ -and DQ/DQ -, respectively. A qualitative agreement is obtained even if the MDFT results are overestimated by about 10 % with respect to the ab initio-based simulation results.

Electron transfer of benzoquinone derivatives

We now take advantage of the computational efficiency of MDFT to carry on an extensive study of the impact of functionalization on the electron transfer properties of quinones. More precisely, we determine the effect of chemical substitutions by studying the six derivatives of BQ displayed in Figure 7.6. These are obtained by systematically substituting the H atoms with a functional group. For simplicity, we only consider BQ derivatives substituted with a single type of functional group. The chemical functions considered in this work are methyl (-CH 3 ), ethyl (-CH 2 CH 3 ), methoxy (-OCH 3 ), amino (-NH 2 ), hydroxy (-OH), fluoro (-F), chloro (-Cl), thiol (-SH), and carboxyl (-COOH).

We first calculate the reorganization free energies of the oxidized and reduced states. As for the benchmark molecules, they differ by less than 5 % for the whole set of molecules so that all the redox couples can be considered as being well described by Marcus theory and attributed a single λ value. The corresponding reorganization free energies are displayed in Figure 7.7. Firstly, we notice that all the functionalized BQ have a reorganization free energy lower than that of the parent molecule. There is also a clear trend of the evolution of the reorganization free energy with the number of functional groups: the more the molecule is functionalized the lower is the reorganization free energy. Finally, there is little difference between the different di-functional quinones, but we can see that molecules where the functional groups In an attempt to rationalize this behaviour we estimated the volume of each functional group as the difference of cavity volume between BQ and the monofunctionalized BQ derivatives. Cavity volume is defined as the volume of the space where the solvent density is below 10 -2 × n 0 . Examples of such volumes for AQ and BQ are displayed in Figure 7.4. The variation of λ with respect to the volume of functional groups is displayed in the bottom panel of Figure 7.8. A strong correlation between the two quantities is observed for the whole range of functional groups, with the exception of the amino one. Recall that the activation free energy, and thus the rate constant of electron transfer reaction, is connected to the reorganization free energy. This could be a way to fine tune the kinetics of the electron transfer reaction. However, we emphasize that this approach cannot capture the full picture of the electron transfer reaction since internal reorganization of the solute, in particular due to electronic effects, is neglected. The latter may be important in the case of complex chemical entities and affect the kinetics aspect as well.

Since we have access to the reorganization free energies and the solvation structure for the full set of functionalized molecules, we finally test the relation 4.6 that links λ to the solute radius. This is a difficult test since this relation was derived for spherical solutes, which is clearly not the case here as evidenced by the molecular structure displayed in Figure 7.6 and by the solvent density isosurface around BQ and AQ displayed in Figure 7.4. To do so, we define the cavity radius r c as the radius of a sphere of identical volume to the cavity created by the solute. The reorganization free energy is plotted as a function of the inverse of the cavity radius 4.6, is 424 kJ Å mol -1 . The value of ϵ op is chosen to be consistent with our treatment of the electronic polarizability by scaling the the excess charge [136] according to q eff = q √ ϵop . We find an excellent agreement between our numerical data and the theoretical predictions.

If we take a closer look into the difference between functional groups, we can see that carboxylic, amino, methoxy and hydroxy groups lead to λ values slightly below the linear regression yet chloro and fluoro groups have a reverse effect. Functional group with larger dipole moment are thus lying below the linear regression. However, this effect remains weak compared to the general trend with respect to the molecular radius.

We could argue that the decrease of the reorganization free energy when the solute radius increases is simply a consequence of the reduced charge density of the solute. A good way to assess this assumption is to split the reorganization free energy into its different contributions using Equation 3.82. The ideal, external and excess part of the reorganization free energies of BQ and its monofunctionalized derivatives are displayed in Figure 7.9 as a function of the inverse of the cavity radius. All contributions exhibit variations much larger than the one of λ. Thus, the linear evolution of the reorganization free energy with the cavity radius cannot be explained by simple argument on the solute-solvent interaction but is rather due to non trivial compensation between this term, solvent-solvent contributions and entropic term.

Chapter conclusions

In this work we carry a systematic study of the influence of functionalization on the reorganization free energies of BQ derivatives. We first validate our methodology by carrying on a detailed comparison of the solvent structural properties around AQ predicted using MDFT and classical MD. The force fields used in the two methods differ since (i) molecules are rigid in MDFT while they are flexible in MD (ii) the solvent model is coarse grained in MDFT while it is all-atom in MD. Despite these differences the agreement between the two approaches is very satisfactory. The density peaks in the first solvation shell are well reproduced, and more importantly the polarization maps are highly similar indicating that the solvent orientation properties are also well captured by MDFT. The agreement between the two methods is not limited to the solvation structure since the free energy profiles computed for both redox states with MD and MDFT are almost identical. We found that the AQ/AQ -electron transfer half-reaction is in good agreement with Marcus' prediction since the free energy profiles of the two species are well fitted by parabolas of identical curvature. This establishes MDFT as an appropriate method to study the electron transfer reaction of quinone derivatives while it requires 10 3 less computational times than MD to study the same system. To be fair, using a non-polarizable force field in MD would reduce the computational cost by one order of magnitude, but the gap would remain very large. This computational efficiency allows to systematically study the influence of functionalization on the reorganization free energy of electron transfer half reaction involving BQ derivatives. Here again, we found that the whole set of molecules considered in this study follows the Marcus picture with identical values for the reorganization free energies of the oxidized and reduced states. The values of the reorganization free energy of each functionalized molecule is lower than the value for the non-functionalized parent BQ. The more the molecule is functionalized, the lower is its reorganization free energy. Finally, we found a correlation between the value of the reorganization free energy and the volume of the functional group. This lead us to check the linear scaling of the reorganization free energy with the inverse of the cavity radius. This relation is indeed verified but with a slope that differs from the theoretical values, the difference being attributed to the non-sphericity of the cavity created by the solute.

Overall, in this chapter, we show that MDFT is a suitable method to account for the functional effects on the reorganization free energies of electron transfer reaction. In the future, systematic studies could thus be performed in order to understand the rate constant of the redox reactions. Of course having a full picture of the reactions free energy profile would require to account for internal degrees of freedom using electronic structure methods. A promising approach to fully address this problem could be the quantum mechanics/MDFT hybrid approach that our group recently proposed [142]. In this approach the solute is dealt with at the quantum level using eDFT while the solvent degrees of freedom are accounted at the classical level using MDFT. We believe that this method could be the appropriate compromise between precision and computational cost to make the calculation of electron transfer reaction rate constant in solution feasible in the future.

Chapter 8

Lithium ion adsorption on graphite surface The content of this chapter is an adapted version of "Multi-scale simulation of the adsorption of lithium ion on graphite surface: From quantum Monte Carlo to molecular density functional theory" published in The Journal of Chemical Physics [START_REF] Ruggeri | Multi-scale simulation of the adsorption of lithium ion on graphite surface: From quantum monte carlo to molecular density functional theory[END_REF] and has been reproduced here with the permission of the copyright holder.

Carbon materials play a very important role in electrochemistry. In energy storage applications, carbon is used as an electrode in the form of graphite for Liion batteries [143], of hard carbon for Na-ion batteries [144], and of nanoporous materials for supercapacitors [145]. In all these examples, the carbon materials are put in contact with a charged electrolyte, and the interfacial structure and dynamics play a crucial role in the operation of the devices [146,147]. In order to optimize the performance, molecular simulations can play an important role [148], but they need to accurately account for the interactions between the ions and the carbon surface.

Over the years, a large number of electronic density functional theory (eDFT) and molecular dynamics (MD) simulations were devoted to the study of electrochemical interfaces between carbon and ions or liquid electrolytes. Although they provide qualitatively similar results, some discrepancies may be observed when analyzing quantitative properties, such as adsorption energies, preferential "binding" distances, or adsorption profiles. For example, as reported by Valencia et al. [149], in the case of bare lithium adsorption on graphite surfaces, eDFT calculations always display binding energies larger than 1 eV with a preference for the center of the C 6 hexagonal rings (hollow sites) with respect to adsorption above C atoms (top sites) or C-C bonds (bridge sites) [149]. The lithium atom loses a valence electron, which is entirely transferred to the carbon surface. Yet, the adsorption energy may differ by more than 0.5 eV when changing the exchange-correlation (XC) functional [150].

In particular, the inclusion of the dispersion effect further increases the binding energy of the lithium ion [START_REF] Fan | Adsorption of single li and the formation of small li clusters on graphene for the anode of lithium-ion batteries[END_REF].

Similar problems arise in classical molecular dynamics (classical MD) simulations due to the choice of different interaction potentials. On one hand, the choice of the electrostatic model for the electrode atoms, which can be treated using fixed-charges or the constant potential method [113], will impact the localization of the charge induced by the ions on the carbon. On the other hand, several force fields may be used for the short-range repulsion and dispersion interactions.

Important discrepancies were already reported in the case of the adsorption of water molecules on carbon surfaces [147]. In two studies on carbon nanotubes [START_REF] Al-Hamdani | How strongly do hydrogen and water molecules stick to carbon nanomaterials?[END_REF] and other carbon nanostructures [153], Michaelides and co-workers compared the results of eDFT calculations using a large set of XC functionals with diffusion Monte Carlo (DMC) calculations. DMC is a Quantum Monte Carlo (QMC) method, which explicitly accounts for electronic correlation and exchange. Although for fermions the fixed node error prevents DMC from being exact, it is more accurate than eDFT and belongs to the class of variational methods that can be systematically improved. This method is much more costly than conventional eDFT in terms of computational time, but it is expected to capture the intermolecular interactions with a high accuracy. It can thus be viewed as a good reference for benchmarking purposes. Based on these results [START_REF] Al-Hamdani | How strongly do hydrogen and water molecules stick to carbon nanomaterials?[END_REF]153], it appears that the choice of the XC approximation in eDFT is crucial, but also that it is difficult to predict the accuracy since the adsorption energy results from a subtle balance of the interactions, especially at a medium range.

Our collaborators followed a similar approach to study the adsorption of lithium on graphite. The choice of the system was made based on its relevance for energy applications and because the lithium ions have a small number of valence electrons, thus allowing the cost of the reference DMC calculations. The adsorption profiles on hollow, top and bridge sites were computed by varying the distance between the lithium and the carbon material. DMC results are then compared with several XC approximations. The comparison between the adsorption profiles on hollow site calculated with DMC and several XC functionals are shown in Figure 8.1.

Yet, knowing the adsorption energy of a single ion may not be enough. Indeed, the solvent molecules within the electrolyte will also interact strongly both with the surface and with the ions, which can result in large variations of the adsorption properties under realistic conditions [START_REF] Merlet | Influence of solvation on the structural and capacitive properties of electrical double layer capacitors[END_REF]. The main quantity to be determined is then the free energy of the adsorption profile, which should be obtained by sampling the whole liquid degrees of freedom. This quantity is not accessible to electronic structure methods due to their large computational cost, and it is necessary to resort to classical mechanics based methods instead. Classical MD is usually the method of choice for such purposes, but it suffers from high inefficiency for systems studied under large dilution conditions. In such cases, it may conveniently be replaced by molecular density functional theory (MDFT).

In Chapter 5, we introduced the constant potential electrode model and the theoretical background for the incorporation of the model in MDFT. In this chapter, we build upon these developments to study the effects of water solvent on the adsorption of lithium on graphite surface. The classical lithium-carbon interaction is parameterized by our collaborators using the eDFT calculations with the best choice of XC functional benchmarked on DMC. The present study is a preliminary step towards the comprehension of realistic electrochemical electron transfer reactions using MDFT.

Computational details

We mentioned in Section 3.3 that it is essential to apply correction to the hyper netted chain (HNC) functional for good estimates of the solvation free energy. In this work, we use the same bridge functional as in Chapter 6 [67,68] The external potential ϕ(r, Ω) is created by the graphite electrodes and the lithium ion whose interactions are described as the sum of Lennard-Jones and electrostatic interactions. The Lennard-Jones and charge parameters of graphite atoms and lithium ion are collected in Table 8.1. For the carbon atoms, the choice of the force field of Werder et al. was made based on a previous MD study, in which it was shown to provide a good estimate of the energy of adsorption of a water molecule on a graphene surface [START_REF] Al-Hamdani | How strongly do hydrogen and water molecules stick to carbon nanomaterials?[END_REF]. The ensemble averaged electrode charges are calculated with the iterative optimization method introduced in Chapter 5 with a Gaussian width of 0.40 Å. Hence the total charge on the carbon is fixed to be equal to -e.

In the first step of the iterative process, here we first perform the functional minimization without electrostatic interactions. Then, carbon charges are optimized in the presence of the inhomogeneous water charge density and the lithium cation. The functional is minimized again but in the presence of lithium charge and the previously determined electrode charges. The process is repeated until it converges, with a convergence criterion of 5 × 10 -4 on the relative change in solvation free energy between two consecutive steps.

MDFT calculations were performed using an in-house Fortran code and electrode charges are optimized using the constant potential MD software MetalWalls [START_REF] Coretti | MetalWalls: simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes[END_REF][START_REF] Marin-Laflèche | Met-alWalls: a classical molecular dynamics software dedicated to the simulation of electrochemical systems[END_REF]. We use a 24.672 × 21.366 × 40 Å

3 simulation box with 74 × 64 × 120 grid nodes and an angular grid of 196 orientations per grid node. We run calculations for a distance z between the electrode plane and the lithium varying between z = 1.0 and z = 10 Å with an increment of 0.2 Å between 1 and 6 Å and of 0.5 Å otherwise.

Fitting the carbon-lithium potential

In order to incorporate solvent effects, it is necessary to develop accurate classical interaction potentials. It is not possible to fit them directly on the QMC calculations due to the limited amount of data as shown in Figure 8.1. Instead, we pick the most accurate XC functional, vdW-DF-C09, and calculate the Li-graphite binding energy for a large number of distances. The intermolecular interaction should, in principle, account for four different effects: electrostatics, polarization, short-range repulsion, and dispersion. In our electrostatic model, the two former effects are explicitly introduced through the use of a +1 point charge on the lithium and of the calculation of partial (Gaussian) charges on the carbon atoms. These partial charges are calculated for each lithium-carbon distance using the same methodology as described in Chapter 5.

Concerning the short-range repulsion and the dispersion effects, the two main potentials used in the literature are the Lennard-Jones and the Born-Huggins-Mayer (BHM) ones. However, it appears that the electrostatic interaction was sufficient to account for the attractive part of the binding energy, as shown in Figure 8.2. The fitted potential should, therefore, add very few, if no contribution for the dispersion interaction, which agrees with the previous observation on the use of dispersioncorrected functionals. A well-known drawback of the Lennard-Jones potential is that it is not possible to fit the short-range repulsion and the dispersion term separately since they both involve the same parameters. We have, therefore, chosen a BHM potential instead, the analytical form of which is

V BHM (d) = A exp (-bd) - C 6 d 6 (8.1)
where A and b are the parameters describing the intensity and the range of the repulsion interaction, while C 6 is the dipole-dipole dispersion interactions, and d is the inter-particle distance. In principle, higher order terms could be included for dispersion, but as discussed above, this term is almost negligible in the case of the lithium ion. The fitted potential reproduces with a high accuracy the vdW-DF-C09 curves for the three types of absorption sites as shown in Figure 8.2. It also matches very well with the QMC results (on which it was not fitted) adsorption curves. The corresponding parameters are A = 91.17, b = 2.518, and C 6 = 1.107 (all numbers are given in atomic units).

Adsorption of the lithium ion in the presence of water

The fitted potential can directly be used in any classical molecular simulation, such as MD. Since we focus here on the adsorption free energy of the lithium ion on the carbon surface, we prefer to use MDFT, which is a computationally more efficient alternative. The solvation free energy of a single system can be computed within a few minutes on a single central processing unit (CPU), while it would require tens of CPU hours with MD. The free energy profiles obtained for the three adsorption sites in the presence of liquid water are shown in Figure 8.3A. The profiles are very different from the gas phase results. The minimum at ≈ 1.8 Å completely disappears and is replaced by a strongly repulsive wall. This shows that there is no preferential adsorption of lithium on the graphite surface in the aqueous phase.

This result is in qualitative agreement with a recent MD study on the adsorption of hydrated ions on graphene [START_REF] Loche | Breakdown of linear dielectric theory for the interaction between hydrated ions and graphene[END_REF], which provided a repulsive free energy profile over the whole range of considered distances. Yet, the latter study did not include any Coulombic interaction between the ion and the carbon surface, which is the main driving force for adsorption in the gas phase as discussed above. It is thus interesting to examine the various contributions to the total free energy, which are provided in Figure 8.3B. We observe that the electrostatic attraction between Li and C is, in fact, counterbalanced by the solvation free energy. The latter contains the electrostatic interactions of the water molecules with both the lithium and the graphite surface, which results in a strong screening effect. The extent of this screening effect was studied in detail in a recent study focused on gold surfaces [156]: the presence of the water molecules strongly impacts the polarization of the surface. Consequently, the total free energy is almost equal to the BHM contribution over the whole range of distances, except between 2 and 5 Å, where the solvation free energy overcomes the ion-surface Coulombic interaction, resulting in a more repulsive potential.

The effect of the solvent can be further analyzed by plotting the density profiles for various distances between the ion and the surface (Figure 8.4). At z = 8 Å, two regions with larger densities emerge, corresponding to the surface adsorbed water molecules at a distance of 3 Å from the surface [START_REF] Jeanmairet | Study of a water-graphene capacitor with molecular density functional theory[END_REF] on one hand and to the lithium ion first solvation shell on the other hand. At a distance of 5 Å, the solvation shell starts to overlap with the adsorbed layer at close contact to the electrode, which results in small depletion zones in the latter. These depletion zones remain observable at smaller distances, but the impact on the free energy becomes negligible with respect to the large short-range repulsion between the carbon and the lithium.

Chapter conclusions

In this study, we obtained classical interaction potentials which reproduce very well the energies calculated with electronic calculations. We then used the parameterized force field to compute the free energy profile for the adsorption of a lithium ion at the graphite surface in the presence of water as a solvent. MDFT was preferred over classical MD for this calculation due to its much lower computational cost. We showed that the low energy minimum obtained in the gas phase completely vanishes, resulting in an overall repulsive profile over the whole range of studied lithium-carbon distances. This is due to the screening of the attractive Coulombic term by the water molecules on the one hand and on the other hand to the interferences between the densities corresponding to the first adsorbed layer on carbon and to the first solvation shell of the cation.

The proposed multi-scale approach may be extended to a large variety of systems in the future. In particular, it would be interesting to study how adsorption varies within the alkali and alkali-earth cationic series. The study of larger organic ions would be very useful for the scope of supercapacitor devices, but this would require additional work to account accurately for the flexibility of the ion. Finally, the study of more polarizable electrode materials such as gold or platinum would be relevant in the context of catalysis. Changing metal would lead to very different adsorption properties for the solvent so that one can expect large differences in the adsorption free energy profiles of ions at electrodes.

Part III Conclusion and Appendices

Chapter 9

Conclusion and Perspectives

General conclusion

Electrochemical energy storage devices provide solutions to rectify the intermittent renewable energy production in various application fields. Specifically, organic redox flow battery is a promising candidate for stationary applications. The power density and energy density of a redox flow battery are closely related to the redox reactions taking place at the electrode/liquid electrolyte interfaces. Our group are particularly interested in the kinetics aspect of the redox reaction which, together with the transport properties of the electrolyte, determines the power output of the device.

Marcus theory introduced a theoretical framework that helps us understand the reaction kinetics of the redox reactions involving electron transfer. Moreover, molecular simulation methods that complement experimental works have been developed based on Marcus theory. These methods provide a perspective at molecular scale which may explain those electron transfer processes that Marcus' approximation fails to describe. Molecular simulations can also be applied to predict electron transfer properties of novel systems, such as heterogeneous electron transfer of quinones. Molecular dynamics (MD) is until now the simulation method of choice for the studies of electron transfer. However, the heavy computational cost prohibits MD from being used for systematic studies and high throughput screening. Such obstacles may be overcome by using molecular density functional theory (MDFT) with which we calculate properties of a system more efficiently via the variational procedure.

The main objective of this thesis is to develop MDFT to enable fast computational studies of electron transfer taking place at the electrochemical interfaces. Before the commencement of the thesis, MDFT at the standard level of approximation (homogeneous reference fluid approximation) and the two pressure correction methods that improve the predicted solvation free energy have been developed to study bulk solutions. Also, the framework that uses MDFT to study electron transfer and the strategy to incorporate the constant potential electrode model in MDFT have both been proposed and tested for simple systems.

During the thesis, the detailed assessment of the correctness of the two pressure correction methods on the study of electron transfer has been performed. On analysing the free energy profiles and the reorganization free energies calculated 9.2. PERSPECTIVES without correction and with two different correction methods, we concluded that the a posteriori correction strategy is not suitable for the study of electron transfer. Only the bridge functional correction strategy which respects the variational principle can be used.

On the other hand, we shown that the solvation structure and free energy profiles of anthraquinone (AQ) solvated in acetonitrile (MeCN) calculated by MD simulations using polarizable force fields can be well reproduced by MDFT calculations using generic force fields. We then make use of the numerical efficiency of MDFT to calculate systematically the reorganization free energies of a series of 1,4benzoquinone (BQ) derivatives. We found that the electron transfer processes of all these quinones in MeCN solvent can be well described by Marcus' approximation and the reorganization free energy decreases as the volume of functional group increases. The combined results of the above two works provide a guideline for high throughput screening of redox active molecules for redox flow battery application based on their electron transfer properties.

In order to study electrochemical interfaces, it is necessary to build model electrodes in MDFT. During the thesis, a new version of Fortran subroutines that incorporate the constant potential electrode model in MDFT are implemented. A Python interface is also built for improving the numerical efficiency of MDFT calculations of electrochemical interfaces. We then relied on these developments to study the effects of water solvent on the adsorption of lithium on graphite electrode surface. We found that the low energy minimum obtained in the gas phase completely vanishes, resulting in an overall repulsive adsorption free energy profile over the whole range of studied lithium-carbon distances. This is due to the screening of the attractive Coulombic term by the water molecules on the one hand and on the other hand to the interferences between the densities corresponding to the first adsorbed layer on carbon and to the first solvation shell of the cation.

Perspectives

In Chapter 6, we have shown that although the weighted density approximation (WDA) bridge functional correction conforms to the variational principle, the improvement made by the correction is limited. This might be due to the angular independent nature of the WDA bridge functional that we used. The impact on the orientational polarization of the solvent, which is the dominant effect for ions, is low. If accurate predictions of free energy profiles and reorganization free energies are expected, more sophisticated, angular dependent bridge functionals should be developed and applied.

In Chapter 7, we have seen that both structural and energetic properties of quinones solvated in MeCN calculated by MDFT without corrections are good as compared to MD simulations. Nevertheless, if some bridge functionals are developed for MeCN solvent, they should be able to further improve the quantitative agreement of the results. Another possible improvement route, especially when comparing to the results given by ab-initio based methods, may be developing quantum mechanics/MDFT hybrid approaches that can take solute's internal degrees of freedom into CHAPTER 9. CONCLUSION AND PERSPECTIVES account. As an example, our group has developed an approach in which the solute is dealt with at the quantum level using electronic density functional theory (eDFT) while the solvent degrees of freedom are accounted at the classical level using MDFT. With such approach, it is possible to have a full picture of the reaction free energy profiles. We believe that this method could be the appropriate compromise between precision and computational cost to make the calculation of electron transfer rate constant in solution feasible in the future.

Finally, the methodology established in Chapter 8 may be applied to a large variety of systems involving different redox active species, solvents, electrode materials, and surface morphology. In particular, it would be interesting to study how adsorption varies within the alkali and alkali-earth cationic series. The study of larger organic ions, such as quinones, would be very useful for the scope of energy storage devices, but this would require additional work mentioned above to account accurately for the solute's internal degrees of freedom. Furthermore, the study of more polarizable electrode materials such as gold or platinum would be relevant in the context of catalysis. Changing metal would lead to very different adsorption properties for the solvent so that one can expect large differences in the adsorption free energy profiles of ions at electrodes. These studies concerning adsorption would be preliminary steps towards the comprehension of realistic electrochemical electron transfer reactions. The framework reported above may be combined with this multiscale approach to study electron transfer at realistic electrochemical interfaces.

Appendix B Functional derivative

A functional is a natural extension of the concept of a function. Let f (v) be a function of v where v = [v 1 , v 2 , ..., v N ] is a N-dimensional vector. f (v) maps points in the N -dimensional space spanned by v to scalars (we will limit ourselves to real scalars). A functional F [f (v)] then further takes the function f (v) as argument and return again scalars. Let's take Ω id [n(r)] defined in Equation 3.12 as an example. where Ω id [n(r)] is a functional of particle number density n(r) which is itself a function of the coordinate r in the 3-dimensional space. The value of Ω id [n(r)] is related to the value of particle number density over its entire domain of definition, here the 3-dimensional space. The variation of f (v) due to an infinitesimal change of v is written as

df (v) = f (v + dv) -f (v) = N i=1 ∂f (v) ∂v i dv i . (B.2)
Similarly, if F is a functional of f (v), i.e., it depends on the value of f (v) on an interval of v, a ≤ v ≤ b, the variation of the functional can be expressed as cell potential The measure of the potential difference between the two redox half reactions in an electrochemical cell. [START_REF] Gentil | Aqueous organic and redox-mediated redox flow batteries: a review[END_REF][START_REF] Yuan | Membranes in non-aqueous redox flow battery: a review[END_REF][START_REF] Park | Investigation of charge transfer kinetics at carbon/hydroquinone interfaces for redox-active-electrolyte supercapacitors[END_REF]65 diabatic At the limit of weak electronic coupling between reactants and products states, the potential energy surfaces and free energy profiles are called diabatic. The occurrence of the electronic transition is restricted to the transition state located at the crossing point of the two free energy profiles. ix, 44, 46, 48-50, 122, 123 electron transfer A process by which an electron moves from one atom or molecule to another.

δF [f (v)] = F [f (v) + δf (v)] -F [f (v)] = b a δF [f (v)] δf (v) δf (
x , xv, 8-11, 13, 20, 36, 41-48, 51, 56, 57, 61, 66, 71, 72, 74, 75, 77, 80, 83, 84, 88-91, 95, 100, 109-111, 117 energy density The amount of energy stored in an energy storage device divided by its own volume. [START_REF] Gentil | Aqueous organic and redox-mediated redox flow batteries: a review[END_REF][START_REF] Yuan | Membranes in non-aqueous redox flow battery: a review[END_REF][START_REF] Park | Investigation of charge transfer kinetics at carbon/hydroquinone interfaces for redox-active-electrolyte supercapacitors[END_REF]109 force field The analytical expression for interaction potential. xiii, xv, 36, 56, 57, 84-86, 95, 98, 100, 103, 110 free energy profile Free energy of a system as a function of a reaction coordinate. ix, x, xv, 11, 42, 45-60, 71, 73, 74, 77, 80, 84, 88-90, 95, 101, 103, 105, 109-111, 122, 125, 126, 128 inner-sphere The solvent molecules in the first coordination shell of the redox species. Upon electron transfer, the inner-sphere contributes to the energy change through the changes in the bond length and bond angle between the redox species and the molecules in the coordination shell. [START_REF] Marcus | Electrostatic free energy and other properties of states having nonequilibrium polarization. I[END_REF][START_REF] Marcus | On the theory of oxidation-reduction reactions involving electron transfer. II. applications to data on the rates of isotopic exchange reactions[END_REF]57,[START_REF] Sergiievskyi | Fast computation of solvation free energies with molecular density functional theory: thermodynamic-ensemble partial molar volume corrections[END_REF]125 macrostate A thermodynamic equilibrium state characterized by a combination of three fixed thermodynamic variables that together determine the ensemble. 13-15

microstate A microscopic state of the system characterized by the position and momentum of all the constituting particles. 13-15, 17, 44, 49, 50, 53, 56, 123 outer-sphere The solvent molecules outside the inner-sphere. Marcus assumed that the outersphere behaves like an unsaturated dielectric continuum and, under the linear response approximation, respond linearly to the change of electronic state of the redox species. 46-48, 57, 61, 122 potential energy surface The potential energy of the system as a function of the coordinates of all the atoms in the system. ix, 44-46, 49, 50, 52, 53, 56, 122, 123 power density The amount of power an energy storage device can deliver divided by its own volume. [START_REF] Park | Investigation of charge transfer kinetics at carbon/hydroquinone interfaces for redox-active-electrolyte supercapacitors[END_REF]61,[START_REF] Zhao | Molecular density functional theory of solvation: From polar solvents to water[END_REF]109 quinones A family of redox active organic molecules. The main structure of quinones are the aromatic rings with a pair of hydrogen atoms substituted by oxygen atoms. xv, 7-11, 48, 71, 83-85, 89, 109-111 reorganization free energy The free energy cost to reorganize the solvent configurations corresponding to the minimum of the free energy profile of the reactants to those corresponding to the minimum of the free energy profile of the products while staying on the diabatic potential energy surface of the reactants. ix, xi, xiii, xv, 45-50, 53-55, 59, 61, 71, 75-80, 84, 85, 88, 89, 91-95, 109, 110, 121-123 trajectory The evolution of the system at equilibrium over time. 50-53, 55, 66, 83, 85, 86, 88 vertical energy gap The potential energy difference between the system containing reactants plus solvent and the system containing products plus solvent evaluated at the same atomic coordinates of all the atoms. ix, x, 42, 44, 50, 52-57, 88, 90, 122, 123, 125- The parameter describes the range of the repulsion interaction of the BHM potential B(r N , p N ) Some physical quantity of a system that can be expressed as a function of the coordinate and momentum of all the particles constituting the system B(r N , p N ; N ) Some physical quantity of a system that may depend on the number of particles

1. 1

 1 Schematic representation of a rechargeable lithium ion battery during discharge process. Reprinted from Ref. [1]. . . . . . . . . . . . . . . . 1.2 Schematic representation of a vanadium redox flow battery during a discharge process. Reprinted with permission from [3]. Copyright 2011 American Chemical Society. . . . . . . . . . . . . . . . . . . . . 1.3 Molecular structures of (A) 1,2-benzoquinone and (B) 1,4-benzoquinone (BQ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Molecular structures of (A) naphtoquinone and (B) anthraquinone (AQ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Molecular structures of hydroquinone . . . . . . . . . . . . . . . . . . 3.1 The value of x ln x and (x -1) as functions of x . . . . . . . . . . . . 3.2 Schematic representation of the 3 Euler angles (θ, ϕ, ψ) . . . . . . . . 4.1 Schematic representation of the reactants (R), products (P), and the surrounding solvent molecules in reaction 4.1. The arrows denote the water molecular dipole. The longer Fe-OH 2 bond length in the cation is indicated by the larger ionic radius [102]. . . . . . . . . . . . . . . . 4.2 Schematic representation of diabatic (solid) and adiabatic (dashed) potential energy surfaces of systems containing reactants plus solvent (R) and systems containing products plus solvent (P). The vertical energy gap defined in Equation 4.18 is also depicted. . . . . . . . . . 4.3 Schematic representation of the two-state model. Free energy profiles of the reactants plus solvent (R) and the products plus solvent (P) are plotted as functions of the global reaction coordinate. The reaction free energy, activation free energies, and reorganization free energies are also depicted. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Schematic representation of the inverted region [102]. . . . . . . . . . 5.1 A set of point charges (illustrated on the left hand side) can be considered as a set of point charges screened by a set of Gaussian charge distributions and compensated by the same set of Gaussian charge distributions in reverse sign (illustrated on the right hand side) . . . . 5.2 A snapshot of a molecular dynamics (MD) simulation of an electric double layer capacitor composed of 322 EMIM-TFSI ions pairs and 2 graphite electrodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF FIGURES 5.3 The flowchart illustrating the procedure of finding ensemble averaged properties of an electrochemical cell using molecular density functional theory (MDFT). . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.1 Free energy profiles of the (A) Cl → Cl + and (B) Cl -→ Cl electron transfer half reactions computed using the hyper netted chain (HNC) functional (full) without and (dashed) with the partial molar volume (PMV) correction and (dotted) using the weighted density approximation (WDA) bridge functional. The MD data obtained by Hartnig and Koper are the circles [121]. . . . . . . . . . . . . . . . . . . . . . 74 6.2 Free energy of Cl + and Cl -computed using the HNC functional (full)

80 7. 1

 801 Difference between functional values calculated with HNC and WDA functional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.5 Equilibrium radial polarization density around Cl -and Cl + computed with HNC and WDA functionals. . . . . . . . . . . . . . . . . . . . . Molecular structures of duroquinone (DQ) . . . . . . . . . . . . . . . 84 7.

  82 and 4.59. . . . . . . 94 8.1 (A) Comparison of the adsorption energies obtained with various electronic density functional theory (eDFT) functionals and diffusion Monte Carlo (DMC) for the adsorption of the lithium on the hollow site of graphite. (B) Same as (A) but subtracting the adsorption energy at a distance of 3 Å. . . . . . . . . . . . . . . . . . . . . . . . 99 8.2 Comparison of the fitted potential with the Quantum Monte Carlo (QMC) and vdW-DF-C09 energies for the three different absorption sites. Panel (A): hollow site; panel (B): bridge site; panel (C): top site.102 8.3 (A) Adsorption free energy for a lithium ion on the graphite surface in the presence of water, computed using MDFT, for the three adsorption sites. The energy variation in the absence of water is also shown for comparison. (B) Contributions to the total free energy for the hollow adsorption site in the MDFT calculation. . . . . . . . . . . 104 8.4 Projection of the solvent densities computed using MDFT for various lithium-carbon distances. . . . . . . . . . . . . . . . . . . . . . . . . 105 xi List of Tables 3.1 The partial charges and Lennard-Jones parameters of the oxygen (O) and hydrogen (H) of the extended simple point charge (SPC/E) water model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The partial charges and Lennard-Jones parameters of the united atom acetonitrile (MeCN) model [93]. . . . . . . . . . . . . . . . . . . . . . 6.1 Reorganization free energies (in kJ/mol) calculated with the HNC functional, with the PMV correction using equations 4.59 and 4.61, with the PMV correction using equations 4.58 and 4.60 and with the WDA bridge functional. . . . . . . . . . . . . . . . . . . . . . . . . . 8.1 Force field parameters used in the MDFT simulations. Mixed parameters are computed using the Lorentz-Berthelot rules (except for the C-Li interaction, which does not affect the MDFT results). . . . . . . xiii

  energy storage devices . . . . . . . . . . 1.1.1 Rechargeable batteries . . . . . . . . . . . . . . . . . . . . 1.1.2 Redox flow batteries . . . . . . . . . . . . . . . . . . . . . 1.2 Organic molecules for energy storage . . . . . . . . . . . . 1.2.1 Organic redox flow battery . . . . . . . . . . . . . . . . . 1.2.2 Quinones as organic redox active materials . . . . . . . . . 1.2.3 Reaction kinetics . . . . . . . . . . . . . . . . . . . . . . . 1.3 Theoretical studies of electron transfer . . . . . . . . . . 1.3.1 Marcus theory of electron transfer . . . . . . . . . . . . . 1.3.2 Molecular simulations on electron transfer . . . . . . . . . 1.3.3 Molecular density functional theory (MDFT) . . . . . . . 1.3.4 The application of MDFT to electron transfer . . . . . . . 1.3.5 Pressure correction to MDFT . . . . . . . . . . . . . . . . 1.3.6 Constant potential electrode in MDFT . . . . . . . . . . .

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: Molecular structures of (A) 1,2-benzoquinone and (B) 1,4-benzoquinone (BQ)

1. 2 .Figure 1 . 5 :

 215 Figure 1.5: Molecular structures of hydroquinone

Contents 2 . 1

 21 Entropy and temperature . . . . . . . . . . . . . . . . . . 13 2.2 Ensemble theory and micro-canonical ensemble . . . . . 15 2.3 Canonical ensemble . . . . . . . . . . . . . . . . . . . . . . 16 2.4 Grand canonical ensemble . . . . . . . . . . . . . . . . . . 17

3. 1 .Figure 3 . 1 :

 131 Figure 3.1: The value of x ln x and (x -1) as functions of x

Figure 3 . 2 :

 32 Figure 3.2: Schematic representation of the 3 Euler angles (θ, ϕ, ψ)

Figure 4 . 1 :

 41 Figure 4.1: Schematic representation of the reactants (R), products (P), and the surrounding solvent molecules in reaction 4.1. The arrows denote the water molecular dipole. The longer Fe-OH 2 bond length in the cation is indicated by the larger ionic radius [102].

Figure 4 . 3 :

 43 Figure 4.3: Schematic representation of the two-state model. Free energy profiles of the reactants plus solvent (R) and the products plus solvent (P) are plotted as functions of the global reaction coordinate. The reaction free energy, activation free energies, and reorganization free energies are also depicted.

Figure 4 . 4 :

 44 Figure 4.4: Schematic representation of the inverted region [102].

Contents 5 . 1

 51 Electrode and electrolyte models . . . . . . . . . . . . . . 62 5.2 Periodic boundary condition (PBC) . . . . . . . . . . . . 62 5.2.1 Electrostatic interactions in three dimensional PBC . . . . 63 5.2.2 Electrostatic interactions in two dimensional PBC . . . . 64 5.2.3 Constant potential and electroneutrality constraints . . . 65 5.3 Constant potential electrode model in MDFT . . . . . . 66 5.3.1 Approaches to the ensemble average . . . . . . . . . . . . 66 5.3.2 Practical aspects . . . . . . . . . . . . . . . . . . . . . . . 68

Figure 5 .

 5 Figure 5.1: A set of point charges (illustrated on the left hand side) can be considered as a set of point charges screened by a set of Gaussian charge distributions and compensated by the same set of Gaussian charge distributions in reverse sign (illustrated on the right hand side)

Figure 5 . 2 :

 52 Figure 5.2: A snapshot of a MD simulation of an electric double layer capacitor composed of 322 EMIM-TFSI ions pairs and 2 graphite electrodes.

Figure 5 . 3 :

 53 Figure 5.3: The flowchart illustrating the procedure of finding ensemble averaged properties of an electrochemical cell using MDFT.

Contents 6 . 1

 61 Selected correction methods . . . . . . . . . . . . . . . . . 72 6.2 Computational details . . . . . . . . . . . . . . . . . . . . . 72 6.3 Results and discussion . . . . . . . . . . . . . . . . . . . . 73 6.3.1 Assessment of partial molar volume correction . . . . . . . 73 6.3.2 Assessment of bridge functional correction . . . . . . . . . 77 6.4 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . 80

Figure 6 . 1 :

 61 Figure 6.1: Free energy profiles of the (A) Cl → Cl + and (B) Cl -→ Cl electron transfer half reactions computed using the HNC functional (full) without and (dashed) with the PMV correction and (dotted) using the WDA bridge functional. The molecular dynamics (MD) data obtained by Hartnig and Koper are the circles [121].

Figure 6 . 3 :

 63 Figure 6.3: Average vertical energy gap computed using the HNC (black solid line)and WDA functional (red solid line) as a function of the atomic charge q. The dotted and dashed lines are the tangent to the HNC curve at q = -1 and q = 1, respectively.

Figure 6 . 4 :

 64 Figure 6.4: Difference between functional values calculated with HNC and WDA functional.

Figure 6 . 5 :

 65 Figure 6.5: Equilibrium radial polarization density around Cl -and Cl + computed with HNC and WDA functionals.

Contents 7 . 1

 71 Computational details . . . . . . . . . . . . . . . . . . . . . 85 7.2 Results and discussion . . . . . . . . . . . . . . . . . . . . 86 7.2.1 Benchmark of MDFT results . . . . . . . . . . . . . . . . 86 7.2.2 Electron transfer of benzoquinone derivatives . . . . . . . 89 7.3 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . 95

Figure 7 . 1 :Figure 7 . 2 :

 7172 Figure 7.1: Molecular structures of duroquinone (DQ)

Figure 7 . 3 :

 73 Figure 7.3: Two dimensional maps of acetonitrile (MeCN) density and polarization maps in the aromatic plane of AQ computed with AIMD, classical MD and MDFT. Grey, blue, and red circles represent carbon, hydrogen, and oxygen atoms, respectively. Polarization is represented by arrows oriented in the direction of the projection of P and whose length is proportional to ||P ||

Figure 7 . 4 :

 74 Figure 7.4: Solvent density isosurfaces of (A) AQ and (B) 1,4-benzoquinone (BQ) at 10 -2 × n 0 calculated by MDFT. Oxygen, hydrogen, and carbon atoms are displayed in red, cyan and white respectively.

Figure 7 . 5 :

 75 Figure 7.5: Free energy profiles of electron transfer half reactions of (A) AQ and anthraquinone anion (AQ -) calculated by classical MD and MDFT, (B) BQ and 1,4-benzoquinone anion (BQ -), (C) DQ and duroquinone anion (DQ -) calculated by MDFT. As described in the text, the vertical energy gaps were shifted by the same constant value, which was chosen in order to align the minima for the AQ/AQ -redox couple. Solid lines are fit of MDFT results using Marcus theory.

Figure 7 . 6 :

 76 Figure 7.6: Primary molecular structure of (A) 2-R-BQ, (B) 2,3-R-BQ, (C) 2,6-R-BQ, (D) 2,5-R-BQ, (E) 2,3,5-R-BQ, and (F) 2,3,5,6-R-BQ derivatives with R's represent the functional groups. R's in the figure represent the sites to be substituted by functional groups.

Figure 7 . 7 :Figure 7 . 8 :

 7778 Figure 7.7: Reorganization free energies of functionalized BQ derivatives as a function of substitution positions. The notation of primary molecular structures follows those in Figure 7.6. The black dashed line is the reorganization free energy of BQ.

Figure 7 . 9 :

 79 Figure 7.9: Ideal (black), external (blue), and excess (red) contributions to the total (green) reorganization free energy of BQ and its monofunctionalized derivatives computed usingEquations 3.82 and 4.59.

Contents 8 . 1

 81 Computational details . . . . . . . . . . . . . . . . . . . . . 100 8.2 Fitting the carbon-lithium potential . . . . . . . . . . . . 101 8.3 Adsorption of the lithium ion in the presence of water . 101 8.4 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . 103

Figure 8 . 1 :

 81 Figure 8.1: (A) Comparison of the adsorption energies obtained with various electronic density functional theory (eDFT) functionals and diffusion Monte Carlo (DMC) for the adsorption of the lithium on the hollow site of graphite. (B) Same as (A) but subtracting the adsorption energy at a distance of 3 Å.

Figure 8 . 3 :

 83 Figure 8.3: (A) Adsorption free energy for a lithium ion on the graphite surface in the presence of water, computed using MDFT, for the three adsorption sites. The energy variation in the absence of water is also shown for comparison. (B) Contributions to the total free energy for the hollow adsorption site in the MDFT calculation.

Figure 8 . 4 :

 84 Figure 8.4: Projection of the solvent densities computed using MDFT for various lithium-carbon distances.

Ω 3 n

 3 id [n(r)] = k B T R 3 n(r)(ln [n(r)Λ 3 ] -1)drµ R

Fadiabatic 46 anolyte

 46 [f (v)] = a(v)f (v)dv, δF [f (v)] = a(v)δf (v)dv.-linear functional in the following formF [f (v)] = f (v) ln f (v)dv (B.6) for which δF [f (v)] = [δf (v) ln f (v) + f (v)δ ln f (v)]dv = [ln f (v) + 1]δf (v)dv.When the electronic coupling between reactants and products states is non-negligible, the potential energy surfaces and free energy profiles are called adiabatic. ix, 44, The electrolyte tank of a redox flow battery containing reductant. 5, 8 average vertical energy gap The free energy change upon changing the external potential from reactants state to products state while keeping the solvent density field unchanged. x, 57-59,73, 75-77, 80, 88 catholyte The electrolyte tank of a redox flow battery containing oxidant. 5, 8

2 M

 2 127 ∆F elec M (η) The electrostatic part of external functional difference ∆e The charge transferred from one reactant to the other ∆E(r N ) Vertical energy gap of diabatic potential energy surface ∆F The reaction free energy change ∆F (x) The difference between the two free energy profiles of the reactants and the products state ∆G Reaction free energy ∆G ‡ The activation free energy ∆G ‡ M,ϵ The activation free energy of the forward (M = R) or the reverse (M = P) reaction defined with respect to vertical energy gap ∆G ‡ M The activation free energy for the forward (M = R) or the reverse (M = P) reaction ∆G ϵ The reaction free energy defined with respect to vertical energy gap ∆V Partial molar volume of the solute δ(rr i ) The Dirac delta distribution δ 3 (rr k ) The three dimensional Dirac delta distribution centered at point r k ϵ A given value of vertical energy gap ϵ 0 The vacuum permittivity ϵ LJ The depth of the Lennard-Jones potential well ϵ M The mean value of the Gaussian of the reactants (M = R) or the products (M = P) state ϵ op , ϵ s Optical and static frequency dielectric permittivity of the solvent ϵ R , ϵ P The values of vertical energy gap corresponding to the minimum of the free energy profiles of the reactants (R) or the products (P) state η The interpolation coefficient γ(r) The integral equation describes the indirect correlation ĉ(k = 0) Fourier transform of c (2) b (r -r ′ ) taken at k = 0 n(r) The microscopic particle density The reorganization free energy of the forward (M = R) or the reverse (M = P) reaction defined with respect to vertical energy gap λ M The reorganization free energy of the forward (M = R) or the reverse (M = P) reaction λ o Outer-sphere reorganization free energy λ HNC M The reorganization free energies obtained with the hyper netted chain (HNC) functional for the oxidation state M Ω Grand potential Ω[f ] A functional of probability density f (r N , p N ; N ) Ω b Grand potential of the homogeneous reference fluid having density n 0 Ω ϕ [n(r)] The full functional Ω id Grand potential of ideal gas Ω id [n(r)] A functional of arbitrary particle density n(r) of ideal gas in grand-canonical ensemble whose minimum is the grand potential of ideal gas Ω ± The regions occupied by the positive or negative electrode ϕ(r) Externally applied potential on the system Φ(r N ) The total external energy of a system due to externally applied potential ϕ(r) ϕ ′ (r) An external potential other than ϕ(r) that gives rise to the same n eq (r) ϕ elec (r, Ω) The electrostatic part of the external potential energy ϕ LJ (d) The Lennard-Jones part of the external potential energy ϕ η (r, Ω) The biased external potentials ϕ M (r) The external potential generated by the solute molecules in reactant states (M = R) or products state (M = P) Φ M (r N ) The external energy of the reactant states (M = R) or products state (M = P) ψ(r) The intrinsic chemical potential Ψ ± The value of uniform electric potential imposed in the regions occupied by the positive or negative electrode, Ω ± , with r being a point in the regions Ω ± ρ Particle density ρ(r, Ω) An arbitrary particle density field that depends on both position r and orientation Ω of the solvent molecules ρ b Equilibrium particle density of the homogeneous reference bulk fluid ρ η (r, Ω) The equilibrium solvent density fields given by the biased external potential ϕ η (r, Ω) ρ c (r) The charge density of the solute molecule ρ eq (r, Ω) The equilibrium particle density that depends on both position and orientation ρ P (r, Ω) The equilibrium solvent density of products state ρ R (r, Ω) The equilibrium solvent density of reactants state ρ s (r) The inhomogeneous charge distribution generated by the inhomogeneous solvent density ρ eq (r, Ω) ρ j (r) The Gaussian charge distribution assigned to an electrode atom, labelled j ρ k (r) The charge distribution of a point charge σ LJ The distance at which the Lennard-Jones potential energy is zero σ The variance of the Gaussian of the reactants (M = R) or the products (M = P) state σ w The Gaussian width of the Gaussian weighting function w(r) σ Q The Gaussian width that can be tuned to account for the metallicity of the electrode ρ s (r, Ω) The charge density of a single solvent molecule located at the origin with an orientation Ω n(r) The weighted density using a Gaussian weighting function w(r) Ξ Grand partition function Ξ id Grand-partition function of ideal gas A The parameter describes the intensity of the repulsion interaction of the Born-Huggins-Mayer (BHM) potential a A parameter for angular independent WDA bridge functional uniquely defined by imposing zero compressibility route pressure a 1 , a 2 The ionic radii of reactants 1 and 2 including the inner-sphere b

b ( 2 ) 6 F ( 0 )

 260 (r, r ′ ) Bridge functionCPre-exponential factor c(1) (r; n eq (r)) One body direct correlation function c(2) (r, r ′ ; n eq (r)) Two body direct correlation function c(2)b (r -r ′ ) The two body direct correlation function of homogeneous fluid C The dipole-dipole dispersion interactions of the BHM potential The free energy at the crossing point of the two free energy profiles defined with respect to vertical energy gap f (r N , p N ; N ) The probability density of finding a system in grand-canonical ensemble having N particles and energy H at thermodynamic equilibriumf (v) A function of v F (x ‡ )The free energy at the crossing point of the two free energy profilesF [f (v)] A functional of function f (v)f ′ (r N , p N ; N ) Probability densities other than f (r N , p N ; N ) with the only condition that they are normalized, i.e. Tr cl f ′ = 1.f ′ eqThe equilibrium probability density associated to ϕ ′ (r)f N (r N , p N )The probability density of a micro-canonical ensemblef N 0 (r N , p N )The probability density of finding a system in canonical ensemble having energy H at thermodynamic equilibriumF MFree energy of the reactants plus solvent (M = R) or the products plus solvent (M = P)F M (ϵ)The free energy profile of the reactants (M = R) or the products (M = P) state with vertical energy gap being the reaction coordinateF M (x)Free energy profile of the reactants plus solvent (M = R) or the products plus solvent (M = P) with respect to the reaction coordinate x F exc The excess Helmholtz free energyF exc [n(r)]The excess part of full functionalF extThe external Helmholtz free energy F id Helmholtz free energy of ideal gasF id [n(r)] A functional of arbitrary particle density n(r) of ideal gas in grand-canonical ensemble which reduces to the Helmholtz free energy of ideal gas for equilibrium particle density of ideal gasF intThe intrinsic Helmholtz free energy F int [n(r)] The intrinsic part of full functional g (2) (r, r ′ ) Pair distribution function h The Planck constant h (2) (r, r ′ ) Two body total correlation function k Reaction rate constant k B Boltzmann constant l x , l y , l z The length of an orthorhombic simulation box in the x, y, z directions, An arbitrary particle number density field n 0 (r) Equilibrium particle density of the reference state N 1 , N 2 Number of particles in the subsystems 1 and 2 N Q Number of electrode atoms N q Number of electrolyte atoms N c The number of atoms of the solute molecule

Table 3 .

 3 2:The partial charges and Lennard-Jones parameters of the united atom MeCN model[START_REF] Edwards | A computer simulation study of the dielectric properties of a model of methyl cyanide: I. the rigid dipole case[END_REF].

	site ϵ LJ [kJ/mol] σ LJ [Å]	q [e]
	O	0.65	3.166	-0.8476
	H	0	0	0.4238
	Table 3.1: The partial charges and Lennard-Jones parameters of the oxygen (O)
	and hydrogen (H) of the SPC/E water model.
	site ϵ LJ [kJ/mol] σ LJ [Å]	q [e]
	N	0.416	3.3	-0.398
	C	0.416	3.4	0.129
	CH 3	1.59	3.6	0.269

use only rigid molecule models and consider only non-bonded interactions. Throughout the thesis, all the repulsion and dispersion interactions between the solute and solvent are modelled by the Lennard-Jones potential energy with Lorentz-Berthelot mixing rule. The expression for the Lennard-Jones potential energy used in the in-house developed code of MDFT is

  3 3 Å 3 . As for the orientations of solvent molecules, a common choice is 196 orientations per grid node.

	Chapter 4
	Electron transfer
	4.7.3 Constructing free energy profiles with MDFT . . . . . . .
	4.7.4 Definitions of numerous free energies in MDFT . . . . . .

Contents 4.1 Fundamentals of electron transfer . . . . . . . . . . . . . . 4.1.1 Electron transfer mechanism . . . . . . . . . . . . . . . . 4.1.2 Global reaction coordinate . . . . . . . . . . . . . . . . . . 4.2 Linear response approximation . . . . . . . . . . . . . . . 4.2.1 Inner-sphere and outer-sphere . . . . . . . . . . . . . . . . 4.2.2 Reorganization free energy . . . . . . . . . . . . . . . . . . 4.2.3 Electron transfer at electrochemical interface . . . . . . . 4.3 Statistical physics for electron transfer . . . . . . . . . . . 4.3.1 Constructing free energy profiles . . . . . . . . . . . . . . 4.3.2 The ideal global reaction coordinate . . . . . . . . . . . . 4.4 Molecular dynamics on electron transfer . . . . . . . . . 4.4.1 Working principles of molecular dynamics . . . . . . . . . 4.4.2 Application to electron transfer . . . . . . . . . . . . . . . 4.5 Vertical energy gap . . . . . . . . . . . . . . . . . . . . . . 4.6 Microscopic linear response approximation . . . . . . . . 4.7 MDFT on electron transfer . . . . . . . . . . . . . . . . . 4.7.1 Reexamination of the statistical physics . . . . . . . . . . 4.7.2 Average vertical energy gap . . . . . . . . . . . . . . . . .

  theory. Readers interested in models beyond the linear response approximation are referred to "Chapter 4 Electron Transfer Reactions in Enzymes: Seven Things that Might Break Down in Vanilla Marcus Theory and How to Fix Them if They Do" written by Lande et al. [58].

Table 6 .

 6 1: Reorganization free energies (in kJ/mol) calculated with the HNC functional, with the PMV correction using equations 4.59 and 4.61, with the PMV correction using equations 4.58 and 4.60 and with the WDA bridge functional.

Table 8 .

 8 1: Force field parameters used in the MDFT simulations. Mixed parameters are computed using the Lorentz-Berthelot rules (except for the C-Li interaction, which does not affect the MDFT results).

			8.1. COMPUTATIONAL DETAILS
	Atoms ϵ LJ [kJ/mol] σ LJ [Å]	q [e]
	C	0.2364	3.214	1
	Li	0.07648	2.216	Fluctuating

STUDY OF ELECTRON TRANSFER

the linear response approximation of Marcus theory, while Cl -→ Cl deviates from it. Overall, the WDA functional does not modify significantly the results obtained without correction [60]. This might be because it is an angular independent correction having a low impact on the polarization of the solvent, which is the dominant effect for ions.

List of Symbols

⟨B⟩

The ensemble average of some physical quantity B

⟨N ⟩

The average number of particles of a grand-canonical ensemble

Free energy functional of particle density ρ(r, Ω)

F[ρ eq (r, Ω)] The free energy functional of the equilibrium particle density

The intrinsic part of out-of-equilibrium solvation free energy

The solvation free energy calculated with HNC functional and partial molar volume (PMV) correction

The solvation free energy calculated with functional including angular independent WDA bridge functional The cross electrostatic interaction energy between Gaussian distributed charges and point charges

U qq

The electrostatic interaction energy between point charges

LIST OF SYMBOLS

x Specific value of the physical quantity that serve as the reaction coordinate X(r N ) The physical quantity chosen to serve as the reaction coordinate

x ‡

The location of the crossing point of the two free energy profiles

x M The values of the reaction coordinate corresponding to the minimum of the free energy profiles of the reactants (M = R) or the products (M = P) state z Activity

Z N

The configuration integral A 1 , A 2 Subsystems 1 and 2 that divide a system in the micro-canonical ensemble Ox, Red Redox active species in its oxdiant and reductant form