Ionic Liquid and Au

To gain a deeper insight into the nature of the elementary act, we model the same ... Our estimates of solvent reorganization energy, work terms, and ...
5 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Ferrocene/Ferrocenium Redox Couple at Au(111)/Ionic Liquid and Au(111)/Acetonitrile Interfaces: A Molecular-Level View at the Elementary Act Victoria A. Nikitina,† Sergey A. Kislenko,‡,# Renat R. Nazmutdinov,*,§ Michael D. Bronshtein,§ and Galina A. Tsirlina† †

M.V. Lomonosov Moscow State University, Leninskie Gory 1/3, 199991 Moscow, Russia Joint Institute for High Temperatures of RAS, Izhorskaya 13/19, 125412 Moscow, Russia § Kazan National Research Technological University, Kazan, 420015 Republic of Tatarstan, Russia # Institute of Problems of Chemical Physics of RAS, Academician Semenov avenue 1, 142432 Chernogolovka, Moscow, Russia ‡

S Supporting Information *

ABSTRACT: A quantum mechanical theory is employed to describe heterogeneous ferrocene (Fc)/ferrocenium (Fc+) electron transfer across Au(111)/[bmim][BF4] and Au(111)/acetonitrile interfaces. Classical molecular dynamics simulations were performed to calculate the potential of mean force for Fc and Fc+ and to estimate the solvent reorganization energy. The structure of the reaction layer and the solvent shell of Fc and Fc+ as derived from molecular dynamics are thoroughly investigated. The molecular structure of ferrocene and ferrocenium species, as well as the reactant−electrode orbital overlap are addressed on the basis of a quantum chemical approach. The experimental dielectric spectra for both types of solvents are used for quantum corrections of the outer-sphere reorganization energy as well as for estimations of the effective frequency factor in the limit of strong and weak electronic coupling. The dependence of the electronic transmission coefficient on the electrode−reactant distance is calculated for several orientations of the ferrocenium cation relative to the electrode surface which was represented by a cluster. Emphasis is put on the molecular nature of the elementary act and its qualitatively interesting features for both interfaces. The electron transfer rate constants are calculated and discussed in the viewpoint of available experimental data.

1. INTRODUCTION

employed at initial stages but is hardly meaningful for the media composed of ions rather than dipoles. Qualitatively, fast ET in ILs can be explained by either anomalously low solvent reorganization energies, and/or (for diabatic ET) by some specific structure of reaction layer, and/ or (for adiabatic ET) by contributions of some fast solvent relaxation modes. The former aspect was addressed in refs 9−12 by means of MD simulations for hypothetic spherical reactants. The main conclusion of these studies is the absence of the significant difference of medium reorganization energies in ILs and molecular solvents. This points formally to the failure of dielectric continuum models for ILs, as the Pekar factors for ILs are noticeably lower. However, no reorganization energy calculations have been performed for commonly applied real reactants considering their molecular features. Furthermore, the calculated reorganization energy values were not corrected for the effect of electronic polarization, which was

Despite the growing number of studies on electrochemical kinetics in room-temperature ionic liquids (ILs), attempts to model electron transfer (ET) in these media are still noticeably lacking. Although some attempts to address ET in IL in terms of the Marcus theory and using molecular dynamics (MD) simulations were made recently, the redox processes at an electrode/IL interface (to the best of our knowledge) were never investigated in the framework of modern quantum mechanical theories combined with quantum chemical modeling. The key challenge for theory is to explain observed high rate constants for ET in ILs, always exceeding 10−4 cm·s−1 at bare metals1−7 given very slow relaxations of this type of solvents. In the framework of traditional treatment based on the data of dielectric spectroscopy,8 the agreement of calculated rate constants for the Fc+/Fc redox reaction with experimental values for ILs cannot be achieved without rather speculative assumptions (which are to be checked). The major drawback of this analysis is a very straightforward continuum approach to the calculation of solvent reorganization energy, which can be © 2014 American Chemical Society

Received: July 20, 2013 Revised: December 6, 2013 Published: January 16, 2014 6151

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

The most reliable information is obtained when the ratio kIL/ kAN is analyzed rather than absolute rate constant values in AN (kAN) and IL (kIL). This paper is organized as follows. In Section 2, all computational details are described. In Section 3, general relationships for the ET rate constant are considered from the viewpoint of quantum-mechanical theory. In Section 4, key contributions to the rate constant and this quantity itself are calculated and discussed in comparison with available experimental data. Some concluding remarks can be found in the Section 5.

shown to be significant (see, for example, refs 13−15). Adequate theoretical models for the description of electrostatic interactions between a charged solute and the components of ionic solvent are only starting to appear16,17 and require comprehensive testing (from both experiment and atomistic simulations). The influence of solvent dynamics (friction) on the rate of ET in IL is another challenging issue. The authors11,12 investigated a homogeneous ET between two solute molecules in [emim][PF6] and AN combining the Kramers and Grote− Hynes (GH) theories with molecular dynamics simulations. They have demonstrated that solvent dynamics when crossing the barrier in the elementary act can be adequately described in terms of the GH theory. This theory presumes a timedependent solvent friction (i.e., non-Markovian effects). In this study we focus on making estimations of rate constants in IL 1-butyl-3-methylimidazolium tetrafluoroborate [bmim][BF4] on the basis of quantum mechanical theory of ET in polar solvents for a ferrocene/ferrocenium (Fc/Fc+) redox couple and compare the results with the available experimental information. To gain a deeper insight into the nature of the elementary act, we model the same system in acetonitrile (AN) in the framework of the same approach. The Fc/Fc+ redox system is one of the generally adopted model reactants for the studies of electrode kinetics in ILs. Standard rate constants reported for Fc/Fc+ at solid electrodes (Au, Pt, glassy carbon) are rather high, which induces serious complications for quantitative experiments. The reported values unavoidably scatter in a wide interval of 10−4−10−1 cm·s−1, and this scatter demonstrates no correlation with electrode material nature, mostly with the measurement technique.1,3−7 The accuracy of rate constant values is so low because they are at the edge of voltammetry capability even when microelectrodes are applied. As typical ohmic distortions of experimental voltammograms often result in too low rate constant values, it seems to be reasonable to limit the range of experimental rate constants to 10−3−10−1 cm·s−1, thus neglecting the obviously underestimated values. In acetonitrile ET is also very fast, but the scatter in the rate constant literature values is significantly smaller because of the less complicated arrangement of experiments: 1−10 cm·s−1 (see ref 18). Some more precise kinetic data on long-range Fc oxidation in ILs at alkanethiol-modified surfaces are available now as well,19,20 which can be helpful for a quantitative ET modeling. These lower values can be obtained with much better accuracy, but all technical problems of IL electrochemistry (sensitivity to air and moisture, high ohmic drops, absence of reliable reference electrodes) still can be the reason for a noticeable scatter. Key information going from long-range ET data is the distance dependence of rate constants and their components. We consider [bmim][BF4] to be the most defined IL medium for quantitative data treatment and modeling. Dielectric spectroscopy data21 and reliable all-atom force field parameters22 are available for this IL, making it one of the most “well characterized”. Our estimates of solvent reorganization energy, work terms, and reaction layer structure are obtained from molecular dynamics simulations of Fc+/Fc at Au(111) in [bmim][BF4] and in AN. We employed modern quantum mechanical theory of ET23 to address the effect of orbital overlap, as well as to correct the reorganization energy for solvent quantum modes. Using the same computational scheme for ionic and molecular solvent allows us to highlight the differences in various factors affecting ET in these polar media.

2. METHODS AND DETAILS OF MODEL CALCULATIONS 2.1. Quantum Chemical Calculations. The quantum chemical calculations of Au clusters and of ferrocene in oxidized and reduced states were performed at the density functional theory (DFT) level with the hybrid functional B3LYP as implemented in the Gaussian 09 program suite.24 A basis set of DZ quality was employed to describe the valence electrons of Au and Fe atoms. The effect of inner electrons was included in a relativistic effective core potential (LanL2) developed by Hay and Wadt.25 The standard 6-311++g(d,p) basis set was used to describe the electrons of C and H atoms. The geometry of the Fc molecule was fully optimized in the gas phase for the oxidized and reduced states without symmetry restrictions and then fixed at certain orientations relative to the Au cluster. The main results were obtained with a two-layer Au54 cluster; test calculations were performed also for two-layer Au20, Au30, and Au40 and three-layer Au66 clusters. The geometry parameters of the metal clusters were fixed at crystallographic values of gold crystalline. The solvation free energies for the Fc+ cation were calculated in terms of the self-consistent reaction field (SCRF) with the help of the Polarized Continuum Model (PCM)26−28 as implemented in Gaussian 09 code. All calculations were run on dual CPU Pentium IV workstations. 2.2. Molecular Dynamics Simulations. All MD simulations were performed using the DL_POLY classic package.29 The calculations were run on the supercomputers MVS-100 K (Joint Supercomputer Center of the Russian Academy of Sciences), SKIF MSU “Chebyshev” and “Lomonosov” (Moscow State University). The force field of Fc was taken from the work of Lopes et al.,30 while the force field parameters of Fc+ were derived by ourselves. By analogy with the Fc model, ferrocenium was treated as a rigid molecule with only one internal degree of freedomrotation of cyclopentadienyl rings around their centroids. The geometrical parameters of Fc+ molecule were taken from the above-mentioned DFT calculations. The partial atomic charges were calculated with the help of the CHELPG method31 and Boltzmann averaging over eclipsed and staggered equilibrium conformations. The van der Waals interactions are extracted from the work of Lopes et al.30 with no changes. The parameters are summarized in Table S1 (Supporting Information). Solvents were described by all-atoms models taken from the literature. Acetonitrile was described by the model of Nikitin et al.32 The Liu force field22 was used for [bmim][BF4]. Force fields for both solvents are nonpolarizable; i.e., electronic polarization of the medium is neglected, and therefore, the optical dielectric constant is equal to 1. To address the electronic polarization effect in the energy values obtained from the MD simulations, we used a scaling factor13,33,34 (see Section 4.2.2). 6152

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

MD simulations of bulk solutions were performed in the NPT ensemble under a pressure of 1 bar and at temperature of 300 K kept constant using a Nose−Hoover thermostat.35,36 We used periodic cubic boxes containing single solute species (Fc or Fc+) dissolved in 150 ion pairs of IL or 576 acetonitrile molecules. The Ewald method with a real space cutoff value of 15 Å was used to calculate the electrostatic interaction. The equations of motion were solved using the Verlet leapfrog integration algorithm with a time step of 1 fs. The cutoff radius of the van der Waals interaction was 15 Å. When performing MD simulations we started from different initial configurations for all electrolytes investigated. The final results were obtained by averaging over two configurations and over 3 ns of dynamics per configuration. To generate an equilibrated system at 300 K it was first annealed at 400 and 350 K during 0.5 ns of simulation at each temperature. Then, the temperature was reduced to 300 K, and the systems were further equilibrated for 1 ns before collecting data. To simulate the Au(111)/[bmim][BF4] and Au(111)/AN interfaces we used the simulation box shown in Figure S1 (Supporting Information). This box represents two parallel Au(111) single-crystal surfaces 35.16 × 35.53 Å2 in area and 9.57 Å in width, arranged at a distance of 62.0 Å; the space between the surfaces was filled with the solvent. The number of [bmim][BF4] ion pairs in the solvent region was 224, which corresponds to a density of 1.147 g·cm−3 in the central part of the box, obtained by an independent NPT simulation of bulk [bmim][BF4] at 350 K and 1 atm. For AN the number of solvent molecules was 851, and the density at 300 K in the central region of the box was 0.786 g·cm−3. 3D periodic boundary conditions were employed. To reduce the interaction between the simulation box and its images, the periodicity in the direction perpendicular to the surfaces was elongated to 25 nm. Thus, the periodic images were separated by a region of vacuum of ∼17 nm thick. We used the GolP force field37,38 to describe the interaction of the [bmim][BF4] IL and AN with the Au(111) surface. The force field takes into account polarization of gold (image charge effects). MD simulations of the interface in IL were performed in the NVT ensemble at a temperature of 350 K to avoid glassy behavior, whereas for AN the temperature was 300 K. Four simulations starting from different initial configurations were performed to reduce a statistical error. The interface structure was obtained by averaging over four configurations and over 3 ns of dynamics per configuration. All starting configurations were equilibrated during 1 ns of simulation at 350 K before collecting data. Some details on the calculations of potential of mean force are given in Appendix A.

equation for the rate constant in the case of strong electronic coupling (adiabatic ET)23,39−41 k=

1 δx exp(−σ )exp( −Wi /kBT ) τeff exp[− (λ + Wf − Wi − e0η)2 /4λkBT ]

(1)

where Wi and Wf are the work terms for reactant and product; λ is the total reorganization energy; τeff is the effective time scale which characterizes the solvent relaxation times; exp(−σ) is tunneling factor (the physical meaning will be explained below); δx is reaction volume; and η is overvoltage. In the case of weak orbital overlap (heterogeneous ET in nonadiabatic and intermediate limits) and for a given reactant− electrode separation in a fixed orientation, it is convenient to recast the rate constant first in the similar form ω k = eff κeδx exp(−σ )exp( −Wi /kBT )exp( −ΔEa /kBT ) 2π (2)

where ωeff is the classic polarization frequency; the averaged activation barrier is ΔEa; and the electronic transmission coefficient is κe. Pertinent details of the calculation of ΔEa and κe are presented in Appendix B. At zero overvoltage ΔEa ≈ (λ + Wf − Wi)2/4λ, but it is not true for higher overvoltage. It is important to stress that in the ET nonadiabatic limit ωeff in eq 2 is reduced by the ℏωeff quantity in the Landau−Zener factor (see Appendix B). The resulting equation for the rate constant does not depend on the polarization frequency, and we need the ωeff values only to calculate the electronic transmission coefficient. Nevertheless, we will use further eq 2 for the sake of convenience. Since the inner sphere contribution to the total reorganization energy is negligibly small for Fc+/Fc (ca. 0.5 kJ·mol−1),42 λ in eqs 1 and 2 results from the medium response solely. The outer-sphere reorganization energy was computed with the help of MD simulations (see the next sections) and corrected for the quantum solvent modes, i.e., reduced by a factor of ξ which is below 1. The meaning of such a correction is cutting off the contribution of high-frequency solvent degrees of freedom from the classical activation barrier which makes it smaller. This factor, in turn, can be calculated using experimental complex dielectric spectra, ε(ω)23 ξ=

2 πC

∫0

ω*

Im ε(ω) dω ω |ε(ω)|2

(3)

where C is the Pekar factor; C = 1/ε∞ − 1/ε0; and ε∞ and ε0 are optical (fast) and static (slow) dielectric constants. In our calculations we used the dielectric spectra described in refs 21 and 43, and the fitting functions proposed by the authors were adopted. ε∞ values were obtained from refractivity measurements (2.02 for IL44 and 1.79 for AN45) rather than from the fitting of experimental dielectric spectra because the accuracy of fast mode extraction in subpicsecond regions is not too high. Another effect of the correction to quantum solvent modes arises from the tunneling factor exp(−σ) in eqs 1 and 2

3. GENERAL RELATIONSHIPS According to the ET theory, there are several basic quantities responsible for the rate constant value, namely, the work terms of reactant and product, reorganization energy, and preexponential factor. The latter includes electronic transmission coefficient and (in general) effective frequency. All these quantities are interrelated and crucially dependent on molecular structure of the reaction layer, which requires special modeling efforts. The key quantities depend on both metal and solvent properties, and the rate constant calculation requires selfconsistent models of both constituents. If the linear response from both the medium and the reactant outer sphere is assumed, it is convenient to write the following

σ=

2λs πC



∫ω*

Im ε(ω) dω ℏω 2 |ε(ω)|2

(4)

where λs is the outer-sphere reorganization energy which is not corrected to factor ξ. 6153

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

The integration limit ω* in eqs 3 and 4 roughly separates classic and quantum frequency regions. We found the ω* value (ca. 150 ps−1) by constructing two characteristic plots of the course-value function F(θ, ω)23 F(θ , ω) =

ch{β ℏω/2} − ch{(1 − 2θ )β ℏω/2} sh{β ℏω/2}

IL and AN as contacting zero-charged Au(111) using a force field that accounts for Au polarization, to include solvent− metal specific interactions and allow comparison with available experimental information (mostly reported for gold electrodes). In general, our results are in qualitative agreement with the simulation at other solid surfaces. The mass and charge densities of the ionic liquid near an uncharged Au(111) surface exhibit damped oscillations (with a period of ca. 5 Å) up to the distance 20 Å (Figure 1a). At

(5)

where θ is the transfer coefficient (close to 0.5 at low overvoltage) and β = 1/kBT. The course-value function results from a statistical averaging of an ensemble of oscillators and tends to 1 in a quantum limit. This ω* value is in a good agreement with previous estimates known from the literature.23 The factor ξ is responsible for decreasing the medium reorganization energy due to the cutting out of solvent quantum modes from the ET activation barrier, whereas the factor exp(−σ) addresses the probability of tunneling of atoms through the Franck−Condon barrier. The ξ and σ values are collected in Table 1. Due to well-known Table 1. Model Parameters Contributing to the Rate Constant of Fc+/Fc Reduction/Oxidation at the Au(111)/IL and Au(111)/AN Interfaces Calculated for Adiabatic (a) and Nonadiabatic (na) Regimes [bmim+][BF4−] parameters

AN

a

na

a

na

x0, Å Wi, eV Wf, eV ξ σ λbulk, eVc

4.5 0.16 0.

8.5 0.07 0.

5 −0.1 −0.1

9 0. 0.

λs, eVd −1 τ−1 eff (ωeff), ps κe

0.45 0.9 1

0.7 0.38 0.80a 0.68b

0.95 0.18 1.15a 1.16b 0.62 24.5 0.01

0.69 79.1 1

Figure 1. Mass and charge densities of the [bmim][BF4] (a) and AN (b) along the normal to the Au(111) surface.

0.91 28.9 0.01

longer distances the density recovers its bulk value. Three dense layers can be clearly observed near the uncharged surface (each of ca. 4−5 Å thickness). In AN, the amplitudes of mass and charge oscillations appear to be significantly lower, but some structuring near the metal surface is also found (Figure 1b). The concentration profiles of imidazolium rings and [BF4]− anions along the normal to the Au(111) surface show that the centers of the closest anions and cations are both located at ca. 3.5 Å from the metal surface, and no charge alternation is observed between the highest mass density layers (Figure 2). These layers are composed of both cations and anions, although the distance of closest approach is a little shorter for cations, which results in a first positive peak in

a

Calculated on the basis of the MD simulations. bCalculated using a continuum model. cSee Section 4.2.2. dCorrected to image term and quantum solvent modes (See Section 4.2.2).

experimental difficulties, reliable data on dielectric spectra are not available for a high-frequency region (for 2−3 order interval preceding ω* value), and the accuracy of calculations using eqs 3 and 4 becomes sensitive to some detail of fitting the spectra. This decreases, of course, the level of quantitative predictions. As emphasis in the present work is put on a comparative study, we can reconcile with such a consequence.

4. RESULTS AND DISCUSSION 4.1. Structure of the Reaction Layer. The number of simulations of ILs at solid interfaces is considerably high. Typically, graphite (or graphene) or the uniformly charged wall are treated in atomistic simulations of ILs to study the ion rearrangements when changing the charge of the electrode.46−51 However, experimental atomic force microscopy (AFM), scanning tunneling microscopy (STM), and electrochemical studies which provide the most valuable information regarding the structure and dynamics of the electrode/IL interface are typically related to metal single-crystal electrodes (typically gold).52−57 Capacitance data for gold electrodes are also available for comparison with the results of simulations and theoretical estimates.56,58−60 In the present work we simulated

Figure 2. Concentration profiles of imidazolium rings and anions [BF4]− along the normal to the Au(111) surface. 6154

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

charge density profile. In AN, the charge in the liquid is also oscillating near the uncharged metal surface due to specific interactions of Au with the positively charged CH3 group of AN. The results for IL are in agreement with previous simulations of ILs near an uncharged solid surface46−49 and experimental AFM measurements.52,55−57 In all of these experiments and simulations the layering of IL is observed even near the uncharged surface, and it becomes more significant when the electrode charge is increased. The next question to be addressed is what is the free charge on the electrode at zero overvoltage for the Fc+/Fc reaction at Au(111) in AN and [bmim][BF4]. The potential of zero charge (p.z.c.) of Au(111) in AN is ca. +0.11 V versus Fc+/Fc;61−63 therefore, both cathodic and anodic reactions under study take place in the region of negative surface charges, at least at low overvoltage. As the specific adsorption of AN is considered to be relatively weak at the formal potential of the Fc+/Fc couple,64,65 no significant additional structuring near the interface is expected at negative surface charges as compared to zero charge. The p.z.c. in [bmim][BF4] at Au(111) is still uncertain, but it is known to be in the range from −0.23 to +0.08 V versus AgCl/Ag, as follows from differential capacitance measurements,58 i.e., in the range from −0.39 to −0.08 V versus Fc+/Fc. Both anodic and cathodic reactions in IL proceed, therefore, in the region of not too high positive surface charges if low overvoltage is assumed. According to ref 49 the charge on the electrode is screened by a layer of ions, which tends to overcompensate it. Clearly, charging the Au(111) electrode positively would result in a significant rearrangement of ions. However, as it was pointed out in ref 49, the changes refer to the amplitudes of the density oscillations, not the positions of the maxima and minima. As we do not know exactly what charges correspond to the region of Fc+/Fc redox potential, we calculate the rate constant for the uncharged surface, bearing in mind that for a charged surface the barriers for the ion to approach the metal surface can be different to some extent (in ref 49 the difference was pointed out to be less than 5 kJ· mol−1). Since reliable models for a charged metal electrode/IL interface to calculate the electrostatic work terms are absent so far, the latter were determined from the PMF calculations. Due to the limitations in accuracy of MD energy computations, the results obtained will be regarded more as qualitative rather than quantitative. Due to the charge density and mass density oscillations, PMF profiles in IL and AN also show several minima and maxima (Figure 3). The PMF profiles in IL exhibit two minima: at 4.5 and 8.5 Å (Figure 3a). The transition of the reacting species from the minimum at 8.5 Å to the closer to surface minimum at 4.5 Å is associated with an energy barrier. The most striking feature is that this energy barrier is by 15 kJ·mol−1 higher for Fc+ than for Fc (ca. 30 and 15 kJ·mol−1, respectively). We have found that for the Fc+ ion in IL equilibration takes noticeably longer time as compared with a neutral Fc (as well as with both species in AN); this was taken into account when calculating the PMF. As the model surface is uncharged and the position of the first minimum does not correspond to the contact approach of Fc to the electrode, we conclude that solvent−solute interactions mostly determine the barriers for the reactant and product approach, and the contribution of the reactant− electrode interactions is less significant. Note that the mean

Figure 3. Potential of mean force for Fc and Fc+ in [bmim][BF4] (a) and AN (b) and shematic representation of Fc near the Au surface (c).

force potential is a free energy profile, and it is difficult to judge about the interplay between entropy and enthalpy terms without a special analysis. When charging the electrode we expect the increase in the barrier height, as was demonstrated in simulations of spherical charged probes in the vicinity of charged and uncharged graphite.49 In the study of Lynden-Bell et al.49 the difference in the PMF barriers impeding the positively charged and neutral probe from approaching the wall is ca. 5 kJ·mol−1 for the uncharged wall in [dmim][Cl] IL and 2 kJ·mol−1 for the positively charged wall. When the charge on the positive wall is further increased, the barrier for the neutral probe becomes larger than that for the positive probe. Resting on the results presented in ref 49, we do not consider the difference in barriers near the charged and uncharged interface, as we do not know the exact charges on the metal at the potential of Fc+/Fc redox reaction. As the positions of the minima do not change significantly upon changing the charge on the electrode, the effect is hardly essential. We use the computed profiles to obtain approximate distances of Fc+(Fc) closest approach x0 and an estimate of the work terms Wi (Fc+) and Wf (Fc). The PMF profiles for Fc+ and Fc interacting with AN almost coincide (Figure 3b). The two minima in the PMF curves are located almost at the same distances x0 as for IL: 5 and 9 Å. The barrier separating these minima is also close to that in IL (30 kJ·mol−1). There is no difference, however, in the barrier heights for Fc+ and Fc. In what follows we assume equilibrium between two minima at the PMF when calculating work terms. 6155

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

4.2. Solvation of the Reactants. 4.2.1. Equilibrium Solvation. To get molecular level insight into the solvation of the reactants we performed MD simulations for Fc and Fc+ in the bulk of [bmim][BF4] and AN. Figure 4 shows radial

(6)

significantly. The Fc+ solvation shell is still constituted by a mixture of cations and anions, the latter being in a slight excess. Within the 8.0 Å from the Fe+ there are formally 6.6 anions and 6.2 cations. As we can see from these coordination numbers, the screening of charge is complete in the first solvation shell of the ferrocenium cation, which is formed by both the anions and cations, and at distances larger than 10 Å the solvent is rather structureless. The first solvation shell of Fc+ is more dense compared to that for the molecule. The distance of the anions’ closest approach to reacting species is reduced by ca. 1 Å, while the positions of the cations remain almost unchanged. However, the shell of Fc+ appears to be less structured, as the heights of the second and third maxima are less than those for neutral Fc. In AN, rearrangement of the solvation shell upon oxidation of Fc manifests itself in the reorientation of AN dipoles. When a positive charge appears at the Fc molecule, the negatively charged nitrogen atoms approach closer to the Fe center. This rearrangement induces small changes in the peak positions in g(r) functions (see Figure S4, Supporting Information). Our simulations of the equilibrium solvation shells confirm the difference in reorganization scheme in two solvents: in ILs the solvation shell can be restructured via translational motion of ions, while in AN the solvent dipoles should reorient. 4.2.2. Solvent Reorganization. We have to stress from the very beginning that the linear response theory (LRT) and, therefore, the solvent reorganization energy concept can be used in our case as the crude approximation. In the vicinity of a layered metal/IL interface the LRT works poorly, as nongauss fluctuations of solvent molecules are highly probable. Redox reactions at the aqueous interface of an active protein might be another example of systems where nonlinear solvent effects take place.66 Deviations from the LRT were thoroughly analyzed in the literature.67,68 The direct calculation of the ET free energy surfaces which address a nonlinear solvent response (for example from MD simulations) is, however, an exceedingly complicated computational issue. That is why we discuss the activation barriers of the reaction under study in terms of the solvent reorganization energy, which is enough to obtain qualitative interesting results in a comparative study. To investigate the solvent reorganization, we have used the approach developed by other authors.34,69 (see Appendix C for details). To obtain the activation energy of the ET reaction the free energy surfaces were fitted by parabolas in the vicinity of their minima (Figure 5). The distribution of fluctuations in

where the upper limit of integration was 8.0 Å (first minimum in the g(r), corresponding to the first solvation shell radius) and ρ is the bulk number density of IL (the number of ions per volume). Within the distance of 8 Å the solvation shell of the neutral Fc molecule comprises ca. 6 anions and 5.5 cations; i.e., it appears to be essentially neutral, with a slight excess of anions (this excess is more probably caused by a somewhat arbitrary choice of the solvation shell limit than by physically meaningful factors). The structure of the first solvation shell of the Fc+ cation differs to some extent from the solvent distribution around the neutral ferrocene molecule. First of all, the distance of the anions’ closest approach is reduced to 5 Å, while the centers of imidazolium rings are still located within 6 Å from the Fe center. The position of the first minimum does not change

Figure 5. Free energy surfaces along the reaction coordinate calculated for Fc and Fc+ in bulk ionic liquid [bmim][BF4] (T = 400 K) and bulk AN (T = 300 K).

Figure 4. Radial distribution functions g(r) for the Fe atom and the centers of mass of the anion and imidazolium ring around Fc and Fc+ in [bmim][BF4]. Coordination numbers ncoord for cations and anions are shown by dashed lines.

distribution functions g(r) of cations and anions of ILs around the neutral Fc (Figure 4a) and Fc+ cation (Figure 4b). The snapshots of the first solvation shells in IL and g(r) in AN can be found in the SI (Figures S3 and S4). From the g(r) in IL one can conclude that no drastic solvation shell rearrangement is observed upon charging Fc molecule. For the neutral Fc molecule centers of the anions are located at the distance of ca. 6 Å from the Fe atom, and the distance to the centers of imidazolium rings is only slightly different from 6 Å, indicating that there is no preferential Fc−anion or Fc−cation interaction. The coordination numbers ncoord were calculated as follows ncoord = ρ

∫o

R1

g (r )4πr 2dr

6156

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

were made and λbulk values in ILs were found to be nearly equal to or higher than those in AN. The resulting corrected for image terms λs values (as well as uncorrected λbulk values) are listed in Table 1 for the distances corresponding to two minima in the PMF: 5 and 9 Å. The minimum closest to the electrode surface is associated with adiabatic ET (a), while the minimum at 9 Å relates to ET in the nonadiabatic (na) regime. The justification of such a consideration is provided in the next sections. 4.3. Effective Frequency Factor. In the quantum mechanical theory of ET the effective frequency factor is determined by the classical part of the medium polarization

solvent order parameter can be found in the SI (Figure S5). It follows from our analysis that the main contribution to the solvent reorganization originates from the nearest coordination shell of Fc(Fc+), which is most likely caused by very efficient electrostatic screening. A strong asymmetry of the free energy surfaces is observed; this conspicuous feature calls for at least qualitative explanations which are also given in Appendix C. Similar asymmetric features were observed in simulations of reaction involving neutral and charged halogen atoms70 and O2−/O2 reaction in water.67 At the same time Linden-Bell9 found no noticeable reorganization asymmetry in MD simulations of simple spherical solutes interacting with IL. In general, the reorganization energy asymmetry has to be taken into account when computing the rate constant values. Rate constants computed in such a way can be compared with the experimental rate constant values measured in a wide overvoltage range. Since the only experimental information we possess refers to zero overvoltage, and the accuracy of measured k values is very low, we will limit ourselves to the calculation of an averaged reorganization energy, using the value of the barrier in Figure 5 to estimate λs/4 values. The barrier heights in AN and IL are nearly equal (ca. 55 kJ· mol−1). The similarity of the reorganization energies of ILs and molecular solvents was reported in a number of studies. For a spherical small probe bearing +1 charge the reorganization energy was found to be ca. 300 kJ·mol−1 in both [dmim][PF6] IL and AN.9 This estimate is in satisfactory agreement with our values (220 kJ·mol−1) taking into account that the probe in ref 9 is smaller than Fc. In simulations of Shim and Kim12 the barrier height in IL was found to be by 10−15 kJ·mol−1 higher as compared with that in AN. It was shown that application of nonpolarizable force fields results in substantial overestimation of the λs value.13−15 Following this approach we used a scaling factor (1/ε∞) to account for the difference in the optical dielectric constants of IL (see Section 2.2). The λs value in AN computed from MD was scaled in the same way. Remember that ε∞ values for both solvents were taken from refractivity data. In this work, the homogeneous redox reaction Fc ↔ Fc+ + e in the IL [bmim][BF 4 ] and acetonitrile was initially investigated. At this stage we considered only the solute species (Fc or Fc+) in solvent bulk. For heterogeneous ET, the interaction between the solute and a metal electrode was not taken into account explicitly due to the computational difficulties (see the top of this section). When we go to calculations of the ET rate constant at the Au(111) interface (Section 4, see below), the reorganization energy obtained from bulk simulation is corrected by the image term λs = λbulk − C(e20/4x), where x is the distance from the solute center to the image plane. The bulk solvent reorganization energy value estimated from MD λbulk may be compared with the predictions of continuum 2 theory. The corresponding λcont bulk value was estimated as e0C/ + 2reff, where reff is the effective radius of Fc calculated from the value of electrostatic solvation energy obtained by using the PCM (see ref 71 for details). For AN, λbulk and λcont bulk values are very close (ca. 1.15 eV). For IL the difference in reorganization energies is more pronounced (λcont bulk is ca. 0.1 eV lower than λbulk). Both from MD and the continuum model we obtain a similar trend: the reorganization energy in IL is significantly lower than that in AN. In this respect our results differ from the λs estimates reported in refs 9−12, where no corrections for the nonpolarizability of the force field and solvent quantum modes

ωeff2 =

2 πC

∫0

ω*

ω Im ε(ω)dω || ε(ω)||2

(7)

where the integration limit ω* has the same meaning as was already discussed above (see explanations to eq 5). We employed this equation to estimate the electronic transmission coefficient (see the next section) and (formally) the rate constant k in the case of weak electronic coupling. The calculated values of ωeff are shown in Table 1. As can be seen from the table, for the IL the fast modes contribute mainly to the frequency factor. Moreover, eq 7 predicts close ωeff values for both solvents. According to MD simulations performed by the authors11,12 the frequency factors obtained for the RTIL and AN are also very close. In the case of strong electron coupling (adiabatic ET limit), solvent dynamics effects may play an important role. It was shown earlier in refs 39−41 and 72 that when solvent dynamics affects the ET rate the effective frequency factor is governed by either slow or fast solvent modes. As was already mentioned above, the Grote−Hynes theory can be used as well. It is more convenient, however, to employ the formalism developed by the authors in ref 73 for the adiabatic ET. In terms of this approach the frequency factor is calculated as a reciprocal solvent time scale function, τ(q), which depends on the solvent correlation function M(τ) (see details in Appendix D) τ(q) = exp( −q2 /2)

∫0





⎧ ⎫ ⎡ q 2 M (τ ) ⎤ ⎪ ⎪ 1 ⎨ ⎥ − 1⎬ exp⎢ ⎪ ⎪ 2 ⎣ 1 + M (τ ) ⎦ ⎩ 1 − M (τ ) ⎭

(8)

where at zero electrode overvoltage q = λs/2kBT. It can be readily shown that using the experimental dielectric spectrum43 and the inverse Laplace transform technique M(τ) for AN (“two-Debye” spectrum) is exactly expanded into the sum 2

2

M(τ ) =

∑ δi exp(−τ /τi*) i=1

(9)

where δi are contributions to the solvent reorganization energy 2 from the ith solvent mode (∑i = 1 δi = 1) and τ*i are their characteristic times (see Table 2). For [bmim][BF4] (dielectric spectrum with “non-Debye” contributions21) M(τ) was calculated numerically and then fitted by three exponents 3

M(τ ) ≈

∑ δi exp(−τ /τi*) i=1

6157

(10)

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

orientation only one frontier orbital (the first or the second one) gives the most significant contribution to the κe, and the other contribution can be neglected. Model lg κe versus x dependencies are displayed in Figure 7; scatters of points are to

Table 2. Parameters of the Correlation Function M(τ) Calculated for AN and [bmim][BF4] AN IL

τi*, ps

δi

0.08; 0.4 0.35; 13; 185

0.73; 0.27 0.55; 0.25; 0.2

The corresponding parameters are listed in Table 2. Note that for this case the normalizing of δi values is not attained automatically like for exact expansion, and we have to do this additionally. Some examples of solvent correlation functions constructed for different solutes by using the coarse-grained ionic liquid model and molecular dynamics simulations are presented in ref 11, where the authors employed another type of fitting. The characteristic times of M(τ) obtained for [bmim][BF4] agree qualitatively with those reported by the authors74 on the basis of MD simulations performed for 1,3dimethylimidazolium chloride. The model calculations predict a noticeable difference in ωeff values for [bmim][BF4] and AN. One can conclude that for the IL the middle time (τ2*) contributes basically to the resulting frequency factor in the adiabatic region. 4.4. Electronic Transmission Coefficient. The model dependence of κe versus the Au(111) surface−Fc separation was fitted by exponential decay κe = κ0 exp( −βx)

Figure 7. Lg κe vs z dependencies calculated for different orientations of Fc+ relative to the model Au cluster; linear approximations are shown with dotted line.

be attributed to the surface corrugation. Four lg κe vs x dependencies were averaged (dotted line in Figure 7). The dependence of the average between the orientations κ̃e vs x shows two distinct regions: at x < 6 Å, κ̃e is nearly constant, and the exponential approximation is inaccurate. This corresponds to the adiabatic regime of ET. At x > 6 Å, κ̃e vs x dependence is exponential and can be associated with diabatic ET. The β value obtained from the slope at x > 6 Å is ca. 1.2 Å−1. This is in good agreement with results for long-distance ET at the thiolmodified Au(111)/IL interface.20 4.5. Calculated versus Experimental Rate Constants. Comparing the model and experimental rate constants is a complicated issue, and any pertinent attempt leads to a number of approximations and assumptions. Fortunately, some approximations from the computational side (e.g., the ignored asymmetry of reorganization energy and nonzero electrode charge) are less crucial just for the experimental data available only for zero overvoltage and not too high electrode charge. Again, our model estimates of the rate constants are related formally to reduction because the electrode−Fc+ overlap was computed. Note, however, that the widely applied Nicholson approach75 to evaluate the rate constant from cyclic voltammetry data results in the apparent rate constant kapp. This quantity has a sense of the exchange current density, where oxidation and reduction rate constants (kox and kred) are combined and “screened”. The contribution from work terms is included implicitly in this effective rate constant. As can be readily seen from eqs 1 and 2 applied to kox and kred, when combining work terms (Wi and Wf) in two exponents, contributions of these terms do not induce any difference between rate constant values. It should also be kept in mind that ET taking place in the course of reduction and oxidation is described by the same transmission coefficient. Actually, in the vicinity of the crossing point of the reaction free energy surfaces the resonance splitting is the same for oxidation and reduction processes. These circumstances enable a very straightforward comparison with experiment; additional assumptions are not needed. In addition to the kinetic parameters considered above the reaction volume δx is necessary to calculate the rate constant. To obtain an estimate of this quantity for the both ET regimes, we computed k(x) using the distance-dependent values of λs, κe, Wi, and Wf. As the integration interval of distance-dependent rate constant includes an intermediate area (separating

(11)

Calculations were performed for four different orientations of Fc relative to the Au54 cluster (Figure 6; orientation considered

Figure 6. Different orientations of Fc+ relative to the model Au cluster. Tilt angle is calculated between the axis of Fc+ (connecting two rings) and the Au surface.

above corresponds to Figure 6a) using an original code written with the help of the Mathematica 7 program suite. A ferrocene molecule is known to have two degenerated frontier molecular orbitals (see Figure S6 in SI), and meaningful estimations require considering the contribution of both orbitals to the kinetic parameters. We found, however, that for each 6158

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

• less noticeable difference of the rate constants in two different ET limits resulting in part from a high effective frequency in the nonadiabatic regime; • relatively high solvent reorganization energy as compared with predictions based on continuum models (for example, using the Marcus formula) but still lower than that in molecular solvent, at least AN. These features are probably more important than the seeming differences in solvent dynamics of ILs and molecular solvent (roughly treated as viscosity effect and related to the slowest relaxation time). As was already stressed, according to our calculations the effective IL frequency in the adiabatic limit is significantly higher than that for its slowest relaxation mode as measured in dielectric spectra. When comparing the ET reactions occurring under adiabatic mode in ILs and molecular solvents, we see that the rate constants in molecular solvent like AN are noticeably higher, which might be attributed to some extent to solvent dynamics effect (τ−1 eff factor, see eq 1). If, however, nonadiabatic reactions in both media are considered, the differences between IL and molecular media become less significant, and for certain distances of the closest approach the rate constants in ILs can become even higher than those in molecular solvent. We argue that the literature data on rather high rate constants in ILs for various reactions (as compared with the same redox processes in molecular liquids) are probably related to ET proceeding in the intermediate or nonadiabatic regime. All the stated conclusions can be maintained only on the basis of a combined consideration of all factors and parameters affecting the rate constants, which is possible only in the framework of quantum mechanical theory. A reliable test of our predictions might be possible with the appearance of high precision kinetic data in a wide overvoltage region for the same reactants in IL and AN at bare metal surfaces, which are absent so far. Alternatively, modeling of the metal/alkanethiol/IL interface is of interest in the future, as the experimental data for slower ET in the presence of thiol artificial barriers are more quantitative and available for nonzero overvoltage. We also expect that new experimental data on the dielectric spectra of solvent will be available in the nearest future which minimizes the uncertainty window in high-frequency region. In this work we combined molecular dynamics simulations, quantum chemical approaches, and DFT calculations with modern theory of ET, which plays the leading role. Although such a way is in general very powerful and attractive and makes it possible to gain a deeper insight into the elementary act at a proper molecular level, the quantitatively reliable prediction of observable quantities still cannot be made. That is why we are focused mostly on qualitatively interesting trends; some of them are compelling albeit call for additional support. It would be very tempting in the future to address solvent dynamics effects on ET kinetics in IL in more detail. Another challenging issue is the possible breakdown of the linear response theory at a metal electrode/ionic liquid interface. With the advent of more powerful computer facilities it becomes possible to perform ab initio MD simulations at charged interfaces which are so promising for the microscopic modeling of interfacial ET kinetics.

adiabatic and nonadiabatic regions), we have to deal with an interpolation scheme (an example of an interpolation scheme for homogeneous ET processes can be found in ref 76). In our calculations a simple linear interpolation of the rate constants in the intermediate region was used; we found that more sophisticated schemes do not affect the general observed trend. x bulk The integrals ∫ k(x)dx were computed for x0 values in x0

adiabatic and nonadiabatic limits, while the xbulk value amounts to 20 Å. This enables us to evaluate the reaction volumes δx (see Table 3). The resulting δx values demonstrated a minor dependence on the width of the intermediate region. Then the rate constants kIL and kAN were calculated according to eqs 3 and 5. Table 3. Rate Constants (cm·s−1) and Reaction Volumes (δx, Å) Calculated for Ferrocenium Reduction at the Au(111) Surface ET regime

δxIL

δxAN

kIL

kAN

adiabatic nonadiabatic experiment

1.6 0.9

0.3 1.5

0.1 0.05 (10−3−10−1)1,3−7

1000 0.15 (1−10)18

The integration of rate constant over x in the adiabatic limit for AN results in the reaction volume of ca. 0.3 Å, which seems to be only slightly lower than the conventional value for molecular liquids. The computed rate constant in AN in adiabatic limit was found to be remarkably high, ca. 1000 cm· s−1. This agrees qualitatively with the experimental observations of very fast electrode kinetics. Although the experimental rate constant value ranges in the interval 1−10 cm·s−1, taking into account the limits of the experimental method, as well as a number of model assumptions (such as surface charge neglecting), a satisfactory agreement with experiment can be argued. At the distance of 9 Å (nonadiabatic limit) the rate constant is 4 orders of magnitude lower (0.15 cm·s−1), and the reaction volume was calculated to be ca. 1.5 Å. The experimental data in AN might be attributed, therefore, to the adiabatic ET regime. In IL the δx values amount to 1.6 and 0.9 Å in adiabatic and nonadiabatic limits, respectively. The calculated rate constants are 0.1 cm·s−1 (adiabatic) and 0.05 cm·s−1 (nonadiabatic). Both values are close to the upper limit of rate constants evaluated from experiment. Because of the limited accuracy of estimations, as well as a scatter in the experimental values, it becomes difficult in this case to distinguish between adiabatic and nonadiabatic regimes. Formally our results do not contradict experimental work2 where the authors have maintained the adiabatic nature of ET for a Fc+/Fc redox couple at the Au(111)/IL interface. Anyway, more evidence is needed to corroborate the computational prediction and to assign the experimental information more reliably. Despite of all the mentioned uncertainties the ratio kAN/kIL was estimated to be 104, which falls into the interval of experimentally observed ratios.

5. CONCLUDING REMARKS Resting on our calculations of rate constants we can list the following main specific features of ET in ionic liquid ([bmim][BF4]): • higher barrier impeding charged reactants to approach the electrode (in particular high Wi for Fc+);



APPENDIX A. POTENTIAL OF MEAN FORCE To compute the potential of mean force (PMF) for Fc and Fc+ in [bmim][BF4] and AN as a function of a distance to the 6159

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

base square did not exceed 8−10 Å2. The number of such layers is an input parameter (ca. 20−30 in this work). When performing calculations the ferrocenium cation in a given orientation was scanned parallel to each layer with some step (the typical number of scan points ranged from 10 to 15 for one layer). Due to the surface corrugation, the rate constants calculated for one (y, z) layer noticeably differ, and on the second step we select the maximal k value within each layer (kmax), so that finally we have a reduced array of the rate constants for different reactant−electrode separations. Then these rate constants can be integrated over distance x and discussed in terms of the reaction volume δx. It should be noted that each point of the grid was calculated using eq B1 with summation over all electron energy levels in the metal cluster, N (for the Au54 cluster N amounts to 503). According to our results, at zero overvoltage only several electron levels lying in the vicinity of the HOMO of the metal cluster contribute to the resulting rate constant. At the final step we recast the results obtained previously in the form like eq 3. To extract the averaged distance-dependent electronic transmission coefficients, for each kmax we constructed the Arrhenius plot value in a narrow temperature interval. Of course, this should be regarded more as a computational expedient than modeling a real temperature effect on ET kinetics. Another way to average the “partial” electronic transmission coefficient was used by Gosavi and Marcus.80,81 The authors80,81 treated an electrode as a solid and integrated the squared electronic matrix element in the space of k-vectors. The “partial” electronic transmission coefficients (κ(j) e , see eq B1) were treated in the framework of the Landau−Zener theory23

Au(111) surface we employed one of the variants of restrained MD simulation and thermodynamic integration.78,79 The PMF was calculated by integrating the average force in the direction perpendicular to the surface ⟨f x(x)⟩, acting on the whole solute (Fc or Fc+) W (x ) = −

∫x

x

⟨fx (x)⟩dx

(A1)

bulk

where xbulk is any position in the bulk ionic liquid. To estimate the average force, we ran a series of MD simulations where the solute molecule is restrained at different distances from the surface with the step of 0.5 Å. Simulations were performed for the fixed molecule orientation, in which cyclopentadienyl rings were arranged perpendicular to the Au surface. The orientation and distance of Fc (Fc+) from the surface were fixed with the harmonic potentials U(x) = k(x − xi)2 (k = 3200 kJ·mol−1·Å−2), acting on virtual atoms placed in centroids of the rings (Figure S2 in the Supporting Information). On average, the net force ⟨f x(x)⟩ on the solute due to the solvent and the surface is balanced by that due to the harmonic res restraints ⟨f res x (x)⟩, i.e., ⟨f x(x)⟩ + ⟨f x (x)⟩ = 0. So, we can obtain ⟨f x(x)⟩ indirectly using the equation ⟨fx (x)⟩ = −⟨f xres (x)⟩

(A2)

⟨f res x (x)⟩

At each distance x from the surface the force was averaged over 1 ns after equilibration during 0.5 ns. For Fc+ in IL the equilibration period was extended to 1.5 ns due to the slow IL relaxation.



APPENDIX B. CALCULATION OF THE RATE CONSTANT OF HETEROGENEOUS ET IN THE CASE OF WEAK COUPLING Since we employ the cluster model of a metal surface, it is convenient to use the relation for the ET rate constant as a finite sum

κe(j) ≈ 1 − exp{−2πγe(εj)}

The Landau−Zener factor γe was calculated as follows γe(εj) =

N

k=

∑ j=1

κe(j)

(B2)

2 |Vif (εj)|2 ℏωeff

π λskBT

(B3)

where Vif(εj) = ∫ ψi(εj)Vψ̂ fdΩ is the resonance integral which was computed on the basis of perturbation theory (ψi(εj) is the j-th molecular orbital of the model cluster corresponding to the eigen value εj; ψf is the acceptor molecular orbital of reactant; and V̂ is the perturbation operator); ωeff is the effective frequency calculated by using eq 7; and λs is the distancedependent solvent reorganization energy, corrected for image terms and quantum solvent modes. A factor 2 in eq B3 results from the doubly occupied energy levels of model clusters with odd number of electrons. The perturbation operator V̂ is defined as a molecular potential induced by the ferrocenium ion (Fc+) in oxidized state:

2

exp[−(λ + εF − εj − e0η) /4λkBT ] (B1)

where εj is the j-th electronic energy level of a metal electrode; εF is the Fermi level of metal, i.e., the highest occupied molecular orbital (HOMO) of the cluster; κ(j) e is the electronic transmission coefficient resulting from the overlap of the j-th orbital of the metal with the acceptor molecular orbital of a reactant; and N is the number of the electron energy levels in the metal cluster. Usually a more familiar relation is employed for the rate constant in the nonadiabatic limit,23 which results from the integration over all the electronic energy levels in a metal; it is convenient to use (the density of electronic states of the metal electrode in this case). The following scheme was employed to calculate the ET rate constant in the nonadiabatic limit. As the cluster model of the metal electrode we employed enables us to address the surface corrugation, at the first step we computed an array of the k values, which form a grid in some elongated parallelepiped. The parallelepiped is located in a liquid layer closest to the surface, with the base parallel to the first layer of the metal cluster; the base center coincides with the center of the cluster first layer. For the sake of convenience, the parallelepiped is divided into sections, i.e., (y, z) layers, which are parallel to its base and aligned along its height. In praxis the

M

V ̂ ( r ⃗) =

∑ m=1

qm |R⃗ m − r |⃗

(B4)

where qm are the atomic charges calculated using the CHELPG scheme; the radius-vector R⃗ m describes the position of the atoms; and M is the number of atoms. Although in general the perturbation V̂ (r)⃗ can be readily corrected to the effect of screening of an electric field by fast solvent modes using a simple scaling factor (1/ε∞), we did not take it into account in our comparative study. Some pertinent details of the above-mentioned model can be found in ref 82 as well. There is, however, an important difference between the 6160

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

two model approaches. The authors of ref 82 constructed the effective wave function of a metal electrode using the distancedependent total electronic density obtained from periodical DFT calculations. In our work we use a two-layer Au54(27 + 27) cluster which mimics the (111) gold face, to calculate the molecular orbitals of the metal surface. This makes it possible to describe the electronic structure of a metal in more detail. In part, the overlap of the acceptor orbital of reactant with the s-, p-, and d-orbitals of the metal is addressed exactly; the effect of corrugation and surface defects can also be considered in a proper way. Test calculations of the electronic transmission coefficients were performed for other gold clusters, two-layer Au20, Au30, and Au40 and three-layer Au66. It was found that the size of the Au54 cluster is enough to obtain reliable results. Note that the accuracy of our model decreases with decreasing the cluster−reactant separation. This results from the failure of the perturbation theory at short interatomic distances where we can use the unperturbed electronic wave functions in initial and final states no longer.

It is evident from eq C4 that some difference between both reorganization energies can originate purely from a difference in the spring constants ki and kf. It was already mentioned that the nearest coordination numbers of anions and cations in the Fc(Fc+) solvation shells do not differ significantly for Fc and Fc+. It is easy to see, on the other hand, that g(r) functions in IL describing the coordination shells of the reactant and product (Figure 4) demonstrate a more structured (ordered) character of the ferrocene solvent shell as compared with ferrocenium. In AN the solvation shell of Fc+ is also less structured than that of the neutral molecule, as can be seen from the shapes and positions of minima and maxima in g(r) functions. Such a feature gives some evidence in favor of nonequality of ki and kf (or more exactly to say, ki > kf). Just this observation is the most important for a qualitative explanation of the origin of the reorganization asymmetry.



APPENDIX D. DEFINITION OF THE SOLVENT CORRELATION FUNCTION The solvent correlation function M(τ) describes the fluctuating reactant energy level in solution and can be defined as follows73



APPENDIX C. CONSTRUCTING OF THE REACTION FREE ENERGY SURFACES AND QUALITATIVE EXPLANATION OF THE REORGANIZATION ASYMMETRY Free energy surfaces of an ET reaction are calculated as a function of the generalized solvent coordinate X, which is defined as the difference in the potential energy between the oxidized and reduced states of the solute for the same configuration of the solvent X = Vo − Vr

M (τ ) =

Q (τ ) =

M (τ ) =



kf ki (xf − xi)2 and λf = (xf − xi)2 2 2





0

(D2)

⟨δ ΔE(τ ), δ ΔE(0)⟩ ⟨(δ ΔE(τ ))2 ⟩

(D3)

ASSOCIATED CONTENT

S Supporting Information *

Snapshot of the simulation box, visualization of solvation shells of Fc and Fc+ in IL, g(r) in AN, probabilities of the solvent order parameter fluctuations, and molecular orbitals for Fc+. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

(C3)

Notes

The authors declare no competing financial interest.

where ki and kf are the spring constants of oscillators, and xi and xf are the their equilibrium positions. Then the solvent reorganization energies for reduction (λi) and oxidation (λf) are defined as follows λi =



∫−∞ exp(−iωτ)⎢⎣ ε(1ω) − ε1 ⎥⎦ dωω

where δΔE(τ) = ΔE(τ) − ΔEeq(τ) and ΔE(τ) is the Coulomb part of the nonequilibrium solute−solvent interaction energy, and ΔEeq(τ) is the Coulomb energy of the equilibrium solute− solvent interaction.

(C2)

where kB is the Boltzmann constant; T is temperature; and G0 is the difference in free energy between the two minima. In the case of electrochemical reactions G0 is equal to the electrode overvoltage, η. When calculating the free energy surfaces we restricted them to the case of zero η because there are no experimental data related to various η values for comparison (at least in the absence of artificial barrier layers). A simple analysis is suggested below to explain qualitatively the origin of the strong reorganization asymmetry we observed. Let us represent the potential energy surfaces along the polarization (solvent) coordinate x in terms of effective oscillators kf ki (x − xi)2 and Uf (x) = (x − xf )2 2 2

i 2π

with ε0 being the static (low-frequency) dielectric constant (the dependence 1/ε(ω) − 1/ε0 plays the role of a generalized Pekar factor). The function M(τ) can be calculated directly from MD simulations74

(C1)

Go(x) = −kBT ln(Po(X ))

Ui(x) =

(D1)

where Q(τ) is a function of the dielectric spectrum ε(ω)

During the MD simulation the value of X fluctuates. So, it is possible to calculate the probability distributions for oxidized and reduced states Po(X) and Pr(X) of finding various values of X(t). Then, the free energy surfaces can be calculated according to the equations

Gr(x) = −kBT ln(Pr(X )) + G0

Q (τ ) Q (0)



ACKNOWLEDGMENTS This work was supported in part by the RFBR (projects No. 14-03-00935a, 12-03-31748). The computational resources of the Joint Supercomputer Center of RAS and Supercomputing Center of Lomonosov Moscow State University77 are highly

(C4) 6161

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

Electron Transfer Reactions: High Precision Steady-State Measurements of the Standard Electrochemical Rate Constant for Ferrocene Derivatives in Alkyl Cyanide Solvents. J. Electroanal. Chem. 2005, 580, 78−86. (19) Dolidze, T. D.; Khoshtariya, D. E.; Illner, P.; Eldik, R. v. Heterogeneous Electron Transfer at Au/SAM Junctions in a RoomTemperature Ionic Liquid Under Pressure. Chem. Commun. 2008, 2112−2114. (20) Khoshtariya, D. E.; Dolidze, T. D.; Eldik, R. v. Multiple Mechanisms for Electron Transfer at Metal/Self-Assembled Monolayer/Room-Temperature Ionic Liquid Junctions: Dynamical Arrest versus Frictional Control and Non-Adiabaticity. Chem.Eur. J. 2009, 15, 5254−5262. (21) Stoppa, A.; Hunger, J.; Buchner, R.; Hefter, G.; Thoman, A.; Helm, H. Interactions and Dynamics in Ionic Liquids. J. Phys. Chem. B 2008, 112, 4854−4858. (22) Liu, Z.; Huang, S.; Wang, W. A Refined Force Field for Molecular Simulation of Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2004, 108, 12978−12989. (23) Kuznetsov, A. M.; Ulstrup, J. Electron Transfer in Chemistry and Biology; Wiley: Chichester, U.K., 1999. (24) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision A.1; Gaussian, Inc.: Wallingford CT, 2009. (25) Hay, P. J.; Wadt, W. R. Ab Initio Effective Core Potentials for Molecular Calculations - Potentials for K to Au Including the Outermost Core Orbitals. J. Chem. Phys. 1985, 82, 299−310. (26) Cancès, E.; Mennucci, B.; Tomasi, J. A New Integral Equation Formalism for the Polarizable Continuum Model: Theoretical Background and Applications to Isotropic and Anistropic Dielectrics. J. Chem. Phys. 1997, 107, 3032−3041. (27) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Ab Initio Study of Solvated Molecules: A New Implementation of the Polarizable Continuum Model. Chem. Phys. Lett. 1996, 255, 327−335. (28) Cossi, M.; Scalmani, G.; Rega, N.; Barone, V. New Developments in the Polarizable Continuum Model for Quantum Mechanical and Classical Calculations on Molecules in Solution. J. Chem. Phys. 2002, 117, 43−54. (29) Smith, W.; Forester, T. R.; Todorov, I. T. The DL POLY Classic User Manual, D. L., UK. The DLPOLY website is http://www.ccp5.ac. uk/DL_POLY_CLASSIC/. (30) Lopes, J. N. C.; Couto, P. C. d.; Piedade, M. E. M. d. An AllAtom Force Field for Metallocenes. J. Phys. Chem. A 2006, 110, 13850−13856. (31) Breneman, C. M.; Wiberg, K. B. Determining Atom-Centered Monopoles from Molecular Electrostatic Potentials. The Need for High Sampling Density in Formamide Conformational Analysis. J. Comput. Chem. 1990, 11, 361−373. (32) Nikitin, A. M.; Lyubartsev, A. P. New Six-Site Acetonitrile Model for Simulations of Liquid Acetonitrile and its Aqueous Mixtures. J. Comput. Chem. 2007, 28, 2020−2026. (33) Tanaka, S.; Sengoku, Y. Nuclear Quantum Effects on Electron Transfer Reactions in DNA Hairpins. Phys. Rev. E 2003, 68, 031905− 031910. (34) Warshel, A. Dynamics of Reactions in Polar Solvents. Semiclassical Trajectory Studies of Electron-Transfer and ProtonTransfer Reactions. J. Phys. Chem. A 1982, 86, 2218−2224. (35) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (36) Nose, S. A Unified Formulation of the Constant Temperature Molecular-Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (37) Iori, F.; Corni, S. Including Image Charge Effects in the Molecular Dynamics Simulations of Molecules on Metal Surfaces. J. Comput. Chem. 2008, 29, 1656−1666. (38) Iori, F.; Felice, R. D.; Molinari, E.; Corni, S. GolP: An Atomistic Force-Field to Describe the Interaction of Proteins with Au(111) Surfaces in Water. J. Comput. Chem. 2009, 30, 1465−1476.

appreciated. We also thank Tamara Zinkicheva for the help with calculations.



REFERENCES

(1) Barnes, A. S.; Rogers, E. I.; Streeter, I.; Aldous, L.; Hardacre, C.; Compton, R. G. Extraction of Electrode Kinetic Parameters from Microdisc Voltammetric Data Measured under Transport Conditions Intermediate between Steady-State Convergent and Transient Linear Diffusion As Typically Applies to Room Temperature Ionic Liquids. J. Phys. Chem. B 2008, 112, 7560−7565. (2) Dolidze, T. D.; Khoshtariya, D. E.; Illner, P.; Kulisiewicz, L.; Delgado, A.; van Eldik, R. High-Pressure Testing of Heterogeneous Charge Transfer in a Room-Temperature Ionic Liquid: Evidence for Solvent Dynamic Control. J. Phys. Chem. B 2008, 112, 3085−3100. (3) Fietkau, N.; Clegg, A. D.; Evans, R. G.; Villagron, C.; Hardacre, C.; Compton, R. G. Electrochemical Rate Constants in Room Temperature Ionic Liquids: The Oxidation of a Series of Ferrocene Derivatives. ChemPhysChem 2006, 7, 1041−1045. (4) Fontaine, O.; Lagrost, C.; Ghilane, J.; Martin, P.; Trippé, G.; Fave, C.; Lacroix, J.-C.; Hapiot, P.; Randriamahazaka, H. N. Mass Transport and Heterogeneous Electron Transfer of a Ferrocene Derivative in a Room-Temperature Ionic Liquid. J. Electroanal. Chem. 2009, 632, 88−96. (5) Lovelock, K. R. J.; Cowling, F. N.; Taylor, A. W.; Licence, P.; Walsh, D. A. Effect of Viscosity on Steady-State Voltammetry and Scanning Electrochemical Microscopy in Room Temperature Ionic Liquids. J. Phys. Chem. B 2010, 114, 4442−4450. (6) Siraj, N.; Grampp, G.; Landgraf, S.; Punyain, K. Cyclic Voltammetric Study of Heterogeneous Electron Transfer Rate Constants of Various Organic Compounds in Ionic Liquids: Measurements at Room Temperature. Z. Phys. Chem. 2013, 227, 105−119. (7) Zhang, J.; Bond, A. Conditions Required To Achieve the Apparent Equivalence of Adhered Solid- and Solution-Phase Voltammetry for Ferrocene and Other Redox-Active Solids in Ionic Liquids. Anal. Chem. 2003, 75, 2694−2702. (8) Fawcett, W. R.; Gaál, A.; Misicak, D. Estimation of the Rate Constant for Electron Transfer in Room Temperature Ionic Liquids. J. Electroanal. Chem. 2011, 660, 230−233. (9) Lynden-Bell, R. M. Can Marcus Theory Be Applied to Redox Processes in Ionic Liquids? A Comparative Simulation Study of Dimethylimidazolium Liquids and Acetonitrile. J. Phys. Chem. B 2007, 111, 10800−10806. (10) Roy, D.; Maroncelli, M. Simulation of Solvation and Solvation Dynamics in an Idealized Ionic Liquid Model. J. Phys. Chem. B 2012, 116, 5951−5970. (11) Shim, Y.; Kim, H. J. Free Energy and Dynamics of ElectronTransfer Reactions in a Room Temperature Ionic Liquid. J. Phys. Chem. B 2007, 111, 4510−4519. (12) Shim, Y.; Kim, H. J. Adiabatic Electron Transfer in a RoomTemperature Ionic Liquid: Reaction Dynamics and Kinetics. J. Phys. Chem. B 2009, 113, 12964−12972. (13) Ando, K. Solvent Nuclear Quantum Effects in Electron Transfer Reactions. III. Metal Ions in Water. Solute Size and Ligand Effects. J. Chem. Phys. 2001, 114, 9470−9477. (14) King, G.; Warshel, A. Investigation of the Free Energy Functions for Electron Transfer Reactions. J. Chem. Phys. 1990, 93, 8682−8692. (15) Vladimirov, E.; Ivanova, A.; Rösch, N. Effect of Solvent Polarization on the Reorganization Energy of Electron Transfer from Molecular Dynamics Simulations. J. Chem. Phys. 2008, 129, 194515− 194524. (16) Kobrak, M. N.; Li, H. Electrostatic Interactions in Ionic Liquids: the Dangers of Dipole and Dielectric Descriptions. Phys. Chem. Chem. Phys. 2010, 12, 1922−1932. (17) Xiao, T.; Song, X. Reorganization Energy of Electron Transfer Processes in Ionic Fluids: A Molecular Debye-Hückel Approach. J. Chem. Phys. 2013, 138, 114105−114117. (18) Clegg, A. D.; Rees, N. V.; Klymenko, O. V.; Coles, B. A.; Compton, R. G. Marcus Theory of Outer-Sphere Heterogeneous 6162

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

(39) Zusman, L. D. Calculating of the Rate Constant for an Elementary Act of Electron Transfer with an Arbitrary Dielectric-Loss Spectrum for the Solvent. Theor. Exp. Chem. 1984, 19, 381−386. (40) Zusman, L. D. The Theory of Electron Transfer Reactions in Solvents with Two Characteristic Relaxation Times. Chem. Phys. 1988, 119, 51−61. (41) Zusman, L. D. Electron Transfer Reactions at an Electrode in Solvents with Two Characteristic Relaxation Times. Sov. Electrochem. 1989, 24, 1125−1132. (42) Baranski, A. S.; Winkler, K.; Fawcett, W. R. New Experimental Evidence Concerning the Magnitude of the Activation Parameters for Fast Heterogeneous Electron Transfer Reactions. J. Electroanal. Chem. 1991, 313, 367−315. (43) Asaki, M. L. T.; Redondo, A.; Zawodzinski, T. A.; Taylor, A. J. Dielectric Relaxation and Underlying Dynamics of Acetonitrile and 1Ethyl-3-Methylimidazolium Triflate Mixtures Using THz Transmission Spectroscopy. J. Chem. Phys. 2002, 116, 10377−10385. (44) Kumar, A. Estimates of Internal Pressure and Molar Refraction of Imidazolium Based Ionic Liquids as a Function of Temperature. J. Solution Chem. 2008, 37, 203−214. (45) Lide, D. R. CRC Handbook of Chemistry and Physics, 86th ed.; CRC Press: Boca Raton, 2005. (46) Fedorov, M. V.; Lynden-Bell, R. M. Probing the Neutral Graphene-Ionic Liquid Interface: Insights from Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2012, 14, 2552−2556. (47) Feng, G.; Zhang, J. S.; Qiao, R. Microstructure and Capacitance of the Electrical Double Layers at the Interface of Ionic Liquids and Planar Electrodes. J. Phys. Chem. C 2009, 113, 4549−4559. (48) Kislenko, S. A.; Samoylov, I. S.; Amirov, R. H. Molecular Dynamics Simulation of the Electrochemical Interface Between a Graphite Surface and the Ionic Liquid [BMIM][PF6]. Phys. Chem. Chem. Phys. 2009, 11, 5584−5590. (49) Lynden-Bell, R. M.; Frolov, A. I.; Fedorov, M. V. Electrode Screening by Ionic Liquids. Phys. Chem. Chem. Phys. 2012, 14, 2693− 2701. (50) Vatamanu, J.; Borodin, O.; Smith, G. D. Molecular Insights into the Potential and Temperature Dependences of the Differential Capacitance of a Room-Temperature Ionic Liquid at Graphite Electrodes. J. Am. Chem. Soc. 2010, 132, 14825−14833. (51) Vatamanu, J.; Borodin, O.; Smith, G. D. Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-N-propylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes. J. Phys. Chem. B 2011, 115, 3073−3084. (52) Atkin, R.; Borisenko, N.; Druschler, M.; Abedin, S. Z. E.; Endres, F.; Hayes, R.; Huber, B.; Roling, B. An In Situ STM/AFM and Impedance Spectroscopy Study of the Extremely Pure 1-Butyl-1methylpyrrolidinium Tris(pentafluoroethyl)trifluorophosphate/ Au(111) Interface: Potential Dependent Solvation Layers and the Herringbone Reconstruction. Phys. Chem. Chem. Phys. 2011, 13, 6849−6857. (53) Hayes, R.; Borisenko, N.; Tam, M. K.; Howlett, P. C.; Endres, F.; Atkin, R. Surface Structure at the Ionic Liquid-Electrified Metal Interface. Acc. Chem. Res. 2008, 41, 421−431. (54) Hayes, R.; Borisenko, N.; Tam, M. K.; Howlett, P. C.; Endres, F.; Atkin, R. Double Layer Structure of Ionic Liquids at the Au(111) Electrode Interface: An Atomic Force Microscopy Investigation. J. Phys. Chem. C 2011, 115, 6855−6863. (55) Su, Y.; Yan, J.; Li, M.; Zhang, M.; Mao, B. Electric Double Layer of Au(100)/Imidazolium-Based Ionic Liquids Interface: Effect of Cation Size. J. Phys. Chem. C 2013, 117, 205−212. (56) Su, Y.-Z.; Fu, Y.-C.; Yan, J.-W.; Chen, Z.-B.; Mao, B.-W. Double Layer of Au(100)/Ionic Liquid Interface and Its Stability in Imidazolium-Based Ionic Liquids. Angew. Chem., Int. Ed. 2009, 48, 5148−5151. (57) Zhang, X.; Zhong, Y.-X.; Yan, J.-W.; Su, Y.-Z.; Zhang, M.; Mao, B.-W. Probing Double Layer Structures of Au (111)-BMIPF6 Ionic Liquid Interfaces from Potential-Dependent AFM Force Curves. Chem. Commun. 2012, 48, 582−584.

(58) Alam, M. T.; Masud, J.; Islam, M. M.; Okajima, T.; Ohsaka, T. Differential Capacitance at Au(111) in 1-Alkyl-3-methylimidazolium Tetrafluoroborate Based Room-Temperature Ionic Liquids. J. Phys. Chem. C 2011, 115, 19797−19804. (59) Gnahm, M.; Pajkossy, T.; Kolb, D. M. The Interface Between Au(111) and an Ionic Liquid. Electrochim. Acta 2010, 55, 6212−6217. (60) Lin, L. G.; Wang, Y.; Yan, J. W.; Yuan, Y. Z.; Xiang, J.; Mao, B. W. An In Situ STM Study on the Long-Range Surface Restructuring of Au(111) in a Non-Chloroaluminumated Ionic Liquid. Electrochem. Commun. 2003, 5, 995−999. (61) Panzram, E.; Baumgartel, H.; Roelfs, B.; Schroter, C.; Solomun, T. A Capacitance and Infrared Study of the Electrical Double Layer Structure at Single Crystal Gold Electrodes in Acetonitrile. Ber. Bunsenges. Phys. Chem. 1995, 99, 827−837. (62) Pavlishchuk, V. V.; Addison, A. W. Conversion Constants for Redox Potentials Measured versus Different Reference Electrodes in Acetonitrile Solutions at 25 °C. Inorg. Chim. Acta 2000, 298, 97−102. (63) The Potential of Zero Charge; Trasatti, S., Lust, E., Eds.; Kluwer Academic/Plenum Publishers: New York, 1999; Vol. 33. (64) Davidson, T.; Pons, B. S.; Bewick, A.; Schmidt, P. P. Vibrational Spectroscopy of the Electrode/Electrolyte Interface. Use of Fourier Transform Infrared Spectroscopy. J. Electroanal. Chem. 1981, 125, 237−241. (65) Pons, S.; Davidson, T.; Bewick, A. Vibrational Spectroscopy of the Electrode-Solution Interface: Part III. Use of Fourier Transform Spectroscopy for Observing Double Layer Reorganization. J. Electroanal. Chem. 1982, 140, 211−216. (66) Martin, D. R.; Matyushov, D. V. Non-gaussian Statistics and Nanosecond Dynamics of Electrostatic Fluctuations Affecting Optical Transitions in Proteins. J. Phys. Chem. B 2012, 116, 10294−10300. (67) Hartnig, C.; Koper, M. T. M. Molecular Dynamics Simulation of the First Electron Transfer Step in the Oxygen Reduction Reaction. J. Electroanal. Chem. 2002, 532, 165−170. (68) Small, D. W.; Matyushov, D. V.; Voth, G. A. The Theory of Electron Transfer Reactions: What May be Missing? J. Am. Chem. Soc. 2003, 125, 7470−7478. (69) Kuharskii, R. A.; Bader, J. S.; Chandler, D.; Sprick, M.; Klein, M. L.; Impey, R. W. Molecular Model for Aqueous Ferrous-Ferric Electron Transfer. J. Chem. Phys. 1988, 89, 3248−3257. (70) Hartnig, C.; Koper, M. T. M. Molecular Dynamics Simulations of Solvent Reorganization in Electron-Transfer Reactions. J. Chem. Phys. 2001, 115, 8540−8546. (71) Nikitina, V. A.; Nazmutdinov, R. R.; Tsirlina, G. A. Quinones Electrochemistry in Room-Temperature Ionic Liquids. J. Phys. Chem. B 2011, 115, 668−677. (72) Belosov, A. A.; Kuznetsov, A. M. Simple Approximative Method in the Kramers Theory of Chemical Reactions. Sov. Electrochem. 1988, 23, 1564−1567. (73) Sparpaglione, M.; Mukamel, S. Dielectric Friction and the Transition from Adiabatic to Nonadiabatic Electron Transfer in Condensed Phases. II. Application to Non-Debye Solvents. J. Chem. Phys. 1988, 88, 4300−4311. (74) Bhargava, B. L.; Balasubramanian, S. Dynamics in a RoomTemperature Ionic Liquid: A Computer Simulation Study of 1, 3Dimethylimidazolium Chloride. J. Chem. Phys. 2005, 123, 144505− 144512. (75) Nicholson, R. S. Theory and Application of Cyclic Voltammetry for Measurement of Electrode Reaction Kinetics. Anal. Chem. 1965, 37, 1351−1355. (76) Zusman, L. D. Outer-Sphere Electron Transfer in Polar Solvents. Chem. Phys. 1980, 49, 295−304. (77) Sadovnichy, V.; Tikhonravov, A.; Voevodin, V.; Opanasenko, V. ″Lomonosov″: Supercomputing at Moscow State University. In Contemporary High Performance Computing: From Petascale toward Exascale; CRC Press: Boca Raton, USA, 2013. (78) Tolokh, I. S.; Vivcharuk, V.; Tomberli, B.; Gray, C. G. Binding Free Energy and Counterion Release for Adsorption of the Antimicrobial Peptide Lactoferricin B on a POPG Membrane. Phys. Rev. E 2009, 80, 031911−031923. 6163

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164

The Journal of Physical Chemistry C

Article

(79) Vivcharuk, V.; Tomberli, B.; Tolokh, I. S.; Gray, C. G. Prediction of Binding Free Energy for Adsorption of Antimicrobial Peptide Lactoferricin B on a POPC Membrane. Phys. Rev. E 2008, 77, 031913−031924. (80) Gosavi, S.; Marcus, R. A. Nonadiabatic Electron Transfer at Metal Surfaces. J. Phys. Chem. B 2000, 104, 2067−2072. (81) Gosavi, S.; Marcus, R. A. Temperature Dependence of the Electronic Factor in the Nonadiabatic Electron Transfer at Metal and Semiconductor Electrode. J. Electroanal. Chem. 2001, 500, 71−77. (82) Nazmutdinov, R. R.; Berezin, S. A.; Soldano, G.; Schmickler, W. Orbital Overlap Effects in Electron Transfer Reactions across a Metal Nanowire/Electrolyte Solution Interface. J. Phys. Chem. C 2013, 117, 13021−13027.

6164

dx.doi.org/10.1021/jp4072108 | J. Phys. Chem. C 2014, 118, 6151−6164