10746
J. Phys. Chem. 1996, 100, 10746-10753
Electron Transfer Across the Electrode/Electrolyte Interface: Influence of Redox Ion Mobility and Counterions August Calhoun and Gregory A. Voth* Department of Chemistry, UniVersity of PennsylVania, Philadelphia, PennsylVania 19104-6323 ReceiVed: February 27, 1996; In Final Form: April 16, 1996X
A microscopic model for electron transfer (ET) across an electrolyte/metal electrode interface is studied by molecular dynamics (MD) simulation. The Anderson-Newns Hamiltonian is used to describe the adiabatic ET for an Fe2+/Fe3+ ion in water at a Pt(111) electrode surface using an implicit representation of the electron transfer event. The diabatic and adiabatic solvent free energy curves for the electron transfer are calculated; however, the present study significantly expands upon our earlier simulations [Straus, J. B.; Calhoun, A.; Voth, G. A. J. Chem. Phys. 1995, 102, 529)]. In particular, the motion of the redox ion is now allowed, as well its corresponding image in the metal surface arising from the electrode polarization. Importantly, a counterion is also added to the simulations in order to study its effect on the ET free energy curves. In all cases, the effects are found to be significant.
I. Introduction Electron transfer (ET) processes lie at the heart of many chemical, biological, and technological processes. For example, the microscopic steps in, e.g., photosynthesis, redox reactions in batteries, and corrosion, hinge on the underlying ET dynamics. In this work, we report further progress in our ongoing effort to simulate heterogeneous ET at the microscopic level, specifically across the electrolyte/electrode (metal) interface.1 In our earlier molecular dynamics (MD) simulation,1 an aqueous ferric (Fe3+) ion was studied which was fixed 5.1 Å above a Pt(111) electrode surface. The electronic population of the lowest unoccupied atomic orbital of the ferric ion at the surface was described by the Anderson-Newns Hamiltonian.2,3 Accordingly, this Hamiltonian provides the adiabatic potential energy function through which the forces on the water nuclei can be derived for a molecular dynamics simulation. By virtue of our simulation approach, both the diabatic and adiabatic solvent free energy curves for the ET could be calculated, and it was demonstrated that the diabatic curves are reasonably parabolic as a function of the solvent reaction coordinate, ∆E. Furthermore an approximate solution4 to the Anderson-Newns Hamiltonian based on the parabolic assumption was found to adequately reproduce the exact adiabatic free energy curve (and hence the solvent reorganization barrier). Interestingly, a path integral MD simulation for the same system discovered that the quantization of the water hydrogen nuclei at room temperature significantly affects the barrier and thermodynamic driving force for the ET reaction.1 In the present paper, the Anderson-Newns-based MD methodology is extended to incorporate the effects of both a mobile redox ion and counterions. This extension of the model leads to a re-definition of the reaction coordinate for the ET event and allows for a quantification of the effects of redox ion mobility and counterion dynamics on the ET free energy curve. Moreover, the ion pair structure and dynamics can be studied within the context of the electron transfer reaction at the electrode, and this leads to some interesting perspectives for future research. The subsequent sections are organized as follows: In section II, the generalizations of the simulation model to include redox X
Abstract published in AdVance ACS Abstracts, June 1, 1996.
S0022-3654(96)00603-X CCC: $12.00
ion mobility and counterion effects are outlined. In section III, results from simulations are then presented and discussed which focus on these effects. Concluding remarks are given in section IV. II. Model A. Hamiltonian. The simulation of heterogeneous ET across the electrode/electrolyte interface represents a considerable challenge due to the many-electron nature of the problem, as well as the complexity of the interactions in such systems. As a first step, ET between a redox ion in solution and a metal electrode was simulated1 using the Anderson-Newns Hamiltonian2,3,5,6 as adapted for electrochemical systems.4,7-9 In these simulations, a single Fe3+ redox ion was studied near a Pt(111) electrode, and the ion was assumed to be fixed in position. In general, the Anderson-Newns Hamiltonian for a single redox ion in the electrochemical context can be written as
H ) Hsolution + Hel
(2.1)
where Hsolution is the Hamiltonian for the electrolyte solution (including the solvent/solvent, solvent/surface, and all other interactions which do not explicitly involve the redox ion orbital). Here, Hel is the“ET” part of the Hamiltonian and is given in second quantized form as
Hel ) (a + ∆E)na + ∑(knk + Vakc†ack + Vkac†k ca)
(2.2)
k
where a is the vacuum energy level of the redox ion electronic orbital |a〉 which is involved in the ET. The ion orbital |a〉 has associated with it a number operator na, which is equal to either 0 or 1, depending on its occupancy, and ∆E is the instantaneous energy shift of the ion electronic state |a〉 due to its interaction with the fluctuating electrochemical environment (i.e., it is the collective “solvent” coordinate). The precise microscopic definition of the latter variable depends on the system one is studying, and this will be discussed in detail below. The sum over k in eq 2.2 is over the electronic states in the semiinfinite metal electrode. The latter summation includes contributions arising from both the occupancy of the metal electronic states and interaction terms leading to the charge transfer between the metal states and the ion orbital (i.e., c†i and ci are electron © 1996 American Chemical Society
ET Across the Electrode/Electrolyte Interface
J. Phys. Chem., Vol. 100, No. 25, 1996 10747
creation and annihilation operators, respectively, and Vak is the electronic resonance integral which couples the redox ion to the metal electronic levels). In much of the following discussion, the Fe2+/Fe3+ heterogeneous ET at a Pt electrode studied in ref 1 is the paradigmatic case study. Within the adiabatic assumption of separability between the solvent and electronic degrees of freedom, one can use resolvent operator techniques to explicitly simplify the electronic Hamiltonian Hel in eq 2.2. In the case of a single redox ion, the resulting ground state electronic energy E0(∆E) can be analytically obtained as a function of four physically important quantities: (1) the Fermi level of the metal electrode f, (2) the vacuum energy level a of the redox ion orbital, (3) the solvation or redox ion level shift variable ∆E, and (4) the effective energy broadening ∆ of the ion orbital |a〉 due to its interaction with the metal electronic states. The explicit expression for E0(∆E) is given by
1 1 E0(∆E) ) ∆E + (a + ∆E - f) × 2 π f - (a + ∆E) ∆ + tan-1 ln[(a + ∆E - f)2 + ∆2] (2.3) ∆ 2π
[
]
where the level broadening parameter ∆ is given in this equation by
∆() ) π∑|Vak|2δ( - k)
(2.4)
k
In principle, the parameters ∆ and f depend on the instantaneous configurations of the electrolyte molecules and the redox ion. However, in ref 1 it was assumed that these were constant. In particular, the quantity ∆() is insensitive to the energy for energies in the range of interest and is evaluated at ) f for convenience. If, however, the redox ion is free to move as was not the case in ref 1, ∆ is not independent of the distance of the redox ion from the surface. This issue will be revisited later. The explicit specification of the redox ion level shift variable ∆E in the Hamiltonian requires some thought in order to carry out explicit MD simulations.1 In principle, one must assign ∆E to certain microscopic potential energy terms, e.g., the terms that depend on the Coulombic interaction between the solvent and the redox ion. In the case of the Fe2+/Fe3+ system, these terms can be written with respect to the ferric ion, Fe3+, with a positive core charge of Z ()3), and the occupancy of the electronic redox ion orbital is described by the number operator na. If the ion orbital is empty, then na ) 0, and the Hamiltonian effectively reduces to that of the Fe3+ ion solvated near the electrode surface. If the ion orbital is filled, however, then na ) 1, and the Hamiltonian effectively describes the solvated ferrous ion, Fe2+, near the surface. In the case of a fixed redox ion and a perfectly polarizable metal electrode, the overall Coulombic interaction involving the redox ion is found to be given by1
VCoulombic ) (Z - na)(Vsolvent-ion - Vsolvent-ion image)
(2.5)
Here, Vsolvent-ion indicates the potential energy between the solvent and a single positive charge at the ion’s position, while Vsolvent-ion image indicates the potential energy between the solvent and the charge of opposite sign at the position of the ion's image. For the fixed redox ion, in order to make the connection with the term ∆E from the adiabatic solution of the AndersonNewns Hamiltonian (eq 2.2), one specifies it by inspection to be1
∆E ) -(Vsolvent-ion - Vsolvent-ion image)
(2.6)
The Z -dependent ionic “core” terms arising from eq 2.5 can be included with the rest of the interactions in the solution Hamiltonian Hsolution. As shown in ref 1, once a reasonable microscopic potential has been specified that describes the interactions between the solvent molecules, the solvent and the electrode, and the ion and the solvent, one can carry out MD simulations based on the overall adiabatic Hamiltonian
H ) Hsolution + E0(∆E)
(2.7)
In the expression for E0(∆E) given in eq 2.3, the Fermi level f of the electrode, the redox ion orbital vacuum energy level a, and its broadening ∆ due to interaction with the electrode are specified as input parameters. In the MD simulations of a fixed redox ion based on this approach, the solvent “reaction coordinate” for the interfacial ET process is clearly given by ∆E in eq 2.6.1 A general comment is in order here regarding the population operator, na. It is noteworthy that na does not appear explicitly within the Hamiltonian that is used to run MD trajectories, eq 2.7. Instead, for the discussion of a MD simulation, a useful and illustrative quantity is the average occupancy of the orbital on the ion, 〈na(∆E)〉,5,6 rather than the population operator. This quantity varies continuously as a function of ∆E from 0 to 1, i.e.,
() (
)
f - (a + ∆E) 1 1 tan-1 〈na(∆E)〉 ) + 2 π ∆
(2.8)
Whereas the population operator is a binary quantity, the average orbital occupancy provides a more intuitive picture of the electron transfer within an MD simulation since it is a continuously varying quantity. The value of 〈na(∆E)〉 is also directly accessible from the value of ∆E at each molecular dynamics time step. Equation 2.8 stands as a useful summary of the resulting behavior of the Anderson-Newns Hamiltonian.1 B. Free Energy Curves. One of the most important quantities in heterogeneous ET systems that can be determined from the present simulation approach is the ET activation free energy and thermodynamic driving force. In general, the ET free energy curve is given up to an additive constant by
F(∆e) ) -kBT ln[∫dq δ[∆e - ∆E(q)]e-βV(q)] (2.9) where ∆e is a particular value of the collective energy gap coordinate ∆E(q), and the variable q denotes all microscopic system coordinates which contribute to the instantaneous value of ∆E. For simulations of the heterogeneous ET adiabatic free energy curves,1 V(q) is the total adiabatic potential energy function implicit in eq 2.7. On the other hand, for the reactant and product diabatic free energy curves, the ion with an empty and filled redox ion orbital |a〉, respectively, is simulated and the appropriate V(q) is used in eq 2.9.1,10-14 In either case, appropriately modified1,10 umbrella sampling techniques15,16 are used to ensure adequate sampling. If the diabatic free energy curves for the reaction coordinate ∆E are reasonably parabolic, then it can be shown that the adiabatic ET free energy curve F(∆e) is given by4,8
F(∆e) )
1 (∆e - ∆emin)2 + E0(∆e) 4λ
(2.10)
where λ is the diabatic reorganization free energy, ∆emin is the minimum of the diabatic free energy curve for the empty redox ion orbital state, and ∆e is the instantaneous value of the reaction
10748 J. Phys. Chem., Vol. 100, No. 25, 1996
Calhoun and Voth
coordinate, ∆E. In ref 1, the above form was found to fit the exact simulation data for the adiabatic free energy curves reasonably well. The Fermi level f can be treated as a continuously adjustable parameter in the simulation because one is free to adjust the overpotential η of an electrode. The Fermi level for which eq 2.10 will yield a symmetric double well (zero overpotential) can be determined by the relation4
η ) a - f + ∆emin - λ
(2.11)
C. Modifications for a Mobile Redox Ion. In the case of a mobile redox ion, the specification of the term ∆E in eq 2.2 becomes more complicated because the ion interacts dynamically with the electrode surface as well as with the solvent molecules and their images. In this case, one has the Coulombic interactions of the ion with the solvent, surface, and solvent images, given by
VCoulombic ) (Z - na)(Vsolvent-ion - Vsolvent-ion image Vsolvent images-ion + (Z - na)Vmetal-ion) (2.12) where Vsolvent-ion is the Coulombic interaction of the solvent with a singly charged ion, Vsolvent-ion image is the Coulombic interaction of the solvent with a singly charged image of the ion, Vsolvent images-ion is the Coulombic interaction of the images of the solvent with the ion, and Vmetal-ion is the Coulombic image-based interaction of the singly charged ion with the metal electrode surface (see eq 2.20). Equation 2.12 can be expanded, and the core ionic terms (those with only Z dependence) can be separated from those involving the redox ion orbital occupation number na. From this analysis, ∆E in the AndersonNewns Hamiltonian can be written by inspection as
∆E ) -(Vsolvent-ion - Vsolvent-ion image - Vsolvent images-ion + (2Z - 1)Vmetal-ion) (2.13) For the case of the fixed redox ion discussed earlier and in ref 1, only the electrostatic solvent fluctuations and their images composed the reaction coordinate (eq 2.6). From the above expression for ∆E, one can see that this is no longer true in the case of a mobile redox ion. Indeed, the redox ion level energy is also influenced by its interaction with the metal electrode and the solvent images. These terms can be significant in their own right, and the ionic motions can be coupled to the solvent polarization fluctuations in ways which are not present in the fixed redox ion case. The mobile redox ion scenario also introduces another important factor into the simulation. That is, the coupling between the redox ion level and the electrode electronic states is strongly distance dependent,7 so this motion can influence the adiabatic ET barrier through eq 2.3. To be more specific, the parameter ∆ is explicitly given by
∆(q, ) ) π∑|Vak(q)|2δ( - k)
(2.14)
k
where Vak(q) is now the distance-dependent matrix element which couples the ion orbital to the kth metal electronic state. The parameter ∆ can be calculated explicitly by ab initio methods, or, as in the present simulations, it can be treated empirically by the “universal” form,17
∆(zion) ) ∆0e-a(zion-z0)
(2.15)
where ∆(zion) is a function only of the vertical distance, zion, of the ion from the surface. In the present case, the parameters R
and z0 were chosen so that ∆(zion) matches the value used in our earlier work of 0.05 eV at 5.1 Å, with a contact value of 4.0 eV. In general, the value of ∆0 is not known without an explicit ab initio calculation, but the value used in the present work is thought to be within the reasonable range.1,7,18 From the above expression, one sees that fluctuations in the ion distance away from the surface will exponentially decrease the coupling term ∆(zion) which, in turn, increases the adiabatic ET barrier height. Within the context of variational transition state theory, one is free to pick other definitions of the reaction coordinate than the ones defined above. For example, in the spirit of our earlier calculations,1 one could pick a reaction coordinate involving just the solvent polarization fluctuations and their images as in eq 2.6. A third possibility for a reaction coordinate might be to include also the Vsolvent images-ion term in eq 2.6, i.e.,
∆E ) -(Vsolvent-ion - Vsolvent-ion image - Vsolvent images-ion) (2.16) The three possibilities for reaction coordinates discussed here will be referred to as the full, solvent and solvent + image reaction coordinates, respectively. Variational transition state theory states that the choice of reaction coordinate giving the highest activation barrier in the resulting adiabatic ET free energy curve, and thus the lowest ET rate, can be considered to be the “optimal” choice of reaction coordinate. This idea will be explored in section IIIC. Before proceeding to the next subsection, it should be noted that another approach for incorporating redox ion mobility could be employed (see, e.g., ref 14). In such an approach, one assumes an adiabatic separation of the redox ion motion from the electron transfer event. That is, one assumes the electron transfer rate between the electrode and the ion is always much faster than the motion of the ion, irrespective of the distance from the electrode. Within this assumption, one calculates an ET rate kET(zion) for fixed values of zion, as was done in ref 1 for the adiabatic limit or as in ref 14 for the nonadiabatic limit, and then averages that rate over the properly normalized equilibrium probability P(zion) for finding the ion a distance zion from the electrode surface, i.e.,
kET = ∫dzion Pion(zion) kET(zion)
(2.17)
In the present paper, however, this approach is not used because we are specifically interested in the influence of the redox ion motions on the adiabatic ET free energy curves and do not wish to invoke additional assumptions. D. Modifications for Counterions. Another step toward increasing the realism of the interfacial ET simulations is to add counterions to the simulation cell. For example, such an extension will provide insight into how the presence of counterions influences the adiabatic ET curves and hence the ET rate. The generalization of the Hamiltonian required to include the effect of counterions is more straightforward than that to include the mobility of the redox ion. The reaction coordinate simply inherits several extra terms with the addition of counterions due to the Coulombic interaction of the counterions with the redox ion, the image of the redox ion, and the image of the counterions with the redox ion. The reaction coordinate can thus be generalized to be
∆E ) -(Vsolvent-ion - Vsolvent-ion image - Vsolven timages-ion + (2Z - 1)Vmetal-ion + Vion-counter - Vion image-counter Vion-counter images) (2.18) The additional terms involving the interaction of the counterions
ET Across the Electrode/Electrolyte Interface
J. Phys. Chem., Vol. 100, No. 25, 1996 10749
with the water molecules, the redox ion core, and the electrode surface are included in Hsolution in eq 2.7. With this approach, one can consistently define an adiabatic free energy curve by sampling the values of the reaction coordinate as defined in eq 2.18. Since the relative motion between the redox ion and a counterion should be very slow compared to the fluctuations of the solvent (water) dipoles, one can adopt an adiabatic approach to define the overall ET rate. In this case, the ET rate kET(rrel) is calculated for a fixed relative distance rrel between the redox ion and counterion. Then, the overall ET rate kET is given by
kET = ∫drrel P(rrel) kET(rrel)
(2.19)
where P(rrel) is the properly normalized probability distribution for the redox ion and counterion to be a distance rrel apart. In the simulations described below, the value of kET(rrel) is reported for the Fe2+/Fe3+-Pt(111) electrode reaction with a Clcounterion at several fixed distances rrel. It should be noted that ion pairing is an experimental observation,19 so this is another clear reason why it is important to expand the heterogeneous ET investigations to include the effects of counterions. E. Details of the Simulation Model. The simulations for the single redox ion consisted of 671 water molecules in a box of dimensions 41.5 × 24.0 × 22.5 Å (x, y, z directions) with the Fe3+ ion core, periodic boundary conditions in the x and y directions, and the Anderson-Newns1 description of the electron transfer from the electrode to the redox ion orbital. The bottom of the box was bounded by a corrugated Pt(111) potential,20,21 while the top was bounded by a featureless slab.1,22 The SPC/F model was used to describe the water/water interactions, and it incorporates the internal motions of the water molecules through a harmonic intramolecular potential.23,24 The iron/water25 and water/platinum potentials20,21 were the same as in our previous work.1 A simulation temperature of 300 K was maintained by attaching a Nose´ thermostat to every particle.26,27 The polarizability of the metal electrons is extremely important for an ion in solution at a metal electrode surface, and these interactions were approximately treated using appropriate image charges.1 An image of the iron ion was thus included and the solvent allowed to interact with it. The iron ion also interacts with an image of itself in the electrode as one expects from a polarizable medium. This approach is reasonable since the redox ion in this case remains 5 Å or greater from the surface. The latter potential, appearing as Vmetal-ion in eqs 2.12, 2.13, and 2.18 includes screening effects from the surface electrons and is given by28
Vmetal-ion )
1 4(zion + se-3s/2ao)
(2.20)
where s2 ) πa0/4kf, kf is the Fermi wave vector, and a0 is a Bohr radius. The images of the solvent were included as well, and the iron ion interacted with these, but the solvent did not interact with its own image, as this interaction was implicitly included in the solvent/surface potential.20 The advantage of using that latter potential instead of a perfectly polarizable image plane is that is also includes the effects of surface corrugation in the water/electrode interactions. (The surface corrugation is less important for the ion/surface interaction since the ion remains over 5 Å from the surface.) The Coulombic potential and forces were smoothly cut off after 10.3 Å to a value of zero at 12.0 Å. It should be noted that the water interaction
Figure 1. Water oxygen density for the water/electrode system in the z direction. A well-defined double layer is evident at the platinum surface, with a bulklike density in the center of the box.
with the corrugated platinum potential results in an ice-like water structure20 and a well-defined double layer at the water/platinum interface (cf. Figure 1). As stated before, a chloride counterion was used in some of the simulations to study its effect on the ET free energy curve. This anion has been studied by others at a metal interface28,29 and as an ion pair with a sodium cation.30 The chloride ion interacts with the water solvent through the potential given by30
VH2O-Cl- )
A2 r12
-
C2 r6
+∑ i
qiqClri
(2.21)
where A2 × 10-3 ) 4050 (kcal Å12)/mol and C2 ) 1480 (kcal Å6)/mol.31 The iron ion interacted with the chloride ion through a normal Coulombic potential and the same repulsive term as in the iron/oxygen interaction for the Fe3+/H2O potential. (We note that at contact a chloride ion can actually form a bound complex with a ferric ion and may catalyze the electrochemical charge transfer process through a bridging mechanism.32 However, at the distances studied in the present paper the chloride ion is always constrained to be outside the first water solvent shell of the iron ion in order to focus on the purely electrostatic counterion effect on the ET free energy curve.) The chloride ion interacted with the platinum surface through an image potential of the same form as the ion/platinum image potential. The simulation with a chloride ion was carried out in a larger simulation cell of dimensions 41.6 × 33.6 × 21.0 Å (x, y, z) with 1008 water molecules to minimize additional cutoff force effects. The Coulombic forces were smoothly cut off after 15.0 Å to a value of zero at 16.5 Å. A note is in order regarding the errors associated with the free energy curves calculated from the simulations. The free energy as a function of the reaction coordinate was in error by approximately 0.6 kcal/mol, and hence the barrier heights and driving forces experience this error as well. The error was calculated by dividing the free energy averaging into 10 intervals, and the standard deviation of the free energy at each point on the curve was computed. The error in the reorganization energy, λ, is on the order of 0.8 kcal/mol, and this error was calculated from the parabolic least squares fit to the free energy curve and a propagation of errors. To explore the possibility of systematic error in the sampling, the free energy sampling windows for each different type of simulation were compared to cases where the simulation temperature was heated to 500 K for 10 ps, then cooled down, and re-equilibrated, and the free energy sampling carried out again. The results from these calculations were that they were statistically indistinguishable from those which did not employ the heating-cooling cycle.
10750 J. Phys. Chem., Vol. 100, No. 25, 1996
Figure 2. Diabatic free energy curves for the mobile Fe3+ ion. The solid line is the result as calculated from the simulation, while the dashed line is the parabolic fit. The ferrous curve (dotted line) was generated from the ferric curve via eq 3.1. There is substantial deviation from the parabolic behavior, as assumed in Marcus theory, away from the minimum.
III. Results and Discussion A. Ion Image Effects. An MD simulation was first carried out for the diabatic state of the system with the ferric ion, water solvent, electrode surface, and the appropriate images. For these simulations, the water solvent was first equilibrated with the iron ion fixed 5.1 Å above the surface and then the iron ion released to move freely within the simulation cell. The inclusion of the solvent images interacting with the iron ion was observed to affect the equilibrium height of the ion from the platinum surface. Without the solvent images, the equilibrium height of the ferric ion from the surface was approximately 5.0 Å, but with the solvent images included it was found to be approximately 5.5 Å. This result is expected since the solvent images add an additional screening effect. The physical environment of the ion was found to be very different depending upon the distance from the surface. The water oxygen density profile shown in Figure 1 indicates that near the surface the density of water molecules changes rapidly with distance from the surface. This result suggests that to have a proper modeling of the redox ion near the surface, the images of the solvent need to be included, otherwise the ion has a strong attraction to its own image screened only by the waters between the ion and the surface, and thus it stays too close to the surface. The inclusion of the solvent images further screens the image attraction between the ion and its image, leading to a larger equilibrium distance from the surface. B. Diabatic Free Energy CurvessMobile Redox Ion. The first step toward determining the adiabatic ET free energy is to determine the diabatic free energy curves, the equilibrium value of the reaction coordinate ∆E in those states, and the reorganization free energy. The curvature and the position of the minimum of the free energy as a function of the electronic reaction coordinate can next be used in the calculation of the overpotential from eq 2.11. The overpotential can then be adjusted as desired by varying the value of the Fermi level in the Anderson-Newns Hamiltonian. The diabatic state for the ferric ion was first equilibrated for 12.5 ps with no umbrella potential. Harmonic umbrella windows were then applied to the reaction coordinate,1 ∆E, and the free energy curve was mapped out with 60 such windows. Each window was equilibrated for 5 ps if the initial configuration was chosen from a neighboring window or for 12.5 ps if the initial configuration was not from a neighboring window. The reaction coordinate probability distribution was then sampled for 12.5 ps in each window. The diabatic free energy curves from the simulation are shown in Figure 2. The free energy sampling was carried out only for the ferric case, and the ferrous free energy curve was generated by the relation4
Calhoun and Voth
Figure 3. Adiabatic free energy curve for the Fe2+ f Fe3++ ereaction as a function of the full reaction coordinate in eq 2.13 as calculated from simulation (solid line) and by eq 2.10 (dashed line).
F2+ (∆e) ) F3+ (∆e) + ∆e
(3.1)
A parabolic fit to the diabatic free energy for the ferric ion is shown in Figure 2 in which the fitting was carried out around the minimum. In particular, the curve was fit to the parabolic form given by the first term on the right hand side of eq 2.10. The value of ∆emin was specified directly from the MD simulation, while the reorganization free energy λ was determined from the fit to the MD data in a symmetric fashion. By using this procedure, the minimum was found to be at 526.7 kcal/mol, and the reorganization free energy, λ, given by 56.0 kcal/mol. There is a substantial deviation from a parabolic approximation to the free energy away from the minimum.10 C. Adiabatic Free Energy CurvessMobile Redox Ion. The adiabatic free energy curves were next calculated with the full Hamiltonian in eq 2.7. A group of 21 umbrella sampling windows was used to calculate the free energy curve, each window being equilibrated independently for 12.5 ps, and free energy statistics were then collected for between 12.5 and 40 ps, depending on the convergence. The adiabatic free energy curve for the full reaction coordinate defined in eq 2.13 and a Fermi level of f ) -10.2 eV was calculated and is shown in Figure 3. The curve has the form of a nearly symmetric double well, the left well corresponding to the Fe2+ state and the right well corresponding to the Fe3+ state. Also shown is the predicted adiabatic curve using the analytic expression from eq 2.10 and the values of λ and ∆emin determined from the diabatic Fe3+ simulation. The agreement is semiquantitative, though there is a clear difference in the thermodynamic driving force for the (Fe2+ f Fe3++ e- ) direction. For the mobile redox ion case, the overall free energy curve is shifted by approximately 80 kcal/mol for the full reaction coordinate (eq 2.13) relative to our earlier result,1 having a fixed iron ion and the corresponding solvent reaction coordinate given in eq 2.6. The position of the barrier along the reaction coordinate is correspondingly shifted. For the mobile redox ion, the barrier for the oxidation reaction, Fe2+ f Fe3++ e-, is 12.3 kcal/mol, while the barrier for the reverse (reduction) reaction is 14.1 kcal/mol. The forward barrier is in good agreement with the experimental result32 of 13.5 ( 0.4 kcal/ mol, even though the value of the electronic coupling between the redox ion and the electrode was only estimated in the present study and there was a slight (1.8 kcal/mol) overpotential. (The activation barrier in the simulation comes primarily from the solvent reorganization and only weakly depends on the strength of the electronic coupling for the parameters chosen.). These values are to be compared with the results from simulations[1] where the iron ion was fixed above the surface which gave a barrier of 11.3 and 12.8 kcal/mol for the oxidation and reduction direction, respectively. Allowing the iron ion to move and naturally find its equilibrium height, which is slightly higher than that used in the fixed ion calculations, results in a slightly
ET Across the Electrode/Electrolyte Interface
Figure 4. Adiabatic free energy curve as a function of the full reaction coordinate, eq 2.13 (solid line), and as a function of the solvent reaction coordinate eq 2.6 (dotted-dashed line). The latter curve is shifted by the exclusion of the ion/surface interaction term.
higher barrier for the ET reaction, as one might expect from the corresponding decrease in the coupling to the electrode ∆(zion) with increasing distance. The qualitative features of the fixed ion calculation, however, are very similar, indicating that the solvent polarization fluctuations, as defined by the solvent reaction coordinate in eq 2.6, provide the dominant contribution to the barrier. To examine the latter idea more carefully, one can study other definitions of the reaction coordinate for the case of the heterogeneous ET. For example, one might not include the potential energy of the ion interacting with the electrode surface in the definition of the reaction coordinate. The same is true for the term arising from the images of the solvent interacting with the ion in eq 2.12. The resulting reaction coordinate is then just given by eq 2.6, i.e., the solvent reaction coordinate appropriate for the fixed ion simulations.1 To explore this idea in the simulation, the free energy was calculated for the case of the mobile redox ion, but with an umbrella sampling of the solvent reaction coordinate in eq 2.6. In Figure 4, the result of this calculation is compared to the result from Figure 3 which employed a sampling of the full reaction coordinate in eq 2.12. The results are similar to each other, although one is shifted from the other by a constant amount roughly equal to the average ion/electrode interaction energy. The result in Figure 3 confirms that the major factor in determining the free energy curve for the mobile ion simulation is indeed the solvent polarization fluctuations. On the other hand, when the curves for the different reaction coordinates are considered within the context of variational transition state theory,33 the full reaction coordinate in eq 2.12 represents a better choice of dividing surface because the oxidation free energy barrier calculated with this reaction coordinate is the highest, 12.3 kcal/mol versus 11.1 kcal/mol. In this picture, the redox ion fluctuations account for approximately a factor of 5 reduction in the ET rate. (As in the previous calculations, 25 umbrella windows were used for the solvent free energy curve. Each free energy window was equilibrated for 12.5 ps, and free energy averaging done for between 12.5 to 25 ps, depending upon the convergence of the statistics.) The effect of the fluctuations of the iron ion relative to the electrode surface were further analyzed through simulation by fixing the iron ion at its equilibrium distance of 5.5 Å above the surface and then calculating the ET free energy curve (cf. Figure 5). In this case, the free energy curve for the fixed ion was calculated as a function of the full reaction coordinate (eq 2.12), but the term for the ion/electrode interaction contributed a simple constant, since the ion could not fluctuate. The free energy barrier for ET is 13.4 kcal/mol for the oxidation direction, Fe2+f Fe3++ e- (compared to 12.3 kcal/mol for the mobile ion), and 13.1 kcal/mol for the reduction direction (compared to 14.1 kcal/mol for the mobile ion). This change in the free
J. Phys. Chem., Vol. 100, No. 25, 1996 10751
Figure 5. Adiabatic free energy curve for electron transfer with a fixed redox ion at its equilibrium distance from the electrode for the Fe3+ state (solid line). This curve was calculated with the same Fermi level as the mobile ion simulation (also shown), and the same reaction coordinate (eq 2.13) was sampled in each case for consistency.
Figure 6. Adiabatic free energy curve for the iron ion fixed at two different distances from the platinum surface.
energy curve can be directly attributed to the fluctuations of the redox ion. While the change in the curve would appear to be small, the fluctuations in the ion-electrode distance actually reverse the direction of the spontaneous ET reaction to the oxidation direction and affect the equillibrium constant by a factor of approximately 20. In order to explore the dependence of the adiabatic free energy curves on the distance of the redox ion from the surface, they were calculated for the iron ion at two different heights. The results are shown in Figure 6. The barrier to the Fe2+f Fe3++ e- oxidation reaction is significantly higher in the case where the iron ion is fixed at 7.5 Å above the surface than when the ion is fixed 5.1 Å. This is because the ion/electrode coupling is a function of the ion distance from the surface as specified by eq 2.15. The barrier in the ferric reduction direction is 19.8 kcal/mol when the ion is further from the surface, versus 12.8 kcal/mol when the ion is closer to the surface. This difference in barrier height translates into 6 orders of magnitude in the transition state theory rate constant. Such a large difference in rate could allow for the iron ion diffusion to compete with the solvent fluctuations as the rate limiting step in the ET reaction, at least for large distances from the surface. On the other hand, for large distances from the surface the question of whether the ET rate is adiabatic or nonadiabatic arises. Based on the ion/electrode coupling used in the present study (eq 2.15), the value of ∆(zion) is 1.7 kBT for zion ) 5.1 Å and 0.2 kBT for zion ) 7.5 Å. The latter value would seem to indicate the ET process is nonadiabatic at that distance, resulting in a further reduction in the rate. It should be noted, however, that a recent ab initio study has found that the electronic density from the substrate surface states may significantly extend into the first water layer through a coupling to the unfilled water electronic orbitals.34 Thus, the effective value of ∆ in real systems could be larger than that predicted from solid/vacuum interface considerations alone. Furthermore, a theoretical study35 has suggested that the effective value of ∆ should
10752 J. Phys. Chem., Vol. 100, No. 25, 1996
Calhoun and Voth
TABLE 1: Properties of the ET Free Energy Curves for Various Ion Pair Distances ion pair distance (Å)
forward barrier (kcal/mol)
reverse barrier (kcal/mol)
∆G0 (kcal/mol)
λ (kcal/mol)
∆emin (kcal/mol)
4.8 7.5 10.0 no counterion
41.7 14.1 12.0 12.3
2.1 16.5 18.0 14.1
+38.0 -2.4 -6.0 -2.0
70 60 60 57.4
469 488 502 531.8
actually be ∆ + πkBT, further increasing the likelihood of adiabatic ET. D. Free Energy Curves-Counterion Effects. The ET free energy curves were recalculated with a chloride counterion included in the simulation, as outlined in section IID. The local structure of the solvent and the relative position of the counterion to the iron ion are important in determining the free energy curves. The chloride ion interacts directly with the iron ion core through a Coulombic interaction and also changes the solvent structure and fluctuations around the ion. The calculations of the ET free energy curves can be carried out for fixed ion-counterion distances because the time scale of the ionic motion is slower than the solvent fluctuations relevant to ET (cf. eq 2.19). The calculations in the previous subsection show the systematic dependence of the ET free energy curves on the ionelectrode distance. This effect is important here because the lone ion obtains a different equilibrium height than when a counterion is nearby. There is also the question of how the chloride ion affects the ET reorganization free energy for the reaction. The latter effects can arise from a change in the solvent fluctuations near the redox ion, as well as from the direct Coulomb interaction of the counterion with the redox ion orbital. In order to systematically address the above issues through simulation, the ion pair distance was holonomically constrained, and the ET free energy curves were then calculated. Since the ion pair distance was not allowed to fluctuate, the reaction coordinate could be chosen to be the same as the reaction coordinate for the mobile iron ion alone, eq 2.13, or the more complicated reaction coordinate given in eq 2.18. The former choice facilitates a direct comparison to the free energy curves calculated without a counterion present, so that is the option chosen here. The image terms arising from the counterion are, however, still included in the Hamiltonian for the molecular dynamics simulation. To further facilitate comparison with the single mobile redox ion case, the electrode Fermi level used in the simulations was chosen to be the same value as the in the latter, i.e., f ) -10.2 eV. For each simulation, the ions were first fixed 5.1 Å above the electrode and then the solvent was equilibrated for 20 ps. The ion pair was next allowed to move, but constrained to remain a fixed distance apart. The free energy curves were then calculated for the 4 interionic distances, and each curve was sampled using 20 umbrella windows. Each umbrella window was first equilibrated for 10 ps and then sampled for 10 ps. The results of the simulations are shown in Figure 7. The ET free energy curves are clearly dependent upon the ion pair separation. One feature of the free energy curve for the smallest ion pair separation sampled, 4.8 Å, is a significant shift such that there is a large driving force for reduction to the ferrous state. This effect may arise because at this ion pair separation the chloride ion is on the edge of the first solvation shell of the iron ion and hence disrupts the solvation. There is a distinct second solvation shell,25 because of the high charge on the iron, and the chloride ion displaces part of that second shell as well. At a distance of 7.5 Å, the chloride ion is in the third solvation shell of the iron ion and hence has a smaller effect on the free energy curve relative to the single ion case, but it still clearly changes the ET free energy curve. At a separation of 10 Å,
Figure 7. Adiabatic free energy curves for electron transfer with the iron-chloride separation holonomically constrained at several distances, and with no counterion present. In each case, the Fermi level of the electrode is chosen to be that for the no counterion case, i.e., f ) -10.2 eV, and the reaction coordinate in eq 2.13 was sampled.
the curve seems to show further, but slow, convergence to the case of no counterion. In general, both the barrier height and the thermodynamic driving force are nonmonotonically shifted as the chloride ion is moved further away from the iron ion. It is unclear how to predict these shifts using simple principles, so this is clearly a challenge for future research. A summary of the estimated properties of the ET free energy curves for various ion pair separations is given in Table 1. The solvent reorganization energy is significantly shifted to 70.0 kcal/ mol when the chloride is 4.8 Å from the iron ion, compared to that for the other ion pair separations listed in Table 1, which are approximately given by 60 kcal/mol. Furthermore, the minimum of the curves in the Fe3+ state experiences a large shift as well. These changes in turn affect the overpotential as defined by eq 2.11, especially for an ion pair separation of 4.8 Å. Clearly, the counterion effects have significant implications for the ET rate across the water/electrode interface, and they illustrate the true microscopic complexity of the problem. IV. Concluding Remarks In this paper, the MD methodology for simulating adiabatic ET across an electrode/electrolyte interface has been significantly extended to include the effects of redox ion mobility and counterions. Within the underlying Anderson-Newns model, the definition of the ET reaction coordinate was accordingly extended to include the fluctuating interaction of the mobile ion with the electrode, with the images of the system due to the polarizability of the metal electrons, and with the counterions and their images. While a variational TST argument confirmed that the solvent polarization modes still provide the dominant contribution to the ET barrier and reaction free energy, the effects of the redox ion mobility were found to be important as well. These effects were manifested by tangible changes in the ET barrier height (which reduce the ET rate by a factor of ∼5) and thermodynamic driving force (which change the equilibrium constant by a factor of ∼20). The redox ion motion also introduces a source of nonlinearity into the system which is not completely reproduced by the analytic theory4,7 which assumes a linear response picture for the reaction coordinate fluctuations. The appropriate modification of the theory to include these effects will be the subject of future research.
ET Across the Electrode/Electrolyte Interface The MD simulations performed with a chloride counterion suggest that the microscopic details of the ET reaction are significantly more complicated in the presence of counterions. When the ET free energy curve was calculated with the counterion constrained to be at various distances from the redox ion, two important effects were observed. First, the iron ion was found to have a greater equilibrium height above the electrode, and this in turn raises the barrier to the ET reaction while influencing the reaction free energy. In fact, for the closest interionic distance (4.8 Å), the counterion disturbs the second solvation shell, clearly affecting both the reorganization free energy and redox ion solvation. As the interionic constraint is increased, the ET free energy curves change in a nonmonotonic fashion while slowly converging to the isolated redox ion limit at large separations. The present simulations have provided us with valuable insight into the various microscopic effects which can influence electron transfer across an electrode/electrolyte interface. Indeed, this study should be viewed as one more step toward our goal of providing a “complete” microscopic computational approach to study heterogeneous ET in electrochemical systems. Future studies from this research group and elsewhere will seek to further generalize the underlying simulation model to include multiple redox ions, a larger concentration of counterions, a more accurateseven first-principlessdescription of the electronic coupling between the redox ion(s) and the electrode, increasingly complex electrode surfaces, a more accurate and efficient treatment of the long-range forces, and explicit quantum dynamical studies of the ET dynamics beyond the TST approximation. While this is surely a challenging list of research targets, it seems worthwhile for us to continue the effort, one step at a time, particularly considering the central importance of electrochemical processes in the chemical, materials, and biological sciences. Acknowledgment. This research was supported by the Office of Naval Research and a Department of Defense AASERT Award. The computations reported in this study were supported in part by grants of computer time from the National Science Foundation (NSF) Supercomputer Centers under Grant No. MCA94P017 and from the Department of Defense High Performance Computer Centers. Some of the computations were also carried out on an IBM SP2 parallel computer purchased in part through an NSF Academic Research Infrastructure Grant (DMR94-13373).
J. Phys. Chem., Vol. 100, No. 25, 1996 10753 References and Notes (1) Straus, J. B.; Calhoun, A. W.; Voth, G. A. J. Chem. Phys. 1995, 102, 529. (2) Anderson, P. W. Phys. ReV. 1961, 124, 41. (3) Newns, D. M. Phys. ReV. 1969, 178, 1123. (4) Schmickler, W. J. Electroanal. Chem. 1986, 204, 31. (5) Grimley, T. B. Progress in Surface and Membrane Science; Academic: San Francisco, 1975; Vol. 9. (6) Muscat, J. P.; Newns, D. M. Prog. Surf. Sci. 1978, 9, 1. (7) Schmickler, W. Chem. Phys. Lett. 1995, 237, 152. (8) Sebastian, K. L. J. Chem. Phys. 1989, 90, 5056. (9) Smith, B.; Hynes, J. T. J. Chem. Phys. 1993, 99, 6517. (10) Straus, J. B.; Voth, G. A. J. Phys. Chem. 1993, 97, 7388. (11) Rose, D. A.; Benjamin, I. J. Chem. Phys. 1994, 100, 3545. (12) Rose, D. A.; Benjamin, I. Chem. Phys. Lett. 1995, 234, 209. (13) Xia, X.; Berkowitz, M. L. Chem.Phys. Lett. 1994, 227, 561. (14) Smith, B. B.; Halley, J. W. J. Chem. Phys. 1994, 101, 10915. (15) A. King, G.; Warshel J. Chem. Phys. 1990, 93, 8682. (16) Warshel, A.; Parson, W. W. Annu. ReV. Phys. Chem. 1991, 42, 279. (17) Moser, C. C.; Keske, J. M.; Warncke, K.; Farid, R. S.; Dutton, P. L. Nature 1992, 355, 796. (18) Nazmutdinov, R. R.; Spohr, E. J. Phys. Chem 1994, 98, 5956. (19) Spohn, P. D.; Brill, T. R. J. Phys. Chem. 1989, 93, 6224. (20) Raghavan, K.; Foster, K.; Motakabbir, K.; Berkowitz, M. J. Chem. Phys. 1991, 94, 2110. (21) Spohr, E.; Heinzinger, K. Ber. Bunsen-Ges. Phys. Chem. 1988, 92, 1358. (22) Hautman, J.; Halley, J. W.; Rhee, Y. J. J. Chem. Phys. 1989, 91, 467. (23) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (24) Lobaugh, J. H. Ph.D. thesis, Department of Chemistry, University of Pennsylvania, Philadelphia, PA, 1995. (25) Kuharski, R. A.; Bader, J. S.; Chandler, D.; Sprik, M.; Klein, M. L J. Chem. Phys. 1988, 89, 3248. (26) Nose´, S. J. Chem. Phys. 1984, 81, 511. (27) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2695. (28) Rose, D. A.; Benjamin, I. J. Chem. Phys. 1991, 95, 6856. (29) Rose, D. A.; Benjamin, I. J. Chem. Phys. 1992, 98, 2283. (30) Dang, L. X. J. Chem. Phys. 1992, 97, 1919. (31) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1984, 106, 903. (32) Curtiss, L. A.; Halley, J. W.; Hautman, J.; Hung, N. C.; Nagy, Z.; Rhee, Y. J.; Yonco, R. M. J. Electrochem. Soc. 1991, 138, 2032. (33) Truhlar, D. G.; Hase, W. L.; Hynes, J. T. J. Phys. Chem. 1983, 87, 2664 (5523E). (34) Ursenbach, C. P.; Voth, G. A. J. Chem. Phys. 1995, 103, 7569. (35) Gorodyskii, A. V.; Karasevskii, A. I.; Matyushov, D. V. J. Electroanal. Chem. 1991, 315, 9.
JP960603Q
