Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Guillaume
Jeanmairet
*^{ab},
Benjamin
Rotenberg
^{ab},
Maximilien
Levesque
^{c},
Daniel
Borgis
^{cd} and
Mathieu
Salanne
^{abd}
^{a}Sorbonne Université, CNRS, Physico-Chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, F-75005 Paris, France. E-mail: guillaume.jeanmairet@sorbonne-universite.fr
^{b}Réseau sur le Stockage Électrochimique de l'Énergie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France
^{c}PASTEUR, Département de chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 75005 Paris, France
^{d}Maison de la Simulation, CEA, CNRS, Université Paris-Sud, UVSQ, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

Received
10th October 2018
, Accepted 11th December 2018

First published on 12th December 2018

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.

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 fluctuations of this quantity are Gaussian. Such Gaussian statistics give rise to the famous parabola picture of Marcus for the free energy profile. A strict Gaussian behaviour is equivalent to a linear response of the solvent to the field generated by the solute; it also implies that the two free energy parabolas have the same curvature because the solvent fluctuations are identical for the two states.^{7,8}

This validity of the Gaussian assumption has been verified 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–16} and observed since in classical^{17–19} and ab initio^{20,21} simulations. Several extensions of Marcus theory have been put forward to take into account the various origins of non-linearity.^{21–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 configurations 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–28} This typically implies to run simulations on half a dozen fictitious 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 fluid in the presence of a perturbation by introducing a functional of the fluid density. Minimization of the functional yields the grand potential at equilibrium fluid 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 field. 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, after 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 first 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 influence 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.

Values of the order parameter differing from x_{0} correspond to solvent configurations that are not in equilibrium with state 0. The more the solvent configuration 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, ΔW. Two others key quantities appear in Fig. 1: the reorganization free energies λ_{0} (resp. λ_{1}) which represent the cost in free energy to solvate state 0 (resp. 1) in a solvent in equilibrium with the other 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 field 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 λ = λ_{0} = λ_{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.

MDFT is designed to study solvation problems. The fluid perturbed by the presence of one or several solutes is described by its density field. In the following we will always consider liquids, referred to as the solvent. It is advantageous to define a new functional F, as the difference between the grand potential functional of the solvent in the presence of the solute Θ and that of the homogeneous solvent,

F[ρ] = Θ[ρ] − Θ_{H}. | (1) |

With this definition, 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,Ω) where r is the position in cartesian coordinates and Ω the orientation described by three Euler angles (θ,ϕ,ψ). 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,

F[ρ] = F_{id}[ρ] + F_{ext}[ρ] + F^{exact}_{exc}[ρ]. | (2) |

In eqn (2) the ideal term corresponds to the entropy of the non-interacting fluid, which reads

(3) |

(4) |

(5) |

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

F^{exact}_{exc}[ρ] = F_{exc}[ρ] + F_{b}[ρ], | (6) |

(7) |

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 ≈10 cpu minutes with our lab-developed program while it requires ≈100 cpu hours to compute the same quantity with commercially available MD codes.

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 ρ_{0} and ρ_{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

where η = 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 μ is the chemical potential of the solvent and (X^{N},p^{N}) is a point in phase space X denotes the couple (r,Ω) describing the CM and orientation of a solvent molecule with momentum p. Ξ_{η} 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 defined as linear combinations of V_{0} and V_{1}

with {R} denoting the whole set of coordinates of solvent molecules, but they are actually related by

where on the left hand side 〈..〉_{η} denotes the thermodynamic average on the potential energy surface η. 〈ΔE〉_{η} 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.

which is equivalent to the linear response formula often used in molecular simulations. We also note that, as for the usual variable ΔE in MD, the exact relation introduced by Warshel is satisfied:

ΔW = F_{1}[ρ_{1}] − F_{0}[ρ_{0}] + ΔE_{0}. | (8) |

The first two terms represent the solvent contribution to the free energy while ΔE_{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 field generated by the solute. In the MDFT framework the density field contains all the structural equilibrium information of the solvent, including its polarization. We can consider MDFT as a more general field theory than that used by Marcus. Nevertheless, the density field itself remains a complicated object. To facilitate our understanding it is useful to define a solvent reaction coordinate i.e. a scalar quantity which is uniquely defined by the density field. 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

H_{η} = K + U + V_{η}. | (9) |

f_{η}(X^{N},p^{N}) = Ξ_{η}^{−1}exp[−β(H_{η}(X^{N},p^{N}) − μN)], | (10) |

Ξ_{η} = Tr[exp(−β(H_{η} − μN))], | (11) |

(12) |

V_{η} = V_{0} + η(V_{1} − V_{0}) with η ∈ [0,1]. | (13) |

This defines the corresponding set of Hamiltonians (eqn (9)), probability distributions (eqn (10)) and grand partition functions (eqn (11)) for any value of η. Since for physically relevant cases V_{0} and V_{1} differ by more than a constant, any value of η defines a unique potential V_{η} (up to an irrelevant constant). Because of this uniqueness of the potential, a unique equilibrium solvent density ρ_{η} is associated with any value of η. 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 define the average vertical energy gap, related to an equilibrium density ρ_{η} by

(14) |

This quantity represents the energy difference between states 1 and 0 solvated in the solvent of density ρ_{η}. We prove in Appendix A that, for the family of potentials in eqn (13), 〈ΔE〉_{η} is an adequate order parameter since it uniquely defines ρ_{η}. Thus, the free energy of any state is a function of 〈ΔE〉_{η}. For instance, for state 0 it reduces to

F_{0}(〈ΔE〉_{η}) ≡ F_{0}[ρ_{η}]. | (15) |

Note that the average free energy gap defined in eqn (14) differs from the microscopic version, ΔE, used in MD;

ΔE({R}) = E_{1}({R}) − E_{0}({R}), | (16) |

〈ΔE({R})〉_{η} = 〈ΔE〉_{η}, | (17) |

We can now express the reorganization free energies displayed in Fig. 1 as

λ_{0} = F_{0}(〈ΔE〉_{1}) − F_{0}(〈ΔE〉_{0}) = F_{0}[ρ_{1}] − F_{0}[ρ_{0}] = ΔW − 〈ΔE〉_{1}, | (18) |

λ_{1} = F_{1}(〈ΔE〉_{0}) − F_{1}(〈ΔE〉_{1}) = F_{1}[ρ_{0}] − F_{1}[ρ_{1}] = −ΔW + 〈ΔE〉_{2}. | (19) |

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. λ = λ_{0} = λ_{1}, eqn (18) and (19) reduce to:

(20) |

F_{1}(〈ΔE〉_{η}) = F_{0}(〈ΔE〉_{η}) + 〈ΔE〉_{η}. | (21) |

This is a corollary of eqn (34) with η = 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_{η} 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 field 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).

with A = 0,1. To construct the FEC as in Fig. 1 we first minimize the functional of eqn (2) for several values of η to obtain ρ_{η}, next compute the value of the average vertical energy gap, and finally evaluate F_{0} and F_{1} for the different ρ_{η}.

As shown above, the free energies corresponding to these values of the energy gaps are

F_{A}(〈ΔE〉_{η}) = F_{A}[ρ_{η}] | (22) |

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.

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 σ = 4.404 Å and ε = 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–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 shift 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.

Fig. 2 Pairs of free energy curves of (a) Cl^{0}/Cl^{−} and (b) Cl^{0}/Cl^{+} as a function of the average vertical energy gap. The black solid line and the dashed red line correspond to the MDFT results for the atom and the ions, respectively. Those results are compared to Hartnig's work^{19} reported as a function of the absolute vertical energy gap. The black circles correspond to the atom and the red squares to the ions. |

The agreement between MD and MDFT is satisfactory, the main difference being a shift of the MDFT curves towards negative values of 〈ΔE〉_{η} 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 fitting the data within 90 kJ mol^{−1} from the minimum with the following expression

(23) |

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 fitting 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.

Species |
λ
_{MDFT} (kJ mol^{−1}) |
λ
_{MD} (kJ mol^{−1}) |
---|---|---|

Cl^{0} |
233 | 153 (132) |

Cl^{−} |
297 | 263 (252) |

Cl^{+} |
216 | 165 (177) |

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 defined in Fig. 1, because the meaningful physical quantity is a free energy, not the curvature of a fitting 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 λ parameter defining 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 significantly reduced compared to those listed in Table 1, another consequence of the smaller curvature of state 0.

Species | Cl^{0}/Cl^{−} (kJ mol^{−1}) |
Cl^{0}/Cl^{+} (kJ mol^{−1}) |
---|---|---|

Cl^{0} |
297 | 214 |

Cl^{−} |
264 | N/A |

Cl^{+} |
N/A | 218 |

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 〈ΔE〉_{η} with the parameter η. As mentioned earlier, such a curve is linear in Marcus theory. The evolution of the average vertical energy gap with η is presented in Fig. 3. For the neutral to positive charge transfer, the vertical energy gap does vary linearly with the coupling parameter η. 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 〈ΔE〉_{η} with respect to η 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 first 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 first 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

(24) |

(25) |

If we first consider the neutral and positive solutes, we observe in Fig. 4 that the preferential orientation of water in the first solvation shell reverses between the neutral solute and the cation. Around the cation, the water molecules in the first 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 first 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 first maximum originating from hydrogen is located at 4.3 Å. There is a shrinking of the first solvation shell for the anion, in agreement with Hartnig et al.^{19} It is confirmed 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.

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 fixed in the MDFT calculation it is not necessary to resort to biasing techniques to constrain its position and no fluctuations blur the reorganization free energies. We consider the Cl^{0} → Cl^{+} ET with the parameters introduced in Section 2.1 and study the influence 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 σ = 3.37 Å and ε = 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.

Fig. 5 Snapshot of system under consideration: the flat wall is shown in grey, the solute in blue. The solute is moved along the z direction perpendicular to the wall. |

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 profile around the neutral (left 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 confinement 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 λ_{1} into its components F_{1}[ρ_{1}] and F_{1}[ρ_{0}] in the lower panel of Fig. 6. While F_{1}[ρ_{0}] exhibits a marked maximum around 5 Å, F_{1}[ρ_{1}] has a maximum around 3.4 Å which is flatter. Their difference consequently gives rise to the oscillatory behavior of λ_{1} around 5.5 Å.

The first solvation shell of the neutral solute at 5.5 Å is in contact with the first fluid 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}[ρ_{0}]. Considering the density for the charged solute, the “contact” between the solvation shell of the solute and the fluid 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}[ρ_{1}] is rather flat compared to F_{0}[ρ_{1}] and why the position of the maximum is shifted to the left. 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 λ_{1} and λ_{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 satisfied at all distances. This is clearly not the case for the Cl^{0} → Cl^{−} ET. Indeed, far from the wall λ_{0} and λ_{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 Å and z = 6.7 Å correspond to the first 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 shifted towards positive values, which is consistent with the above-mentioned 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 λ 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):

(26) |

(27) |

Fig. 9 shows the various contributions to the reorganization free energy for the neutral and the charged solutes. To our knowledge, this is the first time that such a decomposition of the reorganization free energy is reported.

A first 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. λ_{0} = λ_{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 λ_{0} and λ_{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.

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 confirmed 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 finally 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 first solvation shell, thereby reducing the free energy cost.

This work is a first 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 first 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 fixed 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.

(28) |

By inverting the primed and unprimed quantities we get

(29) |

If we now sum eqn (28) and (29) we arrive at

(30) |

The integral on the r.h.s of eqn (30) vanishes as a consequence of eqn (14) and the assumption that 〈ΔE〉_{η′} = 〈ΔE〉_{η} leading to a contradiction. Consequently, for this family of external potential V_{η} there is a unique 〈ΔE〉_{η} which corresponds to a given probability distribution f_{η}. We hence have a one to one mapping between all the following quantities

η ↔ V_{η} ↔ f_{η} ↔ ρ_{η} ↔ 〈ΔE〉_{η} | (31) |

Because there is a bijection between 〈ΔE〉_{η} and a probability distribution, then the free energy of any state uniquely depends on 〈ΔE〉_{η}. To express the free energy as a function of 〈ΔE〉_{η}, it is sufficient to take advantage of the one to one mapping between ρ_{η} and 〈ΔE〉_{η}, to obtain the expression of eqn (15).

It is worth noticing that we can actually define F_{0}(〈ΔE〉_{η}) as the Legendre transform of Θ_{η} with respect to η, as 〈ΔE〉_{η} is the conjugate variable of η:

(32) |

Moreover,

(33) |

Therefore the second derivative of Θ_{η} is negative, so that Θ_{η} is convex and its Legendre transform exists. We can hence define the Legendre transform of F_{η} by

(34) |

V^{s}_{η} = V_{0} + s(η)(V_{1} − V_{0}) | (35) |

(36) |

Again the δ and γ indexes can be interchanged to show a one-to-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,

(37) |

The mapping between 〈ΔE〉^{s}_{η} and ρ^{s}_{η} leads to:

F_{0}(〈ΔE〉^{s}_{η}) = F_{0}[ρ^{s}_{η}]. | (38) |

This relation does not depend on the choice of the parametrization s, but the values of 〈ΔE〉^{s}_{η} and F_{0}[ρ^{s}_{η}] 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(η) ∈ [0,1], there exists a unique α ∈ [0,1] such that

s(η) = α |

V^{s}_{η} = V_{α} |

Since the potential uniquely defines the functional this implies the same equality between all the properties in eqn (37), i.e.

f^{s}_{η} = f_{α} |

ρ^{s}_{η} = ρ_{α} |

〈ΔE〉^{s}_{η} = 〈ΔE〉_{α} |

F_{0}[〈ΔE〉^{s}_{η}] = F_{0}[〈ΔE〉_{α}] |

This completes the proof that any parametrization of the intermediate potential leads to the same FEC.

W^{v}_{η} = 〈V_{0} − V_{η}〉_{η}. | (39) |

Fig. 10 Schematic representation of the thermodynamic cycle used to recover eqn (40). 0 corresponds to the state under consideration, while η corresponds to a fictitious solute which interacts with the solvent via an external potential V_{η}. The solvation states in equilibrium with states 0 and η are respectively denoted by S_{0} and S_{η}. To compute the FEC we need to compute , the free energy cost to modify the solvent configuration around state 0 from S_{0} to S_{η}. |

The final quantity required to close the cycle is the free energy cost to solvate 0 into its equilibrium solvent configuration. This corresponds to the solvation free energy of state 0, ΔF_{0}. By closing the cycle, Chong and Hirata obtain the following formula:

(40) |

If we replace the solvation free energy by the functional of the present work, eqn (40) becomes

(41) |

They do not mention the use of any finite 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 shifted all the FEC such that the minimum of each curve is equal to 0. As a consequence their curves do not cross for 〈ΔE〉_{η} = 0, as should be the case by definition. Because they did not report the values of the solvation free energy, or equivalently the values of the shift applied to each curve, we decided to freely shift 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 shifted vertically to fulfill the zero gap condition.

- R. A. Marcus, J. Electroanal. Chem., 1997, 438, 251–259 CrossRef CAS .
- R. A. Marcus, J. Chem. Phys., 1956, 24, 979–989 CrossRef CAS .
- R. A. Marcus, J. Chem. Phys., 1956, 24, 966–978 CrossRef CAS .
- R. A. Marcus, J. Phys. Chem., 1990, 94, 1050–1055 CrossRef CAS .
- J. R. Miller, L. T. Calcaterra and G. L. Closs, J. Am. Chem. Soc., 1984, 106, 3047–3049 CrossRef CAS .
- J. K. Hwang and A. Warshel, J. Am. Chem. Soc., 1987, 109, 715–720 CrossRef CAS .
- M. Tachiya, J. Phys. Chem., 1989, 93, 7050–7052 CrossRef CAS .
- M. Tachiya, J. Phys. Chem., 1993, 97, 5911–5916 CrossRef CAS .
- R. A. Kuharski, J. S. Bader, D. Chandler, M. Sprik, M. L. Klein and R. W. Impey, J. Chem. Phys., 1988, 89, 3248–3257 CrossRef CAS .
- J. Blumberger and M. Sprik, Theor. Chem. Acc., 2006, 115, 113–126 Search PubMed .
- T. Simonson, Proc. Natl. Acad. Sci., 2002, 99, 6544–6549 CrossRef CAS PubMed .
- F. Sterpone, M. Ceccarelli and M. Marchi, J. Phys. Chem. B, 2003, 107, 11208–11215 CrossRef CAS .
- A. d. l. Lande, F. Cailliez and D. R. Salahub, in Simulating Enzyme Reactivity, 2016, pp. 89–149 Search PubMed .
- T. Kakitani and N. Mataga, J. Phys. Chem., 1985, 89, 8–10 CrossRef CAS .
- T. Kakitani and N. Mataga, J. Phys. Chem., 1986, 90, 993–995 CrossRef CAS .
- T. Kakitani and N. Mataga, J. Phys. Chem., 1987, 91, 6277–6285 CrossRef CAS .
- E. A. Carter and J. T. Hynes, J. Phys. Chem., 1989, 93, 2184–2187 CrossRef CAS .
- Z. Li, G. Jeanmairet, T. Méndez-Morales, M. Burbano, M. Haefele and M. Salanne, J. Phys. Chem. Lett., 2017, 1925–1931 CrossRef PubMed .
- C. Hartnig and M. T. M. Koper, J. Chem. Phys., 2001, 115, 8540–8546 CrossRef CAS .
- J. Blumberger, Phys. Chem. Chem. Phys., 2008, 10, 5651–5667 RSC .
- R. Vuilleumier, K. A. Tay, G. Jeanmairet, D. Borgis and A. Boutin, J. Am. Chem. Soc., 2012, 134, 2067–2074 CrossRef CAS PubMed .
- D. V. Matyushov and G. A. Voth, J. Chem. Phys., 2000, 113, 5413–5424 CrossRef CAS .
- D. W. Small, D. V. Matyushov and G. A. Voth, J. Am. Chem. Soc., 2003, 125, 7470–7478 CrossRef CAS PubMed .
- G. Jeanmairet, D. Borgis, A. Boutin and R. Vuilleumier, in CHAPTER 18: Extension of Marcus Rate Theory to Electron Transfer Reactions with Large Solvation Changes, 2013 Search PubMed .
- G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199 CrossRef .
- A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1989, 63, 1195–1198 CrossRef CAS PubMed .
- M. R. Shirts and J. D. Chodera, J. Chem. Phys., 2008, 129, 124105 CrossRef PubMed .
- Z. Tan, E. Gallicchio, M. Lapelosa and R. M. Levy, J. Chem. Phys., 2012, 136, 144102 CrossRef PubMed .
- P. H. Fries and G. N. Patey, J. Chem. Phys., 1985, 82, 429–440 CrossRef CAS .
- D. Chandler and H. C. Andersen, J. Chem. Phys., 1972, 57, 1930–1937 CrossRef CAS .
- F. Hirata and P. J. Rossky, Chem. Phys. Lett., 1981, 83, 329–334 CrossRef CAS .
- A. Kovalenko and F. Hirata, Chem. Phys. Lett., 1998, 290, 237–244 CrossRef CAS .
- T. Imai, A. Kovalenko and F. Hirata, Mol. Simul., 2006, 32, 817–824 CrossRef CAS .
- N. D. Mermin, Phys. Rev., 1965, 137, A1441–A1443 CrossRef .
- R. Evans, Adv. Phys., 1979, 28, 143 CrossRef CAS .
- R. Ramirez, R. Gebauer, M. Mareschal and D. Borgis, Phys. Rev. E, 2002, 66, 031206 CrossRef PubMed .
- R. Ramirez and D. Borgis, J. Phys. Chem. B, 2005, 109, 6754–6763 CrossRef CAS PubMed .
- M. Levesque, V. Marry, B. Rotenberg, G. Jeanmairet, R. Vuilleumier and D. Borgis, J. Chem. Phys., 2012, 137, 224107 CrossRef PubMed .
- D. Borgis, L. Gendre and R. Ramirez, J. Phys. Chem. B, 2012, 116, 2504–2512 CrossRef CAS PubMed .
- R. Ramirez, M. Mareschal and D. Borgis, Chem. Phys., 2005, 319, 261–272 CrossRef CAS .
- G. Jeanmairet, M. Levesque, R. Vuilleumier and D. Borgis, J. Phys. Chem. Lett., 2013, 4, 619–624 CrossRef CAS PubMed .
- G. Jeanmairet, N. Levy, M. Levesque and D. Borgis, J. Phys.: Condens. Matter, 2016, 28, 244005 CrossRef PubMed .
- L. Ding, M. Levesque, D. Borgis and L. Belloni, J. Chem. Phys., 2017, 147, 094107 CrossRef PubMed .
- S. Luukkonen, L. Belloni, D. Borgis and M. Levesque, 2018, arXiv:1806.03118 [physics].
- P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864 CrossRef .
- J.-P. Hansen and I. McDonald, Theory of Simple Liquids, Academic Press, 3rd edn, 2006 Search PubMed .
- M. Frigo and S. Johnson, Proc. IEEE, 2005, 93, 216–231 Search PubMed .
- L. Blum and A. J. Torruella, J. Chem. Phys., 1972, 56, 303–310 CrossRef CAS .
- L. Blum, J. Chem. Phys., 1972, 57, 1862–1869 CrossRef CAS .
- L. Blum, J. Chem. Phys., 1973, 58, 3295 CrossRef CAS .
- M. Levesque, R. Vuilleumier and D. Borgis, J. Chem. Phys., 2012, 137, 034115 CrossRef PubMed .
- G. Jeanmairet, M. Levesque and D. Borgis, J. Chem. Phys., 2013, 139, 154101 CrossRef PubMed .
- G. Jeanmairet, M. Levesque, V. Sergiievskyi and D. Borgis, J. Chem. Phys., 2015, 142, 154112 CrossRef PubMed .
- R. A. Marcus, Discuss. Faraday Soc., 1960, 29, 21–31 RSC .
- S.-H. Chong and F. Hirata, Mol. Simul., 1996, 16, 3–17 CrossRef CAS .
- H. Sato, Y. Kobori, S. Tero-Kubota and F. Hirata, J. Chem. Phys., 2003, 119, 2753–2760 CrossRef CAS .
- M. A. Kastenholz and P. H. Hünenberger, J. Chem. Phys., 2006, 124, 124106 CrossRef CAS PubMed .
- M. A. Kastenholz and P. H. Hünenberger, J. Chem. Phys., 2006, 124, 224501 CrossRef PubMed .
- J. Puibasset and L. Belloni, J. Chem. Phys., 2012, 136, 154503 CrossRef PubMed .
- L. Belloni and I. Chikina, Mol. Phys., 2014, 112, 1246–1256 CrossRef CAS .
- L. Belloni, J. Chem. Phys., 2017, 147, 164121 CrossRef PubMed .
- R. C. Remsing, I. G. McKendry, D. R. Strongin, M. L. Klein and M. J. Zdilla, J. Phys. Chem. Lett., 2015, 6, 4804–4808 CrossRef CAS PubMed .
- M. W. Cole and J. R. Klein, Surf. Sci., 1983, 124, 547–554 CrossRef CAS .

## Footnote |

† Note that this is true for any V_{η} = V_{0} + s(η) (V_{1} − V_{0}), as long as s is a strictly increasing continuous function with s(0) = 0 and s(1) = 1. |

This journal is © The Royal Society of Chemistry 2019 |