- Email: [email protected]
Large-scale computer simulation of an electrochemical bond-breaking reaction
14 May 1999
Chemical Physics Letters 305 Ž1999. 94–100
Large-scale computer simulation of an electrochemical bond-breaking reaction August Calhoun a , Marc T.M. Koper b, Gregory A. Voth a b
a,)
Department of Chemistry and Henry Eyring Center for Theoretical Chemistry, UniÕersity of Utah, Salt Lake City, UT 84112, USA Schuit Institute of Catalysis, Laboratory of Inorganic Chemistry and Catalysis, EindhoÕen UniÕersity of Technology, NL-5600 MB EindhoÕen, The Netherlands Received 11 January 1999; in final form 12 March 1999
Abstract A novel Hamiltonian is employed to explicitly simulate an electrochemical bond-breaking reaction in which an electron-transfer reaction is directly coupled to the dissociation of a molecular species. The free energy surface as a function of both the collective solvation coordinate of the electron transfer and the intramolecular bond length of a CH 3 Cl molecule is computed by virtue of a classical molecular dynamics ŽMD. simulation. The method is also easily generalized to treat a variety of electrochemically catalyzed phenomenon. The simulation data show very significant deviations from the predictions of standard analytical theory. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction The enormous advances in the capabilities of modern-day computers promise to revolutionize our understanding of molecular-scale processes, including electron- and charge-transfer reactions in the liquid phase. In the traditional viewpoint of the Marcus–Levich–Dogonadze theory, the subsystem exhibiting electron or charge exchange is coupled linearly to the solvent polarization fluctuations, which are in turn modeled in the linear response approximation. This solvent fluctuation model has been successfully applied and extended to a variety of solution-phase reactions, and computer simulations have now made it possible to test more carefully –
)
Corresponding author. E-mail: [email protected]
and possibly improve – its fundamental underlying assumptions. A particularly interesting and technologically important class of solution-phase reactions are electrochemical reactions occurring at metalrsolution interfaces w1–4x. In recent years, molecular dynamics computer simulations have aided significantly in building a more molecular picture of the structure and dynamics at electrified metalrsolution interfaces w5–12x. Molecular dynamics ŽMD. simulations w9,10x of a typical outer-sphere electron-transfer ŽET. reaction such as the Fe 3qrFe 2q couple at the platinumrwater interface have shown that the Marcus solvent picture is reasonably accurate for such reactions, but they have also made it clear that for reactions in which the initial or final state is adsorbed onto the metal surface quite a different picture is needed in which the ion has to lose part of its
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 3 5 3 - X
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
solvation shell and to displace adsorbed water molecules from the electrode surface. So far, there is no clear picture of electrochemical reactions involving bond breaking at metal electrodes. Clearly, such reactions are of great importance in electrochemical synthesis and electrocatalysis. Saveant has proposed an analytical model for ´ concerted bond-breaking and electron-transfer reactions in which he adopted the Marcus linear response solvent picture w1–4x. In this model, the electron transfer from the electrode to the molecule leads to the concerted dissociation of the molecule into two fragments. Saveant and his co-workers have pub´ lished extensive experimental studies of such reactions, of which the reductive cleavage of alkyl halides is the prototypical example. Recently, we have proposed a generalization of Saveant’s model by explic´ itly incorporating the metal electronic levels and coupling to the electrode through an explicit Hamiltonian approach w13x. The underlying Hamiltonian has been used extensively in the modeling of various types of electrode reactions. However, in both Saveant’s and our model of Ref. w13x, the solvent is ´ characterized by a single reorganization energy l, the central parameter in the Marcus solvent model, with little evidence for the accuracy of the approximation for such a complex kinetic process Ži.e., one involving both solvent fluctuation-driven ET and bond breaking.. The Hamiltonian approach of Ref. w13x offers two important advantages over the original Saveant ´ model. First, the strength of the electronic interaction is explicitly included which has allowed us recently to extend the Saveant ´ model to treat real electrocatalytic reactions such as dissociative adsorption w14x. Secondly, as we have shown in earlier work w9–11x, the solvent part of the Hamiltonian can be described by an explicit molecular dynamics ŽMD. model, and in this way the role of the solvent in electron-transfer reactions can be treated essentially without approximations Žapart from the approximations involved in model solvent potentials.. In this Letter, we report the first computer simulation of an electrochemical bond-breaking reaction, specifically the reductive cleavage of methyl chloride at a platinum electrode in water .
CH 3 Cl q ey ™ CH 3 qCly .
Ž 1.
95
Fig. 1. A schematic diagram of the electrochemical BBET process. This diagram shows the electron being transferred to the CH 3 Cl redox species followed by a dissociation of the radical. The reaction occurs in a concerted manner.
A pictorial view of the reaction is given in Fig. 1. ŽNote that, in accordance with the conclusions from the experimental work of Saveant ´ and co-workers, we assume the methyl chloride is not adsorbed onto the metal electrode.. The present study consists of mapping out the free energy surface for the above reaction as a function of two coupled coordinates: the traditional generalized solvent fluctuation coordinate from the Marcus theory, and the distance r between the dissociating methyl and chloride fragments. The resulting two-dimensional free energy surface is unique as it is, to our knowledge, the first of its kind in the literature on electrochemical charge-transfer reactions. More importantly, our simulations illustrate very clearly the importance of including explicit molecular level information in order to have a quantitative physical picture of bond-breaking electrode reactions. It will be shown that the solvent and solute molecularity induces new free energy barriers not present in the previous models, and the real solvent response to reaction Eq. Ž1. is nonlinear, leading to an ‘effective’ solvent reorganization energy that is solvent coordinate dependent. This last conclusion contrasts importantly with the findings from MD simulations of outer-sphere electron transfer w9x, for which the Marcus picture was found to give quite a reasonable description of the simulation results.
2. Theoretical approach The Hamiltonian in the present simulation model is similar to that of Koper and Voth w13x, consisting of three parts: H s Hsolvent q Helect q H bb ,
Ž 2.
96
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
where Hsolvent is the standard classical MD Hamiltonian for the solvent and solvent–electrode interactions w9x, Helect is the electronic part, and H bb is an additional term that characterizes the effect of the donation of an electron from the electrode into a molecular antibonding orbital on the redox species Žthe CH 3 Cl in this case.. The electronic term has its origin in the the Anderson–Newns Hamiltonian w15– 18x which has been extended and applied to a number electrochemical systems w9–11,19–23x. To be more specific, the electronic part of the Hamiltonian, Helect , describes a single acceptor orbital on a redox species in solution coupled to the continuum of metal electrode states. Within the Born–Oppenheimer approximation, it is possible to find an analytical expression for the electronic energy of the system, E0 , as a function of D e, the collective Marcus solvation coordinate w11,24x. The resulting expression can be readily incorporated into classical MD simulations, because it is a function of only the nuclear positions of the system w9–11x. Also, the average electronic population of the redox species orbital, ² n a :, can be written as a function of these nuclear positions w17,19x. It is a continuous function of the solvent coordinates that can be used to compute the average population of the redox orbital at any point in the MD simulation. The electronic energy and average population is thus a unique function of the solvent coordinates, which is defined in the traditional Marcus picture as being a collective reaction coordinate consisting of the sum of all of the charge interactions with a unit electronic charge at the redox site, and the redox image Žreflected through the metal surface in the perfect conductor approximation for the metal. w10x. This can also be thought of as a difference potential, where the only change between the two electronic states is the ionic charge. MD simulations can then be performed with this Hamiltonian. The bond-breaking term in the Hamiltonian in Eq. Ž2. was proposed by Koper and Voth w13x. If the electrode states are coupled to an antibonding orbital in the redox molecule, then the ET event can be followed by a molecular dissociation Žsee Fig. 1.. Experimental and theoretical investigations have shown that this is indeed the case, and that the bond-breaking electron transfer ŽBBET. proceeds in a concerted single-step manner w1–4x. Thus, a natural
second quantized operator to describe this effect can be constructed from the antibonding orbital population operator, n a , given by w13x H bb s Ž 1 y n a . V bonding q n aVantibonding .
Ž 3.
This Hamiltonian switches off a bond Ž V bonding . in the redox species as an antibonding orbital is populated Žby the electrochemical ET.. Another potential term, in this case an antibonding function Ž Vantibonding ., is switched on as the antibonding orbital is occupied. A more general formalism can be utilized where a function like Eq. Ž3. couples two generic physical states of the redox molecule. The H bb term in the present model is given by Hb b s Ž 1 y n a . w V Morse q Vsolvent x q n a Vrep ,
Ž 4.
where the V Mo rse and Vrep are Koper and Voth’s terms w13x. The parameters for the Morse functions are taken from Ref. w4x. Vsolvent is the reactant methyl chloride–solvent Coulombic interaction. It should be noted that the Koper–Voth analytical model of Ref. w13x does not contain the latter term Ži.e., it is assumed that the reactant is non-polar.. It is not straightforward to incorporate the dipolar solvation into the analytic model, but it can be done approximately by modeling the reactant solvation through a linear response model and then neglecting the remaining dependence of the redox orbital population, ² n a :, on the term Vsolvent . This term is much smaller than the electrostatic interaction with the charged CH 3 Cly metastable species.
3. Simulation details and results The specific system in the present study was a CH 3 Cl molecule near a PtŽ111. surface in water. This electrode and solvent were chosen because reasonably accurate interaction potentials are available. Experimentally, such bond-breaking reactions are often studied at gold in non-aqueous solvents, but this should not distract from the fundamental results to be reported below. An antibonding orbital on the redox molecule was coupled to the metal electronic states with the Anderson–Newns Hamiltonian, and the population of the redox orbital was in turn coupled to the CH 3 Cl intramolecular bond with the bond breaking part of the total Hamiltonian.
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
In this simulation model, the CH 3 Cl molecule was represented as a two-site diatomic molecule, so the methyl fragment was treated in the united atom approximation. The TIPS parameters for liquid CH 3 Cl was used to describe the water–MeCl interactions in the MD simulation w25x. A total of 320 SPCrF2 w11x water molecules were included in the molecular dynamics cell which had dimensions of ˚ = 19.16 A˚ = 42.33 A˚ Ž x–y–z .. The water 24.93 A molecules were restricted to be in the top half of the molecular dynamics cell by the metal surface, and the image of the MeCl and of the electronic population on the redox site were reflected through the metal surface in the perfect conductor approximation for the electrode w11x. The solvent, therefore, interacts with the MeCl and the MeCl image, both of which have a variable charge. The images of the water molecules were not included in the present model because this effect is accounted by in a static way by the water–metal interaction w11x. The entire molecular dynamics cell was replicated in three dimensions via the Ewald summation, which, including the ion image, results in a charge neutral cell w11x. Periodic boundary conditions were used, and the z box size was adjusted to yield the bulk density of water in the center. The x and y dimensions were chosen to be commensurate with an integer number of PtŽ111. unit cells w11x. The equations of motion were integrated with a multiple time-step algorithm based on the velocity Verlet method w11x, and a Nose–Hoover chain of length 2 was attached to each particle to maintain temperature control and ergodic sampling w11x. A small time-step of 14 fs was used to integrate the intramolecular modes while a large time-step of 1 fs was used to integrate the intermolecular modes in the equations of motion. The electronic coupling was chosen to be 55 meV, because this value is in a reasonable range for a typical ET reaction. Also, this is in accord with ab initio calculations w11x. The vacuum ET acceptor energy level of the MeCl molecule was taken to be that of Cl Žy3.617 eV., and the Fermi level of the electrode was chosen to be 7.98 eV in order to give a zero electrode overpotential in the Koper–Voth analytical model. The free energy surface was then computed using the MD model described above in which the two reaction coordinates Žthe collective solvent coordinate and the methyl chloride bond length.
97
define a two-dimensional free energy surface for the reaction. In the Marcus–Saveant ´ picture, it is assumed that the solvent free energy responds linearly to changes in the charge distribution, which is equivalent to a quadratic dependence on the solvent polarization coordinate D e. Therefore, this solvent free energy is usually written as Fsolvent s FMeCl , min q
1 4l
2 Ž D e y D e0 . ,
Ž 5.
where FMeCl, min is the solvation free energy of the MeCl molecule near the interface, and l is the solvent reorganization energy, with D e 0 being the minimum of the solvent free energy curve. If the Marcus–Saveant solvent picture applies, these pa´ rameters can be computed from a parabolic fit to an MD simulation of the diabatic undissociated methyl chloride state, and can in turn be used to calculate the entire potential energy surface as prescribed by a modified version of the model of Koper and Voth w13x. The actual free energy surface computed from an explicit MD simulation with the BBET Hamiltonian is shown in Fig. 2. For comparison, the free energy surface computed from the analytic model of Ref. w13x using the linear response model for the solvent free energy ŽEq. Ž5.. is also shown. The free energy surface computed by computer simulation is a result of over 2700 umbrella sampling calculations. The parts of the free energy with the most error are near the edges of the surface. The important part of the surface, the barrier and well regions, were sampled the most. The errors here are ; 1 kcalrmol Žcomputed from the standard deviation for different averaging intervals of the free energy calculation.. The simulation free energy surface shown in the upper part of Fig. 2 represents 114,000 node hours of computation on parallel computers. This is equivalent to running a single present-day high-end workstation for 13 years. From the results in Fig. 2, it is seen that the gross qualitative features of the analytic model are similar to the actual free energy surface computed from the simulation. First, there is a bound well where the CH 3 Cl is stable and the electron has not been transferred from the electrode to the CH 3 Cl. There is also an exit channel, because the CH 3 Cly radical is not
98
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
Fig. 2. A contour plot of the adiabatic free energy surface computed from explicit MD simulation Žtop panel. and from the linear response model for the solvent Žbottom panel.. The well Župper left. and exit channel Žlower right. are clearly visible in both surfaces. The well corresponds to the neutral CH 3 Cl molecule, and the exit channel to the unstable charged radical. A solvent-induced barrier in the exit channel is also seen in the simulation data, occurring where the fragments are still in contact.
stable. Thus, the analytical free energy surface contains the essential physics of the BBET model: a stable CH 3 Cl species and an unstable CH 3 Cly radical. There is, however, at least one key qualitative difference between the free energy surfaces computed via computer simulation and the analytic model. In particular, there is a well in the exit channel of the simulation result where the methyl fragment and the chloride fragment are still in contact. This well, which is ; 7 kcalrmol in depth, is a result of explicit microscopic interactions in the MD model. When the fragments are in contact, or separated by a single water molecule, the result is a configuration with a lower energy than when the fragments are separated, but not separated enough to allow more water molecules in between the fragments. There is an energy ‘cost’ associated with making this hole in the solvent Žthe intermediate state., therefore the contact fragment pair and solvent
separated fragment pair are distinguished by a free energy barrier between them. Another key qualitative difference is a positive free energy of reaction Žoverpotential. in the simulation data which is not present in the analytic model for the given set of parameters. The contour plots in Fig. 2 also reveal important quantitatiÕe differences between the simulation results and the linear response model predictions. The curvature, minima, and reaction free energy Žoverpotential. differ considerably between the linear response solvent and the MD simulation data. This is due to the non-linear response of the solvent in the MD simulation as the CH 3 Cl bond is broken, which is in turn related to the change in the charge on the redox center from y0.25 in the CH 3 Cl molecule to y1 for the chloride ion, and to the change in the volume exclusion effect from the Me fragment at different values for r. The primary quantitative difference between the simulation results and the analytical prediction is a 10 kcalrmol positive overpotential for the former in forming first the metastable MeCly ion and then subsequently the Me and Cly fragments. The analytical model thus predicts a steady-state concentration of the metastable ionic intermediate Žand ultimately the products. which is too large by seÕen orders of magnitude. Interestingly, the rate constant Žactivation free energy. for the ET event from the electrode to the MeCl is predicted to be the roughly same by both simulation and theory. However, the thermodynamics Žequilibrium constant, etc.. of the reaction at the specified electrode conditions are found to differ greatly. It should be noted that the above comparison of the analytic model with the simulation results is based on the usage of a common electrode Fermi level in the underlying Hamiltonian ŽEq. Ž2.. which appears in both approaches. One might argue that a better strategy is to compare the prediction for the rate constant from the analytical model with that from the simulation at a common overpotential, i.e., to adopt a more phenomenological approach which has its footing in experimentally specified quantities. Along these lines, the expression for the activation free energy D Fact in the analytical model is given to a good approximation by w1,13x D Fact s
lqD 4
2
h
ž
1q
lqD
/
,
Ž 6.
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
99
Fig. 3. A plot of l, the reorganization energy, as a function of r, the CH 3 Cl dissociation coordinate. The reorganization energy is determined from global harmonic fits to the free energy surface at fixed values of r.
where h is the overpotential, l is the reorganization free energy, and D is the gas-phase dissociation energy of the MeCl bond. ŽIt should be noted that the sum l q D acts like an overall reorganization energy in this theory.. For the case of a 10 kcalrmol positive overpotential as observed in the simulation, the above expression predicts a kinetic barrier to ET of 42 kcalrmol, while the simulation barrier is found to be 37 kcalrmol. This difference in barrier corresponds to a difference of more than three orders of magnitude in the ET rate constant between the analytical prediction and the simulation result. Clearly, the fundamental quantitative differences between the simulation results and the analytical predictions do not disappear when viewed in this alternative light. When the MD simulation data are analyzed in some detail, it is found that in reality the effective reorganization energy l depends on the charge state of the system, which is different for the molecule and the fragment state, and on the bond length or inter-fragment distance r. Furthermore, a l determined from a fit locally about a minimum often gives a poor fit of the entire free energy curve as obtained from the umbrella sampling procedure. Therefore, an alternative method of extracting the
effective l from the MD simulations is to fit the entire free energy curve of one state, at least up to the transition state, to a parabola. The l obtained in this way is close to another definition of the solvent reorganization energy, namely as the vertical energy shift from the diabatic state under consideration to the minimum of the other diabatic state Žat equilibrium., as it effectively includes deviations from nonparabolicity of the overall diabatic energy curve. The solvent reorganization energies determined in this way, and their dependence of the bond length or inter-fragment distance r, for both the molecule and the fragment states, are shown in Fig. 3. The l as determined from a local fit about the methyl chloride minimum is also shown as a dashed line in Fig. 3 and is seen to give a poor description of the real solvent response.
4. Concluding remarks In conclusion, for the first time very large-scale MD simulations have been carried out for a bondbreaking electrochemical reaction in terms of two coupled reaction coordinates, the solvent fluctuation
100
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
and molecular dissociation coordinates. The simulations were compared to an analytical free energy surface calculated using the linear response solvent model. Although the latter model is a convenient tool for estimating trends and gross features of the reaction, it fails to reproduce the molecular details that are found to be very important for bond-breaking reactions. Specifically, the simulations reveal a free energy well for the contact fragment pair due to the molecularity of the solvent Žthe solvent ‘cage effect’., and the solvent polarization fluctuations driving the bond-breaking process are seen to be significantly non-linear, as manifested by an effective solvent reorganization energy which is reaction coordinate dependent. The latter effects contribute to large quantitative differences between the simulation results and the predictions of the existing analytical theory. The present results thereby provide a clear and unique illustration of the indispensable role accurate and detailed computer simulations can play in obtaining a proper understanding of complicated electrochemical reactions.
Acknowledgements This research was supported by the United States Office of Naval Research and the Royal Netherlands Academy of Arts and Sciences ŽKNAW.. Computations were carried out at the Aeronautical Systems Center Major Shared Resource Center at WrightPatterson Air Force Base through a Department of Defense Grand Challenge Computing Grant, and at the National Center for Supercomputing Applications ŽNCSA. through NSF Metacenter grant number
MCA94P017P. We also gratefully acknowledge the University Utah Center for High Performance Computing for a generous grant of Silicon Graphics Origin 2000 time. References w1x J. Saveant, Acc. Chem. Res. 26 Ž1993. 455. ´ w2x J. Bertran, I. Gallardo, M. Moreno, J. Saveant, J. Am. Chem. ´ Soc. 114 Ž1992. 9576. w3x C.P. Andrieux, A.L. Gorande, J. Saveant, J. Am. Chem. Soc. ´ 114 Ž1992. 6892. w4x J. Saveant, J. Am. Chem. Soc. 109 Ž1987. 6788. ´ w5x X. Xia, M.L. Berkowitz, Chem. Phys. Lett. 227 Ž1994. 561. w6x D.A. Rose, I. Benjamin, Chem. Phys. Lett. 234 Ž1995. 209. w7x D.L. Price, J.W. Halley, J. Chem. Phys. 102 Ž1995. 6603. w8x J.B. Straus, G.A. Voth, J. Phys. Chem. 97 Ž1993. 7388. w9x J.B. Straus, A. Calhoun, G.A. Voth, J. Chem. Phys. 102 Ž1995. 529. w10x A. Calhoun, G.A. Voth, J. Phys. Chem. 100 Ž1996. 10746. w11x A. Calhoun, G.A. Voth, J. Electroanal. Chem. 450 Ž1998. 253. w12x A. Calhoun, G.A. Voth, J. Chem. Phys. 109 Ž1998. 4569. w13x M.T.M. Koper, G.A. Voth, Chem. Phys. Lett. 282 Ž1998. 100. w14x M.T.M. Koper, G.A. Voth, J. Chem. Phys. 109 Ž1998. 1991. w15x P.W. Anderson, Phys. Rev. 124 Ž1961. 41. w16x D.M. Newns, Phys. Rev. 178 Ž1969. 1123. w17x W. Schmickler, J. Electroanal. Chem. 204 Ž1986. 31. w18x A. Kornyshev, W. Schmickler, J. Electroanal. Chem. 185 Ž1985. 253. w19x K.L. Sebastian, J. Chem. Phys. 90 Ž1989. 5056. w20x M.T.M. Koper, W. Schmickler, Chem. Phys. 211 Ž1996. 123. w21x M.T.M. Koper, J. Phys. Chem. 101 Ž1997. 3168. w22x M.T.M. Koper, J.H. Mohr, W. Schmickler, Chem. Phys. 220 Ž1997. 95. w23x B. Smith, J.T. Hynes, J. Chem. Phys. 99 Ž1993. 6517. w24x R.A. Marcus, J. Chem. Phys. 43 Ž1965. 679. w25x W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 Ž1983. 926.
Chemical Physics Letters 305 Ž1999. 94–100
Large-scale computer simulation of an electrochemical bond-breaking reaction August Calhoun a , Marc T.M. Koper b, Gregory A. Voth a b
a,)
Department of Chemistry and Henry Eyring Center for Theoretical Chemistry, UniÕersity of Utah, Salt Lake City, UT 84112, USA Schuit Institute of Catalysis, Laboratory of Inorganic Chemistry and Catalysis, EindhoÕen UniÕersity of Technology, NL-5600 MB EindhoÕen, The Netherlands Received 11 January 1999; in final form 12 March 1999
Abstract A novel Hamiltonian is employed to explicitly simulate an electrochemical bond-breaking reaction in which an electron-transfer reaction is directly coupled to the dissociation of a molecular species. The free energy surface as a function of both the collective solvation coordinate of the electron transfer and the intramolecular bond length of a CH 3 Cl molecule is computed by virtue of a classical molecular dynamics ŽMD. simulation. The method is also easily generalized to treat a variety of electrochemically catalyzed phenomenon. The simulation data show very significant deviations from the predictions of standard analytical theory. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction The enormous advances in the capabilities of modern-day computers promise to revolutionize our understanding of molecular-scale processes, including electron- and charge-transfer reactions in the liquid phase. In the traditional viewpoint of the Marcus–Levich–Dogonadze theory, the subsystem exhibiting electron or charge exchange is coupled linearly to the solvent polarization fluctuations, which are in turn modeled in the linear response approximation. This solvent fluctuation model has been successfully applied and extended to a variety of solution-phase reactions, and computer simulations have now made it possible to test more carefully –
)
Corresponding author. E-mail: [email protected]
and possibly improve – its fundamental underlying assumptions. A particularly interesting and technologically important class of solution-phase reactions are electrochemical reactions occurring at metalrsolution interfaces w1–4x. In recent years, molecular dynamics computer simulations have aided significantly in building a more molecular picture of the structure and dynamics at electrified metalrsolution interfaces w5–12x. Molecular dynamics ŽMD. simulations w9,10x of a typical outer-sphere electron-transfer ŽET. reaction such as the Fe 3qrFe 2q couple at the platinumrwater interface have shown that the Marcus solvent picture is reasonably accurate for such reactions, but they have also made it clear that for reactions in which the initial or final state is adsorbed onto the metal surface quite a different picture is needed in which the ion has to lose part of its
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 3 5 3 - X
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
solvation shell and to displace adsorbed water molecules from the electrode surface. So far, there is no clear picture of electrochemical reactions involving bond breaking at metal electrodes. Clearly, such reactions are of great importance in electrochemical synthesis and electrocatalysis. Saveant has proposed an analytical model for ´ concerted bond-breaking and electron-transfer reactions in which he adopted the Marcus linear response solvent picture w1–4x. In this model, the electron transfer from the electrode to the molecule leads to the concerted dissociation of the molecule into two fragments. Saveant and his co-workers have pub´ lished extensive experimental studies of such reactions, of which the reductive cleavage of alkyl halides is the prototypical example. Recently, we have proposed a generalization of Saveant’s model by explic´ itly incorporating the metal electronic levels and coupling to the electrode through an explicit Hamiltonian approach w13x. The underlying Hamiltonian has been used extensively in the modeling of various types of electrode reactions. However, in both Saveant’s and our model of Ref. w13x, the solvent is ´ characterized by a single reorganization energy l, the central parameter in the Marcus solvent model, with little evidence for the accuracy of the approximation for such a complex kinetic process Ži.e., one involving both solvent fluctuation-driven ET and bond breaking.. The Hamiltonian approach of Ref. w13x offers two important advantages over the original Saveant ´ model. First, the strength of the electronic interaction is explicitly included which has allowed us recently to extend the Saveant ´ model to treat real electrocatalytic reactions such as dissociative adsorption w14x. Secondly, as we have shown in earlier work w9–11x, the solvent part of the Hamiltonian can be described by an explicit molecular dynamics ŽMD. model, and in this way the role of the solvent in electron-transfer reactions can be treated essentially without approximations Žapart from the approximations involved in model solvent potentials.. In this Letter, we report the first computer simulation of an electrochemical bond-breaking reaction, specifically the reductive cleavage of methyl chloride at a platinum electrode in water .
CH 3 Cl q ey ™ CH 3 qCly .
Ž 1.
95
Fig. 1. A schematic diagram of the electrochemical BBET process. This diagram shows the electron being transferred to the CH 3 Cl redox species followed by a dissociation of the radical. The reaction occurs in a concerted manner.
A pictorial view of the reaction is given in Fig. 1. ŽNote that, in accordance with the conclusions from the experimental work of Saveant ´ and co-workers, we assume the methyl chloride is not adsorbed onto the metal electrode.. The present study consists of mapping out the free energy surface for the above reaction as a function of two coupled coordinates: the traditional generalized solvent fluctuation coordinate from the Marcus theory, and the distance r between the dissociating methyl and chloride fragments. The resulting two-dimensional free energy surface is unique as it is, to our knowledge, the first of its kind in the literature on electrochemical charge-transfer reactions. More importantly, our simulations illustrate very clearly the importance of including explicit molecular level information in order to have a quantitative physical picture of bond-breaking electrode reactions. It will be shown that the solvent and solute molecularity induces new free energy barriers not present in the previous models, and the real solvent response to reaction Eq. Ž1. is nonlinear, leading to an ‘effective’ solvent reorganization energy that is solvent coordinate dependent. This last conclusion contrasts importantly with the findings from MD simulations of outer-sphere electron transfer w9x, for which the Marcus picture was found to give quite a reasonable description of the simulation results.
2. Theoretical approach The Hamiltonian in the present simulation model is similar to that of Koper and Voth w13x, consisting of three parts: H s Hsolvent q Helect q H bb ,
Ž 2.
96
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
where Hsolvent is the standard classical MD Hamiltonian for the solvent and solvent–electrode interactions w9x, Helect is the electronic part, and H bb is an additional term that characterizes the effect of the donation of an electron from the electrode into a molecular antibonding orbital on the redox species Žthe CH 3 Cl in this case.. The electronic term has its origin in the the Anderson–Newns Hamiltonian w15– 18x which has been extended and applied to a number electrochemical systems w9–11,19–23x. To be more specific, the electronic part of the Hamiltonian, Helect , describes a single acceptor orbital on a redox species in solution coupled to the continuum of metal electrode states. Within the Born–Oppenheimer approximation, it is possible to find an analytical expression for the electronic energy of the system, E0 , as a function of D e, the collective Marcus solvation coordinate w11,24x. The resulting expression can be readily incorporated into classical MD simulations, because it is a function of only the nuclear positions of the system w9–11x. Also, the average electronic population of the redox species orbital, ² n a :, can be written as a function of these nuclear positions w17,19x. It is a continuous function of the solvent coordinates that can be used to compute the average population of the redox orbital at any point in the MD simulation. The electronic energy and average population is thus a unique function of the solvent coordinates, which is defined in the traditional Marcus picture as being a collective reaction coordinate consisting of the sum of all of the charge interactions with a unit electronic charge at the redox site, and the redox image Žreflected through the metal surface in the perfect conductor approximation for the metal. w10x. This can also be thought of as a difference potential, where the only change between the two electronic states is the ionic charge. MD simulations can then be performed with this Hamiltonian. The bond-breaking term in the Hamiltonian in Eq. Ž2. was proposed by Koper and Voth w13x. If the electrode states are coupled to an antibonding orbital in the redox molecule, then the ET event can be followed by a molecular dissociation Žsee Fig. 1.. Experimental and theoretical investigations have shown that this is indeed the case, and that the bond-breaking electron transfer ŽBBET. proceeds in a concerted single-step manner w1–4x. Thus, a natural
second quantized operator to describe this effect can be constructed from the antibonding orbital population operator, n a , given by w13x H bb s Ž 1 y n a . V bonding q n aVantibonding .
Ž 3.
This Hamiltonian switches off a bond Ž V bonding . in the redox species as an antibonding orbital is populated Žby the electrochemical ET.. Another potential term, in this case an antibonding function Ž Vantibonding ., is switched on as the antibonding orbital is occupied. A more general formalism can be utilized where a function like Eq. Ž3. couples two generic physical states of the redox molecule. The H bb term in the present model is given by Hb b s Ž 1 y n a . w V Morse q Vsolvent x q n a Vrep ,
Ž 4.
where the V Mo rse and Vrep are Koper and Voth’s terms w13x. The parameters for the Morse functions are taken from Ref. w4x. Vsolvent is the reactant methyl chloride–solvent Coulombic interaction. It should be noted that the Koper–Voth analytical model of Ref. w13x does not contain the latter term Ži.e., it is assumed that the reactant is non-polar.. It is not straightforward to incorporate the dipolar solvation into the analytic model, but it can be done approximately by modeling the reactant solvation through a linear response model and then neglecting the remaining dependence of the redox orbital population, ² n a :, on the term Vsolvent . This term is much smaller than the electrostatic interaction with the charged CH 3 Cly metastable species.
3. Simulation details and results The specific system in the present study was a CH 3 Cl molecule near a PtŽ111. surface in water. This electrode and solvent were chosen because reasonably accurate interaction potentials are available. Experimentally, such bond-breaking reactions are often studied at gold in non-aqueous solvents, but this should not distract from the fundamental results to be reported below. An antibonding orbital on the redox molecule was coupled to the metal electronic states with the Anderson–Newns Hamiltonian, and the population of the redox orbital was in turn coupled to the CH 3 Cl intramolecular bond with the bond breaking part of the total Hamiltonian.
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
In this simulation model, the CH 3 Cl molecule was represented as a two-site diatomic molecule, so the methyl fragment was treated in the united atom approximation. The TIPS parameters for liquid CH 3 Cl was used to describe the water–MeCl interactions in the MD simulation w25x. A total of 320 SPCrF2 w11x water molecules were included in the molecular dynamics cell which had dimensions of ˚ = 19.16 A˚ = 42.33 A˚ Ž x–y–z .. The water 24.93 A molecules were restricted to be in the top half of the molecular dynamics cell by the metal surface, and the image of the MeCl and of the electronic population on the redox site were reflected through the metal surface in the perfect conductor approximation for the electrode w11x. The solvent, therefore, interacts with the MeCl and the MeCl image, both of which have a variable charge. The images of the water molecules were not included in the present model because this effect is accounted by in a static way by the water–metal interaction w11x. The entire molecular dynamics cell was replicated in three dimensions via the Ewald summation, which, including the ion image, results in a charge neutral cell w11x. Periodic boundary conditions were used, and the z box size was adjusted to yield the bulk density of water in the center. The x and y dimensions were chosen to be commensurate with an integer number of PtŽ111. unit cells w11x. The equations of motion were integrated with a multiple time-step algorithm based on the velocity Verlet method w11x, and a Nose–Hoover chain of length 2 was attached to each particle to maintain temperature control and ergodic sampling w11x. A small time-step of 14 fs was used to integrate the intramolecular modes while a large time-step of 1 fs was used to integrate the intermolecular modes in the equations of motion. The electronic coupling was chosen to be 55 meV, because this value is in a reasonable range for a typical ET reaction. Also, this is in accord with ab initio calculations w11x. The vacuum ET acceptor energy level of the MeCl molecule was taken to be that of Cl Žy3.617 eV., and the Fermi level of the electrode was chosen to be 7.98 eV in order to give a zero electrode overpotential in the Koper–Voth analytical model. The free energy surface was then computed using the MD model described above in which the two reaction coordinates Žthe collective solvent coordinate and the methyl chloride bond length.
97
define a two-dimensional free energy surface for the reaction. In the Marcus–Saveant ´ picture, it is assumed that the solvent free energy responds linearly to changes in the charge distribution, which is equivalent to a quadratic dependence on the solvent polarization coordinate D e. Therefore, this solvent free energy is usually written as Fsolvent s FMeCl , min q
1 4l
2 Ž D e y D e0 . ,
Ž 5.
where FMeCl, min is the solvation free energy of the MeCl molecule near the interface, and l is the solvent reorganization energy, with D e 0 being the minimum of the solvent free energy curve. If the Marcus–Saveant solvent picture applies, these pa´ rameters can be computed from a parabolic fit to an MD simulation of the diabatic undissociated methyl chloride state, and can in turn be used to calculate the entire potential energy surface as prescribed by a modified version of the model of Koper and Voth w13x. The actual free energy surface computed from an explicit MD simulation with the BBET Hamiltonian is shown in Fig. 2. For comparison, the free energy surface computed from the analytic model of Ref. w13x using the linear response model for the solvent free energy ŽEq. Ž5.. is also shown. The free energy surface computed by computer simulation is a result of over 2700 umbrella sampling calculations. The parts of the free energy with the most error are near the edges of the surface. The important part of the surface, the barrier and well regions, were sampled the most. The errors here are ; 1 kcalrmol Žcomputed from the standard deviation for different averaging intervals of the free energy calculation.. The simulation free energy surface shown in the upper part of Fig. 2 represents 114,000 node hours of computation on parallel computers. This is equivalent to running a single present-day high-end workstation for 13 years. From the results in Fig. 2, it is seen that the gross qualitative features of the analytic model are similar to the actual free energy surface computed from the simulation. First, there is a bound well where the CH 3 Cl is stable and the electron has not been transferred from the electrode to the CH 3 Cl. There is also an exit channel, because the CH 3 Cly radical is not
98
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
Fig. 2. A contour plot of the adiabatic free energy surface computed from explicit MD simulation Žtop panel. and from the linear response model for the solvent Žbottom panel.. The well Župper left. and exit channel Žlower right. are clearly visible in both surfaces. The well corresponds to the neutral CH 3 Cl molecule, and the exit channel to the unstable charged radical. A solvent-induced barrier in the exit channel is also seen in the simulation data, occurring where the fragments are still in contact.
stable. Thus, the analytical free energy surface contains the essential physics of the BBET model: a stable CH 3 Cl species and an unstable CH 3 Cly radical. There is, however, at least one key qualitative difference between the free energy surfaces computed via computer simulation and the analytic model. In particular, there is a well in the exit channel of the simulation result where the methyl fragment and the chloride fragment are still in contact. This well, which is ; 7 kcalrmol in depth, is a result of explicit microscopic interactions in the MD model. When the fragments are in contact, or separated by a single water molecule, the result is a configuration with a lower energy than when the fragments are separated, but not separated enough to allow more water molecules in between the fragments. There is an energy ‘cost’ associated with making this hole in the solvent Žthe intermediate state., therefore the contact fragment pair and solvent
separated fragment pair are distinguished by a free energy barrier between them. Another key qualitative difference is a positive free energy of reaction Žoverpotential. in the simulation data which is not present in the analytic model for the given set of parameters. The contour plots in Fig. 2 also reveal important quantitatiÕe differences between the simulation results and the linear response model predictions. The curvature, minima, and reaction free energy Žoverpotential. differ considerably between the linear response solvent and the MD simulation data. This is due to the non-linear response of the solvent in the MD simulation as the CH 3 Cl bond is broken, which is in turn related to the change in the charge on the redox center from y0.25 in the CH 3 Cl molecule to y1 for the chloride ion, and to the change in the volume exclusion effect from the Me fragment at different values for r. The primary quantitative difference between the simulation results and the analytical prediction is a 10 kcalrmol positive overpotential for the former in forming first the metastable MeCly ion and then subsequently the Me and Cly fragments. The analytical model thus predicts a steady-state concentration of the metastable ionic intermediate Žand ultimately the products. which is too large by seÕen orders of magnitude. Interestingly, the rate constant Žactivation free energy. for the ET event from the electrode to the MeCl is predicted to be the roughly same by both simulation and theory. However, the thermodynamics Žequilibrium constant, etc.. of the reaction at the specified electrode conditions are found to differ greatly. It should be noted that the above comparison of the analytic model with the simulation results is based on the usage of a common electrode Fermi level in the underlying Hamiltonian ŽEq. Ž2.. which appears in both approaches. One might argue that a better strategy is to compare the prediction for the rate constant from the analytical model with that from the simulation at a common overpotential, i.e., to adopt a more phenomenological approach which has its footing in experimentally specified quantities. Along these lines, the expression for the activation free energy D Fact in the analytical model is given to a good approximation by w1,13x D Fact s
lqD 4
2
h
ž
1q
lqD
/
,
Ž 6.
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
99
Fig. 3. A plot of l, the reorganization energy, as a function of r, the CH 3 Cl dissociation coordinate. The reorganization energy is determined from global harmonic fits to the free energy surface at fixed values of r.
where h is the overpotential, l is the reorganization free energy, and D is the gas-phase dissociation energy of the MeCl bond. ŽIt should be noted that the sum l q D acts like an overall reorganization energy in this theory.. For the case of a 10 kcalrmol positive overpotential as observed in the simulation, the above expression predicts a kinetic barrier to ET of 42 kcalrmol, while the simulation barrier is found to be 37 kcalrmol. This difference in barrier corresponds to a difference of more than three orders of magnitude in the ET rate constant between the analytical prediction and the simulation result. Clearly, the fundamental quantitative differences between the simulation results and the analytical predictions do not disappear when viewed in this alternative light. When the MD simulation data are analyzed in some detail, it is found that in reality the effective reorganization energy l depends on the charge state of the system, which is different for the molecule and the fragment state, and on the bond length or inter-fragment distance r. Furthermore, a l determined from a fit locally about a minimum often gives a poor fit of the entire free energy curve as obtained from the umbrella sampling procedure. Therefore, an alternative method of extracting the
effective l from the MD simulations is to fit the entire free energy curve of one state, at least up to the transition state, to a parabola. The l obtained in this way is close to another definition of the solvent reorganization energy, namely as the vertical energy shift from the diabatic state under consideration to the minimum of the other diabatic state Žat equilibrium., as it effectively includes deviations from nonparabolicity of the overall diabatic energy curve. The solvent reorganization energies determined in this way, and their dependence of the bond length or inter-fragment distance r, for both the molecule and the fragment states, are shown in Fig. 3. The l as determined from a local fit about the methyl chloride minimum is also shown as a dashed line in Fig. 3 and is seen to give a poor description of the real solvent response.
4. Concluding remarks In conclusion, for the first time very large-scale MD simulations have been carried out for a bondbreaking electrochemical reaction in terms of two coupled reaction coordinates, the solvent fluctuation
100
A. Calhoun et al.r Chemical Physics Letters 305 (1999) 94–100
and molecular dissociation coordinates. The simulations were compared to an analytical free energy surface calculated using the linear response solvent model. Although the latter model is a convenient tool for estimating trends and gross features of the reaction, it fails to reproduce the molecular details that are found to be very important for bond-breaking reactions. Specifically, the simulations reveal a free energy well for the contact fragment pair due to the molecularity of the solvent Žthe solvent ‘cage effect’., and the solvent polarization fluctuations driving the bond-breaking process are seen to be significantly non-linear, as manifested by an effective solvent reorganization energy which is reaction coordinate dependent. The latter effects contribute to large quantitative differences between the simulation results and the predictions of the existing analytical theory. The present results thereby provide a clear and unique illustration of the indispensable role accurate and detailed computer simulations can play in obtaining a proper understanding of complicated electrochemical reactions.
Acknowledgements This research was supported by the United States Office of Naval Research and the Royal Netherlands Academy of Arts and Sciences ŽKNAW.. Computations were carried out at the Aeronautical Systems Center Major Shared Resource Center at WrightPatterson Air Force Base through a Department of Defense Grand Challenge Computing Grant, and at the National Center for Supercomputing Applications ŽNCSA. through NSF Metacenter grant number
MCA94P017P. We also gratefully acknowledge the University Utah Center for High Performance Computing for a generous grant of Silicon Graphics Origin 2000 time. References w1x J. Saveant, Acc. Chem. Res. 26 Ž1993. 455. ´ w2x J. Bertran, I. Gallardo, M. Moreno, J. Saveant, J. Am. Chem. ´ Soc. 114 Ž1992. 9576. w3x C.P. Andrieux, A.L. Gorande, J. Saveant, J. Am. Chem. Soc. ´ 114 Ž1992. 6892. w4x J. Saveant, J. Am. Chem. Soc. 109 Ž1987. 6788. ´ w5x X. Xia, M.L. Berkowitz, Chem. Phys. Lett. 227 Ž1994. 561. w6x D.A. Rose, I. Benjamin, Chem. Phys. Lett. 234 Ž1995. 209. w7x D.L. Price, J.W. Halley, J. Chem. Phys. 102 Ž1995. 6603. w8x J.B. Straus, G.A. Voth, J. Phys. Chem. 97 Ž1993. 7388. w9x J.B. Straus, A. Calhoun, G.A. Voth, J. Chem. Phys. 102 Ž1995. 529. w10x A. Calhoun, G.A. Voth, J. Phys. Chem. 100 Ž1996. 10746. w11x A. Calhoun, G.A. Voth, J. Electroanal. Chem. 450 Ž1998. 253. w12x A. Calhoun, G.A. Voth, J. Chem. Phys. 109 Ž1998. 4569. w13x M.T.M. Koper, G.A. Voth, Chem. Phys. Lett. 282 Ž1998. 100. w14x M.T.M. Koper, G.A. Voth, J. Chem. Phys. 109 Ž1998. 1991. w15x P.W. Anderson, Phys. Rev. 124 Ž1961. 41. w16x D.M. Newns, Phys. Rev. 178 Ž1969. 1123. w17x W. Schmickler, J. Electroanal. Chem. 204 Ž1986. 31. w18x A. Kornyshev, W. Schmickler, J. Electroanal. Chem. 185 Ž1985. 253. w19x K.L. Sebastian, J. Chem. Phys. 90 Ž1989. 5056. w20x M.T.M. Koper, W. Schmickler, Chem. Phys. 211 Ž1996. 123. w21x M.T.M. Koper, J. Phys. Chem. 101 Ž1997. 3168. w22x M.T.M. Koper, J.H. Mohr, W. Schmickler, Chem. Phys. 220 Ž1997. 95. w23x B. Smith, J.T. Hynes, J. Chem. Phys. 99 Ž1993. 6517. w24x R.A. Marcus, J. Chem. Phys. 43 Ž1965. 679. w25x W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 Ž1983. 926.