Extension of Marcus Picture for Electron Transfer ... - ACS Publications

Article pubs.acs.org/JACS

Extension of Marcus Picture for Electron Transfer Reactions with Large Solvation Changes Rodolphe Vuilleumier, Kafui A. Tay,† Guillaume Jeanmairet, Daniel Borgis,* and Anne Boutin Ecole Normale Supérieure, Département de Chimie, UMR 8640 ENS-CNRS-UPMC, 24 Rue Lhomond, 75005 Paris, France ABSTRACT: The standard Marcus theory of charge transfer reaction in solution, relying on a Gaussian solvation picture, or, equivalently, on a linear response approximation, and involving two parameters, the reorganization energy and the reaction free-energy parameter, may fail when the solvation has a different character in the reactant and product state. We propose two complementary theoretical extensions of Marcus theory applying to those cases, based either on a two-Gaussian-states solvation picture, or on a nonGaussian solvation picture. As illustration, we show that such situations arise even for simple half oxido-reduction reactions involving the Cu+/Cu2+ or Ag0/Ag+ couples, for which electron transfer free-energy surfaces have been generated using first-principle molecular dynamics simulations. The two theoretical extensions are shown to exhibit the correct nonlinear response behavior and to reproduce the simulation results quantitatively, whereas a simple one-Gaussian-state Marcus description breaks down.

I. INTRODUCTION The Marcus theory of charge transfer reactions in solution has provided a very simple two-chemical state picture, based on two intersecting parabolas, that has made it possible to understand the experimental data, to interpret them quantitatively, and to make predictions.1 In the early versions, the solvent was modeled by either a dielectric continuum or a harmonic phonon bath.1−4 Starting with the pioneering work of Warshel in the early 1980s,5 the Marcus theory has fostered the computational modeling and molecular interpretation of those reactions using force field molecular dynamics (MD) and, more recently, first principle molecular dynamics (FPMD).6,7 One of the main contributions of Warshel was to exhibit the so-called energy-gap coordinate, ΔE, as the relevant microscopic reaction coordinate, as suggested originally by Marcus,8 and to show that, to a very good approximation, this quantity obeys Gaussian statistics.9 When translated to free energies, this property gives rise to Marcus’ two parabola picture. Since the study of Kuharski et al.10 for ferrous−ferric ion exchange in water, this Gaussian property has been verified many times, either for electron or proton transfers in solution or in complex environment such as proteins,6,11−18 and with force-field or first principle MD.19−23 It also has been recognized that the strict Gaussian assumption is equivalent to a linear response approximation.24,25 It implies identical solvent fluctuations in the two states and thus identical free-energy curvatures for those states.26 Departure from linear response can be estimated through free-energy integration methods, using ΔE as the relevant reaction coordinate, and it even can be used to design an efficient algorithm to compute the reaction free-energy and reaction barrier.24 It has been noticed in the literature, although this fact is sometimes overlooked, that the Gaussian assumption (or linear © 2011 American Chemical Society

response) is, by essence, incompatible with the existence of different solvation states in reactants and products and, thus, of different curvatures of the free-energy wells, or two different Gaussian widths.26 Such a situation can occur when the solvations of the species have different characters in reactants and products, as postulated theoretically by Kakitani and Mataga27−29 and illustrated by Carter and Hynes for the MD simulation of a neutral to ionic pair internal conversion in a polar solvent: The computed free-energy curvature corresponding to solvent fluctuations around the neutral pair or the ion pair do differ by a factor 1.6.30 More recently, Blumberger encountered a similar deviation from the linear response regime underlying the Marcus theory in the case of the Cu2+/Cu+ oxido-reduction reaction, which he studied using FPMD simulations.31 This nonlinear effect was attributed to the chemical specificity of copper ions with respect to water and a strong coordination change from Cu+ to Cu2+, going from a dihydrate to a 5-fold distorted pyramidal structure. We will show in this paper a similar breakdown of Marcus oxidation theory for the Ag0/Ag+ redox couple, where the reaction goes from a neutral species to a charged species and involves a drastic change in the solvation structure and dynamics. Note that copper and silver in various oxidation states have been studied by FPMD at several other occasions.32−34 The purpose of the present paper is to present extensions of the Marcus theory to cases such as those described above, where the solvation properties are markedly different in the reactant and product chemical states. We thus propose a global and coherent picture that reconciles the observation of two Gaussians of different widths with Marcus theory. We point out Received: July 24, 2011 Published: December 7, 2011 2067

dx.doi.org/10.1021/ja2069104 | J. Am. Chem.Soc. 2012, 134, 2067−2074

Journal of the American Chemical Society

Article

that a clear distinction should be made between a chemical state (with given electronic distribution) and a solvation state (with given solvent structural and fluctuations properties). Our first extension will thus be a direct generalization of Marcus theory based on a two chemical state/two solvation states picture rather than the two-chemical-state/one-solvation-state picture underlying Marcus theory. The second, separate, although complementary extension will be a non-Gaussian theory. Note that there have been several extensions of Marcus theory proposed in the literature, either relying on limiting cases of the harmonic bath Kubo−Toyozawa Hamiltonian35,36 or focused on corrections to the dielectric continuum solvent model.37 We propose here a statistical mechanics approach, involving a minimum number of relevant physical parameters, that is not committed to a specific underlying Hamiltonian or to a (generalized) dielectric solvent description. Although of general scope for charge transfer reactions in solution, our theoretical considerations will be focused on oxidoreduction half reactions, such as Cu+ → Cu2+ + e− and Ag0 → Ag+ + e− and they will be applied to those reactions in bulk water. For this purpose, we will use the published FPMD results of Blumberger for Cu+/Cu2+ in water and an original set of data that we have generated with a similar method for Ag0/ Ag+. Our approach should be valid also for the case of electron transfer full reactions. This is a challenging situation for FPMD, because of the presence of both donor and acceptor in the same simulation box, thus demanding careful treatment of selfinteraction.38,39 This paper is organized as follows: section II describes the theoretical foundation of Marcus theory, formulated in a molecular rather than dielectric continuum solvent framework, and its extensions. The theory is compared to FPMD freeenergy calculations in section III, and a conclusion is proposed in section IV.

pη (ε) = ⟨δ(ΔE({RI}) − ε)⟩η ,

(4)

The logarithm of this probability distribution determines, upon a constant, the Landau free-energy profile for the energy gap ε:

Wη(ε) = A η − kBT ln pη (ε)

(5)

where the constant Aη is the full free energy of the state η. Because the integral of the probability is unity, the Landau free energy satisfies the following:

A η = − kBT ln

∫ dε exp(−βWη(ε))

(6)

with β = 1/kBT. It can be shown that the Landau free energies of the reduced and oxidized states are related by9

W1(ε) = W0(ε) + ε

(7)

This is an exact and key relation: The difference between the reduced and oxidized Landau free energies is exactly the vertical ionization energy. Marcus also made the fundamental assumption that the probability distributions (pη) are Gaussians

⎛ (ε − ΔE )2 ⎞ 1 η ⎟ pη (ε) = exp⎜⎜ − ⎟ 4kBT λ η ⎠ 4πkBT λ η ⎝

(8)

and the Landau free energies are then parabolas

Wη(ε) = A η +

(ε − ΔE η)2 4λ η

k T + B ln(4πkBT λ η) 2

(9)

where we have introduced ΔEη, the equilibrium value of the energy gap, and λη, the reorganization energy, in the product (η = 1) or reactant (η = 0) states. Because of eq 7, the curvatures of W1 and W0 must be equal, that is,

λ1 = λ 0 = λ

(10)

Furthermore, the following two relations among ΔE1, ΔE0, λ, and ΔA hold

II. THEORY The model discussed here can be applied to a general charge transfer reaction of the form (1) D− + A → D + A −

λ=

(D and A being possibly charged molecular species). To fix ideas and stick to our applications, we consider an oxidoreduction half reaction: (2) R → O + e−

1 (ΔE0 − ΔE1) 2

ΔA =

(11)

1 (ΔE0 + ΔE1) 2

(12)

which also justifies calling λ = ΔE0 − ΔA the reorganization energy, the difference between the vertical ionization energy at equilibrium and the free energy difference between the product and reactant states. If we use λ and ΔA as parameters, the Landau free energies for the reactant and product state are written as follows:

where R and O are the reduced and oxidized species. If, for the sake of generality, we denote the reactant and product states as 0 and 1, we will keep in mind the correspondence R ≡ 0 and O ≡ 1. We start by recalling the basis of Marcus theory, formulated with the microscopic energy gap variable ΔE instead of the macroscopic solvent polarization variable used originally by Marcus.1 The derivation is standard but makes it possible to introduce the various quantities of interest, and the relations to be generalized. Marcus’s Theory: A Gaussian Solvation Model. We first introduce the energy gap coordinate, here also the vertical ionization energy, as the difference between the oxidized and reduced state energies for a given solvent configuration, RI (denoting the position of all atoms):

ΔE({RI }) = E1({RI }) − E0({RI })

η = 0 or 1

W0(ε) = A 0 +

k T (ε − λ − ΔA)2 + B ln(4πkBT λ) 4λ 2

W1(ε) = A 0 + ΔA +

(13)

k T (ε + λ − ΔA)2 + B 4λ 2

ln(4πkBT λ)

(14)

The free energy difference ΔA = A1 − A0 can be defined (and computed for a given molecular model) through a thermodynamic integration formula:

(3)

In the “molecular” Marcus theory presented here, it plays the crucial role of the order parameter. The probability distribution of observing a given value ε of the vertical ionization energy in the reduced and oxidized states is expressed as

ΔA =

2068

1

∫0

dη

dE η({RI}) dη

= η

1

∫0

dη⟨ΔE({RI})⟩η (15)



Article

where the subscripted brackets indicate a canonical average with the Hamiltonian H(η) defined by the intermediate PES

pη (ε) = (

E η({RI }) = (1 − η)E0({RI}) + ηE1({RI})

0 I

× δ(ΔE({RI}) − ε) +

= E0({RI }) + ηΔE({RI })

(16)

Wη(ε) = W0(ε) + ηε

/( +

so that, from the Gaussian approximation defined above, we have

∫S ∏ dRI exp(−βEη({RI})) ∫S ∏ dRI exp(−βEη({RI}))) 1 I

k T (ε − λ + 2ηλ − ΔA)2 Wη(ε) = A 0 + + B 4λ 2

(23)

where the indexes S0 and S1 of the integrals indicate integration over the phase space regions associated to the two solvation structures. Introducing

(18)

ΔA η = ηΔA + ηλ − η2 λ

1 I

0 I

(17)

ln(4πkBT λ) + ΔA η

∫S ∏ dRI

exp( −βE η({RI})) × δ(ΔE({RI }) − ε))

Defining, as in eq 4, the probability on this intermediate PES defined for a value of η between 0 and 1, the associated Landau free energy, eq 5, also fulfills the exact relation

AS , η = − kBT ln

(19)

∫S ∏ dRI exp(−βEη({RI}) I

It follows that

(24)

the free energy associated to the solvation structure S and

⟨ΔE({RI })⟩η = ΔA + λ − 2ηλ

(20)

WS , η(ε) = AS , η − kBT ln pS , η (ε)

is a straight line with a slope of −2λ and is equal to ΔA at η = 0.5. To connect to some arguments given in section I, we note that the latter formula amounts to a linear response approximation 2

⟨ΔE({RI })⟩η = ⟨ΔE({RI })⟩0 − ηβ⟨δΔE({RI }) ⟩0

= − kBT ln

0

∫S ∏ dRI exp(−βEη({RI})) I

× δ(ΔE({RI}) − ε)

(21)

(25)

the Landau free energy in solvation state S related to the conditional probability that the vertical ionization energy is equal to ε knowing that the solvation state is state S, we have, according to the Gaussian property postulated for each solvation state, and eq 18 and 19

Assuming from the start that such a linear response relation applies, and invoking the exact thermodynamic perturbation formula

⎛ ηΔE({RI }) ⎞ ΔA η = − kBT ln exp⎜ − ⎟ kBT ⎝ ⎠

∫S ∏ dRI exp(−βEη({RI}))

(22)

WS , η(ε) = AS ,0 +

it is easy to show that eq 21 boils down to an exact second-order cumulant expansion of the average, so that ΔE must be a Gaussian variable, and all the derivations above follow. Gaussian approximation or linear response are thus equivalent qualifiers.25 Note that extensions of the one-Gaussian approximation have been proposed very early by Marcus himself to incorporate molecular vibrations coupled to the electron transfer (and the so-called inner sphere)3,40 and, later on, for treating proton transfer reactions,4 reduction induced dissociation,41,42 or proton coupled electron transfer (PCET) reactions.43−45 Those extensions consist in introducing extra coordinates in addition to the solvent, generally Gaussian variables, that are coupled to the charge transfer and lead to reorganization energies that add to the solvent one. In this work, we rather focus on the ΔE collective coordinate itself and investigate how the (nonlinear) coupling to all other coordinates results in a nonGaussian behavior for ΔE, described below in two different ways. Extension to a Two-Gaussian Solvation (TGS) Model. We now assume that the solvent can be partitioned into two solvation states, labeled S0 and S1, with different fluctuation properties. For each of them, the Landau free energy, as a function of the vertical ionization energy, is supposed strictly Gaussian. The parameters describing the Landau free energy in each solvation state are AS,0, λS and ΔAS, S = S0,S1. The solvation states S0 and S1 are assumed to be favored in the reactant and product states, respectively. The free energy difference between the S1 and S0 solvation states is ΔSA0 = AS1,0 − AS0,0 when the system is in the reactant state. We should then have ΔSA0 > 0 for S0 to be indeed favored in the reactant state. ΔSA1 = ΔSA0 + ΔAS1 − ΔAS0 is then the free energy difference between the two solvation structures in the product state. It should be on the contrary negative, S1 being favored in the product state. Partitioning the solvent configuration space in two regions corresponding to the solvation states S0 and S1, the probability distribution pη(ε) for a vertical energy gap ε in the intermediate PES indexed by η, defined by eq 4, is given by

(ε − λS + 2ηλS − ΔAS )2 4λS

k T + B ln 4πkBT λS + ΔAS , η 2

ΔAS , η = ηΔAS + ηλS − η2 λS

(26) (27)

The vertical energy gap probability of eq 23 can be written as

pη (ε) =

exp(−βWS0 , η(ε)) + exp(−βWS1, η(ε)) exp(−βA S0 , η) + exp(−βA S1, η)

= exp(βA η) exp( −βWη(ε))

(28)

which defines the system free energy and the ε-dependent Landau free energy as a function of η

A η = − kBT ln(exp( −βA S0 , η) + exp( −βA S1, η))

(29)

Wη(ε) = − kBT ln(exp( −βWS0 , η(ε)) + exp( −βWS1, η(ε)))

(30)

It is easy to verify that the exact relation eq 17 for Wη(ε) is fulfilled because it is fulfilled separately for each state. The equilibrium vertical ionization energy, obtained as ∂Aη/∂η, is thus defined as follows:

⟨ΔE({RI })⟩η = [exp(−βA S0 , η)(ΔA S0 + λ S0 − 2ηλ S0) + exp( −βA S1, η) (ΔA S1 + λ S1 − 2ηλ S1)] /[exp( −βA S0 , η) + exp( −βA S1, η)] 2069

(31)



Article

using eq 27, which is the canonical average over the average vertical ionization energies in the two solvation states. When the solvation state S0 is favored (AS0,η ≪ AS1,η), the slope of ⟨ΔE({RI})⟩η as a function of η is −2λS0, and it is −2λS1 when the solvation state S1 is favored. Extension to a Non-Gaussian Solvation (NGS) Model. An alternative theoretical approach consists of departing from the Gaussian assumption of eq 8and 9 and postulating from the beginning a nonharmonic form for the Landau free energy in the reduced state 0

W0(ε) = A 0′ +

double-ζ valence plus polarization (DZVP) basis set, which was specially optimized for molecular systems.49 Core electrons have been replaced by the Goedecker−Teter−Hutter (GTH) norm-conserving pseudopotentials.50−52 The cutoff for the PW representation of the electronic density was set to 280 Ry. The wave function optimization was performed using an orbital transformation method.53 The gradient corrected exchangecorrelation functional BLYP54,55 was employed in the DFT calculations, in the local spin density framework. At each time step of the simulation, the electronic state of each oxidation state, Ag+ and Ag0, was optimized, and the forces on the ions were obtained as a linear combination of the Hellmann− Feynman forces for each state. The system consisted of one silver atom or cation in 64 water molecules with periodic boundary conditions. The box length for the cubic simulation box was 12.4085 Å. The time step for the MD simulation was 0.5 fs, and constant temperature conditions were imposed by a Nose−Hoover thermostat56,57 with a target temperature of 350 K. We generated six simulations corresponding to increasing values of η (η = 0, 0.2, 0.35, 0.5, 0.8, 1.0). Each was separated into 1 ps of equilibration and 5 ps of production at constant energy. Reversibility was checked for η = 0.35, the point where the transition between the solvation states occurs and where the energy gap fluctuations are the largest. The simulations make it possible to estimate the reaction free-energy, ΔA, through the computation of ⟨ΔE⟩η in each window and the integration formula, eq 15. The free-energy curves Wη(ε) and the associated probabilities pη(ε) can be reconstructed too, using weighted histograms techniques.58 The ⟨ΔE⟩η versus η curve obtained by simulation is displayed in Figure 1. The estimated statistical errors bars are 0.03−0.1 eV, that is, about the size of the symbols in the figure.

(ε − ΔE0)2 + am(ε − ΔE0)m 4λ S 0

+ an(ε − ΔE0)n

(32)

with n > m and ΔE0 = λS0 + ΔAS0. According to the considerations of eq 2, the Landau free-energy for the oxidized state is defined by W1(ε) = W0(ε) + ε. We impose that the corresponding curve has a minimum for ε = ΔE1 = −λS1 + ΔAS1, and the curvature at this point is 1/2λS1. These conditions impose that

am =

⎛ ⎞ (− 1)m − 1(n − 1) ⎜ λ(λ S0 + (n − 2)λ S1) ⎟ − 1 ⎜ ⎟ m(n − m)(2λ)m − 1 ⎝ (n − 1)λ S0λ S1 ⎠

(33)

and the symmetric expression for an, when permuting m and n. Here,

1 (ΔE0 − ΔE1) 2 1 = (λ S0 + λ S1 + ΔA S0 − ΔA S1) 2

λ=

(34)

is the linear-response solvent reorganization energy defined by eq 10. One does recover the harmonic case, an = am = 0, for λ = λS0 = λS1. The model thus involves four parameters, λS, and ΔAS for each state, which are reminiscent of (although not equivalent to) those in the previous TGS model. The energy parameter A0′ can be considered as an arbitrary free-energy shift. From there, the probability of a given vertical energy gap on an intermediate PES indexed by η, and the associated free energy and FES are defined by eqs 4 and 6 and can be determined from the knowledge of Wη(ε), which was itself deduced from W0(ε) through eq 17. Furthermore the equilibrium vertical ionization energy is readily computed as

⟨ΔE({RI })⟩η = exp(βA η)

+∞

∫−∞

d ε ε exp( −βWη(ε))

(35)

so that all quantities of interest are known through straightforward one-dimensional integrals over the energy gap variable ε. In the following, we limit ourselves to the simplest, generic approximation for the Landau free energy in eq 32, that is, a fourth-order polynomial with m = 3 and n = 4. Figure 1. Average vertical ionization energy for Ag0/Ag+ in aqueous solution, as a function of the coupling parameter η: CP2K simulation results (circles), Marcus theory (dashed line), Q-model of Matyushov and Voth35 (dotted line), TGS model (red line), and NGS model (blue dotted-dashed line).

III. RESULTS−DISCUSSION We begin by discussing the case of the Ag0 → Ag+ + e− reaction for which we could compute the η-dependent energy gap, ⟨ΔE({RI})η⟩, using first-principle molecular dynamics simulations (FPMD) and the methodology described in refs 19−23 and 31. The FPMD simulations used the Density Functional Theory (DFT) and the Born−Oppenheimer method. They were carried out using the freely available program package QUICKSTEP/CP2K,46 based on a hybrid Gaussian planewave (GPW) approach, which combines a Gaussian basis for the wave functions with an auxiliary plane wave (PW) basis set for the density.47 We choose a triple-ζ valence doubly polarized (TZV2P) basis set for oxygen and hydrogen atoms because it has been shown to provide a good compromise between accuracy and computational cost.48 For silver, we employed a

A highly nonlinear sigmoidal profile is obtained, which departs notably from the linear Marcus prediction, even more so than the previous result of Blumberger for Cu+/Cu2+31(see also Figure 3). In the same figure, we present the best fit corresponding to the TGS or NGS models: It can be seen that both are able to reproduce the nonlinear sigmoidal behavior with very good accuracy. The corresponding parameters are given in Table 1. As noted before, the parameters λS0,1 and ΔAS0,1 have similar significance in the two models, but they are not mathematically equivalent. Furthermore, the TGS model 2070



Article

0 + Table 1. Models Parameters for Agaq /Agaq Obtained by Fitting ⟨ΔE({RI})⟩η (in eV)

TGS NGS TGS using NGS parameters

λS0

ΔAS0

λS1

ΔAS1

ΔSA0

1.25 0.74 0.74

−0.61 −0.08 −0.08

0.63 0.41 0.41

−1.36 −1.58 −1.58

0.4 0.59

contains a fifth parameter, ΔSA0, that measures the free-energy difference between the two solvation states in the reduced charge state. For consistency, we have added to the table a best fit for ΔSA0, otherwise keeping the NGS parameters. Using the corresponding parameter set slightly alters the quality of the TGS curve in Figure 1 but keeps it fully consistent with the CP2K results. The applicability of the TGS model corroborates the physical relevance of the two-solvation-state picture in that case: one expects indeed different structural and dynamic properties of the solvent around a neutral atom such as Ag0 (presenting furthermore a strong excitonic character34) or around its charged cationic species Ag+. This assertion is substantiated in Figure 2 where the Ag0/Hw and Ag+/Hw radial

Figure 3. Average vertical ionization energy for Cu+/Cu2+ in aqueous solution as a function of the coupling parameter η: CPMD simulation results of ref 31 (circles), Marcus theory (dashed line), TGS model (red line), and NGS model (blue dotted-dashed line). + 2+ Table 2. Models Parameters for Cuaq /Cuaq Obtained by a Fitting ⟨ΔE({RI})⟩η (in eV)

TGS NGS TGS using NGS parameters a

λS0

ΔAS0

λS1

ΔAS1

ΔSA0

2.07 1.82 1.82

−0.24 0.21 0.21

1.02 0.81 0.81

−0.77 −1.0 −1.0

0.47 0.73

Taken from Ref 31.

interpreted by the author as the interplay between two Marcus linear curves intersecting close to η = 0.5 and corresponding to the transition between two solvent coordination regimes. In our interpretation suggested by application of both the TGS or NGS models, the two solvation-state picture is completely correct but manifests itself by a (slightly) sigmoidal curve, rather than by the intersection of two linear regimes. The question of the presence or not of an inflection point in the curve is subtle by simple examination of the raw data, and it is impaired by the simulation uncertainties. Note again that this feature is very marked and unambiguous in the case of Ag0/Ag+. Blumberger31 first suggested that the departure from linearity may be due to two very different solvation shells for Cu+ and Cu2+. The physical origin of the two solvation shells for Cu+ and Cu2+ is, however, different than for the case of Ag+/Ag0. The solvation shell of Cu+ has been previously studied.31−33 It was found that it forms a very strong Cu(H2O)2+ complex with linear geometry, further water molecules being pushed far away from the copper. This was attributed to a d−s hybridization of the copper orbitals.32 This hybridization does not occur, however, or not as strongly, in the case of Ag+, as a result of a larger energy gap between the d and s orbitals, as revealed from calculated electronic absorption spectra.33 We therefore attribute the very different solvation regime for Ag+/Ag0 to the neutral rather than charged character of Ag0, combined with its exceptional excitonic character generating a very loose hydration structure; see Figure 2 and the associated discussion. Note that, in contrast to the Cu2+/Cu+ or Ag+/Ag0, the Ag2+/ Ag+ couple behaves in very close accordance with the Marcus linear picture.32,33 To better picture the TGS model and understand the interplay between the TGS and NGS descriptions, we have plotted in Figure 4 the oxidized and reduced diabatic free energy curves W0,1(ε) for the Ag0/Ag+ reaction, using the same set of NGS-parameters (see Table 1) for consistency. Within the TGS framework, for a given reduced or oxidized charge

Figure 2. Ag/Hw radial distribution function for η = 0, that is, Ag0 (full black line), and η = 1, corresponding to Ag+ (dashed red line).

distribution functions computed from the FPMD simulations for η = 0 and η = 1, respectively, are displayed. The water structure around Ag+ is typical of a positive ion, with pronounced first and second shells around r = 3 Å and r = 4.5 Å, respectively, whereas the structure around Ag0 appears very smooth and floppy, even more so than for a hydrated rare gas atom. This floppiness can be attributed to the hyperpolarizability of Ag0, which can be regarded as a Ag+/e− system with a mobile electron interacting with the surrounding water.34 In Figure 1, in addition to the straight linear Marcus theory, we also give the results of the Q-model of electron transfer by Voth and Matyushov,35,36 in which deviation from the strict linear response behavior is allowed through the nonlinear coupling of the energy gap coordinate to an harmonic oscillator bath. The presented curve corresponds to the best fit of this three-parameter model to the simulation data. Even varying the parameters, we found it impossible to reproduce the observed sigmoidal shape. The physics introduced in the model, which was able to account for electron transfer reactions involving large polarizability changes,36 simply does not apply to the present case. Figure 3 presents similarly the ⟨ΔE⟩η results for the Cu+/ Cu2+ redox couple and the comparison to the TGS and NGS models (see the corresponding parameters in Table 2). The simulation data are those of Blumberger.31 They were 2071



Article

Figure 6. Diabatic free-energy curves for Ag0/Ag+ in aqueous solution: CP2K simulations (dots) and NGS model (dotted-dashed lines).

Figure 4. Solvent diabatic free-energy curves for the reduced and oxidized species in Ag0/Ag+: NGS model (thick dotted-dashed lines) and TGS model (thick line). The dashed and dotted thin lines correspond to the two solvation states in the TGS model.

state (η = 0 or 1, marked by colors), the solvation free-energy curves corresponding to the two different solvation states, WS0,η(ε) and WS1,η(ε), are illustrated in the figure by different line styles. The actual free-energy curves Wη(ε), defined by eq 28, exhibit a sharp transition between solvation states when the corresponding curves meet, and they follow the minimum freeenergy path. Concerning the comparison between TGS and NGS, apart from the energy tails below −2 eV or above 1 eV where a proper optimization of the TGS parameters would improve agreement between the approaches, the NGS curves do appear to interpolate smoothly and properly the discontinuous TGS curves in the intermediate region between minima, the one that matters for the description of the reaction. The two models thus carry a fully consistent picture. In Figures 5 and 6, we compare the same TGS and NGS diabatic free-energy curves for the Ag0/Ag+ reaction to the ones

Figure 7. Diabatic free-energy curves for Cu+ (blue) and Cu2+ (red) in aqueous solution: CPMD simulations of ref 31 (dots) and TGS model (lines). The dashed and dotted thin lines correspond to the two solvation states in the TGS model.

Figure 8. Diabatic free-energy curves for Cu+ (blue) and Cu2+ (red) in aqueous solution: CPMD simulations of ref 31 (dots) and NGS model (dotted-dashed lines).

Figure 5. Diabatic free-energy curves for Ag0/Ag+ in aqueous solution: CP2K simulations (dots) and TGS model (lines). The dashed and dotted thin lines correspond to the two solvation states in the TGS model.

Table 3. χ2 Values on Solvent Diabatic Free-Energy Curves for Parameters Fitted on Average Vertical Ionization Energy

obtained by FPMD, using η as the thermodynamic perturbation parameter and the WHAM technique to reconstruct the associated free energies.58 The same comparison is presented in Figure 7 and 8 for the Cu+/Cu2+ redox couple. The agreement is excellent in all cases. To be more quantitative, we report in Table 3 the χ2-values, measuring the mean squared distance between the simulation diabatic free-energy curves, obtained by weighted histograms, and the theoretical ones. It can be seen that the quality of the NGS and TGS results is quite close and

0 + /Agaq Agaq + 2+ Cuaq /Cuaq

parabolic

TGS

NGS

0.085 0.019

0.00020 0.0020

0.00019 0.0034

that, with respect to the linear response parabolic fit, both theories improve the χ2-values by 2 orders of magnitude or more, which is quite remarkable. In those calculations, the optimized parameters of Tables 1 and 2 were used for TGS. In Table 4, we also display for both models the values of the effective reorganization energies defined by λ0 = W0(ΔE1) − 2072



Article

oxidation reaction Cu+ → Cu2+ + e−, and we found it amplified in the case of Ag0/Ag+, with a transition between a rather labile coordination around a neutral, although extremely polarizable, atom to a more rigid 4-fold coordination around its charged cation. To account for such a situation, we have introduced a twoGaussian solvation state model that includes the correct physics by partitioning the solvent phase space into two independent solvation states and carrying out the corresponding statistical mechanical treatment based on the energy gap coordinate. The model was shown to reproduce the correct nonlinear response behavior, in particular, the sigmoidal shape of the curve ⟨ΔE⟩η versus η that is found in FPMD simulations, and that appears as a quite unexpected and “non-traditional” effect. Such agreement could not be obtained with the Q-model of Matyushov and Voth,35,36 which extends the Marcus picture beyond linear response but relies on different physical hypothesis that do not seem to apply here. We have also introduced a non-Gaussian solvation model that starts from a more phenomenological assumption, irrespective of any underlying Hamiltonian: given a nonharmonic, generic, polynomial form for the reactant free-energy curve, the peculiar statistical mechanics of electron transfer reactions, formulated in the microscopic energy gap coordinate, can be carried out. The theory relies on two extra parameters with respect to Marcus theory, and it was shown to be relevant for the reactions studied in this work, with a degree of precision similar to that of the two-Gaussian state model. Such a model is expected to remain valid for a wide range of physicochemical situations where linear response theory is possibly violated, leading to various nonlinear behaviors, in particular, for the curve ⟨ΔE⟩η versus η; this aspect will be rationalized in a forthcoming publication.

Table 4. Effective Reorganization Energies Computed with the One-Gaussian Parabolic Model or the TGS and NGS Models (in eV) 0 + Agaq /Agaq

parabolic TGS NGS

+ 2+ Cuaq /Cuaq

λ0

λ1

λ0

λ1

1.29 1.01 1.00

1.29 1.61 1.64

1.85 1.49 1.50

1.85 2.15 2.33

W0(ΔE0) and λ1 = W1(ΔE0) − W1(ΔE1) (to be clearly distinguished from the model parameters λS0 and λS1 of Tables 1 and 2). We also give the (unique) value corresponding to the parabolic linear response fit. Despite the different model parameters for TGS and NGS, it is very satisfactory to find the reorganization energies λ0 and λ1 in such a close agreement for the two models. Note also they are quite separate and on both sides of the parabolic value, indicating a strong violation of linear response. Finally, Table 5 compares for both reactions the values of the reaction free-energies obtained from eq 15, using either the Table 5. Reaction Free Energy for Oxidation Reactions Computed from First Principle Molecular Dynamics Simulation or the Theoretical Models of the Text (in eV) FPMD simulations TGS NGS

0 + Agaq /Agaq

+ 2+ Cuaq /Cuaq

−0.98 −0.97 −0.98

−0.32 −0.31 −0.32

FPMD results for ⟨ΔE⟩η or eqs 31 and 35 for the TGS and NGS models. Again, the numerical agreement appears excellent, either between the two models or compared with direct numerical integration. For the Ag0/Ag+ reaction, we can compute from the data of Table 5 the oxidoreduction potential of this couple with respect to the normal hydrogen electrode (NHE), using the reference shift of the Hartree potential for this setup59,60 and the absolute potential for the NHE, 4.28 V.61 We then obtain E(Ag0/Ag+) = −1.75 VNHE, taking 3.5 V as the shift of the Hartree potential.59 This value compares favorably with the experimental estimates E(Ag0/Ag+) = −1.9 and −2.1 VNHE.62 In this calculation, we neglect the finite size effects. These are estimated to be much smaller than these differences, as they were found to behave as R3/L3 where R is the solvation radius of the studied species and L the box size.59,63 Finite size effects have, however, a much larger influence on the reorganization free energies, which are presumably underestimated here as a result of a missing contribution from the reorganization of further solvation shells.64 These terms, however, are most likely linear and independent of the solvation state (being associated to far distance reorganization) and thus lead simply to a uniform increase in the slopes of both solvation states.

■

AUTHOR INFORMATION

Corresponding Author

[email protected] Present Address †

25 North Colonnade, Canary Wharf, London E14 5HS, United Kingdom

■ ■

ACKNOWLEDGMENTS We are thankful to James T. Hynes for very fruitful discussions. REFERENCES

(1) Marcus, R. A. J. Chem. Phys. 1956, 24, 956. (2) Levich, V. G.; Dogonadze, R. R. Dokl. Akad. Nauk SSSR 1959, 124, 123. (3) Newton, M. D.; Sutin, N. Annu. Rev. Phys. Chem. 1984, 35, 437. (4) Kuznetsov, A. M.; Ulstrup, J. Electron Transfer in Chemistry and Biology: An Introduction to the Theory; Wiley and Sons: Chichester, 1999. (5) Warshel, A. J. Phys. Chem. 1982, 86, 2218. (6) Warshel, A.; parson, W. Annu. Rev. Phys. Chem. 1991, 42, 279. (7) Van Voorhis, T.; Kowalczyk, T.; Kaduk, B.; Wang, L.-P.; Cheng, C.-L.; Wu, Q. Annu. Rev. Phys. Chem. 2010, 61, 149. (8) Marcus, R. A. Discuss. Faraday Soc. 1960, 29, 21. (9) Hwang, J. K.; Warshel, A. J. Am. Chem. Soc. 1987, 109, 715. (10) Kuharsky, R. A.; Bader, J. S.; Chandler, D. J. Chem. Phys. 1988, 89, 3248. (11) BorgisD.HynesJ. T.J. Chem. Phys.1990 (12) Azzouz, H.; Borgis, D. J. Chem. Phys. 1993, 98, 7361. (13) Marchi, M.; Gehlen, J.; Chandler, D.; Newton, M. J. Am. Chem. Soc. 1993, 115, 4178.

IV. CONCLUSIONS This paper attracts attention to the fact that the straight Marcus theory of charge transfer reaction in solution, relying on a Gaussian solvation picture, or, equivalently, on a linear response approximation, can be violated. This is the case when solvation has a quite different character in the reactant and product states, so that the solvent fluctuations are different. Such an effect was evidenced by FPMD simulations for the elementary 2073



Article

(62) Dubois, V.; Archirel, P.; Boutin, A. J. Phys. Chem. B 2001, 105, 9363. (63) Ayala, R.; Sprik, M. J. Phys. Chem. B 2008, 112, 257. (64) Seidel, R.; Faubel, M.; Winter, B.; Blumberger, J. J. Am. Chem. Soc. 2009, 131, 16127.

(14) Simonson, T. Proc. Natl. Acad. Sci. 2002, 99, 6544. (15) Sterpone, F.; Ceccarelli, M.; Marchi, M. J. Phys. Chem. B 2003, 107, 11208. (16) Ceccarelli, M.; Marchi, M. J. Phys. Chem. B 2003, 107, 5630. (17) Blumberger, J. Phys. Chem. Chem. Phys. 2008, 10, 5651. (18) Tipmanee, V.; Oberhofer, H.; Park, M.; Kim, K. S.; Blumberger, J. J. Am. Chem. Soc. 2010, 132, 17032. (19) Blumberger, J.; Bernasconi, L.; Tavernelli, I.; Vuilleumier, R.; Sprik, M. J. Am. Chem. Soc. 2004, 126, 3928. (20) Tateyama, Y.; Blumberger, J.; Sprik, M.; Tavernelli, I. J. Chem. Phys. 2005, 122, 234505. (21) Blumberger, J.; Sprik, M. J. Phys. Chem. B 2005, 109, 6793. (22) Blumberger, J.; Tateyama, Y.; Sprik, M. Comput. Phys. Commun. 2005, 169, 256. (23) VandeVondele, J.; Ayala, R.; Sulpizi, M.; Sprik, M. J. Electroanal. Chem. 2007, 607, 113. (24) Zhou, H.-X.; Szabo, A. J. Chem. Phys. 1995, 103, 3481. (25) Simonson, T.; Archantis, G.; Karplus, M. Acc. Chem. Res. 2002, 35, 430. (26) Tachiya, M. J. Phys. Chem. 1989, 93, 7050. (27) Kakitani, T.; Mataga, N. J. Phys. Chem. 1985, 89, 8. (28) Kakitani, T.; Mataga, N. J. Phys. Chem. 1986, 90, 993. (29) Kakitani, T.; Mataga, N. J. Phys. Chem. 1987, 91, 6277. (30) Carter, E. A.; Hynes, J. T. J. Phys. Chem. 1989, 93, 2184. (31) Blumberger, J. J. Am. Chem. Soc. 2008, 130, 16065. (32) Bernasconi, L.; Blumberger, J.; Sprik, M.; Vuilleumier, R. J. Chem. Phys. 2004, 121, 11885. (33) Blumberger, J.; Tavernelli, I.; Klein, M. L.; Sprik, M. J. Chem. Phys. 2006, 124, 064507. (34) Spezia, R.; Nicolas, C.; Boutin, A.; Vuilleumier, R. Phys. Rev. Lett. 2003, 91, 208304. (35) Matyushov, D. V.; Voth, G. A. J. Chem. Phys. 2000, 113, 5413. (36) Small, D. W.; Matyushov, D. V.; Voth, G. A. J. Am. Chem. Soc. 2003, 125, 7470. (37) Ichiye, T. J. Chem. Phys. 1996, 104, 7561. (38) Sit, P. H.-L.; Cococcioni, M.; Marzari, N. Phys. Rev. Lett. 2006, 97, 028303. (39) Oberhofer, H.; Blumberger, J. J. Chem. Phys. 2009, 131, 064101. (40) Marcus, R. A. J. Chem. Phys. 1965, 43, 679. (41) Savéant, J.-M. J. Am. Chem. Soc. 1987, 109, 6788. (42) Savéant, J.-M. Acc. Chem. Res. 1993, 26, 455. (43) Soudackov, A.; Hammes-Schiffer, S. J. Am. Chem. Soc. 1999, 121, 10598. (44) Soudackov, A.; Hammes-Schiffer, S. J. Chem. Phys. 2000, 113, 2385. (45) Hammes-Schiffer, S.; Stuchebrukhov, A. A. Chem. Rev. 2010, 110, 6939. (46) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Comput. Phys. Commun. 2005, 167, 103. (47) Lippert, G.; Hutter, J.; Parrinello, M. Mol. Phys. 1997, 92, 477. (48) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Chem. Phys. 2005, 122, 014515. (49) VandeVondele, J.; Hutter, J. J. Chem. Phys. 2007, 127, 114105. (50) Goedecker, S.; Teter, M.; Hutter, J. Phys. Rev. B 1996, 54, 1703. (51) Hartwigsen, C.; Goedecker, S.; Hutter, J. Phys. Rev. B 1998, 58, 3641. (52) Krack, M. Theor. Chem. Acc. 2005, 114, 145. (53) VandeVondele, J.; Hutter, J. J. Chem. Phys. 2003, 118, 4365. (54) Becke, A. Phys. Rev. A 1988, 38, 3098. (55) Lee, C.; Yang, W.; Parr, R. Phys. Rev. B 1988, 37, 785. (56) Nosé, S. J. Chem. Phys. 1984, 81, 511. (57) Nosé, S. Mol. Phys. 1984, 52, 255. (58) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1992, 13, 1011. (59) CostanzoF.SulpiziM.ValleR. G. D.SprikM.J. Chem. Phys.2011 (60) Kathmann, S. M.; Kuo, I.-F. W.; Mundy, C. J.; Schenter, G. K. J. Phys. Chem. B 2011, 115, 4369. (61) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2006, 110, 16066.

■

NOTE ADDED AFTER ASAP PUBLICATION Reference 3 was incorrect in the version of this Communication published ASAP January 10, 2012. The corrected version was posted January 17, 2012.

2074


Extension of Marcus Picture for Electron Transfer ... - ACS Publications

Recommend Documents