A Model for Proton Transfer to Metal Electrodes - American Chemical

10814

J. Phys. Chem. C 2008, 112, 10814–10826

A Model for Proton Transfer to Metal Electrodes Florian Wilhelm,† Wolfgang Schmickler,*,† Renat R. Nazmutdinov,‡ and Eckhard Spohr§ Institute of Theoretical Chemistry, UniVersity of Ulm, D-89069 Ulm, Germany, Kazan State Technological UniVersity, 420015 Kazan, Russian Federation, and Lehrstuhl fu¨r Theoretische Chemie, UniVersita¨t Duisburg-Essen, D-45141 Essen, Germany ReceiVed: January 16, 2008; ReVised Manuscript ReceiVed: April 28, 2008

An empirical valence-bond (EVB) model is developed to describe the transfer of a proton from aqueous solution bulk to a metal electrode surface. The results of density functional calculations are used for parametrizing the model for a Pt(111) surface, which was chosen as a model system. Thereby, the metal surface is addressed in terms of the cluster model. The developed EVB model makes it possible to perform large-scale molecular dynamics (MD) simulations for a metal/electrolyte solution interface. Using MD we explore the proton transfer to a Pt(111) electrode surface with different negative surface charge densities. Preliminary conclusions are reported to elucidate some aspects of the proton transfer mechanism. 1. Introduction The transfer of a proton to a metal electrode is a fundamental process in electrochemistry, but despite decades of research, it is still much less understood than the transfer of other ions. The main difficulty is the strong interaction of the proton with water, with which it can form transitory species such as the Eigen complex H9O4+ (a solvated hydronium H3O+ ion)1 or the Zundel ion H5O2+.2 Indeed, in aqueous solutions the proton is not a well-defined species, because it is rapidly exchanged with the hydrogen atoms of the water molecule. To distinguish these forms using, for exampel, spectroscopy techniques is very complicated from the experimental viewpoint (see, e.g., ref 3). For this reason, some evidence in favor of both the Eigen and the Zundel states for a proton in aqueous solution was reported in experimental literature.4–7 The basic idea of the so-called Grotthuss hopping mechanism8 for proton transfer in bulk water has been known for a long time. Recently, an improved and more detailed understanding of the process in bulk water has been gained using state-ofthe-art molecular dynamics (MD) techniques, which have been developed during the past decade.9–20 These simulations provide very valuable insight into the nature of proton transfer and have helped to elucidate the rate-determining step of the reaction, the breaking of a hydrogen bond in the solvation shell of the proton-transferring cluster (Moses mechanism).21 However, several questions still remain open in bulk water,12,18,19,22–24 and even for advanced ab initio MD simulation techniques substantial deviations from experimental values exist, in addition to problems related to the functionals and methods.12,25 What can be said with some confidence is that the Zundel and Eigen species should rather be regarded as limiting forms22 during the proton transfer, and that the barrier for the process itself is insignificantly small.13,22 It might seem attractive to also use ab initio MD for the investigation of proton transfer to an electrode surface. However, a quick estimate shows that this is not a realistic option. Because * Corresponding author phone: +49731/50-31340; fax: +49731/5022819; email: [email protected]. † University of Ulm. ‡ Kazan State Technological University. § Universität Duisburg-Essen.

of their heavy computational demands, such simulations cover, at best, a range of several picoseconds and an electrode area of perhaps a few square nanometer. Even on a good catalyst such as platinum, transfer of a proton to the metal surface—as opposed to within the water layer, compare ref 26—onto such a small area during such a short time is highly unlikely. The few attempts that have been made in these directions (e.g., refs27 and 28) have therefore been performed at comparatively high electric fields. Although the method employed appears promising for describing the electronic properties of the final reaction step, achieving a realistic characterization of the overall reaction dynamics comprising statistical averaging is still far from accomplishable. Another approach to describe electrochemical reactions that has been extensively used during the last few years29–36 makes use of state-of-the-art density functional theory (DFT) techniques but neglects the finite temperatures at which the electrochemical reactions take place. The water involved is thereby either not explicitly treated because it is estimated that it is not crucial for the rate-determining step of the reaction or in some works it is modeled in an ice-like, often bilayer structure. The latter results from geometry optimization and is comparable to corresponding low-temperature structures found in experiment.37 Undoubtedly, for the chosen structure this method gives accurate results as far as the potential energy changes during electrochemical reactions are concerned. But on the other hand, the solvent dynamics has to be addressed to describe such processes in a comprehensive way. Classical MD, which are based on model potentials and are computationally less demanding, normally cannot describe the breaking of bonds and are therefore also unsuitable for the investigation of electrochemical proton transfer, in which the breaking of hydrogen bonds must play a major role. To overcome the shortcomings of classical MD, the empirical valence-bond (EVB) method, based on the ideas of Pauling and Coulson,38 has been developed. This method was pioneered by Warshel (e.g., refs 39–41) and later adopted for proton transfer in water.42,43 Further development of the method toward a multistate description14–17,20,22,44 has enabled the construction of appropriate models13,18,19 and has led to deeper insight into the details of the mechanism of proton transfer.24

10.1021/jp800414f CCC: $40.75  2008 American Chemical Society Published on Web 06/27/2008

A Model for Proton Transfer to Metal Electrodes

J. Phys. Chem. C, Vol. 112, No. 29, 2008 10815

Because of their relative simplicity, EVB methods allow simulations over a much longer time scale than ab initio methods and are therefore attractive for the simulation of electrochemical proton transfer. However, existing implementations have been designed for homogeneous transfer. We have, therefore, adapted the EVB method to the electrochemical situation. The most challenging part is the determination of the terms that couple the proton to the metal, that is, the Pt(111) electrode we have chosen as a first example. These we have determined by quantum chemical methods. In the following, we present our model, and the calculations on which it is based, in detail. By necessity, the larger part of this article is technical. Although we realize that reading all the details of our model may not be the best way to understand the elementary ideas, for the sake of scientific transparency it is mandatory to present them. We advise the reader who just wants to gain a basic understanding to read Sections 2 and 2.1, which present the concepts, and then turn to Section 4, which shows first results: typical trajectories for a proton transfer to a platinum electrode. The results of systematic MD simulations will be presented in a separate communication. 2. Empirical Valence Bond Model for Proton Transfer Before we present our model, we briefly review the basic features of the EVB model as developed by Warshel.39–41 The original valence bond method for the calculation of the groundstate of a molecular system is based on the approach, which is natural for a chemist, to describe the electronic ground state (resonance hybrid) |ψ〉 as a linear combination (with coefficients ci) of several resonance structures |i〉.

|ψ 〉 ) ∑ ci|i〉

(1)

i

Employing the variation principle of quantum mechanics and the appropriate mathematical minimization techniques yields the corresponding coefficient vector c ) (ci) by the solution of a generalized eigenvalue problem, that is,

Hc ) ScE,

(2)

ˆ |j〉 are the matrix elements of the VBwhere Hij ) 〈i|H Hamiltonian, and Sij ) 〈i|j〉 are the overlap matrix elements. Often, the states |i〉 are assumed to be orthonormal, so that Sij ) δij, and the Hamiltonian is assumed to be symmetric. This leads to a simpler eigenvalue problem to be solved (eq 3).

Hc ) Ec.

(3)

Up to this point the method itself in principle allows the exact solution of the quantum-mechanical problem but involves, depending on the nature of the problem, finding possibly numerous states |i〉 corresponding to resonance structures and the calculation of a large number of quantum-mechanical integrals before solving the eigenvalue problem. Therefore, to make it applicable for the calculation of the potential energy surface (PES) during MD simulations, which demands efficient computational algorithms and the availability of derivatives for force calculation, some simplifications and empirical approximations for the quantum mechanical quantities in eq 3 have been introduced. First of all, the resonance structures, the corresponding VB states of which are called diabatic states in current EVB literature, are chosen according to chemical reasoning, that is, only structures that are expected to give a substantial contribution to the resonance hybrid (adiabatic state) are included in the model. Second, the elements of the VB-Hamiltonian are

represented by empirical expressions. A common approach is to represent the diagonal elements Hii by the energy values that an appropriate force field gives for the bonding situation of the corresponding resonance structure. The off-diagonal elements are often described by analytical functions fitted to either experimental information on barrier heights, obtained, for example, in the vacuum or by ab initio calculations for small systems. Consequently, the Hamiltonian can be described as a function just of the nuclear coordinates of the system with, to some extent, an implicit inclusion of the quantum-mechanical nature of the electronic coordinates. The lowest eigenvalue resulting from the solution of eq 3 gives the expectation value for the ground-state (or adiabatic) energy E0. The adiabatic state can then be calculated using the corresponding eigenvector cg ) (cgi ) (eq 4).

|ψ 〉 ) ∑ cig|i〉

(4)

i

This gives the squares of the ground-state eigenvector coefficients, |cgi |2 the meaning of weights by which the diabatic state |i〉 contributes to the adiabatic state. Furthermore, by employing the Hellmann-Feynman theorem, one can calculate all forces necessary for MD. 2.1. EVB Model for Proton Transfer to a Metal Electrode. As already mentioned in the introduction, the EVB method has been, inter alia, employed for simulations of proton transfer in the aqueous phase by several groups. Early approaches using two EVB states42,43 have been extended to multistate EVB models,14–17,20,22 which are capable of a very detailed description of proton transfer in aqueous environments. Walbran and Kornyshev44 presented a two-state model that has the advantage that it allows long simulations for many protons embedded in large systems, as are necessary for, for example, the simulation of proton transport in polymer electrolyte membranes.45 They also presented an extension of the Toukan-Rahman46 water model, which adds polarizability to the water molecules. This model does suffer somewhat from the limited size of its proton transfer cluster, that is, the entity of all atoms that are treated by the EVB method (see discussion in ref 23), but for many applications this shortcoming may not present a serious problem. We have chosen to use the Walbran-Kornyshev model together with their water model as a basis for our extension of proton transfer to the metal electrode mainly for two reasons: the low computational cost compared to the multistate approaches, and the simplicity of the model, which makes it rather easy to choose a candidate hydrogen atom for the transfer to the surface. As is usually done in simulations with the EVB method, we single out a set of atoms within the system that constitute the transfer cluster. These atoms are dynamically assigned during the simulation, which requires a careful choice of the interactions and the cutoff radii to fulfill energy conservation. The transfer cluster thus can (in principle) and does (in practice) change its composition while the proton diffuses through the bulk. This transfer cluster, see Figure 1, consists of two water molecules and the excess proton, which together form either a Zundel complex H5O2+ or a hydronium ion H3O+ and its best water neighbor for proton transfer within the solvation shell in the aqueous phase. In our scheme, the complex in this state can be regarded as an approximation for the Eigen ion H9O4+. Multistate (MS) EVB simulations reveal that, for the latter limiting species, the positive charge is delocalized over the entire proton complex, that is, the charge on the central hydronium ion only amounts to 0.6 to 0.7.19,20,22 In the present model we approximate the real situation by delocalizing the charge over

10816 J. Phys. Chem. C, Vol. 112, No. 29, 2008

Wilhelm et al.

Figure 1. Diabatic states used in our EVB model; note that the geometry is the same for all pictures, just the bonding situation changes.

only two water molecules. These species are subject to a dynamical interchange during the simulation and can be regarded as limiting cases for proton transport by the Grotthuss mechanism. Additionally, seven metal atoms forming a hexagon are included in the special treatment. They and the “candidate” hydrogen atom, that is, the one that is possibly transferred to the surface, are chosen by a distance criterion; out of the four hydrogen atoms that do not form a hydrogen bond within the (formal) Zundel ion, the one (yellow) that is next to the metal surface is chosen. Then its nearest neighbor on the surface is labeled the “pivot” surface atom (green). The latter is then—together with its six next surface neighbors (light blue)—included in the transfer cluster. Here we deliberately strive to develop a simple model for proton transfer to an ideal surface. The surface has no defects or steps, and other oxygen species, such as oxide, adsorbed OH or OH-, which may be present as dissociation products of water in certain regimes of electrode charge, are neglected, although they may play a role in proton transfer (see, e.g., ref 26). Because these species bear negative partial charges, they might even accelerate proton transfer from the solution bulk to the metal surface. Also, neither different single-crystal surfaces nor surface compounds, alloys, or adatoms in the vicinity of the proton transfer sites are included in the model. Although an extension of the model to include some of these features is, in principle, possible, efforts needed for adjusting the transfer cluster composition, model mechanics, and parametrization would be substantial. Consequently, including such effects, which cannot be neglected for a full description of real electrocatalysis, is very clearly beyond the scope of the current work. 2.1.1. Diabatic States. The choice of the diabatic states in our EVB model is motivated by chemical intuition. The bonding

situation for several of these states is shown in Figure 1. The first two states, Figure 1, panels a and b, form the basis for proton transfer within the aqueous phase. The intra- and intermolecular interactions for these two bonding situations, which determine the matrix elements H11 and H22 of the Hamiltonian, are described by analytical potential functions of the type commonly used in MD. Their functional form and the parameters for the Zundel part are explained in detail in ref 44. The interaction of the hydrogen and oxygen atoms with the metal surface is described by the empirical potential reported in ref 47 for atoms belonging to a water molecule in the corresponding state, whereas for the hydronium part, the oxygen-platinum potential has different parameters (eq 5),

VO-Pt ) [13720.874 exp(-1.3061r) 13526.816 exp(-1.2999r)]f(F) + 7240948.8 exp(-5.3568r)[1 - f(F)] (5) with, as in the original paper, f(F) given by eq 5a.

f(F) ) exp(-0.5208F2)

(5a)

Here, r denotes the distance of the oxygen atom to the metal atom, and F is the length of the projection of the distance vector onto the surface plane. All pre-exponential factors are in units of eV, all distances are given in Ångstro¨ms. The parameters are chosen to reproduce the energy values resulting from DFT calculations on Pt clusters, see Section 3.1 below. The hydrogen-platinum potential is the same as for the case of water.47 Far from the surface, where only these first two diabatic states contribute to the adiabatic state, the model shows exactly the same behavior as the one from ref 44.



TABLE 1: Parameters of Potential Functions VH-Mja ε1 ω req ε2 a s db de

VH-Mia

5.12 eV 2464.77 cm -1 1.505 Å 0.09 eV 6.9 Å -1 0.22 1.6 Å 5.0 Å

εMi ωb rMi

b

VH-O c

εO bO rO

5.12 eV 2464.77 cm -1 1.505 Å

(

(

VH-Mj ) ε1 1 - exp -

))

2

µ ω(r - (req + sFH/)) - ε1 + 2ε ∆H(rz)ε2(1 - exp(-aFH/))2 (6)

where, with u ) (rz - db)/(de - db), we have eq 7.

{

rz e db 1 2 2 ∆H(rz) ) (u - 1) db < rz < de rz g de 0

{

εX(1 - exp(-bX(r - rX)))2 r < rX r g rX 0

b

4.75 eV 0.49 Å -1 1.75 Å

bMi in eq 8 is calculated using these values

TABLE 2: Parameters for the Off-diagonal Hamiltonian Elements, i, j ) 3,..., 9 H12 a,b

Λ0 δ0b

-3.1 eV 0.55 Å

Hij Λ1

-0.145 eV

rmax, 1 rlim, 1

3.8 Å 4.5 Å

H2i Λ2 ι1 rmax 2 rlim, 2 γ ι2 ι3 ι4 R0 βang cH cO

-2.71 eV 0.525 5.05 Å 7.9 Å 5.8824 0.68 Å 0.5161 0.091 116.0° 32.5° 1.3 Å 1.3 Å

a With positive sign in the original paper; unlike here, the sign does not matter with the two-state model there. b From ref 44.

(7)

Here, r denotes the interatomic H/-metal distance, FH/ is the length of the projection of the H/-metal vector onto the (111) surface, and rz is the length of the component perpendicular to the surface. The first two terms in eq 6 introduce a Morse potential with its minimum lying at -ε1 for r ) req and FH/ ) 0, that is, in the on-top position. The next term, together with the off-diagonal matrix elements in eq 13, is designed to reproduce the corrugation of the relatively flat PES for a single hydrogen atom on platinum, as can be found in the literature; for a discussion see Section 2.1.2 below. Furthermore, µ denotes the reduced mass, ω is the vibrational frequency perpendicular to the surface, and req is the equilibrium distance. The values of the parameters can be found in Table 1. ε1, ω, and req were adjusted to reproduce the energy values of DFT calculations on Pt clusters as described in Section 3.1. To make the description of the non-Coulomb interatomic forces complete, one also needs to specfiy the interactions of the hydrogen atom H/ with the next oxygen atom and the remaining hydrogen and metal atoms from the transfer cluster. As H/ in these states belongs neither to a hydronium ion nor to a water molecule, these interactions are not covered by the potentials mentioned above. They are assumed to be purely repulsive, consequently they are modeled by the repulsive branches of Morse-type potentials (eq 8),

VH-X )

εH bH rH

11.55 eV 1.285 Å -1 0.98 Å

a Mj, metal atom to which a bond exists in the corresponding state; Mi, other metal atom. according to the formula within eq 6. c Values are from ref 44.

The remaining states, |3〉 to |9〉, describe a bonding situation with the candidate hydrogen atom H/ being bonded to one of the surface atoms plus the two neighboring water molecules. Note that in these states the former excess proton becomes just a normal hydrogen atom bonded within a water molecule. Again, the intrawater, interwater, and water-metal interactions are modeled by potential functions as described in refs 44 and 47, respectively. The bond of the (now neutral) single hydrogen atom H/ to a metal atom Mj corresponding to the state |j〉, with 3 e j e 9, is represented by a modified Morse-type potential function (eq 6),

VH-H

(8)

where r denotes the interatomic distance of H/ to X ) O, H, Mi, i ) 3... 0.9, i + j for state |j〉, and rX is the corresponding equilibrium distance. For the parameter values see Table 1.

Figure 2. Geometries for the one-dimensional model processes.

State |3〉 (Figure 1c) features the bond of the hydrogen atom to the pivot metal atom, each of the remaining states has a bond to one of its neighbors; state |4〉 is shown as an example (Figure 1d) for six equivalent states. Note that the states |3〉 to |9〉 represent neutral species, so not only bond breaking, but also discharge takes place during transfer to the surface. Table 3 shows a comparison of the charges corresponding to the two types of states. It should be noted that the total effective (i.e., given by their contribution to the adiabatic state) net charge of at most +1 carried by the first two states is neutralized by a negative countercharge of equal magnitude, which is distributed onto the bottom metal layer, giving a crude representation of an image charge. This negative countercharge varies during the simulation according to the nature of the instantaneous adiabatic state, so that in our scheme the total charge is conserved during proton transfer. There is no interplatinum potential present in this study. All metal atoms are kept fixed at their ideal lattice positions. Of course, this is an approximation, but one should not expect major movements of the metal atoms because of their high mass compared to those of the solution species; also, no surface reconstruction is expected for platinum.


Wilhelm et al.

We are now able to calculate the energies within the transfer cluster and the nonbonding interactions to the remaining system, which can be split into two contributions: the non-Coulomb interactions, which are the same for all states and again are calculated according to the force fields from refs 44 and 47, and the Coulomb interactions, which depend on the charge distributions in states |1〉 and |2〉 and |3〉 to |9〉, and which have to be calculated separately; for this purpose the Ewald scheme is used. These energy terms go into H11 to H99, the forces for each state are calculated in the usual manner. See Table 3 for charge values and Section 2.1.3 for details. 2.1.2. Off-diagonal Hamiltonian Elements. What is still missing in our EVB-Hamiltonian are the off-diagonal elements or couplings. Many of them are set to zero according to the following rule: all off-diagonal elements Hij that belong to proton transfer events over more than one atom (i.e., interchanges between states i and j which cannot be described by a single proton transfer, compare ref 14), for exampel, between state |4〉 and |6〉, are set to zero. The states are consecutively numbered in a way that physically neighboring metal atoms correspond to bonding situations in neighboring states. This finally gives the following form for our EVB Hamiltonian (eq 9),

(

H11 H12 0 0 H) 0 0 0 0 0

H12 H22 H23 H24 H25 H26 H27 H28 H29

0 H23 H33 H34 H35 H36 H37 H38 H39

0 H24 H34 H44 H45 0 0 0 H49

0 H25 H35 H45 H55 H56 0 0 0

0 H26 H36 0 H56 H66 H67 0 0

0 H27 H37 0 0 H67 H77 H78 0

0 H28 H38 0 0 0 H78 H88 H89

0 H29 H39 H49 0 0 0 H89 H99

)

(9)

where all elements have been assumed to be real. The offdiagonal elements can be divided into three groups according to the processes that they promote. The first group contains H12, which is mainly responsible for proton transfer in the aqueous phase. For the labeling of the involved atoms, we use O1 and O2 for the two oxygen atoms within the transfer cluster, O1 being nearest to the bridging hydrogen atom HB. HB is the hydrogen atom out of the (formal) Zundel ion H5O2+ which mediates the formation of the hydrogen bond between the two (formal) H2O groups. This coupling function is purely symmetry-dependent, and its functional form and parameters are taken unchanged from ref 44, where the details are described. The resulting formula is shown in eq 10,

[( ) ]

H12 ) Λ0

Q0 2 -1 δ0

2

(10)

where

Q0 ) rO1-HB - rO2-HB

(11)

is the symmetry variable, and Λ0 and δ0 are empirical parameters. The other coupling functions are, besides their respective chemical function within the model, designed to be short-range, which allows us to keep the transfer cluster of a finite and constant size, whereas on the other hand dynamical changes of the composition remain possible. The second group couples two different states out of |3〉 to |9〉, and is mainly responsible for hydrogen diffusion on the

metal surface, and is also important for hydrogen adsorption on all nononefold sites (bridge, hollow) on the surface. The atomic labels relevant for this type of function are H/, the hydrogen atom out of the transfer cluster that is the candidate for proton transfer to the metal surface, or in case of a completed transfer, the neutral hydrogen atom adsorbed at the surface; Mi denotes the appropriate metal atom to which the bond from H/ would go to in one of the states |i〉 out of |3〉 to |9〉. The variable used here is

rsum,1 ) rH/-Mi + rH/-Mj

(12)

so we can write the coupling function in the form shown in eq 13,

Hij ) Λ1 · (x2 - 1)2

(13)

where

x(rsum,1) )

rsum,1 - rmax,1 . rlim,1 - rmax,1

(14)

For rsum, 1 > rlim, 1 or rsum, 1 < rmax, 1-(rlim, 1 - rmax, 1) we set Hij ) 0. The empirical parameters used here are the coupling strength (Λ1), the distance rmax, 1 denoting the value of rsum, 1 where the maximum of Hij is reached, and rlim, 1 denoting the “outer bound” of the function. This functional form guarantees not only that the value of Hij reaches zero when passing rlim, 1, but also that all derivatives become zero, thus avoiding force truncation and accompanying problems with energy conservation. There are quite a number of articles in the literature dealing with hydrogen adsorption on Pt(111), for example refs 36 and 48–52. Although there is some agreement that the barrier for atomic hydrogen migration on the surface is quite low ( cH ;

FO/ > cO ;

rsum,2 > rlim,2 ; rlim,2;

rsum,2 < 2rmax,2 -

|Q/| > 1;

|Qang,2,j| > 1

The parameters for the last type of coupling function are Λ2, γ, ι1, ι2, ι3, ι4, rlim, 2, rmax, 2, R0, βang, cH and cO, the respective values of which are given in Table 2. 2.1.3. Charges and Forces. It should be noted that the water and hydronium models of Walbran and Kornyshev mainly

atom

|1〉, |2〉

|3〉 to |9〉

HH3O OH3O+b PH3O+ HH2O OH2Ob PH2O H/ Ptd Pte

0.33 0.67 -0.66 0.33 0.00 -0.66 0.33c 0.00 -1.00f,g

0.33 0.00 -0.66 0.00 0.00 0.00g

+

(15)

is the symmetry variable analogous to eq 11; eq 16 again represents the influence of the magnitude of distances;

Qang,2,j )

TABLE 3: Formal Charges for the Diabatic States in Atomic Unitsa

a Values for water and the hydronium ion taken from ref 44. Value for pure states only, see text. c Here included either in formal H3O+ or H2O. d Belonging to transfer cluster. e Total charge belonging to bottom surface layer; division by the number of atoms per layer (64 for the performed simulations) yields the individual atom charge. f Effective charge value depends on adiabatic state. g Plus additional surface charge for electrode potential other than pzc. b

extend the original flexible models by Toukan and Rahman46 and Schmitt and Voth,16 respectively, by adding polarizability. To achieve this, a polarization site P coupled to the oxygen atom is introduced, which aims solely to provide electronic polarizability, that is, it has no chemical interactions with the remaining system apart from Coulomb interactions. This pseudoatom has a very low atomic mass of 0.2 au and carries the countercharge for the charges of the two hydrogen atoms in the water model, see Table 3. Its interaction with the oxygen atom is described by a (2,4)-potential, for details see ref 44. Within the diabatic states |1〉 and |2〉, charges are treated in principle as in the original work, that is, the delocalized nature of the excess proton is accounted for by distributing the excess charge between the two oxygen atoms of the transfer complex according to the value of the proton transfer coordinate Q0 from eq 11. A detailed discussion can be found in the original paper, but for convenience a short description shall be given here. Using the atom labels defined above, the variable charges q˜ on the oxygen centers are mixed according to

˜ qO1 ) (1 - f(Q0))qOH O+ + f(Q0)qOH O 3

2

˜ qO2 ) (1 - f(Q0))qOH O + f(Q0)qOH O+ (23) 2

3

where the qX values are the pure state values from table 3 and f(Q0) serves as a charge switching function. It is defined by eq 24,

{

0 Q0 1 ζ Q0 5 2ζ Q0 3 +ζ + f(Q0) ) 2 5 χ0 3 χ0 χ0 1

( ) ( ) ( )

Q0 < -χ0 -χ0 e Q0 e χ0 Q0 > χ0 (24)

and f varies smoothly between 0 and 1 for -χ0 e Q0 e χ0. For a plot see the article by Walbran and Kornyshev;44 note the different labeling of variables and parameters in that article. The values of the parameters are χ0 ) 0.5 Å and ζ ) 15/16. As mentioned, the net charge of +1 for the transfer cluster is neutralized by a countercharge of equal magnitude, which is uniformly distributed over the metal atoms of the bottom metal layer. Within the diabatic states |3〉 to |9〉, the electronic state of the transfer cluster is neutral. The remaining point charges are the same for all seven states, their values are given in Table 3.


Wilhelm et al.

After calculating the adiabatic energy E0 and the groundstate eigenvector cg, the adiabatic charge distribution is calculated according to eq 25, 9

qk )

∑ (cig)2qik

(25)

i)1

where k denotes an arbitrary site within the transfer cluster; qik is the diabatic charge for this site in state |i〉 and can be found in Table 3. In the same manner, an effective countercharge is calculated for the metal atoms of the bottom metal layer from the adiabatic excess charge. As already mentioned, the forces acting on the particles can be calculated, using the Hellmann-Feynman theorem, from the derivatives of all diagonal and off-diagonal elements of the EVB Hamiltonian with respect to every nuclear degree of freedom involved. Explicitly, consider xl denoting such a degree of freedom, then we can write, using the notation introduced above;

Fxl ) -

〈

〉

9 ∂E0 ∂Hij ∂H ) - ψ| |ψ ) cigcjg ∂xl ∂xl ∂xl i,j)1

∑

(26)

Atoms and molecules that are not part of the EVB complex interact using the same force field functions as those presented above for the constituents of the complex, that is, the O-Pt, H-Pt, water-water, water-hydronium and intrawater interactions are the same as for the underlying models,44,47 and all Pt atoms are kept in fixed positions. 3. Fitting the Model Parameters A large part of the potential parameters used in this work are taken from standard models and previous works as described in the last section. But for most of the off-diagonal Hamiltonian elements and for some potential functions, there is, to the best of our knowledge, neither previous work to build on nor applicable experimental data nor first-principles calculations available that could be used for parameter fitting. Some DFT-based data on the reaction has been published in literature,53 but the authors of the latter work have chosen their proton transfer pathways based on different premises than in the present work. Hence, we decided to perform density functional theory (DFT) calculations on very small model systems and processes in vacuum, which enable us to determine our parameter values, while the influence of the electrochemical environment, that is, the solution, is included by the various interaction potentials within the MD simulation. Of course, this assumes pairwise additivity of the interactions of the transfer cluster and the surroundings. 3.1. Density Functional Calculations. We have chosen to use a small platinum cluster as a model for the Pt(111) surface, mainly because this setup allows for the treatment of charged systems, without the difficulties that a periodic model would entail, and at a modest computational cost. Due to the fact that treating platinum is still computationally very demanding, we had to restrict ourselves to one-layer clusters consisting of 7 or 12 atoms, although we are of course aware that for a truly precise modeling of the surface one should consider at least three-layer clusters consisting of a substantially larger number of atoms. On the other hand, in recent literature one can find even smaller systems than ours that are regarded as an appropriate model for reactions on a platinum surface, see, for example, ref 54. The energetic levels of the highest occupied molecular orbitals of our clusters (Pt7 ) -5.99 eV, Pt12 ) -6.02 eV) are found to be in reasonable agreement with the experimental value of

TABLE 4: Adsorption energy of a Hydronium Ion (∆Eads) and Equilibrium Platinum-H/ Distance (r(Pt-H/)) Calculated for a Pt10 (7 + 3) Cluster Using Different Combinations of Basis Sets for Platinum, Oxygen, and Hydrogen Atoms (Pt/O,H)a basis set

∆Eads [eV]

r(Pt-H/) [Å]

LanL2DZ/6-31 g(d, p) LanL2DZ/6-31 g+(d, p) LanL2DZ/6-31 g++(d, p) LanL2DZ/6-31 g++(d, p) LanL2MB/6-31 g(d, p) LanL2MB/6-31+g(d, p)

-1.79 -1.91 -1.97 -2.04 -1.66 -1.99

1.563 1.546 1.547 1.545 1.607 1.557

a All geometries are optimized in on-top position with respect to the hydronium part, cmp. Figure 2a.

the work function of 5.7 eV for Pt(111).55 This indicates that such clusters may serve as reasonable models for bulk electrodes. The model processes described in the next sections were chosen according to chemical intuition in order to gain insight into the dependence of the reaction PES on parameters such as the distance from the surface and the orientation of the molecules. The analytical form and the parameters of the offdiagonal coupling elements were then determined to reproduce the main features of this PES when applied within a force field in the absence of the water layer, that is, under the same vacuum conditions. The quantum chemical calculations were performed at the DFT level with the hybrid functional B3LYP56 as implemented in the Gaussian 03 program suite.57 A basis set of DZ quality was employed to describe the outer electronic shell of Pt atom (5s 5p 5d 6s 6p), whereas the effect of inner electrons was included in a relativistic effective core potential (LanL2) developed by Hay and Wadt.58–60 The standard basis set 6-31 g(d, p)61,62 was used to describe the electrons of O and H atoms. A singlet spin state for the platinum clusters was implied, since we aim at the modeling of the metal surface. To ensure an appropriate description of our system, several test calculations were done employing other basis sets. As can be seen from Table 4, more extended basis sets do not result in any noticeable changes of geometry or adsorption energy values. Therefore, in all further calculations we employed the less computer time-consuming combination LanL2DZ/6-31 g(d, p). To find suitable starting geometries for the model processes, a hydronium ion was optimized on a seven-atom platinum cluster starting from several different initial geometries. The most favorable geometry has one hydrogen atom pointing down, see Figure 2a. In the second most favorable geometry the hydronium lies flat on the surface. Depending on whether the hydrogen atoms point to hollow sites, see Figure 4e, or the neighboring metal atoms, see Figure 4, panels a and b, one finds the total energy being higher by 0.43 and 0.52 eV, respectively. All geometries have in common that the adsorption site is ontop; this seems to be the most favorable adsorption site for the hydronium ion. The three optimized geometries were chosen as starting points for the different model processes described in the following. During all calculations, the metal atoms were kept fixed at their bulk equilibrium nearest-neighbor distance of 2.77 Å. 3.2. One-dimensional Model Processes. To investigate the distance dependence of the energy barrier for proton transfer, the following one-dimensional process was introduced: Starting with a hydronium ion featuring one hydrogen atom pointing toward the surface, one follows the energetic path of the



Figure 3. One-dimensional model processes: (a-c) comparison of EVB to DFT energies for different O-Pt distances, and (d) coupling function for “perpendicular” proton transfer in on-top position for different O-Pt distances (in Å).

Figure 4. Geometries for the two-dimensional model processes, the arrows mark the coordinates for the PES plots.

elongation of the bond from the oxygen atom to this hydrogen atom until bond breaking occurs and a bond to the central metal atom is formed. This was done for different oxygen atom-surface distances ranging from 2.5 Å to 9.0 Å. As an example, Figure 2 shows the corresponding initial and final geometries for a distance of 4.0 Å. Along the reaction pathway, the coordinates of the remaining two hydrogen atoms were allowed to relax while keeping the positions of the oxygen atom and the transferring hydrogen atom fixed.

This process corresponds to an optimum orientation for proton transfer from the liquid phase to the metal surface. The excess proton is, according to the Grotthuss shuttling mechanism, transferred from the bulk to the first water layer facing the surface, where the recipient water molecule is already oriented, or then orients itself in such a way that one hydrogen atom points toward the surface. Three representative energy curves are shown in Figure 3 together with the values obtained using the EVB model for the


Wilhelm et al.

Figure 5. Potential energy surfaces for the 2D model processes: (a-c) proton transfer to nearest on-top position, and (d, e) proton transfer to hollow position.

same situations. There are some remarkable features in these plots; Figure 3a nearly corresponds to the equilibrium oxygen atom-surface distance (of 3.08 Å) for this configuration. As one can see, there is no energetic barrier but only a single minimum. For the equilibrium oxygen atom-surface distance, optimization to the most favorable geometry as mentioned above yields a nearest Pt-H/ distance of 1.6 Å and a O/-H/ distance of 1.5 Å. The Mulliken charge analysis furthermore shows a noticeable partial charge transfer to the metal cluster (qmetal ) 0.485 au). This means that an agostic bond is formed between the three atoms corresponding to the common minimum. The situation changes at longer oxygen-platinum distances. Two pronounced minima emerge, separated by an energy barrier. Figure 3b shows this situation for a distance of 4.0 Å. Here, we can clearly see the barrier, and the minimum corresponding to the initial state is shifted up by more than 0.5 eV compared to the final state, the latter also being energetically higher than for 3.0 Å. These trends continue for even longer oxygen-platinum distances, with the barrier increasing to more than 2.5 eV. From approximately 6.0 Å on, the energy level of

the initial minimum corresponding to the adsorbed hydronium ion remains roughly constant. Figure 3c shows an example for a larger distance. As the initial energy minimum and the barrier remain nearly constant, we have chosen a distance of 9.0 Å (not shown) for the energy curve used for fitting some of the parameters of the hydrogen-platinum potential, see above. One can conclude that, for the optimum distance from the surface and the downward orientation investigated in this section, proton transfer is a barrierless process, whereas for larger distances we see a rapid growth of the barrier. The latter trend can be compared to recent results63 for a similar process on mercury, with two notable differences: the initial state used in the present work consists of a pure hydronium ion facing the metal surface, whereas the authors of ref 63 used two additional water molecules to model the solvent, which of course changes the energy level of the initial state. In the present work, solvent effects are included only within the MD simulations. Second, the chemical difference between the two metals can be seen in the absence of a single-well PES at the equilibrium distance, that is, on mercury no agostic interactions were found.

A Model for Proton Transfer to Metal Electrodes 3.3. Two-dimensional Model Processes. The EVB model needs to give a realistic behavior for all possible orientations of the molecules within the water layer facing the surface. As it is not feasible to treat all possible orientations by quantum methods, one has to decide for the most important orientations. Thus, we investigated two more possible orientations, with the hydronium ion lying flat on the surface. Although not ideal for proton transfer, this orientation enables the hydronium ion to build up to three donor hydrogen bonds to neighboring water molecules for stabilization. As already mentioned, there are two flat on-top orientations of slightly different energy: one featuring the hydrogen atoms pointing to the hollow sites, and one turned by 30° with the hydrogen atoms pointing to the neighboring on-top sites. In these cases, it is not a priori clear what the energetic surface for proton transfer will look like and where the minimum pathway is located. Hence, one has to explore more than one variable to build an appropriate PES. We have chosen to use two variables for each proton transfer process: the two Cartesian coordinates in a plane that perpendicularly cuts the metal plane and comprises the assumed proton transfer pathway. The corresponding coordinates are shown in Figure 4, which also exhibits the respective starting and final geometries. As can be seen, a seven-atom platinum cluster is used for proton transfer “to top”, that is, to the top site of one of the surface neighbors of the platinum atom, on which the hydronium is adsorbed, whereas for proton transfer “to hollow”, that is, to one of the nearest hollow sites, we use a 12-atom cluster. In both cases, to build a two-dimensional (2D) PES, the calculation of many points is required. Consequently, the computational effort turned out to be too high for relaxed calculations, so starting from the relaxed geometry, just the transferring hydrogen atom was moved, and then the singlepoint energy calculated for each structure was used for building the PES. The corresponding energy plots are presented in Figure 5 both for the DFT calculations and the fitted EVB model. Some special points have been marked in each of the DFT plots. In the DFT plots the energetic behavior seems quite similar in both cases. There is a pronounced valley containing the minimum energy pathway for proton transfer. Following the latter, only a relatively small energy barrier for breaking the hydrogen-oxygen bond of at most a few tenths of an electronvolt has to be overcome. Interestingly, in both cases breaking this bond does not involve any stretching along the oxygen-hydrogen axis but bending the molecule toward the surface. Furthermore, as especially Figure 5a shows very clearly, there is a minimum situated near the platinum atom to which the remaining water molecule is adsorbed. Further diffusion to the final on-top and hollow position, respectively, occurs in a shallow energy valley. Consequently, the most important contribution does not come from the destination of the proton but from the platinum atom whose top site is occupied by the hydronium ion. Unlike for the one-dimensional (1D) process treated in Section 3.2, here we find a barrier in the situation where we start from the relaxed hydronium structures, although it is not high. This supports the picture that reorientation of the molecules in the first water layer facing the metal surface may be the rate-determining step.64,65 To obtain a complete picture, one has to consider the influence of the electrochemical environment on the dynamics of the process, which shall be done in the EVB model.

J. Phys. Chem. C, Vol. 112, No. 29, 2008 10823 3.4. Parameter Fit and Comparison of EVB and DFT Potential Energy Surfaces. As has already been noted, most potential functions used in the EVB model and the coupling function H12 are based on previous work. The fitting procedure for the platinum-hydrogen potential of eq 6 and the off-diagonal Hamiltonian elements of eq 13, which describe mainly hydrogen diffusion on the surface, has already been described in Section 2.1.2 above. So what still remains is the description of the coupling elements H2, j, which are mainly responsible for the crucial process, the proton transfer to the metal surface. Studying Figure 3, panels a-c, one can see some of the requirements for this coupling function: (1) it must reproduce the rise of the energetic barrier from zero to several electronvolts with growing distance from the surface; (2) it must become broader with increasing distance, and (3) it must account for the energetic lowering near the equilibrium distance caused by the agostic interactions. It is clear that a purely symmetrydependent function cannot fulfill these requirements, so there must be some additional distance-dependence. Therefore, those parameters in eq 18 that belong to factors depending on the distance to the surface have been adjusted to reproduce the results from Section 3.2. Figure 3 shows the comparison of the EVB to the DFT curves. As one can see, the overall agreement is quite good, although there are some deviations. Special attention has been paid to the region around the equilibrium distance, since this should be the most important for the MD simulations, whereas deviations from the exact shape occurring at larger distances, where the barrier exceeds 1 eV, should not play any important role. Figure 3d shows the shape and magnitude of H23 for the process described in Section 3.2, so one can understand the influence of symmetry and distance. Furthermore, eq 18 contains damping parts for deviations from the on-top position and for bending the internal hydronium angles. They contribute to the barriers that appear in Figure 5, panels b and e. Consequently, the pertinent parameters have been adjusted to reproduce the corresponding barrier seen in the DFT calculations. In addition, all parameters were adjusted to enable dynamical reassembly of the transfer cluster in situations where DFT calculations and chemical reason tell us that there should be no substantial contribution from the offdiagonal elements (pure VB states dominating), while simultaneously conserving total energy. Figure 5c illustrates the shape of the off-diagonal element H23 as a function of the coordinates used for the to top proton transfer. Again, also for the 2D surfaces shown in Figure 5, the general agreement and the reproduction of the key features is quite reasonable for the intended purpose of the EVB model, despite some apparent deviations. Additionally, one must bear in mind the qualitative issues arising from the very limited size of the metal clusters employed and other intrinsic restrictions of the DFT calculations. If higher quality DFT calculations were performed, the EVB model could be refined and perhaps even be simplified, as far as the functional form is concerned. 4. Single Trajectories for Proton Transfer As a first application, we present two typical trajectories for proton transfer events to platinum electrodes with different surface charge densities. The simulations were performed within the canonical (NVT) ensemble employing a tetragonal box of 22.2 × 19.23 × 80.00 Å3 containing 256 metal atoms in a slab of 4 layers, 256 water molecules, and the H5O2+ entity. The proton transfer complex (H5O2)+-Pt7 is initially assigned from these atoms and then dynamically evolves during propagation


Wilhelm et al.

Figure 6. Mechanism of proton transfer to the metal surface; the plots of the height of the transferring hydrogen above the surface illustrate the influence of surface charge density.

of the trajectory. Atoms and molecules that are not part of the complex interact with the same force field functions as those for the constituents of the complex (see Section 2.1.3). The choice of parameters for the coupling functions (see Section 3) was, among others, mandated by the requirement that the resulting force field can be used in a microcanonical ensemble. When used in an NVE simulation with this force field, our simulation code does indeed conserve energy. However, the trajectory simulations discussed below were performed at constant temperature in a tetragonal box, with periodic boundary conditions applied in all directions of space. Because of the large size of the box in the z-direction, the vacuum region between the water layer and the next image of the metal slab should be thick enough to minimize all artificial interactions between the images. The coulomb interactions were calculated using the Ewald summation method. A mean temperature of 298.15 K was maintained employing a Berendsen thermostat66 with a time constant of 2.5 ps for temperature control. The velocity Verlet algorithm was used for integration with a time step of 0.25 fs. Figure 6 gives an impression of the spacial course of proton transfer for low and high surface charge densities. The height of the candidate hydrogen atom H/ (which can change its identity several, possibly many, times until the final transfer) over the metal surface is plotted against time. A very typical feature of Figure 6a is that after some diffusion in the bulk the proton transfer cluster reaches the surface and remains there for a substantial fraction of the overall trajectory time at a distance between 2.5 and 3.0 Å; this corresponds to an adsorbed state. The inset shows the moment of proton transfer at approximately 214 ps; the height decreases to a value of approximately 1.5 Å, which is a characteristic value for a hydrogen atom bonded to the metal surface in the on-top position. The trajectory for a high surface charge, Figure 6b, shows very different features; apart from the much shorter overall trajectory time, the adsorbed state with its characteristic distance between 2.5 and 3.0 Å is completely missing. The vertical lines correspond to a proton transfer event toward the surface according to the Grotthuss shuttling mechanism, where the H/ atom changes its identity to a hydrogen atom bonded to a different water molecule that is much closer to the surface. As can be seen from the time period shortly after 11 ps, the last Grotthuss hopping event into the first water layer facing the metal surface is practically simultaneously accompanied by proton transfer of the H/ atom to the surface, which can be

seen from the typical distances around and below 1.5 Å. One more detail: distances notably below 1.5 Å give a hint to noton-top binding of the neutral hydrogen atom. The squares of the ground-state eigenvector coefficients, (cgi )2, can be interpreted as weights by which the diabatic state |i〉 contributes to the adiabatic state. This enables us to analyze which diabatic state dominates the adiabatic state at any given time. For the proton transfer to the surface, we can identify whether a concerted mechanism takes place or a stepwise one. In the concerted mechanism, the final proton transfer to the metal surface starts from a Zundel-like state where both states |1〉 and |2〉 contribute to the adiabatic state. Thus, the bridging hydrogen atom HB is more or less simultaneously transferred to the oxygen atom O/ that is left by H/ on its way to the surface. In the stepwise mechanism, a solvated hydronium-like ion (approximate Eigen ion) is formed as an intermediate before the actual proton transfer to the surface takes place. In this case, the intermediate structure may have time for reorientation. Two illustrations from single trajectories are shown in Figure 7. The first example, Figure 7a, shows the time around the final proton transfer event to the metal surface from a trajectory at a surface charge density of -7.5 µC/cm2. The first sector of the plot, before about 213.8 ps, features a dynamical equilibrium between the first two states, corresponding to a predominantly “Zundel-like” state, and an adiabatic state with one of the first two states dominating, so an “Eigen-like” state of the proton transfer cluster. Then, at the time of proton transfer around 214 ps, state |2〉 gains weight until we again have a more-or-less pure hydronium ion facing the metal surface, which is then very quickly superseded by the fast rising state |3〉, which corresponds to a neutral hydrogen atom bonded to the pivot platinum atom and two remaining water molecules. The second example, Figure 7b, is again taken from a trajectory at a surface charge density of -18.8 µC/cm2. Before 11.2 ps, state |2〉 dominates the adiabatic state, so we again find an Eigen-like state before proton transfer. In contrast to the case shown in Figure 7a, here the pure hydronium ion is replaced by a mixed state consisting of two or more of the “metal” states, that is, |3〉 to |9〉. At approximately 11.5 ps, state |3〉 becomes dominant, followed by a period where, once again, several metal states contribute to the adiabatic state. This corresponds to different bonding situations for the neutral hydrogen atom adsorbed at the surface. If just the metal state |3〉 contributes, it is most probably bonded in the on-top position to the pivot platinum atom, whereas the contribution of several states gives a hint to hollow or bridge adsorption sites.



Figure 7. Mechanism of proton transfer to the metal surface: (a, b) contributions of the diabatic states; and (c) hydrogen atoms changing parts, see text.

The last plot, Figure 7c, is taken from the same trajectory as 7a. It illustrates another indicator for proton transfer: the intramolecular distances within the proton transfer cluster. The blue curve shows the distance between the candidate hydrogen atom (H/) and the oxygen atom (O/) to which it is bonded in the diabatic states |1〉 and |2〉. The green curve shows the distance between the hydrogen HB, which forms the hydrogen bond within the formal Zundel ion, and O/. Before the proton transfer at about 214 ps, the green curve fluctuates more or less around 1.2 Å, a typical value for intra-Zundel hydrogen bonds, and the blue one around 1.0 Å, which is characteristic for a O-H bond. At the moment of transfer, HB takes over the part of H/ and becomes an ordinary hydrogen atom within a water molecule, whereas the H/-O/ distance strongly increases, indicating bond breaking. 5. Conclusions Empirical valence bond methods are attractive because they allow molecular dynamics simulations over a sufficiently long time to observe reactions. We have extended the WalbranKornyshev model for proton exchange in bulk water to proton transfer to an electrode surface. As a model electrode, we have chosen a Pt(111) surface, and determined by DFT calculations the functions that couple proton states in the solution to states on the surface. Naturally, this method can be applied to other metals as well, and our present model can be improved by more extensive calculations or the usage of more elaborate clusters. It should be stressed that the model presented here does not address quantum tunneling effects. These may be neglected due to the fact that the barrier is small and washed out by zeropoint motions (see, e.g. the results of ab initio path integral

simulations performed for an aqueous solution).11 Therefore, a proton resides most likely in a single-well potential near an oxygen atom. Due to fluctuations of the second coordination sheath of the proton, the bottom of this potential gradually moves from one O atom to another. The relevant activation barrier can be defined as well.67 As a first application, we have calculated trajectories for proton transfer at two different charge densities on the platinum electrode. These examples show that at the lower density of -7.6 µC/cm2 the proton transfer cluster first passes into a weakly adsorbed state about 2.5-3 Å from the surface, before the proton is transferred to the surface and discharged. At a higher charge density of -18.8 µC/cm2, the weakly adsorbed state is missing, and the transfer to the electrode is immediate. Of course, to arrive at conclusions, systematic studies employing a statistically significant number of MD trajectories are necessary. These will be reported in a subsequent communication. Acknowledgment. Financial support from the Deutsche Forschungsgemeinschaft (Schm 344/32) is gratefully acknowledged. RRN wishes to thank the Russian foundation for Basic Research (project No. 05-03-32381a). References and Notes (1) Eigen, M. Angew. Chem. 1963, 75, 489–508. (2) Zundel, G.; Metzger, H. Z. Phys. Chem. Neue Folge 1968, 58, 225– 245. (3) Brooker, M. H. In The Chemical Physics of SolVation. Part B; Dogonadze, R. R.; Ka´lmań, E.; Kornyshev, A. A.; Ulstrup, J., Eds.; Elsevier: Amsterdam, 1986; pp 149-185. (4) Kobayashi, C.; Iwahashi, K.; Saito, S.; Ohmine, I. J. Chem. Phys. 1996, 105, 6358–6366. (5) Heinzinger, K.; Weston, R. E. J. Phys. Chem. 1964, 68, 744–751.

10826 J. Phys. Chem. C, Vol. 112, No. 29, 2008 (6) Conway, B. E. Ionic Hydration in Chemistry and Biophysics; Elsevier: Amsterdam, 1981. (7) Zundel, G.; Fritsch, J. In The Chemical Physics of SolVation. Part B; Dogonadze, R. R.; Ka´lmań, E.; Kornyshev, A. A.; Ulstrup, J., Eds.; Elsevier: Amsterdam, 1986; pp 21-96. (8) von Grotthuss, T. C. J. D. Ann. Chim. 1806, LVIII, 54–74. (9) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1995, 103, 150–161. (10) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. Science 1997, 275, 817–820. (11) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601–604. (12) Izvekov, S.; Voth, G. A. J. Chem. Phys. 2005, 123, 044505. (13) Marx, D. ChemPhysChem 2006, 7, 1848–1870. (14) Vuilleumier, R.; Borgis, D. J. Chem. Phys. 1999, 111, 4251–4266. (15) Schmitt, U. W.; Voth, G. A. J. Phys. Chem. B 1998, 102, 5547– 5551. (16) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361– 9381. (17) Brancato, G.; Tuckerman, M. E. J. Chem. Phys. 2005, 122, 224507-1 to 224507-11. (18) Voth, G. A. Acc. Chem. Res. 2006, 39, 143–150. (19) Swanson, J.; Maupin, C.; Chen, H.; Petersen, M.; Xu, J.; Wu, Y.; Voth, G. J. Phys. Chem. B 2007, 111, 4300–4314. (20) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. J. Phys. Chem. B 2008, 112, 467–482. (21) Agmon, N. Chem. Phys. Lett. 1995, 244, 456–462. (22) Day, T. J.; Soudackov, A. V.; Cuma, M.; Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 2002, 117, 5839–5849. (23) Kornyshev, A. A.; Kuznetsov, A. M.; Spohr, E.; Ulstrup, J. J. Phys. Chem. B 2003, 107, 3351–3366. (24) Lapid, H.; Agmon, N.; Petersen, M. K.; Voth, G. A. J. Chem. Phys. 2005, 122, 014506 to 1014506-11. (25) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Chem. Phys. 2005, 122, 014515 to 1014515-6. (26) Meng, S. Surf. Sci. 2005, 575, 300–306. (27) Sugino, O.; Hamada, I.; Otani, M.; Morikawa, Y.; Ikeshoji, T.; Okamoto, Y. Surf. Sci. 2007, 601, 5237–5240. (28) Otani, M.; Hamada, I.; Sugino, O.; Morikawa, Y.; Okamoto, Y.; Ikeshoji, T. J. Phys. Soc. Jpn. 2008, 77, 024802. (29) Jacob, T.; Goddard, W. J. Am. Chem. Soc. 2004, 126, 9360–9368. (30) Taylor, C. D.; Wasileski, S. A.; Filhol, J.-S.; Neurock, M. Phys. ReV. B 2006, 73, 165402. (31) Filhol, J.-S.; Neurock, M. Angew. Chem. Int. Ed. 2006, 45, 402– 406. (32) Hinnemann, B.; Moses, P. G.; Bonde, J.; Jørgensen, K. P.; Nielsen, J. H.; Horch, S.; Chorkendorff, I.; Nørskov, J. K. J. Am. Chem. Soc. 2005, 127, 5308–5309. (33) Greeley, J.; Jaramillo, T. F.; Bonde, J.; Chorkendorff, I.; Nørskov, J. K. Nat. Mater. 2006, 5, 909–913. (34) Nørskov, J. K.; Christensen, C. H. Science 2006, 312, 1322–1323. (35) Rossmeisl, J.; Nørskov, J. K.; Taylor, C. D.; Janik, M. J.; Neurock, M. J. Phys. Chem. B 2006, 110, 21833–21839. (36) Sku´lason, E.; Karlberg, G. S.; Rossmeisl, J.; Bilgaard, T.; Greeley, J.; Jońsson, H.; Nørskov, J. K. Phys. Chem. Chem. Phys. 2007, 9, 3241– 3250. (37) Ogasawara, H.; Brena, B.; Nordlund, D.; Nyberg, M.; Pelmenschikov, A.; Petterson, L.; Nilsson, A. Phys. ReV. Lett. 2002, 89, 276102. (38) Coulson, C.; Danielsson, U. Ark. Fys. 1954, 8, 245–255. (39) Warshel, A.; Weiss, R. M. J. Am. Chem. Soc. 1980, 102, 2523– 2544.

Wilhelm et al. (40) Hwang, J.-K.; King, G.; Creighton, S.; Warshel, A. J. Am. Chem. Soc. 1988, 110, 5297–5311. (41) Åqvist, J.; Warshel, A. Chem. ReV. 1993, 93, 6218–6226. (42) Lobaugh, J.; Voth, G. A. J. Chem. Phys. 1996, 104, 2056–2069. (43) Sagnella, D. E.; Tuckerman, M. E. J. Chem. Phys. 1998, 108, 2073– 2083. (44) Walbran, S.; Kornyshev, A. A. J. Chem. Phys. 2001, 114, 10039– 10048. (45) Spohr, E.; Commer, P.; Kornyshev, A. A. J. Phys. Chem. B 2002, 106, 10560–10569. (46) Toukan, K.; Rahman, A. Phys. ReV. B 1985, 31, 2643–2648. (47) Spohr, E. J. Phys. Chem. 1989, 93, 6171–6180. (48) Olsen, R. A.; Kroes, G. J.; Baerends, E. J. Chem. Phys. 1999, 111, 11155–11163. (49) Papoian, G.; Nørskov, J. K.; Hoffmann, R. J. Am. Chem. Soc. 2000, 122, 4129–4144. (50) Ka¨lleń, G.; Wahnstro¨m, G. Phys. ReV. B 2001, 65, 033406. (51) Baˇdescu, S¸.; Salo, P.; Ala-Nissila, T.; Ying, S.; Jacobi, K.; Wang, Y.; Bedu¨rftig, K.; Ertl, G. Phys. ReV. Lett. 2002, 88, 136101. (52) Jacob, T.; Goddard, W. A., III J. Phys. Chem. B 2004, 108, 8311– 8323. (53) Ohwaki, T.; Yamashita, K. J. Electroanal. Chem. 2001, 504, 71– 77. (54) Cai, Y.; Anderson, A. B. J. Phys. Chem. B 2004, 108, 9829–9833. (55) CRC Handbook of Chemistry and Physics, 69th ed.; Weast, R. C., Astle, M. J., Beyer, W. H., Eds.; CRC Press: Boca Raton, Florida, 1988. (56) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (57) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, ReVision B.04, Gaussian, Inc., Pittsburgh, Pennsylvania, 2003. (58) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 270–283. (59) Wadt, W. R.; Hay, P. J. J. Chem. Phys. 1985, 82, 284–298. (60) Hay, P. J.; Wadt, W. R. J. Chem. Phys. 1985, 82, 299–310. (61) Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257–2261. (62) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213– 222. (63) Nazmutdinov, R. R.; Bronshtein, M. D.; Wilhelm, F.; Kuznetsov, A. M. J. Eletroanal. Chem. 2007, 607, 175–183. (64) Pecina, O.; Schmickler, W. J. Electroanal. Chem. 1997, 431, 4750(4) . (65) Pecina, O.; Schmickler, W. Chem. Phys. 1998, 228, 265277(13) . (66) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J.Chem. Phys. 1984, 81, 3684–3690. (67) Kuznetsov, A. M.; Ulstrup, J. Russ. J. Electrochem. 2004, 40, 1010–1018.

JP800414F

A Model for Proton Transfer to Metal Electrodes - American Chemical

Recommend Documents