MD simulations of heterogeneous reduction of the tert-butyl bromide molecule

Electrochimica Acta 55 (2010) 2442–2450

Contents lists available at ScienceDirect

Electrochimica Acta journal homepage: www.elsevier.com/locate/electacta

MD simulations of heterogeneous reduction of the tert-butyl bromide molecule Anna Ignaczak a,∗ , Wolfgang Schmickler b a b

Department of Theoretical Chemistry, University of Łód´z, ul. Tamka 12, 91-403 Łód´z, Poland Department of Theoretical Chemistry, University of Ulm, D-89069 Ulm, Germany

a r t i c l e

i n f o

Article history: Received 27 October 2009 Received in revised form 3 December 2009 Accepted 4 December 2009 Available online 13 January 2010 Keywords: Electron transfer Simulations Tert-butyl bromide reduction Rate constant Saddle point avoidance

a b s t r a c t This paper reports results of molecular dynamics simulations on the electrochemical reduction of the tertbutyl bromide molecule according to the scheme: (CH3 )3 CBr + e → (CH3 )3 C• + Br− . A model Hamiltonian is proposed, that is based on the Anderson–Newns concept of the electron transfer process and on quantum calculations performed for the (CH3 )3 CBr molecule and the (CH3 )3 CBr·− radical anion. The rate constants k and the transfer coefficients ˛ obtained at different solvent frictions  x , temperatures T and overpotentials , are presented. At  = 1.3 V, which corresponds to the experimental conditions, the ˛ estimates obtained from the simulations with low solvent frictions are very close to the measured value of 0.2; in more viscous solvents they are slightly higher. The ˛ coefficient is found to depend non-linearly on . It is also shown to depend on the temperature: at lower overpotentials the transfer coefficient decreases with T: ˛ (348 K) < ˛ (298 K) < ˛ (248 K), while for high  the reversed order appears. In all simulations strong saddle point avoidance was observed. The choice of the effective reaction paths is shown to be responsible for some trends found in the results. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction The homogeneous and heterogeneous reductions of alkyl halides have been investigated extensively in the past as a typical example of electron transfer reactions accompanied by significant modification of the reactant’s structure. The fundamental experimental data about kinetics of such reaction were provided by the group of Savéant in a series of papers [1–6], who also proposed the simple theoretical model for description of such processes. In this model [3] the Morse potential was used to describe the C–X bond rupture and combined with the Marcus–Hush theory [7–10]. Although the general trends observed in the experiment were reproduced in the theoretical results, several computed values appeared to differ from the measured ones. For example, the transfer coefficients calculated with this model for the butyl bromide molecules are always higher than the experimental value of 0.2. It was suggested [1] that the poor agreement might be due to the peak widening resulting from the overlap of the RX + e → R· + X− and R· + e → R− reduction waves, but it is difficult to estimate the extent of this perturbation. Other theoretical works [11,12] conducted on the electrochemical reduction of the tert-butyl bromide molecule reported transfer coefficient values above 0.3, thus close to estimates obtained with

∗ Corresponding author. E-mail addresses: [email protected] (A. Ignaczak), [email protected] (W. Schmickler). 0013-4686/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.electacta.2009.12.003

the Savéant’s model. The common element of these theoretical approaches is the use of the Morse potential to describe the C–Br bond cleavage, which may be a cause of similarity in the results. It should be mentioned that in [11] another value, 0.111, was also obtained under different assumptions, showing sensitivity of this parameter to the model used. In the present work we propose the description of the C–Br bond dissociation based on the results of quantum calculations. A similar strategy was used in our earlier studies on the electrochemical reduction of the tert-butyl chloride molecule [13,14]. Unfortunately, the lack of experimental data due to the extremely high overpotentials necessary for the latter compound did not allow for a satisfactory comparison and verification of the model. This is not the case for the electrochemical reduction of the tert-butyl bromide molecule, which was analyzed in detail by Savéant’s group in the past [1–3,5,6]. The potential energy surfaces (PES) derived from the quantum calculations performed for the neutral molecule (CH3 )3 CBr and the anion radical (CH3 )3 CBr·− are combined with the model Hamiltonian formulated by one of us [15] to describe the solvent and electrode effects on the kinetics of the adiabatic electron transfer reactions. This yields a two-dimensional potential energy surface that is introduced into the simulation scheme, allowing to study kinetics and mechanism of this reaction. The details of the model are presented in the next section. In the following part the reaction rates and the transfer coefficients are analyzed with respect to various parameters characterizing the liquid side as well as to the molecule–metal interaction strength. The main goal of our

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

2443

work is to test whether the theoretical model based on quantum calculations yields also the transfer coefficient higher than the experimental estimates, as well as to provide a more detailed insight into the reaction kinetics and mechanism. In the last section the most important achievements of our work are summarized. 2. Methods 2.1. Potential The electron transfer considered in the present work is assumed to occur from an electrode to the tert-butyl bromide molecule that is located at some distance from the metal surface. An important step of this reaction is the C–Br bond dissociation, which occurs simultaneously with the electron transfer process, so the final products are the tert-butyl radical and the bromide anion: (CH3 )3 CBr + e → (CH3 )3 C· + Br−

(1)

The quantum method used to obtain the potential energy surfaces for the molecule and for the anion radical as a function of the C–Br bond length was the B1LYP density functional [16–19] combined with the 6-311G* basis set. It was selected after extensive tests performed using different quantum methods, including standard Hartree–Fock with and without Møller–Plesset electron correlation correction, and several basis sets available in the Gaussian 98 package [20]. In the calculations the three C–C bonds and three C–C–C angles were assumed to be equal. The methyl groups were allowed to rotate independently, but the corresponding C–H bond lengths and the H–C–H and H–C–C angles in the three groups were again assumed to be the same. The final structure obtained from the optimization is shown in Fig. 1, and the detailed parameters are listed in Table 1 together with their experimental estimates found in the literature. As can be seen, and in accordance with our earlier works on the tert-butyl halides [13,32], the structure is highly symmetric, with the three methyl groups oriented such, that in each of them two hydrogens point towards the bromine atom and one away from it. The latter bond is somewhat longer than the remaining two, although the difference is not very large, as can be seen in Table 1. The optimized C–C bond length and C–C–C angles agree very well with the experimental data. The C–Br bond is in the range of measured values, although the more recent spectra suggest rather smaller values. However, this parameter is not crucial in our model, in which the relative changes in the carbon–bromide distance are considered rather than the total bond length. Considering the final purpose of these calculations, which is to construct the potential energy sur-

Fig. 1. Structure of the tert-butyl molecule obtained from the quantum calculations.

face to be used in the simulations, the most important parameter for a proper description of the intramolecular changes occurring in the reaction course is the carbon–bromide vibrational frequency, therefore the best reproduction of this property was our main criterion for the selection of the quantum method. For the tert-butyl bromide molecule the B1LYP/6-311G* method gives an excellent agreement with the experimental estimate of 514 cm−1 . Considering relatively wide range of experimental estimates for the C–Br dissociation energy, our computed value appears to be very close to the value 2.82 eV of Savéant [3]. The potential energy surfaces obtained for the neutral molecule and the anion radical are plotted in Fig. 2. The molecule potential has been shifted such that it has the minimum at E = 0 eV, which satisfies the assumption that the energies of the reactant and the products (defined in Eq. (1)) must be equal at equilibrium. Similarly to the tert-butyl chloride case, also here the PES for the reduced form exhibits some long-distance interaction between the departing bromide anion and tert-butyl radical. On the curve this is marked first by some small wave observed at dC–Br = 3.3 Å when the inversion in the tert-butyl radical occurs: from this moment the methyl groups are pointed towards the bromide. The shallow minimum seen around 4 Å results from the weak electrostatic interaction between the anion and the radical, which vanishes as they separate further. The points obtained from quantum calculations we have fitted to analytical potentials of the form similar to that used previously

Table 1 Experimental data available for the (CH3 )3 CBr system and our parameters obtained from the B1LYP/6-311G* quantum calculations. Experimental results

dC–Br /Å

dC–C /Å

[21] [22] [23] [23] [24] [25] [26] [27] [28] [29] [3] [30] [31] This work

2.06 1.92 1.98 1.94 1.92 1.94 1.975 1.974

1.55 1.54a 1.54a 1.54a 1.54 1.53

a b c

dC–H /Å 1.09a 1.093a 1.093a

∠CCC/◦

Edys (C–Br)/eV

ω/cm−1

110 111.3 111.3a 109.28a 111.5 108.9

2.767 2.732 2.82 2.936 2.0358

1.5249

Values assumed in the model. Value for the two C–H bonds in each methyl group directed towards Br. Value for the C–H bond pointing away the bromine atom.

1.0901b 1.0954c

111.98

2.785

514 515.02

2444

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

tronic states on the reactant and on the metal via their energies and occupation numbers, as well as their interaction. The latter can be generally described by the energy broadening parameter

 2 Vk  ı(ε − εk ), where ε is the electronic energy. As in all

=

k

Fig. 2. Potential energy profiles for the C–Br bond elongation in the tert-butyl bromide molecule and the corresponding radical anion obtained from the B1LYP/6311G* calculations (symbols) and with the fitted potentials (lines) given by the Eqs. (2) and (3).

in our studies on the tert-butyl chloride reduction [13]: Vi (y) = A1 exp(−A2 y) + A3 exp(−A4 y) + A5 exp(−A6 y) + A7

(2)

Vf (y) = B1 exp(−B2 y) + B3 exp(−B4 y) + B5 exp(−B6 y) + B7 exp(−B8 (y − yinv )2 )

(3)

In the anion potential the last term has been additionally introduced to the previous form used in [13] to reproduce the irregularity on the potential curve due to the sudden inversion in the tert-butyl part. In the final form used in the simulations the variable y in Eq. (2) is defined as y = dC–Br − d0 , where d0 is the C–Br equilibrium distance in the molecule, while in Eq. (3) it is defined as y = dC–Br − (d1 − d0 ), where d1 = 4 Å – the minimum on the anion PES. The parameters A and B are presented in Table 1, and the potential curves obtained with them are shown in Fig. 2 along with the quantum points. The fitted potentials of the molecule and the radical anion cross at dC–Br = 2.5 Å and predict for the reductive bond dissociation the activation energy of 0.525 eV, which is much smaller than predicted by the Morse potential [3]. Following the strategy introduced by Koper and Voth [33], the intramolecular effects described with these potentials were combined with the model Hamiltonian proposed by one of us for the outer sphere electron transfer reactions at electrochemical interfaces [15]. Since this part of the model was already described in detail in our several earlier works [13,14,34–36] as well as in the original paper of Schmickler [15], we will only summarize here briefly the main concept of it. The Hamiltonian consists of two most important contributions to the electron transfer: He – the electronic part and Hph – the solvent part. The first term includes descriptions of elec-

our earlier works, also here we apply “the wide-band approximation” which assumes that during the reaction the interaction between the reduced particle and the electrode does not change, thus the parameter  is set to be constant. The second term describes the effect of the solvent reorganization. The solvent is modelled as a phonon bath, to which the harmonic approximation is applied. Thus, it is considered as a set of harmonic oscillators, characterized by the dimensionless coordinates q and momenta p and the corresponding frequencies ω. Obviously, the changes in the solvent configuration influence the electronic states of the redox system, and this effect is also included into this part of the Hamiltonian. The interaction between the solvent and the reacting particle is assumed to be linear and is multiplied by the electron occupation number, na , on the antibonding orbital of the reactant. The description of the solvent reorganization was then further simplified by adopting the well-known Marcus concept, in which all changes of the solvent configuration are represented by one generalized reaction coordinate x, and the solvent response – by the solvent reorganization energy . The coordinate x was normalized such, that at the initial state x = 0, while at the final state x = 1. The final form of the Hamiltonian for the bond-breaking electron is the sum of the two contributions described just above and additionally of the third term representing the bond-breaking effect: H = He + Hph + Hbb

(4)

This last term is defined as: Hbb = (1 − na )Hi (y) + na Hf (y), where Hi and Hf are the contributions coming from the reactant (initial state) and the product (final state), respectively. Each of these contributions contains the kinetic and potential energy operators, the latter represented by potentials Vi and Vf , defined in the present work by Eqs. (2)–(3). From the Hamiltonian (4) and with the assumptions described above has been derived the potential energy surface of the following form: E(x, y) = ε˜ a < na > +x2 +

 ln 2



ε˜ 2a + 2



ε2a + 2

+ Vi (y)

where the energy of the reactants antibonding orbital is defined with respect to the Fermi level of the metal as: εa =  − e0 , where  is the overpotential. The first term in Eq. (5) is the electronic energy: ε˜ a (x, y) = εa − 2x + Vf (y) − Vi (y)

(6)

which was obtained by collecting all terms containing the occupation number na , which is defined as: na  = (1/)arc cot(˜εa /).

Table 2 Parameters obtained from fitting the Vi and Vf potentials to the quantum points. Vi parameters

Values [a.u.]

Vf parameters

Values [a.u.]

A1 A2 A3 A4 A5 A6 A7

0.0227462583 2.3445509650 −4.4160009480 × 10−04 3.2880105920 −0.1500799986 0.3424763972 0.1277753404

B1 B2 B3 B4 B5 B6 B7 B8 yinv

2.322762005 × 10−07 2.7759080790 1.1508230220 0.5474199149 −1.1590694950 0.5437036782 3.9194503290 × 10−03 2.9599226680 1.4361922643

Vi – potential for the (CH3 )3 CBr molecule. Vf – potential for the (CH3 )3 CBr·− radical anion.

(5)

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

2445

trochemical reduction of aquocomplexes. The stochastic collision model of Kast et al. describes the motion of the particle 1 of mass m1 , which collides with the bath consisting particles 2 of mass m2 . It is assumed that there is no interaction between the bath particles and that their random motion at the temperature T is governed by the Maxwell–Boltzmann distribution equation f (u) = 1/2 (m2 /2kB T ) exp(−m2 u2 /2kB T ), where kB is the Boltzmann constant and u is the velocity. The position r of the colliding particle 1 in the simulation step n + 1 is then determined from its three previous positions according to: rn+1 = (2 − ˛)rn − rn−1 + ˛rn−2 + 2˛un dt + as (dt)

Fig. 3. Potential energy surface for the electrochemical reduction of (CH3 )3 CBr at equilibrium obtained with Eqs. (5)–(6) and potentials (2)–(3): (a) contour plot; the red line marks the points corresponding to na = 0.99, (b) isosurface. Parameters used:  = 0.01 eV,  = 0.624 eV,  = 0 V. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Fig. 3 presents the contour (Fig. 3a) and isosurface (Fig. 3b) plots of the potential energy surface for the electrochemical reduction of the tert-butyl bromide molecule at equilibrium. Apart from the coefficients listed in Table 2, the following values of parameters characterizing the system were used: the solvent reorganization energy  = 0.624 eV (this value is equal to 25kB T at the room temperature), the reactant–electrode interaction energy  = 0.01 eV, the overpotential  = 0 V. The two regions corresponding to the reactant (R) and the product (P) can be easily distinguished. It should be stressed here, that the product area should be in reality extended to much larger y values, but to provide a more detailed view of the saddle point (SP) vicinity in this figure y variable was limited to a few angstroms. The saddle point is the lowest energy point on the wide energy barrier separating the reactant and product regions and, if the lowest energy passage is assumed, the reaction is expected to proceed via this point. At the equilibrium the saddle point is located at (0.49512,0.46180) and the difference between its energy and the reactant minimum is equal to 0.66 eV. This is significantly lower than the estimates obtained by other authors [3,12]. However, as in our earlier works [14,34–36], also for this system the so-called saddle point avoidance is expected to occur, which means that the reaction not always proceeds via the lowest energy point on the barrier. 2.2. Simulations method The above potential energy surface has been used in the simulations which were performed according to the algorithm of Kast et al. [37,38]. This method was used in our earlier works to study bond-breaking electron transfer reactions as well as the elec-

2

(7)

where as is the particle acceleration, dt is the simulation time-step, and ˛ defines the mass ratio as ˛ = m2 /(m1 + m2 ), which can be easily transformed into the friction parameter  via the relation:  = 2˛/dt. Since our potential energy surface depends on two independent coordinates x and y, the scheme described above was applied in each direction. The mass mx was derived from the solvent reorganization energy under assumption that the frequency ωx = 1 ps−1 (the longitudinal relaxation time of common solvents is usually of the order 1 ps), which, taking into account that the coordinate x is dimensionless, yields mx = 2. The mass my was set to be the reduced mass of the C and Br atoms, which is equal to 10.44 atomic mass units. Due to the fact, that with the assumed ωx the C–Br the frequency ωy is 96 times higher, two different time scales were used in the simulations, with the two time-steps dtx and dty obeying the relation: dtx /dty = ωy /ωx = 96. Thus, one positional change of the solvent corresponds to 96 changes of the C–Br bond. In all simulations the time-step dtx was set to be equal 0.01 ps. As discussed in the literature [39–41] the much faster bond vibrations and the delayed response of the solvent causes that the friction exerted on the bond is much smaller than in the solvent. As in our previous works [14,35,36] the difference in the two frictions was approximated by using the simple relation: ωy /ωx = x /y . A more detailed discussion on the quality of this simple model and justification of its use can be found in [35]. The reaction rate k was obtained as the inverse of the average time required by the system to escape from the reactant well to the product area. The average was taken over all trajectories built for the given reaction conditions; depending on the parameters used, from 1000 up to 4000 trajectories were pursued. Each trajectory started at the bottom of the reactant well (point (0,0) on the potential energy surface shown in Fig. 3) and ended when the occupation probability na was higher than 0.99, and at this point the single simulation time was collected. The relatively steep shape of the barrier in the y direction assures that no re-crossing of the barrier occurs, but it was tested additionally in our simulations and proven to be true. 3. Results 3.1. Reaction rate The rate of the electrochemical reduction of the tert-butyl molecule has been examined with respect to several parameters which can influence the reaction, namely the solvent friction  x , the overpotential , the particle–metal interaction , the temperature T, and the solvent reorganization energy . The results are displayed in Figs. 4–8. The effect of the friction parameter  x was tested for two different overpotentials:  = 1.16 V, which correspond to the activation energy at the saddle point of 0.189 eV, and  = 1.4 V, where Eact (SP) = 0.121 V. As expected, the rates in the second case are much

2446

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

Fig. 4. The rate constant dependence on the solvent friction parameter  x obtained for two different overpotentials.

Fig. 5. The reaction rate as a function of the overpotential  obtained with four different solvent frictions.

Fig. 6. An effect of the metal–reactant interaction  on the rate constant observed at different solvent frictions.

Fig. 7. The rate constant as a function of the solvent reorganization energy.

Fig. 8. Temperature dependence of the reaction rate obtained from the simulations.

higher, but surprisingly the two curves differ significantly in their shape. In the case of higher energy barrier ( = 1.16 V) the results are in a qualitative agreement with Kramers’ turnover theory, namely two regions are distinguished: first, when the process is underdamped (small  x values) the rate increases along with the friction, and the second, when the process is overdamped (large  x ) and the rate decreases when the friction becomes higher. However, the underdamped region covers a relatively wide range of frictions, so the maximum is reached around  x = 70 ps−1 . The turnover region is not sharply peaked, but rather broad and smooth, and when  x is further increased, the rate diminishes slowly. Thus, in a relatively wide range of frictions, between 40 and 110 ps−1 the rate constant does not change very much. The second curve (1.4 V) does not show the maximum at all – the rate increases fast at the region where the friction is extremely low, and the growing tendency is preserved also at higher frictions, although in this region the changes are much smoother. The results presented in Fig. 4 may seem to contradict the experimental findings (for review see [42–46]) reporting that the electron transfer reactions become usually slower in more viscous solvent. However, as pointed out by Weaver [45,46], these conclusions are drawn from measurements performed in different solvents, which means that not only the viscosity factor was changed, but also the activation energetics, while Fig. 4 presents solely the effect of the friction parameter. In fact, the differences seen in the two curves suggest that an interplay of various parameters may be of a great importance for the results obtained. Therefore in our further studies four different frictions were applied to explore the rate dependency on other parameters. An overpotential effect on the rate constant can be seen in Fig. 5. As expected, the rate increases when larger  is applied due to lowering the energy barrier. However, the slopes of the four curves vary for different values of the friction parameter. It is interesting to analyze the steepest curve  x = 120 ps−1 with respect to the other curves, as it crosses them at certain potentials. In agreement with the trends observed in Fig. 4, at  = 1.9 V the reaction is faster when the friction parameter is higher, while at  = 0.9 V the rates obtained with  x = 2 ps−1 and  x = 120 ps−1 are almost the same, and much smaller than in the other two cases. This suggests that when the energy barrier becomes higher, the turnover region is shifted towards lower friction values. This effect is related to a varying reaction mechanism, which will be discussed in Section 3.3. Although the interaction between the reacting molecule and the electrode described with the  parameter also changes the height of the barrier, no curves crossing appears in Fig. 6. The rates order, dictated by  = 1.16 V which was used in these calculations, is preserved for almost all  values. The barrier modifications due to a stronger interaction of the molecule with the metal are much smaller than those caused by overpotentials, as can be seen in

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

2447

Table 3 Activation energies (defined by the position of the saddle point) obtained for different parameters. (a) /V

Eact (SP)/eV

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

0.660 0.611 0.564 0.518 0.474 0.431 0.389 0.350 0.311 0.275 0.240 0.207 0.177 0.148 0.121 0.097 0.074 0.055 0.037 0.023

(b) /eV

Eact (SP)/eV

0.001 0.003 0.005 0.007 0.01 0.03 0.05 0.07 0.1

0.204 0.200 0.197 0.193 0.189 0.164 0.143 0.126 0.103

(c) /eV

Eact (SP)/eV

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.103 0.122 0.142 0.163 0.184 0.205 0.226 0.248

Other parameters used in the calculations: (a)  = 0.01 eV,  = 0.624 eV, (b)  = 1.16 V,  = 0.624 eV, (c)  = 1.16 V,  = 0.01 eV.

Table 3, where activation energies obtained from the saddle point positions are listed for various parameters. Somewhat stronger effect on the barrier height, but also generally on the shape of the potential energy surface, has the solvent reorganization energy (Table 3) and this is again marked in the reaction rates in Fig. 7. The curve  x = 120 ps−1 crosses again all other curves, as happened in Fig. 5, and for the highest  values, above 0.75 eV, this friction value slows down the reaction even more than  x = 2 ps−1 . Obviously, in these conditions the overdamping effect must be much stronger than it was predicted in Fig. 4 for the same overpotential, but smaller  values. The change in the rates order with respect to different friction values is probably an effect of the frictions anisotropy combined with the specific shape of the potential energy surface, but also temperature has some influence on it, as shown in Fig. 8. A closer inspection of the curves reveals that for higher T values the rates k( x = 120 ps−1 ) are above k( x = 40 ps−1 ), while at lower temperatures a reversed order is observed. To explain all these tendencies a detailed analysis of the reaction paths in different conditions is necessary, and this will be performed in Section 3.3. 3.2. Transfer coefficient The transfer coefficient for various conditions of the reaction was obtained from the relation of the rate constants to the overpotential according to: ˛ = kB Td ln(k)/d. We have used several different approaches to calculate this important parameter. In the first method (labelled I) the derivative was calculated by fitting a second-order polynomial to three successive points, and obtaining from it the ˛ value for the middle point. The results obtained with this procedure are shown in Fig. 9a. Notice, that since this method poorly evaluates the boundary points due to the fact that the next neighbouring point is missing, they have been removed from Fig. 9a. In the next approach an attempt was made to eliminate the effect of rate fluctuations due to statistical errors, which might manifest itself in strong perturbations of the results. This was done by using the whole set of points (, ln(k)) belonging to a given series and fitting to them carefully a suitable polynomial, from which the derivatives were next calculated. Two types of functions were tested: third and second-order polynomials. The former function (labelled as method II) reproduces very well the behaviour of orig-

Fig. 9. Transfer coefficient ˛ calculated for different overpotentials  and friction parameters  x : (a) model I and (b) model II (see text and description of Table 4 for details).

inal points, the standard deviation being below 0.05. The transfer coefficient obtained with this method is presented in Fig. 9b. The second-order polynomial, which yields the typical linear relation of ˛ vs. , as predicted in the Marcus theory, was found to give very poor fittings. The standard deviation in most cases was above 0.1 and only the very general growing trend observed in the simulations was mimicked by these curves, therefore these results are not included into our analysis. As predicted, indeed some fluctuations are present in Fig. 9a, but approximately the same tendencies can be found in the two subplots. All curves obtained with different frictions exhibit similar pattern – first ˛ decreases in agreement with a common wisdom, but for extremely high overpotentials it starts to increase again. This naturally corresponds to the shapes of the k vs.  curves. A detailed analysis of Fig. 5 provides an explanation: all curves exhibit a slight change in their slope at very high overpotentials. Our observation of the very poor quality of the parabolic fitting as well as the fact that the curves obtained with different frictions exhibit similar shapes, just shifted along the overpotential axis, seem to suggest that the transfer coefficient for this reaction is not a linear function of the potential. On the other hand, closer inspection of the experimental data provided in the literature (see for example Fig. 12 in [47] or Fig. 2 in [48]) brings some uncertainty about the strict linearity of ˛() as concluded from rather widely scattered points. Independently of the model used in our calculations, several common features can be still found in the subplots of Fig. 9, like the systematic increase of the transfer coefficient with friction. The difference calculated between the highest and lowest curves obtained with the third-order polynomial depends on , but in average it is around 0.05. Notice, that for the conditions considered in this series the highest curve corresponds to the overdamped reaction. As mentioned in the introduction, the transfer coefficient for electrochemical reduction of tert-butyl bromide was found experimentally to be equal 0.2 at the midpoint between the peak and

2448

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

Table 4 Transfer coefficients obtained for  = 1.306 V and four friction parameters at T = 298 K, compared to experimental and other theoretical data. ˛ x

2 ps−1

40 ps−1

70 ps−1

120 ps−1

Theoretical method This work – model I This work – model II

0.2032 0.2125

0.2100 0.2165

0.2180 0.2201

0.2411 0.2534

Other theoretical worksa [1] [3] [5] [11] [12] Experimental data [1] [6]

˛theor 0.31 0.30 0.32 0.33 0.35 0.111 0.334 0.303 0.325 ˛exp 0.20 0.27b

Procedures: I – from a second-order polynomial fitted locally to each three sequential points of ln k vs. , II – from a third-order polynomial fitted to all points. a Two ˛ values in some works come from different parameters assumed in the model. b Low temperature measurement (T < 270 K).

half-peak potential Em = −2.397 [3]. The overpotential to which it ◦ corresponds can be obtained as  = Em − E ◦ − r where E is the standard potential of the reaction (1) and r is the potential difference between the reaction site and the solution. In the work of German et al. [11] it was estimated from the quantities reported in [1] as 1.306 V (since in the present work for simplicity the overpotential sign was omitted, we refer here to the absolute values of ) and we will use the same value to compare our results with experimental and other theoretical estimates. In Table 4 we present values of ˛ calculated at  = 1.306 V using the methods I and II described above and the data found in the literature. As can be seen, our results for all frictions tested are significantly lower than those obtained with models based on the Morse potential and much closer to the experimental values. The higher values of ˛ obtained in those earlier theoretical works may be an indication that the energy barrier was overestimated due to the use of the Morse model for the bond-breaking process. The last experimental ˛ value of 0.27 listed in Table 4 was obtained at low temperatures T < 270 K. We performed two additional series of simulations scanning the dependence of k on  at T = 248 K and T = 348 K, using in all calculations the same friction coefficient  x = 70 ps−1 and from the results the transfer coefficient was again calculated using the method I. In Fig. 10 the results are

Fig. 10. Transfer coefficient as a function of the overpotential obtained with model I at three different temperatures.

presented together with those obtained at T = 298 K. The same nonlinear pattern is observed in all three curves, with the minimum shifted along the overpotential axis towards higher  as the temperature decreases. This results in crossing the curves at a certain point and at the experimentally tested potential  = 1.3 V the ˛(248 K) value appears to be somewhat higher than ˛(298 K). This qualitatively agrees with the experimental results listed in Table 4 for different temperatures, although we must admit that our calculated values are still lower than 0.27. One reason of such a discrepancy can be the relatively low value of 0.624 eV used in the present work for the solvent reorganization energy , while in [5] the value of 1.113 eV was suggested. As shown in Table 3 and Fig. 7 larger  values increase the barrier height and decrease the reaction rate, thus might affect also the slope of the k() curves. 3.3. Reaction mechanism Explanations for the unusual trends observed in the rate constants as well as in the transfer coefficient can be sought in the reaction mechanism. For the systems featuring multidimensional surfaces and frictions anisotropy the so-called saddle point avoidance is likely to occur. In the transition state theory it is assumed that when travelling from the initial to the final state, the system crosses the barrier at the lowest energy point. However, the two-dimensional potential surface provides also other possible reaction paths, corresponding to higher energy points at the barrier. This phenomenon was extensively studied in the past for various systems [49–53], and was also observed in our earlier works [14,34–36]. The saddle point avoidance, if strong, is expected to influence the transfer coefficient, which is supposed to reflect the symmetry of the potential surface at the transition state. This, in turn, should depend on the direction from which the barrier is approached. From our simulations it can be concluded that the preferred pathway is parallel to the y axis, although angular approach to the surface cannot be entirely ignored. As shown in Fig. 11, in the y direction the general barrier shape is the most asymmetric for larger x values, and the symmetry increases as x decreases. Thus, if most trajectories pass the barrier with the solvent configuration adjusted already to the product state (higher x) the ˛ values should be expected to be smaller than in the case when the solvent preserves the reactant state configuration (lower x). However, a more detailed analysis of the barrier profile complicates the matter. It appears, that while at  = 0.9 V the edge of the barrier has almost the same shape for all x values, this changes when higher overpotentials are applied. In the latter case the top of the

Fig. 11. Profiles of the barrier along the C–Br direction y for several values of the solvent coordinate x. To compare the relative changes in slopes all curves were shifted such that the maxima coincide at some (virtual) point y = 0 and E = 0 eV.

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

2449

Fig. 13. Distribution Nf /Ntot (in percentage) of the final coordinates x obtained from Ntot = 1000 trajectories for three different overpotentials.

Fig. 12. (a) Saddle point avoidance observed in the simulations performed to study the friction effect on the reaction rate (shown in Fig. 4). (b) Distribution of the final coordinates x (corresponding to the points when the occupation number na = 0.99) obtained for selected values of the friction parameters. Each distribution was obtained as Nf /Ntot (in percentage) from Ntot = 4000 trajectories, where Nf is defined as a number of runs finished inside each (xi ,xi+1 ) interval of the width of 0.05. The perpendicular line marks the position of the saddle point at the given overpotential.

barrier is slightly smoother and broader for large x than for x = 0.4; going towards smaller x values it becomes less and less broad and its degree of asymmetry does not always follows the pattern found in the general shape of PES. These subtle effects in the barrier shape are observed at a very limited range of y values (±0.001 Å from the barrier top), thus they may play a role only if the system passes the barrier relatively slowly. To examine whether there is any relation between ˛ and the saddle point avoidance, in our simulations the saddle point area was defined to be a rectangle centred at SP and of the size 0.1 x × 0.1 y, where x and y are the differences between the minima on the potential surface in x and y directions, respectively. The number of trajectories not passing this region was collected in each simulation and the results evaluated in percentage with respect to the total number of run will be presented as a qualitative description of this phenomenon. Additionally the final coordinates xf and yf (when the reaction was assumed to occur, i.e. na = 0.99) of all trajectories in each simulation were analyzed with respect to the SP position. For this the whole range of the reaction coordinates in each direction was divided into 60 intervals and the number Nf of trajectories ending in each interval was counted and related to the total number Ntot of trajectories in each simulation. In all runs yf were found to be beyond y(SP) while xf distribution is relatively wide and varies with the reaction conditions. Some additional tests confirmed, that the final coordinate xf differs from the point at which the barrier is crossed by less than 0.01, thus is a good indication how the reaction proceeds. In Fig. 12 the magnitude of the saddle point avoidance as a function of friction is presented. As can be seen, for extremely low

frictions the saddle point avoidance first drops sharply, and after reaching the minimum, systematically increases until for large frictions all trajectories avoid the SP. Fig. 12b provides an information which side is chosen along the x coordinate. For very low  x values most trajectories pass the barrier beyond the saddle point, so the reaction occurs with the solvent configuration corresponding to the reduced state. When the friction increases, the maxima on the distribution curves are shifted closer to the saddle point, so there is larger probability that the reaction will assume the lowest energy path – the distribution of trajectories built with  x = 10 ps−1 is almost symmetrical with respect to the line indicating the SP position. At even higher frictions the saddle point avoidance increases again, but now the system chooses more often to cross the barrier without significant reorganization of the solvent. The shift of the maximum towards lower x values when  x becomes larger explains the trends observed in Fig. 9 and in Table 4, namely larger values of the transfer coefficient for  x = 120 ps−1 . According to our calculations in these conditions the reaction occurs in the region where the barrier is less steep (upper part of the potential surface in Fig. 3a). A similar analysis can be performed for different overpotentials and the results are presented in Fig. 13 for  x = 70 ps−1 . It should be stressed here, that the relation between the transfer coefficient and the saddle point avoidance cannot be as straightforward as for different frictions. Obviously, when the overpotential is increased, the barrier is lowered and it becomes more asymmetric. Thus, the ˛ value will reflect simultaneously the latter effect as well as the different choice of the reaction path upon changed conditions. Indeed, for larger  the saddle point avoidance increases, and this effect is the strongest for the turnover friction of 70 ps−1 , so the system crosses the barrier more often farther away from SP. At the same time, the region in which the reaction is likely to occur becomes narrower – the distribution curves become higher and not as broad as for lower overpotentials. This is important as the transfer coefficient is an average over all possible reaction paths in given conditions and it may be responsible for a slight increase in ˛ at the very high overpotentials. Another important aspect that may play a role here is an increasing (with ) probability that at the activated state not only the bond length is changed but also the solvent reorganization when the overpotentials are higher: for larger  values more trajectories passes the barrier with simultaneous change of the x and y coordinates. In that case the potential energy shape ‘felt’ by the system will be different than that shown in Fig. 11. 4. Summary and conclusions In the present work a stochastic simulation method was applied to study the kinetics of the reduction of the tert-butyl

2450

A. Ignaczak, W. Schmickler / Electrochimica Acta 55 (2010) 2442–2450

bromide molecule at the metal–liquid interface. The main difference between our approach and the earlier theoretical works is the potential describing the C–Br dissociation, which is based on DFT calculations performed for the (CH3 )3 CBr molecule and the (CH3 )3 CBr·− radical anion. The most important achievements and conclusions from our calculations are following: (a) The transfer coefficient ˛ obtained with our model for  = 1.306 V (which corresponds to the conditions of experimental measurements) and T = 298 K was found to be close to the experimental value of 0.2, and much smaller than earlier theoretical predictions suggesting values above 0.3. (b) All series of simulations predict the transfer coefficient to decrease generally with overpotential, but in a non-linear manner. At low overpotential regions it can be assumed to follow locally a linear relation, but in the whole range of  tested it rather exhibits a parabolic pattern. At extremely high overpotentials a slight increase of ˛ is sometimes noted. (c) ˛ appears to increase with the friction parameter. For the conditions given in (a) the increase by about 0.04 is noted when going from  x = 2 ps−1 to  x = 120 ps−1 . (d) The temperature effect on the transfer coefficient is found to depend on the overpotential applied. At lower potentials, at which experiments are usually performed, the following order is found: ˛ (348 K) < ˛ (298 K) < ˛ (248 K). However, at very high overpotentials an opposite order is observed. (e) The activation energies read from the saddle point on the potential energy surfaces are significantly lower than predicted in the earlier works on the subject. Nevertheless, since in all simulations the very strong saddle point avoidance is observed, the apparent activation energies of the reaction are expected to be higher, as the reaction does not follow the lowest energy path. Also, the solvent reorganization energy used in the present work is somewhat lower than in those earlier models. (f) The choice of the effective reaction path is also responsible for some trends observed in the transfer coefficient when the reaction parameters are changed. Certainly, the model used in the present work is relatively simple and cannot take into account all aspects of such a complex phenomenon as the bond-breaking electron transfer. The results, although very promising, might be fortuitous due to the error cancellation and several problems remains unsolved. One important question is how far the interaction between the dissociating molecule and radical ion fragments would be affected by solvent, thus changing the shape of potential energy surfaces obtained in the present work in vacuo. Another problem is the very simplified solvent description. To evaluate these effects further studies with the use of other techniques are necessary. Acknowledgments This work was carried out under the HPC-EUROPA2 project (project number: 228398), with the support of the European Community – Research Infrastructure Action of the FP7.

References [1] C.P. Andrieux, I. Gallardo, J.-M. Savéant, K.-B. Su, J. Am. Chem. Soc. 108 (1986) 638. [2] C.P. Andrieux, J.-M. Savéant, K.-B. Su, J. Phys. Chem. Soc. 90 (1986) 3815. [3] J.-M. Savéant, J. Am. Chem. Soc. 109 (1987) 6788. [4] J. Bertran, I. Gallardo, M. Moreno, J.-M. Savéant, J. Am. Chem. Soc. 114 (1992) 9576. [5] J.-M. Savéant, J. Am. Chem. Soc. 114 (1992) 10595. [6] C.P. Andrieux, J.-M. Savéant, C. Tardy, J. Am. Chem. Soc. 120 (1998) 4167. [7] R.A. Marcus, J. Chem. Phys. 24 (1956) 966. [8] N.S. Hush, J. Chem. Phys. 28 (1958) 962. [9] R.A. Marcus, Annu. Rev. Phys. Chem. 15 (1964) 155. [10] R.A. Marcus, Farraday Discuss. Chem. Soc. 72 (1982) I. [11] E.D. German, A.M. Kuznetsov, V.A. Tikhomirov, J. Electroanal. Chem. 420 (1997) 235. [12] Z.Y. Zhou, Y.M. Xing, J. Mol. Struct. Theochem. 532 (2000) 87. [13] A. Ignaczak, W. Schmickler, J. Electroanal. Chem. 554–555 (2003) 201. [14] A. Ignaczak, Electrochemical reduction of tert-butyl chloride–Computational study, Electrochim. Acta (2009), doi:10.1016/j.electacta.2009.09.015. [15] W. Schmickler, J. Electroanal. Chem. 204 (1986) 31. [16] A.D. Becke, J. Chem. Phys. 104 (1996) 1040. [17] C. Adamo, V. Barone, Chem. Phys. Lett. 274 (1997) 242. [18] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785. [19] B. Miehlich, A. Savin, H. Stoll, H. Preuss, Chem. Phys. Lett. 157 (1989) 200. [20] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, Gaussian 98, Revision A.6, Gaussian, Inc., Pittsburgh, PA, 1998. [21] R.W. Dornte, J. Chem. Phys. 1 (1933) 630. [22] J.W. Beach, D.P. Stevenson, J. Am. Chem. Soc. 60 (1938) 475. [23] J.Q. Williams, W. Gordy, J. Chem. Phys. 18 (1950) 994. [24] P.W. Allen, L.E. Sutton, Acta Cryst. 3 (1950) 46. [25] H.J.M. Bowen, A. Gilchrist, L.E. Sutton, Trans. Faraday Soc. 51 (1955) 1341. [26] A.C. Legon, D.J. Millen, A. Samsam-Baktiari, J. Mol. Struct. 52 (1979) 71. [27] S. Kassi, D. Petitprez, G. Wlodarczak, J. Mol. Struct. 517–518 (2000) 375. [28] C.T. Mortimer, H.O. Pritchard, H.A. Skinner, Trans. Faraday Soc. 48 (1952) 220. [29] J.A. Kerr, Chem. Rev. 66 (1966) 465. [30] D. Griller, J.M. Kanabus-Kaminska, A. Maccoll, J. Mol. Struct. Theochem. 163 (1988) 125. [31] N. Sheppard, Trans. Faraday. Soc. 46 (1950) 527. ´ [32] A. Ignaczak, M. Karpinska, J. Mol. Struct. Theochem. 859 (2008) 98. [33] M.T.M. Koper, G.A. Voth, Chem. Phys. Lett. 282 (1998) 100. [34] A. Ignaczak, W. Schmickler, Chem. Phys. 278 (2002) 147. [35] A. Ignaczak, W. Schmickler, Electrochim. Acta 52 (2007) 5621. [36] A. Ignaczak, Electrochim. Acta 53 (2008) 2619. [37] S.M. Kast, K. Nicklas, H.-J. Bär, J. Brickmann, J. Chem. Phys. 100 (1994) 566. [38] S.M. Kast, J. Brickmann, J. Chem. Phys. 104 (1996) 3732. [39] R.M. Whitnell, K.R. Wilson, J.T. Hynes, J. Chem. Phys. 96 (1992) 5354. [40] R.K. Murarka, S. Bhattacharyya, R. Biswas, B. Bagchi, J. Chem. Phys. 110 (1999) 7365. [41] A. Bagchi, G. Srinivas, K. Miyazaki, J. Chem. Phys. 115 (2001) 4195. [42] W.R. Fawcett, C.A. Foss, Electrochim. Acta 36 (1991) 1767. [43] W.R. Fawcett, M. Opallo, Angew. Chem. Int. Ed. Engl. 33 (1994) 2131. [44] W.R. Fawcett, Electrochim. Acta 42 (1997) 833. [45] M.J. Weaver, G.E. McManis, Accounts Chem. Res. 23 (1990) 294. [46] M.J. Weaver, Chem. Rev. 92 (1992) 463. [47] A. Houmam, Chem. Rev. 108 (2008) 2180. [48] J.-M. Savéant, D. Tessier, Faraday Discuss. Chem. Soc. 74 (1982) 57. [49] S.H. Northrup, J.A. McCammon, J. Chem. Phys. 78 (1983) 987. [50] M. Berkowitz, J.D. Morgan, J.A. McCammon, S.H. Northrup, J. Chem. Phys. 79 (1983) 5563. [51] R.S. Larson, Physica 137A (1986) 295. [52] A.M. Berezhkovskii, L.M. Berezhkovskii, V.Y. Zitserman, Chem. Phys. 130 (1989) 55. [53] A.M. Berezhkovskii, V.Y. Zitserman, S.-Y. Sheu, D.-Y. Yang, J. Kuo, S.H. Lin, J. Chem. Phys. 107 (1997) 10539.