A molecular density functional theory approach to electron transfer reactions

Molecular density functional theory, an efficient computational tool, provides new insights into the study of electron transfer reactions in bulk and interfacial water.

Beyond the dielectric continuum description initiated by Marcus theory, the standard theoretical approach to study electron transfer (ET) reactions in solution or at interfaces is to use classical force field or ab initio molecular dynamics simulations. We present here an alternative method based on liquid-state theory, namely molecular density functional theory, which is numerically much more efficient than simulations while still retaining the molecular nature of the solvent. We begin by reformulating molecular ET theory in a density functional language and show how to compute the various observables characterizing ET reactions from an ensemble of density functional minimizations. In particular, we define within that formulation the relevant order parameter of the reaction, the so-called vertical energy gap, and determine the Marcus free energy curves of both reactant and product states along that coordinate.
Important thermodynamic quantities such as the reaction free energy and the reorganization free energies follow. We assess the validity of the method by studying the model Cl 0 / Cl + and Cl 0 / Cl À ET reactions in bulk water for which molecular dynamics results are available. The anionic case is found to violate the standard Marcus theory. Finally, we take advantage of the computational efficiency of the method to study the influence of a solid-solvent interface on the ET, by investigating the evolution of the reorganization free energy of the Cl 0 / Cl + reaction when the atom approaches an atomistically resolved wall.
Electron transfer (ET) reactions play a central role in a wide range of chemical systems, including energy storage and harvesting in electrochemical devices or biological processes such as aerobic respiration and photosynthesis. This ubiquity can explain the considerable amount of experimental, theoretical and simulation studies that have been devoted to this class of reactions. 1 The widely accepted theory of ET reaction in solution was proposed by Marcus. [2][3][4] It is based on the description of the solvent by a dielectric continuum. The macroscopic uctuations of the solvent are represented by an out-of-equilibrium polarization eld, and the free energy is a functional depending quadratically on this polarization. It eventually provides a simple two-state picture, where the free energy of each state depends quadratically on a reaction coordinate. This famous two-parabola picture has been used with great success to interpret experimental results and to make predictions. 5 However, Marcus theory does not take into account the molecular nature of the solvent which can break the linear assumption of solvent response. In such cases, we must resort to molecular simulation.
The vast majority of simulation studies on ET reaction have been carried out using Molecular Dynamics (MD). For example, the pioneering work of Warshel demonstrated that the vertical energy gap is the appropriate reaction coordinate 6 and that the uctuations of this quantity are Gaussian. Such Gaussian statistics give rise to the famous parabola picture of Marcus for the free energy prole. A strict Gaussian behaviour is equivalent to a linear response of the solvent to the eld generated by the solute; it also implies that the two free energy parabolas have the same curvature because the solvent uctuations are identical for the two states. 7,8 This validity of the Gaussian assumption has been veried in numerous studies since, for ET in solution 9,10 or in complex biological systems, 11,12 using either classical or ab initio MD. However, there is evidence that some systems do not obey the Marcus assumption i.e. the free energy curves of the two states cannot be represented by a pair of identical parabolas. There are several possible origins of such a discrepancy, 13 in particular the fact that reactant and product may have quite different solvation states. This can happen when the ET occurs between neutral and charged states, as predicted by Kakitani and Mataga [14][15][16] and observed since in classical [17][18][19] and ab initio 20,21 simulations. Several extensions of Marcus theory have been put forward to take into account the various origins of nonlinearity. [21][22][23][24] The investigation of ET reactions by MD is quite challenging since it usually requires the computation of the free energy curves as a function of an order parameter which remains a demanding task. Indeed, it necessitates a proper sampling of the solvent congurations around the barrier. If the activation energy is high it requires the use of biases such as umbrella sampling 25 coupled with histogram analysis techniques to reconstruct the unbiased data. [26][27][28] This typically implies to run simulations on half a dozen ctitious intermediate states to study a single system.
To compute free energies, there exist alternative techniques based on statistical theory of liquids, which offer the advantage of keeping a molecular description of the solvent while avoiding to sample the instantaneous microscopic degrees of freedom. Among the different approaches one can mention integral equation theory either in its molecular 29 or multiple sites formulation (RISM) 30,31 and its 3D-RISM version. 32,33 Another method is the classical density functional theory (cDFT) of liquids 34,35 which describes the response of a uid in the presence of a perturbation by introducing a functional of the uid density. Minimization of the functional yields the grand potential at equilibrium uid density. Some of us have previously introduced the molecular density functional theory (MDFT) 36,37 which is able to provide the solvation free energy and the solvation structure of any solute embedded in a molecular solvent described by its inhomogeneous density eld. The solvent density is a function of space coordinates and of the orientation; hence the functional must be minimized on a 6D grid: 3 dimensions for the cartesian coordinates and 3 dimensions for the three Euler angles. This formalism can be used to solvate any simple or complex solutes. 38 We proposed functionals for several solvents 39,40 with particular attention paid to the case of water. 41,42 The most advanced version of the functional is equivalent to the molecular Ornstein-Zernike theory supplemented by the hypernetted-chain closure (HNC) 43 for the solute-solvent correlations and can be minimized efficiently thanks to the use of rotational invariants in an optimal frame. The accuracy on the predictions of solvation free energies is promising as illustrated on the FreeSolv database. 44 We shall take advantage of this accuracy to put forward an efficient way to compute the free energy curves.
An application of the MDFT formalism to ET reactions in acetonitrile was proposed some years ago. 39 In this article we extend this approach and apply it to ET in aqueous solutions. In Section 1, aer recalling some basics of ET theory and giving a very short description of the MDFT framework, we show how to compute the key quantities of ET reactions with MDFT. In particular, we show that the average vertical energy gap is an appropriate order parameter for the ET reaction. We prove that for a given set of external potentials the free energy functional is actually a function of this order parameter. We derive expressions to compute the free energy curves (FEC) and the reorganization free energies.
In Section 2 we rst validate the framework on the simplest solute in water, i.e. a single neutral or charged chlorine atom modeled by a Lennard-Jones site, before studying the inuence of the presence of a solid-solvent interface on the reorganization free energy; to this purpose we investigate the ET of this solute as a function of its distance to an atomistically resolved wall.

Electron transfer reaction
We limit ourselves to the study of ET reactions of solutes which interact with the solvent through a classical force eld. Moreover, the solutes we consider in this article are rigid entities composed of a set of Lennard-Jones sites and point charges. An ET reaction involving two solutes of this type would correspond to an outer-sphere ET because there are no structural changes of the solutes. This implies that the ET reaction is completely controlled by the solvent response, as considered in Marcus' original paper. 3 The physics of the system can be described by the two crossing free energy curves of the system before (0) and aer (1) the ET. A schematic view of the two FEC is presented in Fig. 1 where some of the quantities necessary to describe the process are shown. The order parameter x describes the solvent conguration around the solute, thus the abscissa x 0 of the minimum of the FEC W 0 corresponds to a solvent in equilibrium with state 0. We emphasize that several microscopic solvent congurations correspond to an identical value of the order parameter.
Values of the order parameter differing from x 0 correspond to solvent congurations that are not in equilibrium with state 0. The more the solvent conguration differs from equilibrium, the more the free energy increases. The difference between the minima of the 2 FEC corresponds to the free energy difference between the two states, each surrounded by a solvent in equilibrium with these states, i.e. the reaction free energy, DW. Two others key quantities appear in Fig. 1: the reorganization free energies l 0 (resp. l 1 ) which represent the cost in free energy to solvate state 0 (resp. 1) in a solvent in equilibrium with the other Fig. 1 Schematic representation of a solvent controlled ET reaction. The diabatic free energy curves for state 0 and 1 are represented in solid black and dashed red, respectively. The reorganization free energies l 0 and l 1 for the two states are represented with solid arrows, the reaction free energy with a dotted arrow. state. The difference in free energy between the transition state (the crossing point) and state 0 controls the kinetics of the 0 / 1 reaction.
We emphasize that in the Marcus picture, the solvent is treated as a continuum which responds linearly to the electric eld generated by the solute. This implies that the FEC of the 2 states are identical parabolas. As a consequence, there is a unique reorganization free energy parameter l ¼ l 0 ¼ l 1 . The objective of the present paper is to show how to compute the various quantities appearing in Fig. 1, within MDFT. This could be a way to test the validity of Marcus assumption when the molecular nature of the solvent is taken into account, while taking advantage of the numerical efficiency of MDFT compared to MD.

Molecular density functional theory
We briey recall the fundamentals of MDFT which belong to the more general class of cDFT. Based on the Hohenberg-Kohn ansatz, 45 Mermin introduced the framework of density functional theory (DFT) at nite temperature for the inhomogeneous electron gas. 34 Later, Evans rewrote the theory for a classical system, setting the foundations of cDFT which describes the response of a uid to an external perturbation. 35 MDFT is designed to study solvation problems. The uid perturbed by the presence of one or several solutes is described by its density eld. In the following we will always consider liquids, referred to as the solvent. It is advantageous to dene a new functional F, as the difference between the grand potential functional of the solvent in the presence of the solute Q and that of the homogeneous solvent, With this denition, the functional at its minimum is equal to the solvation free energy. Because both the solute and the solvent are in most cases molecules with several atomic sites, their interactions depend on both the position of the centres of mass (CM) and orientations. Hence, the solvent density will be denoted by r(r,U) where r is the position in cartesian coordinates and U the orientation described by three Euler angles (q,f,j). The density minimizing the functional is the full equilibrium solvent density around the solute which may be integrated to recover the usual radial distribution functions.
The usual strategy to have a workable expression of this functional is to split it into the sum of ideal, excess and external contributions, In eqn (2) the ideal term corresponds to the entropy of the non-interacting uid, which reads The second term in eqn (2) account for the perturbation by the solute. The solute acts on the solvent via an "external" potential V ext , typically the sum of a Lennard-Jones term and of an electrostatic term. In the usual case of pair-wise additive interaction, V ext reads where v ij is the pair potential between site i of solvent and site j of the solute and r iU denotes the position of site i when the solvent molecule has the orientation U. The expression of the external functional is: Finally, the last term corresponds to the solvent-solvent interactions. An exact expression for this term is available 46 but in practice it is rewritten as the sum of two terms where F exc is an approximation of the excess functional. This denes the correction or "bridge" functional F b as the difference between the exact functional and this approximation. To date, the most advanced expression of F exc for water is that recently used by Ding and coworkers which corresponds to the hypernetted chain approximation. 43 where b ¼ (k B T) À1 and c(kr 1 À r 2 k,U 1 ,U 2 ) is the bulk direct correlation function while Dr ¼ r À r H . c depends on the distance between 2 solvent molecules and on their relative orientation dened by six Euler angles instead of 4 in the case of acetonitrile. 39 The calculation of the functional of eqn (7) in this general case is still feasible thanks to the efficient FFT algorithm 47 to handle the spatial convolution and to the use of rotational invariants to handle the angular one. [48][49][50] We have previously developed several approximations for the bridge term F b 51-53 but we shall neglect it in this paper and refer to this functional formulation as the HNC functional. While the numerical problem involved in MD is the sampling of phase space, MDFT involves an optimization which is numerically more efficient. Consequently, to compute solvation free energies of a spherical solute in water, MDFT requires z10 cpu minutes with our lab-developed program while it requires z100 cpu hours to compute the same quantity with commercially available MD codes.
1.3 ET reaction in the MDFT framework 1.3.1 Theory. From a MDFT perspective, the two states 0 and 1 of the ET reaction correspond to two functionals F 0 and F 1 differing only by their external potentials V 0 and V 1 in eqn (5). If we denote by r 0 and r 1 the equilibrium solvent densities of states 0 and 1, obtained by minimization of F 0 and F 1 , the reaction free energy can be expressed as The rst two terms represent the solvent contribution to the free energy while DE 0 is the difference in energy between the 2 solutes in vacuum. In this paper we restrict the study to rigid classical solutes with no intramolecular potentials, so that last term vanishes.
In his original work, Marcus estimated the free energy cost to solvate a solute within a solvent where polarization is not in equilibrium with the electric eld generated by the solute. In the MDFT framework the density eld contains all the structural equilibrium information of the solvent, including its polarization. We can consider MDFT as a more general eld theory than that used by Marcus. Nevertheless, the density eld itself remains a complicated object. To facilitate our understanding it is useful to dene a solvent reaction coordinate i.e. a scalar quantity which is uniquely dened by the density eld. By introducing a class of intermediate potentials interpolating between state 0 and state 1, we show that the average vertical energy gap is an appropriate order parameter. We then derive an expression for the free energies of states 0 and 1 as functions of the average vertical energy gap.
States 0 and 1 are characterized by the following Hamiltonian where h ¼ 0 or 1, K is the kinetic energy and U is the potential energy. For each state, the equilibrium probability distribution in the Grand Canonical ensemble is where m is the chemical potential of the solvent and (X N ,p N ) is a point in phase space X denotes the couple (r,U) describing the CM and orientation of a solvent molecule with momentum p. X h is the corresponding grand partition function: where Tr denotes the classical trace and h is the Planck constant. We now introduce a class of external potentials dened as linear combinations of V 0 and V 1 This denes the corresponding set of Hamiltonians (eqn (9)), probability distributions (eqn (10)) and grand partition functions (eqn (11)) for any value of h. Since for physically relevant cases V 0 and V 1 differ by more than a constant, any value of h denes a unique potential V h (up to an irrelevant constant). Because of this uniqueness of the potential, a unique equilibrium solvent density r h is associated with any value of h. This is a consequence of the cDFT principle 35 which implies a one-to-one mapping between external potential, equilibrium distribution and equilibrium solvent density. † We dene the average vertical energy gap, related to an equilibrium density r h by This quantity represents the energy difference between states 1 and 0 solvated in the solvent of density r h . We prove in Appendix A that, for the family of potentials in eqn (13), hDEi h is an adequate order parameter since it uniquely denes r h . Thus, the free energy of any state is a function of hDEi h . For instance, for state 0 it reduces to Note that the average free energy gap dened in eqn (14) differs from the microscopic version, DE, used in MD; with {R} denoting the whole set of coordinates of solvent molecules, but they are actually related by where on the le hand side h..i h denotes the thermodynamic average on the potential energy surface h. hDEi h is also frequently reported in MD studies of ET since it is another measure of the validity of Marcus theory which predicts that it varies linearly with the coupling parameter. Our approach is closer to Marcus' original work 54 where he mentioned that the "equivalent equilibrium distribution would be obtained in a corresponding equilibrium system in which the charges on the two central ions" are linear combinations of the original ones.
We can now express the reorganization free energies displayed in Fig. 1 as Borgis and coworkers have reported a similar relation to compute the reorganization free energies using MDFT. 39 Under the assumption that Marcus theory is valid -hence that the two reorganizations free energies are equal, i.e. l ¼ l 0 ¼ l 1 , eqn (18) and (19) reduce to: ð ð ½V 1 ðr; UÞ À V 0 ðr; UÞ½r 0 ðr; UÞ À r 1 ðr; UÞdrdU (20) which is equivalent to the linear response formula l ¼ 1 2 ðhDEi 0 À hDEi 1 Þ oen used in molecular simulations. We also note that, as for the usual variable DE in MD, the exact relation introduced by Warshel is satised: This is a corollary of eqn (34) with h ¼ 1. In the next subsection we explain how the average vertical energy gap, reorganization free energies and free energy curves are computed using MDFT.
1.3.2 Computational details. To study a given system, we shall minimize functionals corresponding to different external potentials V h according to eqn (13). We consider only cases for which the Lennard-Jones sites of the solute remain unchanged between state 0 and state 1, so that the energy gap reduces to the difference in the electrostatic potential energy of the solute in the eld generated by the solvent molecules. This can be computed using the electrostatic potential generated by states 0 and 1, while the vertical energy gap can be computed using eqn (14).
As shown above, the free energies corresponding to these values of the energy gaps are with A ¼ 0,1. To construct the FEC as in Fig. 1 we rst minimize the functional of eqn (2) for several values of h to obtain r h , next compute the value of the average vertical energy gap, and nally evaluate F 0 and F 1 for the different r h . An alternative route to compute the FEC was previously proposed by Hirata et al. 55,56 using another implicit solvent method, RISM. We show in Appendix C that the thermodynamic cycle they propose is equivalent to the present scheme, although not expressed in a free energy density functional language. Now that we have shown how MDFT can be used to investigate ET reactions, the following section is dedicated to assess the validity of this approach on simple and complex solutes.

ET between Cl 0 , Cl + and Cl À ions
In this article we focus on the difficult case of aqueous solvation, but calculation for simpler solvents such as acetonitrile or CO 2 are expected to give results of comparable quality. We apply the necessary correction due to periodic boundary conditions to charged solutes 57,58 and an additional correction accounting for the overestimation of the pressure within HNC 53 to both neutral and charged solutes.
To allow comparisons, we chose a system which has been extensively studied using MD by Hartnig et al. 19 This model of chlorine consists of one Lennard-Jones site, with s ¼ 4.404Å and 3 ¼ 0.4190 kJ mol À1 , and a charge equal to À1, 0 or 1 elementary charges e. To compute the FEC of the atom and the 2 ions with a good accuracy we ran MDFT calculations with a solute charge varying in steps of 0.1 elementary charges. We used a 40 Â 40 Â 40Å 3 box with 120 3 spatial grid points and 196 possible orientations per spatial point. The solvent is SPC/E water for which the exact direct correlation function projected on a basis of rotational invariants was obtained by Belloni et al. using a hybrid Monte Carlo plus integral equation approach. [59][60][61] All simulations are carried out at 298.15 K.
The FEC are shown in Fig. 2 which compares the MDFT results (solid lines) to the MD results 19 (symbols). The representation adopted here differs from that used by Hartnig and Koper, since we report the FEC as a function of the vertical energy gap and do not apply an arbitrary vertical shi of the curves. The methodology used to plot the MD data in this representation is described in Appendix D. This representation is better suited to highlight some features of the ET. For instance, we note that both pairs of FEC cross when the average free energy gap is equal to 0, as expected.
The agreement between MD and MDFT is satisfactory, the main difference being a shi of the MDFT curves towards negative values of hDEi h for the cation. The most interesting observation from Fig. 2 is the consistency between the curvature obtained by MD and by MDFT. In particular, Cl 0 and Cl + exhibit a similar curvature while that of Cl À is larger, indicating that the neutral to anion ET does not follow the Marcus picture while the neutral to cation ET does. To be more quantitative, we computed the reorganization free energy associated with the three species, based on the curvature by tting the data within 90 kJ mol À1 from the minimum with the following expression a strategy that was adopted by Hartnig et al. This assumes that the curve can be tted by a parabola. The expression linking the curvature and the parabola parameter in eqn (23) is derived using Marcus theory. The reorganization free energies obtained via MDFT and MD are compared in Table 1. For the MD data, we report the values of the free energies given in the original work, in addition to the one we have recomputed to keep the tting procedure consistent between the two approaches. MDFT overestimates the reorganization free energies compared to MD. However, comparing the values of the reorganization free energies between species leads to conclusions similar to those concerning the curvature of the FEC. The neutral atom and the cation have a similar reorganization free energy, while the anion has a much larger one.
The simple picture emerging from the comparison of curvatures is however misleading. Table 1 suggests that a single reorganization free energy can be associated with each solute, but this does not hold for several reasons: (i) it assumes that the FEC of a solute is a parabola (ii) it neglects the other solutes involved in the ET. We should refer to the reorganization free energy for a given 0 / 1 ET reaction as dened in Fig. 1, because the meaningful physical quantity is a free energy, not the curvature of a tting curve. To illustrate this point, we report in Table 2 the reorganization free energies computed using eqn (18) and (19). These free energies for the Cl 0 / Cl + reaction are almost identical to those reported in Table 1. This is consistent with the Marcus picture: if the two FEC are identical parabolas, there is a unique l parameter dening the curvature and the two free energy differences. However, in the Cl 0 / Cl À ET reaction the reorganization free energy of state 0 is much larger than in the other ET reaction. This is a consequence of the larger curvature of the second state. The reorganization free energy of Cl À is signicantly reduced compared to those listed in Table 1, another consequence of the smaller curvature of state 0.
This difference between the reorganization free energies is a further indication that the 0 / À1 ET does not follow Marcus theory.
Another way to check if the ET reaction follows the Marcus picture is to consider the evolution of the average vertical energy gap hDEi h with the parameter h. As mentioned earlier, such a curve is linear in Marcus theory. The evolution of the average vertical energy gap with h is presented in Fig. 3. For the neutral to positive charge transfer, the vertical energy gap does vary linearly with the coupling parameter h. In contrast, for the neutral to anion ET a non-linear variation is observed indicating a deviation from Marcus theory.
Thus, the curvatures of the FEC, the values of the reorganization free energies and the variation of hDEi h with respect to h consistently indicate a different behavior for the two ET reactions. This has already been noticed by Hartnig et al., who rationalized this observation by arguing that while the distance between the solute and the oxygen of the rst solvation layer remains similar for all oxidation numbers, the hydrogen is much closer to the solute in the case of the anion. This causes a "shrinking" of the rst solvation shell in the case of Cl À which differs considerably from the solvation shell of the neutral and positive solutes. Such a difference in the solvation shells of the two species cannot be properly captured by linear response assumed in Marcus theory.
Since MDFT gives access to the solvent density, we can also investigate the solvation structure. We compute the solvent charge density where s(r,U) is the charge distribution at point r of a single solvent molecule located at the origin, with the orientation U Table 1 Reorganization free energies computed via eqn (23) using data below 90 kJ mol À1 from Fig. 2. For MD, we recomputed the reorganization free energy using points extracted from Fig. 4 where the sum runs over the solvent sites, d is the Dirac distribution, r iU is the position of site i and q i its charge. The spherically averaged one-dimensional solvent charge densities are reported in Fig. 4 for the 3 oxidation numbers as a function of the distance to the solute. For all solutes, we observe a zone of zero charge density for small values of r, i.e. close to the solute, corresponding to the absence of water molecules. Then, alternating regions of positive and negative charge point to a preferential orientation of the solvent in the solvation shells. Finally, zero charge density is reached far from the solute at large r when a bulk behavior without preferential orientation is recovered.
If we rst consider the neutral and positive solutes, we observe in Fig. 4 that the preferential orientation of water in the rst solvation shell reverses between the neutral solute and the cation. Around the cation, the water molecules in the rst solvation shell have their oxygen pointing toward the solute. For the neutral species the hydrogen of water molecules are the closest to the solute. However the positions of the rst extrema are similar, 5.9Å for Cl 0 , 6.3Å for Cl + . This indicates that the two solvation shells essentially differ by the orientation of the water molecules. On the contrary, solvent molecules are much closer to the anion, where the rst maximum originating from hydrogen is located at 4.3Å. There is a shrinking of the rst solvation shell for the anion, in agreement with Hartnig et al. 19 It is conrmed by the comparison of the partial molar volume computed thanks to the equilibrium densities, namely 60Å 3 for Cl 0 and Cl + and 6Å 3 for Cl À . This difference in the solvation shell explains why Marcus theory fails to describe this ET reaction.

Solid/solvent interface
We now turn to the study of the inuence of a solid/solvent interface on the ET reaction. There are only few such studies available due to the computational cost of MD which is to date the only simulation tool used in this context. It is worth mentioning the investigations by Remsing et al. 62 who used MD and by Li et al. 18 based on coarse grained MD. In the former, ions are highly conned between two MnO 2 sheets and connement is kept constant throughout the study. In the ESI of Li's paper, the authors report the evolution of the reorganization free energy when the ion moves towards graphite sheets. They used umbrella sampling to constrain the position of the redox active site in the direction z perpendicular to the surface but no constraint was applied on the lateral coordinates. The dependence of the reorganization free-energy on the distance between solute and electrode was subsequently obtained through binning in the z-direction. In this set-up the position of the solute is not frozen but may uctuate around value of z under consideration, taking all possible values in x and y. The reported reorganization free energy is hence a statistical average.
The computational efficiency of MDFT allows a systematic study of the evolution of the reorganization free energy when the solute carrying the charge approaches an atomistically resolved wall. Because the solute is kept xed in the MDFT calculation it is not necessary to resort to biasing techniques to constrain its position and no uctuations blur the reorganization free energies. We consider the Cl 0 / Cl + ET with the parameters introduced in Section 2.1 and study the inuence of the proximity of a wall made up of 400 atoms arranged as the (100) surface of a fcc crystal. The size of the wall is 40 Â 40Å 2 and the distance between neighbouring atoms is 2Å. Each atom is modeled by a Lennard-Jones site with parameters s ¼ 3.37Å and 3 ¼ 0.23 kJ mol À1 similar to that used to model graphite atoms in previous studies. 63 To study separately the effect of the solvent on the ET we remove direct interactions between the solute and the wall. We used a 40 Â 40 Â 40Å 3 cubic box with 3 grid points perÅ and 196 discrete orientations per grid point.
We move the solute along the z axis perpendicular to the surface as illustrated in Fig. 5, with 175 calculations from z ¼ 2.5Å to z ¼ 20Å in steps of dz ¼ 0.1Å. The reorganization free energies of the Cl 0 / Cl + ET computed using eqn (18) are displayed in panel a) of Fig. 6 in solid black for the charged solute and in dashed red for the neutral solute. The two curves are similar and differ by less than 3 kJ mol À1 . This is a small difference consistent with the result of Li et al. 18 who reported that the ET of an iron atom dissolved in an ionic liquid next to a polarizable planar electrode follows Marcus' scenario.
We observe a decrease in the reorganization free energy as the solute approaches the plane. We can rationalize this observation by realizing that the wall truncates the solvation shell around the solute. This effect is illustrated in Fig. 7 which shows slices of the density prole around the neutral (le column) and charged (right column) solutes for different values of z. As the solute approaches the wall, there are fewer solvent molecules to rearrange when passing from one equilibrium solvation state to the other. This reduces the cost of the reorganization and explains the decrease of the free energy curves for small z. In the limit of total connement the reorganization free energy would vanish. The upper panel of Fig. 6 shows that reorganization free energy of the charged solute exhibits a maximum around 5.5Å. We rationalize this effect by decomposing l 1 into its components F 1 [r 1 ] and F 1 [r 0 ] in the lower panel of Fig. 6. While F 1 [r 0 ] exhibits a marked maximum around 5Å, F 1 [r 1 ] has a maximum around 3.4Å which is atter. Their difference consequently gives rise to the oscillatory behavior of l 1 around 5.5Å.
The rst solvation shell of the neutral solute at 5.5Å is in contact with the rst uid layer adsorbed on the wall. When the solute gets closer, the solvation shell is reduced. This reduces the unfavorable electrostatic term. It also decreases the cavity term which measures the cost of expelling the solvent from the region around the solute. This explains the maximum of F 1 [r 0 ]. Considering the density for the charged solute, the "contact" between the solvation shell of the solute and the uid layer adsorbed on the wall is also found around 5.5Å. However, for the cation the truncation of the solvation shell decreases both the favorable electrostatic term and the unfavorable cavity term. This could explain why F 1 [r 1 ] is rather at compared to F 0 [r 1 ] and why the position of the maximum is shied to the le. For z > 10Å, the reorganization free energies reach a plateau corresponding to the bulk value of Table 2.
The inset of panel (a) shows the difference between l 1 and l 0 for the for the Cl 0 / Cl + ET, it presents a maximum at 5.5Å, i.e.
where the reorganization free energies also have a maximum. The difference never exceeds 6 kJ mol À1 , so that Marcus' hypothesis is satised at all distances. This is clearly not the case for the Cl 0 / Cl À ET. Indeed, far from the wall l 0 and l 1 differ by 34 kJ mol À1 as reported in Table 2. However when the solute approaches the wall, this difference is reduced, indicating a decrease of the deviation the from linear response approximation. Again, this can be rationalized by the truncation of the solvation shell.  To illustrate the numerical efficiency of the method, we also computed the FEC for various positions: z ¼ 3.0Å is in the region where the reorganization free energy decreases, z ¼ 5.5Å Fig. 7 Slices of the solvent density for various values of the distance z from the wall. The neutral solute data are shown is in the left-hand column, while the cation data are shown in the right-hand column. and z ¼ 6.7Å correspond to the rst maximum and subsequent local minimum in Fig. 6. The FEC for the atom and cation are presented in Fig. 8. Each pair of curves cross at a point of vanishing vertical energy gap, as expected. When the solute gets close to the wall the minimum of the cation FEC is shied towards positive values, which is consistent with the abovementioned truncation of the solvation shell. Finally, the parabolas corresponding to z ¼ 3.0Å are wider than those for higher values of z, which is consistent with the smaller value of l close to the wall reported in Fig. 6.
One of the advantages of MDFT is the possibility to split the free energy into entropic, solute-solvent and solvent-solvent contributions according to eqn (2): Fig. 9 shows the various contributions to the reorganization free energy for the neutral and the charged solutes. To our knowledge, this is the rst time that such a decomposition of the reorganization free energy is reported.
A rst conclusion emerging from eqn (3), (6), (26) and (27) is that the ideal and excess contributions are exactly opposite for the neutral and the charged solutes. For both solutes, the ideal term due to the entropic contribution remains quite small and hardly varies with the distance from the electrode.
For the neutral solute the external contribution is small due to the absence of electrostatic interactions and more than 80% of the reorganization free energy is due to the excess term, i.e. the solvent-solvent contribution. In contrast, for the charged solute the main contribution is due to the electrostatic interaction between the solute and the solvent, which is roughly twice in absolute value than the solvent-solvent term. Even if we already know from the previous subsection that the Cl 0 / Cl + transfer does satisfy the linear response approximation, i.e. l 0 ¼ l 1 , it is fascinating to observe the compensation of the three contributing terms resulting in this equality. When the solutes approaches the wall, the linear response approximation gets even better as evidenced in Fig. 6 where the curves of l 0 and l 1 converge. This study also illustrates the interest of MDFT not only to compute the relevant free energies, but also to understand the various contributions to the free energy.

Conclusion
Marcus theory plays a crucial role in the study of ET reactions. This explains why its validity has been investigated extensively using molecular dynamics simulation. However, MD remains computationally very demanding, and has so far been essentially limited to simple systems. Molecular density functional theory has been proposed as an alternative to study solvation because it is computationally much faster, while retaining a molecular description of the solvent. In the present paper, we develop tools to use MDFT to study electron transfer reactions in water using MDFT. We have rst derived how to compute the relevant reaction coordinate: the average vertical energy gap. We Fig. 9 Reorganization free energies and their different contributions for Cl (0) (a) and Cl + (b) as a function of the distance from the wall computed as computed within MDFT. The reorganization free energy is in black, the ideal term is in green, the excess term in blue and the external term in red. Fig. 10 Schematic representation of the thermodynamic cycle used to recover eqn (40). 0 corresponds to the state under consideration, while h corresponds to a fictitious solute which interacts with the solvent via an external potential V h . The solvation states in equilibrium with states 0 and h are respectively denoted by S 0 and S h . To compute the FEC we need to compute DF S0/Sh 0 have also shown how to compute the free energy curves and the reorganization free energies.
We examined the validity of the approach by studying simple solutes, namely the ET reactions between Cl 0 , Cl À and Cl + modeled by a single Lennard Jones site and a point charge. We found a good agreement between the results obtained by MDFT and corresponding MD simulations. We conrmed the effect reported by Hartnig et al., that the ET between neutral and positive solutes is well described by Marcus theory, but not in the case of the transfer between the neutral atom and the anion.
We nally illustrated the potentiality of the method by tackling a more challenging system. We investigated the effect of the presence of a solid/solvent interface on the reorganization free energy, using a model system composed of an atomistically resolved neutral wall which is approached by the solute along the axis perpendicular to the wall. We computed the reorganization free energy for both neutral and charged states and found that they exhibit similar features. The reorganization free energy remains constant when the solute is far from the wall. As it approaches the wall, it exhibits oscillations before decreasing. We rationalized this behavior by considering the evolution of the solvation shell: close to the wall, there is less solvent to reorganize in the rst solvation shell, thereby reducing the free energy cost.
This work is a rst step towards the study of ET reaction in water and at electrode/water interfaces based on MDFT. The solvent effect sometimes called outer-sphere contribution is not the only mechanism playing a role in the ET reaction. The rearrangement of the electron cloud of the solute entering the so-called inner-sphere contribution may also play an important role. This effect is well taken into account in QM/MM calculation. There are mainly two approaches to deal with the MM part in such calculations. The rst one is to use MD, which takes into account the molecular nature of the solvent, but remains computationally costly. The second one is to use PCM-like models in which the solvent is described as a dielectric continuum. This approach neglects the molecular nature of solvent. As a consequence, it always assumes the validity of the linear response approximation and cannot properly describe systems violating Marcus theory. The strength of this method is its numerical efficiency: calculations are almost instantaneous. MDFT is thus a promising alternative to those two approaches to account for solvation in QM calculations: even if it is computationally more demanding than PCM its computational cost remains negligible compared to the cost of the QM calculation while its accuracy is comparable to MD. To that end, we are currently working on coupling MDFT with electronic structure calculations such as electronic density functional theory. We also wish to develop a framework allowing for the description of the polarizability of the wall to describe electrodes at xed electrode potential and study electrochemical reactions. These two objective are currently under investigation, in a attempt to develop a computationally efficient MDFT toolbox to tackle ET reactions.
A Proof that there is a one-to-one mapping between r h and hDEi h A straightforward consequence of eqn (14) is that the average vertical energy gap is uniquely dened by the density eld. Following Mermin and Evans, 34,35 we proceed by reductio ad absurdum to show that the average vertical energy gap uniquely determines the external potential and thus the density. Let us assume there exist two potentials V h and V h 0 with h s h 0 giving rise to the same gap i.e. hDEi h 0 ¼ hDEi h . From the expression of the probability distribution in eqn (10) and as stated in Appendix 1 of Evans's article, 35 From the variational principle of the grand potential we have By inverting the primed and unprimed quantities we get If we now sum eqn (28) and (29) we arrive at The integral on the r.h.s of eqn (30) vanishes as a consequence of eqn (14) and the assumption that hDEi h 0 ¼ hDEi h leading to a contradiction. Consequently, for this family of external potential V h there is a unique hDEi h which corresponds to a given probability distribution f h . We hence have a one to one mapping between all the following quantities where 4 denotes a one-to-one mapping. The bijections between the three quantities V, f and r are always true in the cDFT formalism 35 while the one involving h and V h is true within the class of potentials we have chosen. Because there is a bijection between hDEi h and a probability distribution, then the free energy of any state uniquely depends on hDEi h . To express the free energy as a function of hDEi h , it is sufficient to take advantage of the one to one mapping between r h and hDEi h , to obtain the expression of eqn (15).
It is worth noticing that we can actually dene F 0 (hDEi h ) as the Legendre transform of Q h with respect to h, as hDEi h is the conjugate variable of h: Moreover, Therefore the second derivative of Q h is negative, so that Q h is convex and its Legendre transform exists. We can hence dene the Legendre transform of F h by where we split F h into the sum of its three components as in eqn (2) and use the expression of V h in eqn (13). While the linear parametrization of the external potential is the only one allowing to dene F 0 (hDEi h ) as a Legendre transform, any parametrization leads to the same expression for F 0 (hDEi h ) and to the same FEC as demonstrated in Appendix B.
B Two different parameterizations of V h lead to the same F(hDEi h ) Let us consider the general parametrization for the interpolating potential, where s is a strictly increasing continuous function with s(0) ¼ 0 and s( Again the d and g indexes can be interchanged to show a oneto-one mapping between a value of the coupling parameter, the external potential, the equilibrium probability distribution and the equilibrium density. However this mapping now depends on the chosen parametrization s, The mapping between hDEi s h and r s h leads to: This relation does not depend on the choice of the parametrization s, but the values of hDEi s h and F 0 [r s h ] do. We now show that any parametrization yields the same FEC. Let us consider a strictly increasing continuous function with s(0) ¼ 0 and s(1) ¼ 1. The intermediate value theorem guarantees that s takes all the value between 0 and 1, once only. This is true for all s and in particular for the identity function corresponding to the linear parametrization. This last property implies that for all s(h)˛[0,1], there exists a unique a˛[0,1] such that Since the potential uniquely denes the functional this implies the same equality between all the properties in eqn (37), i.e.
This completes the proof that any parametrization of the intermediate potential leads to the same FEC.

C Thermodynamic cycle proposed by Chong and Hirata
Chong and Hirata proposed the thermodynamic cycle displayed in Fig. 10 where 0 and h are solutes corresponding to external potentials V 0 and V h in eqn (13). 55 The objective is to nd the free energy cost to modify the equilibrium solvent conguration around 0 into a solvent conguration which would be in equilibrium with h, a quantity denoted by DF S0/Sh 0 . Starting from 0 in vacuum, it is transformed into h spending a work W u h . Then, h is solvated in its equilibrium solvent conguration S h . This step corresponds to the solvation free energy DF h . The ctitious solute is transformed into 0 while the solvent conguration is frozen within S h . The free energy cost of this step DF Sh h/0 can be split into the sum of two terms. The rst one is the reversible work to transform the solute in vacuum: it is the opposite of W u h . The second term is the work W v h to transform the solute against the eld created by the solvent conguration S h which can be expressed using our previous notation as: The nal quantity required to close the cycle is the free energy cost to solvate 0 into its equilibrium solvent conguration. This corresponds to the solvation free energy of state 0, DF 0 . By closing the cycle, Chong and Hirata obtain the following formula: If we replace the solvation free energy by the functional of the present work, eqn (40) becomes which is equivalent to our previous nding.
D Modification of Hartnig and Koper's data to plot Fig. 2 In their paper, Hartnig and Koper represented their free energy curves as functions of a generalized order parameter dened as the electrostatic interaction energy between a negative point charge at the site of the solute and the solvent molecules. 19 Because they only considered solutes with a single Lennard-Jones site which is kept unchanged during the ET the vertical energy gap is equal to their order parameter for an anion and to its opposite for a cation. They do not mention the use of any nite size effect corrections while we use that proposed by Hünenberger et al. 57,58 which applies to our case but also to MD simulations with Ewald electrostatics. 57 To convert their order parameter in vertical energy gap we thus (i) multiply it by À1 in the case of the cation, (ii) apply the above mentioned electrostatic corrections with the box length parameter of L ¼ 24.83Å reported in their paper.
Finally, Hartnig and Koper shied all the FEC such that the minimum of each curve is equal to 0. As a consequence their curves do not cross for hDEi h ¼ 0, as should be the case by denition. Because they did not report the values of the solvation free energy, or equivalently the values of the shi applied to each curve, we decided to freely shi vertically the curves corresponding to the atom in order to have the minimum of the curves agree with the value predicted by MDFT. The MD curve for the ion has been subsequently shied vertically to fulll the zero gap condition.

Conflicts of interest
There are no conicts to declare.