A Computationally Efficient Treatment of Polarizable Electrochemical

Article pubs.acs.org/JPCC

A Computationally Efficient Treatment of Polarizable Electrochemical Cells Held at a Constant Potential Matt K. Petersen,† Revati Kumar,† Henry S. White,‡ and Gregory A. Voth*,† †

Department of Chemistry, James Franck Institute and Computation Institute, University of Chicago, 5735 South Ellis Avenue, Chicago, Illinois 60637, United States ‡ Department of Chemistry, University of Utah, 315 South 1400 East, Salt Lake City, Utah 84112, United States ABSTRACT: A computationally efficient method is presented for the treatment of electrostatic interactions between polarizable metallic electrodes held at a constant potential and separated by an electrolyte. The method combines a fluctuating uniform electrode charge with explicit image charges to account for the polarization of the electrode by the electrolyte, and a constant uniform charge added to the fluctuating uniform electrode charge to account for the constant potential condition. The method is then used to calculate electron transport rates using Marcus theory and these rates are incorporated in a Kinetic Monte Carlo - Molecular Dynamics simulation scheme to efficiently model oxidation/reduction reactions in an electrochemical cell.

1. INTRODUCTION Batteries with novel three-dimensional (3D) architectures have been proposed as potential devices with high energy and power density.1 A key feature of these proposed batteries is short interelectrode distances on the order of a few nanometers, resulting in the overlap of the anode and cathode double layers. In order to help understand the processes that dominate ion transport in these nanoscale electrochemical cells, theoretical modeling of these systems is valuable. Although Brownian dynamics simulations have previously been performed,2 atomistic simulation studies are necessary to form a more complete picture of the underlying mechanisms. The computational framework needed to study these systems at an atomistic level is developed in this paper. Ab initio molecular dynamics3 simulations have been used to study electrochemical reactions,4−6 as well as the solution− electrode interface.7,8 However, the computational challenges associated with this method in treating a full electrochemical cell, in particular the large system sizes and time scales needed to simulate an electrode−electrolyte interface, have motivated the use of model potentials for the description of such systems in more empirical classical molecular dynamics (MD) simulations. Several simplifying approximations are made when treating the electrode/electrolyte interface within classical MD simulations. The most straightforward approach, in concept and implementation, has been to model the electrode as a uniformly charged plane.9 More advanced treatments have further included a treatment of the atomically corrugated electrode surface,10,11 or modeled the electrode surface as a collection of atomic sites bearing partial charges contributing to the total electrode charge.12−14 While computationally inexpensive and conceptually simple, these approaches fail to capture the dispersion interactions between the electrolyte and the © 2012 American Chemical Society

electrode, with a potentially significant impact on the structure of the electrical double layer at the electrode.15 One method, the so-called “method of images”, does account for the polarization of the electrode by the electrolyte.11,15,16 In this method, a static image plane is defined  typically coincident with the nuclear plane of the electrode surface. For each charge in the system qi there is a corresponding image charge such that qimg = −qi. In addition to the interactions between all the real charges (i.e., nonimage charges) in the cell, interactions between these real charges and the image charges are included as well, whereas image charge−image charge interactions are excluded. While the method does correctly account for the polarization of the electrode by the electrolyte, the method can only model a charged electrode when the electrolyte has a net charge.17,18 It should be noted that creating a fixed charge on the electrode, either by applying a uniform charge to the electrode or through the image charges of an electrolyte with a net charge, cannot correctly account for the electrode charge due to an imposed constant potential. As will be shown below, the correct charge distribution can be correctly modeled by a combination of the two methods given specific simulation constraints. An approach capable of modeling a polarizable charged electrode adjacent to a non-neutral electrolyte has been formulated by Siepmann and Sprik,19 and applied to bulk simulations of a complete electrochemical cell by Madden and co-workers.20,21 In short, this method models the electrode as a series of Gaussian charge distributions centered about the nuclear positions of the electrode atoms. The charge on each electrode atom at each simulation time-step is obtained by enforcing a constant potential at the electrode atoms. While this method accurately Received: October 25, 2011 Revised: January 12, 2012 Published: January 12, 2012 4903

dx.doi.org/10.1021/jp210252g | J. Phys. Chem. C 2012, 116, 4903−4912

The Journal of Physical Chemistry C

Article

constraint and through the treatment of the long-range electrostatic interactions. From Gauss’ law we know that the electric field integrated over a closed surface is equal to 1/ε0 times the charge enclosed by that surface,23 that is,

reproduces the charge distribution within the electrode, it is computationally expensive. Here, an alternative method is provided that quantitatively reproduces the results of Reed et al.20,22 without the expense of the iterative minimization of the system potential or by otherwise using polarizable atomic potentials. The paper is divided into five sections. Section 2 describes the theoretical model of the electrochemical cell as well as the description of the reaction rates for the oxidation/reduction reactions. The applications of the method to equilibrium MD simulations and nonequilibrium simulations involving redox reactions are presented in section 3. The details of the computations, in particular, the incorporation of the method in a standard MD software package, are discussed in section 4, while the conclusions presented in section 5.

A ε0

∮S E⃗·da ⃗ =

∫ ρ(z) dz

(1)

where E⃗ is the electric field at the surface, A is the area of the electrodes, and ρ(z) is the charge density at z. For a cell composed of discrete charges, q (z*i ) ρ(z*i ) = i AΔz i

(2)

where q(zi*) is the discrete charge at the point zi* on the interval Δzi. When the electric field is integrated over the entire cell, the electric field just inside the second electrode is

2. METHODS 2.1. The Model. Consider an electrochemical cell composed of two parallel electrodes and intervening electrolyte (Figure 1). This cell consists of two electrodes: the left electrode

E=

1 ε0

∫0

D

1 lim ε0A Δz i → 0

ρ(z) dz =

p ⎛ q (z * ) ⎞ i ⎜ i ⎟

∑⎜

⎟Δz i Δz i ⎠ ⎝ i=1 (3)

such that, p

∑ Δzi = D (4)

i=1

for p partitions. Expanding the sum gives 1 (q + qLE + E= A ε0 RE = Figure 1. A schematic diagram of the electrochemical cell. The inner two vertical lines represent the left and right electrodes. The point charge qR within the cell has two images, one corresponding to the left electrode LI and the other to the right electrode RI. The charges qLE and qRE represent the higher order images charges and are distributed evenly over the inner surfaces of the left and right electrodes, respectively.

1 (q + qLE − A ε0 RE

n

n

∑ q R i + 2 ∑ qI i )

i=1 n

i=1

∑ q R i)

(5)

i=1

where qLE and qRE are the higher order image charges induced on the surface of the left and right electrodes by the n charges of the electrolyte, excluding primary images charges. Since the electric field inside a conductor is 0, we have n

∑ qR i = qRE + qLE

(6)

i

(LE) located at z = 0, and the right electrode (RE) located at z = D, both periodic in the x and y dimensions. A charged atom of the bulk electrolyte (R), located at the position (0,0,zR) with charge qR, will induce a charge on the surface of the conducting electrodes. This induced surface charge will create a potential equivalent to one where a charge qI = −qR is located equidistant from the electrode surface but within the electrode.23 That is, the real electrolyte charge R will induce a surface charge in the two electrodes creating a potential equal to that created by a pair of image charges, RI and LI. These image charges will have charge qI = −qR, located at (0,0,−zR) and (0,0,d1 + 2d2), such that d1 = z, d2 is the distance of the real charge from electrode RE, and D = d1 + d2 is the distance between electrodes LE and RE. The above discussion does not account for the possibility of the primary image charges of one electrode inducing secondary image charges in the opposing electrode, which in turn can create further image charges in the first, and so on. As will be shown in the following discussion, these higher order image charges will be accounted for as part of the constant potential

The electrical potential across the electrochemical cell can be found from the electric field,

Φ(z) = −

∫ E (z ) d z

(7)

Integrated over the entire cell, and using the expression for the electric field from eq 3, the potential is Φ(D) = −

=−

1 ε 0A

∫0

D

dz lim

1 lim ε0A Δzj → 0

p ⎛ qi(z*i ) ⎞ ⎟

∑ ⎜⎜

⎟Δz i

Δz i → 0 i = 1 ⎝ Δz i ⎠ p j

∑ ( ∑ qi(z*i ))Δzj j=1 i=1

(8)

Expanding the double sum gives Φ(D) = −

4904

1 ε 0A

n

∑ [ − qR iz R i + (qLEi)(D)]

i=1

(9)



Article

where qLEi is the surface charge induced by the higher order images resulting from the ith charge of the electrolyte on electrode LE such that, n

∑ qLEi

qLE =

(10)

i=1

Noting that the image charge method is derived for a conducting plane held at a potential, V0 = 0, 1 [ ε 0A

n

∑ i=1

− qR z R i + (qLE)(D)] = 0 i

(11)

so, n q zR Ri i

∑

qLE =

D

i=1

(12)

Applying a constant nonzero potential to the cell simply raises the electrode charge uniformly by an amount proportional to the potential. If one expresses that constant potential in terms of the electrode charge, Q0, of an equivalent capacitor, ΔV0 =

Q0 ε 0A

Figure 2. Force experienced by a point charge q1 from a periodically replicated point charge q2 (circles) and from a uniformly charged wall (squares) of equal charge. The force is normalized to q1q2/2Aε0, the force on a point charge q1 due to a uniformly charged plane of total charge q2. The length has been normalized to Lx, the length of the system in the x direction (Note that Lx = Ly).

D (13)

Lx = Ly, and nonperiodic in the z dimension. Replacing the charged wall with a point charge, q2, located at (0,0,0) which is replicated periodically in the x and y directions, the force on the test charge is now given by the circles in Figure 2. The length and force in Figure 2 have been normalized to Lx,y and q1q2/Aε02  the force due to a uniformly charged wall  respectively. At large enough separation, the array of periodic images of the higher order image charges should appear as a uniformly charged wall. From Figure 2 it is clear that this occurs at a length equivalent to one to two periodic images. Approximating the images as a uniformly charged plane has been used to treat primary images at separations where z > min(Lx,Ly).24 Reed et al. have shown that the higher order image charges are uniformly distributed over the entire electrode even at lengths roughly equal to the x and y periodic lengths.20 Since, at their closest, the separation between the real and higher order images charges is never greater than the electrode separation, D, approximating the induced surface charge, qEL and qER, as being uniformly distributed is a reasonable approximation. 2.2. Validation. Initial validation of the model consists of a comparison with the more sophisticated, albeit computationally more intensive, model of Reed et al.20,22 In short, the model of Reed et al. allows for the polarization of the metallic electrode through variable Gaussian charge distributions which are adjusted according to a variational procedure, accounting for the constant potential condition of the electrodes. The key differences between the model of Reed et al. and the one proposed here is the way the electrode charges are adjusted and the use of explicit image charges in the present model. Where Reed et al. meet the constant potential condition by minimizing the total potential of the system through adjusting the atomic charge of the electrode, the present model adjusts the electrode atomic charges according to eqs 14a and 14b. While the model of Reed et al. intrinsically arrives at the electrode surface charge distribution commensurate with that of an image charge, the present model uses explicit image charges, updated according to the position of the corresponding electrolyte charges, to

Using the equality from eq 6, and accounting for the charge due to the constant applied potential, ΔV0, n q zR Ri i q′LE = ∑ + Q0 D (14a) i=1 n

q′RE =

⎛

∑ qR i⎜⎝1 −

i=1

z Ri ⎞ ⎟ − Q0 D ⎠

(14b)

Equations 14a and 14b are the primary result of the analysis. To summarize, assuming that there are two parallel electrodes, an intervening electrolyte of discrete charges, and that primary image charges are induced on each electrode, then the total charge of the higher order images charges are equal to the total of the electrolyte charge. This is seemingly a straightforward conclusion. If it were not so, then the cell would have created charge. More importantly though, eqs 14a and 14b show that the total electrolyte charge is shared between the electrodes in proportion to the position of the center of electrolyte charge. What eqs 14a and 14b show is that the charge is transferred to both electrodes, and the amount of charge transferred to each is proportional to the position of the transferring species within the cell. This position dependence has been demonstrated by Reed et al.,22 albeit as an empirical observation of their model. Here is it shown that this dependence is an a priori consequence of the constant potential boundary condition. Unfortunately, eqs 14a and 14b do not show how this charge is distributed on the surface of the electrode. Considering the charge induced on the surface of the electrodes by the electrolyte was initially assumed to be correctly described by the primary image charges, how then should the charge induced by the opposing electrodes themselves be modeled? Figure 2 shows the force (squares) experienced by a test charge, q1 = e−, located at (0,0,z) due to a charged wall at (x,y,0), composed of a 10 by 10 array of atoms arranged in a simple cubic lattice with a net charge of q2 = −q1, as a function of z. The system is periodic in the x and y dimensions such that 4905



Article

mimic the nonuniform portion of the electrode surface charge distribution. The details of the image charge scheme, treatment of long-range electrostatics, and other issues are discussed in detail in section 4. Two validation systems have been created to investigate the accuracy of the proposed model: System I is composed of two parallel electrodes separated by a distance D = 52.9 Å. The electrodes are composed of three layers of atoms (a total of 2700 atoms) arranged in an FCC lattice with a lattice parameter of 3.92 Å. Each electrode is periodic in the x and y dimensions such that x = 83.1 Å and y = 72.0 Å, with the (111) face exposed to the electrolyte. A single charge, qR = −e is located at (0,0,z) in the space between the two electrodes. System II is composed of two parallel electrodes separated by a distance D = 85.0 Å. The electrodes are composed of three layers of atoms (a total of 294 atoms) arranged in an FCC lattice with a lattice parameter of 3.66 Å. Each electrode is periodic in the x and y dimensions such that x = 25.62 Å and y = 25.62 Å, with the (100) face exposed to the electrolyte. A single charge, qR = −e is located at (0,0,z) in the space between the two electrodes. Figure 3 shows the force experienced by the test charge in System I as it is moved from one electrode to the other for

Figure 4. Energy of System I (refer to section 2.2 in the text) as a function of the test charge position within the cell calculated using the present method (lines) and the method of Reed et al.20 (symbols). The total energy (solid line and diamonds) is decomposed into two terms  the energy due to V0 (dotted line and circles) and the energy from all images (dashed line and squares).

the portion of the energy due to the constant potential drop across the cell (eq 11). Note the distance dependence, that is, the constant slope across the cell. This is the reason the shape of the V0 = 0.0V and V0 = 2.75V force curves in Figure 3 are identical, albeit shifted by a constant force  that force simply being the slope of the V0 line in Figure 4. The dashed line of Figure 4 is the component of the energy due to the interaction of the test charge and all images, including the higher order images induced by the primary images of the opposing electrodes. The solid line is the total energy of the system. The symbols are the result for the method of Reed et al.20 It can be seen that the present method is in excellent quantitative agreement. Most noteworthy is the agreement for the energy due to all images, primary and higher. This is significant because this includes the energy from the higher order image charges that are calculated from eqs 14a and 14b, and further assumed to be uniformly distributed across the electrode surface (as discussed in relation to Figure 2). As discussed above, the net surface charge induced on an electrode from these higher order image charges is a function of the position of the electrolyte charges in the cell. Figure 5 shows the way in which the surface charge due to the higher order image charges varies as the test charge is moved across the cell. By way of an example of the consequences of the variable charge scheme, consider this cell when the test charge is located at z = (1/4)D. From eqs 14a and 14b one sees that qEL = (1/4)qR and qER = (3/4)qR. The total charge on the LE, including the primary image charge (qI = −qR), is then −(3/4) qR, and the total charge on the RE is −(1/4)qR. Note that the total charge of the cell remains neutral for any value of qR. With the details of the present method in place, the focus now moves to extracting electron transfer rate information from the model using the framework of Marcus theory.25,26 2.3. Electron Transfer. The usual Kinetic Monte Carlo (KMC)27,28 method gives the dynamical evolution of the system by utilizing the known rates of processes that can take place. These rates determine the probabilities of events and are used in a Monte Carlo scheme to propagate the system in time. In this work, the simulations of the electron transfer reactions involve regular MD simulations with the oxidation/reduction

Figure 3. Force experienced by the test charge in System I (refer to section 2.2 in the text) as a function of the position in the cell. The present method (lines) and the method of Reed et al.20 (symbols) were calculated at V0 = 0.0V (solid line and circles) and V0 = 2.75V (dashed line and squares).

V0 = 0.0V (solid line) and V0 = 2.75V (dashed line). The circles and squares are for the same system and conditions, respectively, for the method of Reed et al.20 The present method reproduces the force for the cell at the nonzero potential, while the force for the zero applied potential case deviates slightly near the electrode surface. The current approach uses point charges to represent the electrode surface charge, whereas the method of Reed et al. uses Gaussian charge distributions centered about the electrode atom centers. This may explain the difference at close separation. However, the force for the nonzero applied potential should have the same shape as that for the zero applied potential case, albeit shifted downward to more negative values of force. This is the case for the present method, but it does not appear to be so for the method of Reed et al. The reason for this is unclear. Figure 4 gives components of the total energy for System I as the test charge is moved across the cell. The dotted line shows 4906



Article

and A=

2π 2 1 |J | ℏ 4πλkBT

(17)

The reaction coordinate in this case is taken to be the mean vertical energy gap for oxidation. Assuming linear response of the solvent for the electron transfer reaction  that is, the curves of Figure 6 are quadratic  the mean values of the reaction coordinate for the oxidation and reduction reactions, ⟨ΔE1⟩ and ⟨ΔE2⟩, respectively, are related to the difference in the minima of the free energy for the reaction such that, ΔG0 = Figure 5. Electrode charges and interaction energy of System I (refer to section 2.2 in the text) due to the higher order image charges as a function of test charge position. The surface charge on LE (squares), located at z = 0, and charge on RE, located at z = D (circles), and the energy due to the interaction of the test charge with this electrode surface charge (diamonds) are shown.

λ=

(15)

(19)

ΔE1 = δERed → Ox + e− + I + W

(20a)

ΔE2 = − δEOx + e−→ Red + I + W

(20b)

Since the ionization energy and work functions are constants, ⟨ΔE1⟩ and ⟨ΔE2⟩ can be found by averaging δERed→Ox+e− and δEOx+e−→Red over many nuclear configurations of an MD trajectory. Reed et al. have shown22 that, on average,the distance dependence of the energy difference for the redox species is independent of the electrolyte and is simply the energy for creating the same charge in an otherwise empty cell. This is rationalized considering that during the electron transfer process the nuclear coordinates are not allowed to relax and therefore screen the newly created image charge. It is the image charge that is solely responsible for the distance dependence of the average energy difference. Figure 7 shows the average energy of the oxidation/reduction of a single ion, M+ ⇌ M2+ + e−, at V0 = 10.9V for System II described in the validation section. The lines are data calculated from the present method, while the symbols are those of Reed et al.22 Again, the two methods are in excellent agreement. The treatment so far has failed to account for the transfer of the electron from the oxidation/reduction reaction to or from a single electrode. The average energy gap (Figure 7) includes the contribution from the transfer of a partial charge (as described in relation to Figure 5) to both electrodes. While the correct electrode charge is maintained by the partial charge transfer, an adjustment needs to be made to the average energy gap before it can be used to calculate electron transfer rates.22

can be expressed25,26 in terms of the reorganization energy, λ, the free energy of reaction, −ΔG0, and the pre-exponential factor, A, such that, (λ + ΔG0)2 4λ

1 (⟨ΔE1⟩ − ⟨ΔE2⟩) 2

The energy difference between the initial and final states of the reaction can be found using the present method discussed in the preceding section. The energy difference for the oxidation reaction, δERed→Ox+e−, can be found by selecting some Red species in an MD trajectory and calculating the system energy according to the present method. The Red species is then changed to the Ox species (with all nuclear coordinates remaining fixed), and the system energy is again calculated according to the present method. The energy difference for the reduction reaction, ΔEOx+e−→Red is determined in a similar manner. The energy gap coordinate, ΔE, is simply this energy difference combined with the ionization energy for the ion and the work function of the metal electrode,

Figure 6. Representation of the Marcus diabatic energy surfaces for two oxidation states. The reorganization energy λ, the activation free energy (ΔG*), and the free energy difference between the minima of the two states (ΔG0) are shown.

ΔG* =

(18)

and the reorganization energy,

events treated using the KMC selection scheme; that is, when the electroactive species during the course of a regular MD simulation comes close to an electrode, the probability of a reaction taking place is calculated from the known reaction rate and if it is favorable the oxidation or reduction reaction will take place. The motivation for the following discussion in this section is to identify the key values from the preceding method which will allow for the calculation of electrode-reactant distance dependent reaction rates to be used in the KMC like portion of the simulations. Consider two energy surfaces representing the collective coordinates of solvent surrounding an ion with two possible oxidation states as in Figure 6. The rate of electron transfer, given by ket = Ae(−ΔG * /kBT )

1 (⟨ΔE1⟩ + ⟨ΔE2⟩) 2

(16) 4907



Article

Figure 7. ΔE for the oxidation reaction (solid line, circles) and the reduction (dashed line, squares) of the M2+/M+ redox couple in an electrochemical cell with an interelectrode distance of 85 Å and an applied potential of 10.9 V. The data for the symbols are from Reed et al.22

Figure 8. λ (open circles), ΔG0 corrected (open squares) and uncorrected (closed squares) for the electron transfer reaction at the anode for the M2+/M+ redox couple in an electrochemical cell with an interelectrode distance of 85 Å and an applied potential of 10.9 V.

Consider the energy change for the oxidation and reduction reactions at the anode from eqs 20a and 20b. The contribution of the partial charge transfer to the energy gap is simply,

can be made between the rate used by White and White and the Marcus rate discussed above, namely,

ΔV0 (δq+(z) + δq−(z)) 2

⎛ (λ + ΔG0)2 ⎞ ⎟⎟ kanode = A(z) exp⎜⎜ − 4λkBT ⎠ ⎝

(21)

For the system considered in Figure 8, the value of ΔG at zero applied voltage is close to 0. Willard and Chandler have pointed out that the change in ΔG0 with applied voltage is linear in the anode voltage.21 An applied potential of 10.9 V across the electrochemical cell, which gives rise to an anode potential Va of 5.45 V, results in a value of ΔG0 ≈ − 125 kcal/mol or −5.45 eV which is equal to −Va. If the applied voltage, and hence ΔG0, are much lower than the reorganization energy λ, then eq 25 reduces to

where δq+(z) =

(25) 0

qz D

(22) −

−

The charge q is equal to e for oxidation and −e for reduction. Hence the corrected energy gap at the anode is δEanode = δE − ΔV0q+(z)

(23)

⎛ λ − 2VaF ⎞ kanode = A(z) exp⎜ − ⎟ 4kBT ⎠ ⎝

Replacing δE in eqs 20a and 20b with δEanode, it is possible to find the corresponding free energy (eq 18) and reorganization energy (eq 19) for the oxidation reaction at the anode. One can repeat the same Marcus treatment, outlined above, for reduction reactions at the cathode with now the vertical energy gap for reduction as the reaction coordinate. In this case the identities of ΔE1 and ΔE2 are reversed, since the former is now the energy gap for reduction and the latter the energy gap for oxidation. Figure 8 shows the values of ΔG0, both uncorrected and corrected to account for the transfer of a full charge, as well as the reorganization energy, λ, for the oxidation reaction at the anode in System II, described in the validation section, at V0 = 10.9V  the same system and conditions presented in Figure 7. The only term left to evaluate is the pre-exponential term. 2.3.1. Pre-Exponential Term in the Marcus Rate. White and White have carried out Brownian dynamics simulations involving the redox reaction of a single ferrocene molecule between two electrodes at low overpotentials.2 Their rate equation for the oxidation reaction at the anode was given by ⎛ FV ⎞ kanode = k 0 exp( − βz) exp⎜ a ⎟ ⎝ 2kBT ⎠

(26)

We can rewrite the equation as follows ⎛ VF ⎞ kanode = A1(z) exp⎜ a ⎟ ⎝ 2kBT ⎠

(27)

where A1(z) = A(z) exp(−(λ)/(4kBT)). Comparing eq 24 with eq 27, one can arrive at a value of A(z). In the case where ΔG0 has a nonzero value, Ga, at 0 applied voltage, the value at an applied anode voltage Va will be Ga − VaF where F is Faraday’s constant. A1(z) will now be given by A(z) exp(−(λ)/(4kBT)) exp(−(Ga)/(4kBT)). For the reduction reaction at the cathode, where z is now the distance of the redox species from the cathode, a similar equation can be derived: ⎛ VF ⎞ kcathode = A1(z) exp⎜ − c ⎟ ⎝ 2kBT ⎠

(28)

where Vc is the cathode potential and A1(z) = A(z) exp(−(λ)/ (4kBT)) exp((Ga)/(4kBT)). For the system studied in Figure 8, Vc is equal to −5.45 eV. With the free energy, reorganization energy, and appropriately chosen pre-exponential factor (in this case taken from experiments on electron transfer in ferrocene)2,29 for the reaction, it is possible to calculate electron transfer rates (eq 25) for use in a KMC simulation.

(24)

where the constant k0 and β were determined from experiments carried out by Smalley and co-workers,29 Va is the potential applied to the anode and F is the Faraday constant. In order to determine the appropriate pre-exponential factor, a connection 4908



Article

such that Q0 = 0.0e−. At the time origin the constant electrode charge was then set such that Q0 = 0.6e−, and the simulation system was again allowed to equilibrate. The charge on the electrodes quickly reaches an equilibrium value of ±9.8e− within 10 ps. The question is then, given the dielectric response of the SPC/E electrolyte, does this average electrode charge produce the correct voltage drop? Figure 10 shows the potential drop across the simulated capacitor for values of Q0 = 0.0, 0.2, 0.4, and 0.6e−. The value of

3. APPLICATIONS 3.1. Equilibrium Simulation. This section presents results for the method described in the preceding sections applied to an equilibrium all atom simulation of a capacitor. The capacitor was composed of two parallel electrodes and an intervening liquid electrolyte. The two electrodes spanned the periodic dimensions X = 33.63 Å and Y = 33.24 Å, and consisted of three layers of FCC Pt atoms with a lattice constant of 3.92 Å. The (111) faces of the electrodes were exposed to the electrolyte and separated by a distance of 50 Å. The electrolyte was SPC/E water30 at a nominal density of 1 g/cm3. The electrostatic interactions were treated as detailed in section 4, the Lennard−Jones (LJ) parameters for the Pt atoms were σ = 2.41 Å and ε = 0.2013 kcal/mol used with arithmetic mixing. All LJ interactions were cut off at 2.5σ. Four simulation trajectories were created, each with a distinct potential drop between the electrodes of the capacitor. Each simulation was carried out in the NVT ensemble at 298 K for 1 ns. The desired average potential was achieved by setting Q0 = 0.0, 0.2, 0.4, and 0.6 e− (eq 13). From eq 13, a parallel plate capacitor with electrodes where A = (33.63)(33.24) Å2 and D = 50 Å, separated only by a vacuum, will have a potential drop between electrodes of V0 = 0.0, 1.62, 3.23, and 4.85 V. At first consideration it may seem a poor choice to base the electrode charge on the desired potential drop for a capacitor where the electrolyte is vacuum and not the simulation electrolyte of SPC/E water. However, the electrode charge (eqs 14a and 14b) is not only the preset value of Q0, but rather the sum of this preset charge and the charge due to the higher order image charges, which account for the response of the electrolyte to the imposed electric field. Consider the simulation system where the constant charge, Q0, added to the fluctuating electrode charge has been chosen to be 0.6e−. This constant charge quickly polarizes the SPC/E electrolyte adjacent to the electrodes which, in turn, increases the electrode charge through secondary image charges according to eqs 14a and 14b. Figure 9 shows the electrode charge

Figure 10. The potential drop as a function of position in the capacitor described in section 3.1 for applied electrode charges of Q0 = 0.0 (solid), 0.2 (dashed), 0.4 (dashed-dot), and 0.6e− (dotted line). These charges correspond to a capacitor voltage of 0.0 V, 1.62 V, 3.23 V, and 4.85 V, respectively.

the potential drop between the two electrodes from the simulations are V0 = 0.0, 1.62, 3.23, and 4.85 V, respectively, which matches the exact potential drop calculated (using eq 13) from the added constant charge, Q0, for the capacitor with a vacuum for the electrolyte. 3.2. Oxidation/Reduction Reactions. The redox events were treated using a KMC scheme as discussed earlier, which involves the determination of the probability of an oxidation/ reduction event of the electroactive species. These probabilities can be expressed in terms of the Marcus rates, kanode and kcathode, which were discussed in an earlier section. The probability, Poxidation, of an oxidation event taking place at the anode is given by2 Poxidation = 1 − e−kanodeΔt

(29)

where Δt is the time step of the simulation. The reduction event at the cathode is treated in a similar manner. In order to calculate the Marcus curves for the electrochemical cell two sets of equilibrium simulations, at zero applied voltage and with the same electrode set up as in the preceding section, were carried out − one with 1 M NaCl (34 Na+ and 34 Cl−) and the other with 0 M NaCl. The simulations were carried out in the NVT ensemble at 298 K for 1 ns after an initial equilibration of 100 ps. In each case, the cell contained 100 Na0 (the electroactive species) and 1630 water molecules. While Na0 does not exist in water, it is being used in this simulation study to represent the reduced form of the electroactive species in a prototypical test electrochemical cell. A number of configurations were selected at random from the trajectories and the energy differences, δE, were calculated for

Figure 9. Electrode charges for the anode (solid line) and cathode (dashed line) as a function of time for the system described in section 3.1 with an electrode charge of Q0 = ± 0.6e− applied to the electrodes. The starting configuration (which is the time origin for this plot) was taken from an equilibrium simulation with no charge applied to the electrodes, that is, Q0 = 0.0e−. It is seen that the charges quickly reach their equilibrium values of ±9.8e− from an initial value of 0e−.

as a function of simulation time for a simulation system started from a configuration where the electrolyte was equilibrated 4909



Article

the oxidation of a chosen Na0, in a given configuration, to Na+. These energy differences (circles 1 M NaCl and squares 0 M NaCl) are plotted in Figure 11 along with the distance dependence

At steady state, on average, the number of oxidation events is equal to the number of reduction events, and hence the total number of reactions is given by either the number of oxidation or reduction events. In order to study the effect of concentration on the electrochemical cell, 68 water molecules were removed from the cell, keeping all other species the same, and replaced with 34 Na+ and 34 Cl− ions giving rise to a 1 M solution of NaCl. The cell was next allowed to equilibrate at zero applied voltage in a regular MD simulation (with no KMC steps) for 200 ps. The resulting configuration was then used as the starting point in a series of MD-KMC simulations to determine the current as a function of applied voltage. Each simulation was carried out in the NVT ensemble at 298 K for 2 ns. The diffusion limit of the current (Ilim), namely, the current due to mass transfer, is a function of the diffusion coefficient of the electroactive species. The diffusion coefficient (D) of Na0 was determined from simulations of Na0 in bulk water and the limiting current was calculated using31

Ilim =

Figure 11. Scatter plot of δE as a function of the distance of the electroactive species from the left electrode for the equilibrium simulations with the electrochemical systems described in section 3.2. The circles and squares represent the oxidation of Na0 to Na+ at 1 M NaCl and 0 M NaCl, respectively. The diamonds and trianges represent the reduction of Na+ to Na0 at 1 M and 0.5 M, respectively. The distance dependence of the energy of a unit positive test charge (solid line), in the same cell with no electrolyte, as a function of distance from the left electrode is included for comparison.

n·A·D·F ·C d

(30)

where A is the surface area of the electrode, d is the distance between the electrodes, n is the number of electrons involved in the reaction (in this case 1), C is the concentration, and F is the Faraday constant. Figure 12 shows the current density (current per unit area of the electrode) in the electrochemical cell as a function of the

of the energy of a unit positive charge in a cell with no electrolyte. From the figure it is clear that, on average, δE is independent of concentration. In order to determine the reduction curves as a function of concentration the trajectories of two sets of simulations were examined. In the first case, the energy differences for the reduction of the Na+ to Na0 in the 1 M solution were considered. In the second case a simulation was carried out starting from the initial configuration of the 1 M case but with half the Na+ and Cl− ions removed resulting in a 0.5 M solution. The energy differences for the reduction of Na+ to Na0 are plotted in Figure 11 (diamonds 1 M NaCl and triangles 0.5 M NaCl), and it is clear that the reduction curves are independent of concentration as well. These oxidation and reduction curves were used to determine λ and ΔG0 and hence the electron transfer rates as described in section 2.3. The oxidation/reduction probabilities discussed above were utilized for performing a number of simulations at different applied potentials. These KMC-MD simulations incorporated KMC steps for the interconversion between the oxidized and reduced forms of the electroactive species in a regular MD simulation. During the course of the MD simulation, when an electroactive species was within 15 Å or less of an electrode, the probability of the conversion of this species was calculated using the Marcus rates discussed above. If the Monte Carlo step favored conversion, the charge on the species was changed keeping all van der Waals interaction terms the same. The electrochemical cell, with the same electrode configuration as in the previous paragraph, but with 1740 water molecules and 20 electroactive species that can interconvert between Na+ and Na0 states, was studied as a function of applied voltage. The steady state current was calculated as the total number of reactions taking place in a given length of time, with each oxidation or reduction reaction considered as a half reaction.

Figure 12. Plot of the current density (normalized using the limiting current density) − voltage curve from the nonequilibrium simulations of a test electrochemical system described in section 3.2. The circles are the results of the simulation with 0 M NaCl solution, while the straight line gives the best linear fit to the data. The squares give the simulation results for the 1 M case.

applied voltage for the 0 M (circles) and 1 M (squares) NaCl solutions. The applied voltage is in the range of 1−5 V which is within the range of the values used in operating regular batteries. Note that in these simulations water oxidation and reduction are being ignored. Interestingly, the current in the case of the 0 M solution is much higher than the diffusion limit of the current and is fairly linear in the applied voltage. In a typical electrochemical cell where ion transport is not effected by the electrode charge, the current increases rapidly with applied voltage for lower values of V, reaching a plateau or limiting value at higher applied voltages.31 At interelectrode separations of a few nanometers, this 4910



Article

electrode/electrode, and image charge/electrode pairs in the neighbor list. The positions of the electrode and image charge atoms are not integrated, so calculating the pairwise forces is unnecessary. However, excluding these atom pairs from the neighbor list does not eliminate their contribution to the particle− particle particle-mesh (PPPM)34 solver used to treat the electrostatic interaction. The PPPM method is intended for use with systems that are periodic in all dimensions. A typical work around for using PPPM with a slab-like system that is periodic in only two dimensions is to add a large “buffer” of empty space between the two ends of the simulation system in the nonperiodic dimension. The work around is motivated by the large computational cost of two-dimensional electrostatics methods, that is, 2D Ewald.35 Unfortunately, this work around does not fully alleviate the serious spurious effects of using the 3D PPPM method in a 2D simulation.24,36 These spurious effects are due to the net dipole of the system in the nonperiodic dimension. This is of particular concern for simulations of electrochemical cells given the large charge separation created by the charged electrodes. Fortunately there is a simple correction - proportional to the simulation system net dipole moment - that can be used with the 3D PPPM method.36 This correction is included in the LAMMPS base package, and is the method used to calculate all electrostatic interactions for the simulations presented here. 4.2. Redox Reactions. The KMC scheme to treat the oxidation/reduction reactions was implemented in the modified LAMMPS code. Each species in the cell is assigned an identity number in the LAMMPS control file, which in turn allows the specification of the oxidized and reduced forms of the electroactive species. If the distance of the reduced form of the electroactive species from the anode (oxidizing plane) is less than 15 Å, the probability of oxidation is calculated and is used to determine whether a reaction takes place. If a reaction does take place, the identity of the species is changed to that of the oxidized form. A similar scheme is incorporated for the reduction at the cathode. Note that for the sake of simplicity both the oxidized and reduced forms of the electroactive species have the same Lennard−Jones parameters, differing only in their Coulomb interactions due to the difference in charge.

plateauing will no longer take place, resulting in very large currents. This is the onset of “Coulomb transport”32 since the migration of ions now depends on the large Coulomb fields generated by the charges on the surface of the electrodes.1 This Coulomb transport decreases with increasing electrolyte concentration as is clearly seen from the figure. To the best of our knowledge, this is the first atomistic simulation of an electrochemical cell that shows the Coulomb transport regime.

4. COMPUTATIONAL DETAILS 4.1. Image Charges. A modified version of the LAMMPS33 MD package was used for all simulations in this study. The LAMMPS code was chosen because of the number of necessary preexisting features as well as its modular design, which simplifies the addition of features. A single routine consisting of two distinct operations was added to treat the image charge phenomenon. The positions of all atoms in the MD simulation, including atoms representing image charges, are provided to the LAMMPS program through a configuration file. When the program is initiated a control file is read, providing the location of the image planes. These image planes are chosen to coincide with the atomic centers of the first layer of electrode atoms. The first portion of the image charge routine then adjusts the image charge atomic positions consistent with the image charge positions described in section 2.1. Furthermore, the routine verifies that the charge of the image charges are consistent with the charge of the concomitant real charge. The first portion of the routine continues to update the positions of the image charges, relative to the real charges, after each step of the MD simulation, post integration of the real charge positions. If the real charges have variable charge, the routine continues to verify the image charges are consistent with the real charges and alters them if necessary. The second portion of the image charge routine deals with the variable charges of the electrode atoms. If a constant average potential drop between electrodes is desired, an electrode charge, Q0 of eq 13, can be specified in the LAMMPS control file. Then, during the course of an MD simulation, the charges of the surface atoms of each electrode are reassigned according to eqs 13, post integration. As discussed in section 2.1, this surface charge is uniformly distributed across all surface atoms. There are a number of preexisting LAMMPS features that were utilized to implement and increase the efficiency of the image charge model. LAMMPS spatially decomposes the simulation domain to parallelize the simulation over many processors. The only information one processor has about other atoms assigned to another processor concerns only those atoms residing in a “skin” or small overlapping spacial region of adjacent processors. In order to avoid the time-consuming operation of passing information between processors that contain only one atom of an image charge/real charge pair, the image plane routine requires that an image charge/real charge pair reside on the same processor. This is easily accomplished with LAMMPS, which allows for the specification of any arbitrary rectangular processor grid. Any processor grid that spans the direction normal to the electrodes with a single processor will ensure that both atoms of a given image charge/real charge pairs reside on a single spanning processor. LAMMPS allows for the exclusion of specific pairs from the neighbor list build. For the purposes of the present model, it is unnecessary to include image charge/image charge,

5. CONCLUSIONS A computationally efficient atomistic model of an electrochemical cell with an applied voltage has been developed and implemented in the LAMPPS software package. The polarization of the electrode by the electrolyte is incorporated using the method of image charges  a computationally cheaper alternative to explicit polarizable electrodes. The net fluctuating charge on each electrode, due to the higher order image charges of the electrolyte, is distributed uniformly over the surface atoms of the electrode as is the additional charge required to maintain the constant potential difference due to the applied potential. The probability for the oxidation/reduction reaction of a reactant in the electrochemical cell is then obtained from a Marcus Theory treatment of the process. These probabilities were used in simulations to model oxidation/reduction in an electrochemical cell with an interelectrode distance of 5 nm. The current−voltage curve generated from these simulations reveals the lack of a plateau at large applied voltages, a phenomenon that has been termed Coulomb transport. With the computational framework to model such electrochemical cells now in hand, our future research will involve detailed investigations into the nature of Coulomb transport. 4911



■

Article

(31) Bard, A.; Faulkner, L. Electrochemical Methods: Fundamentals and Applications; Wiley: New York, 2001. (32) http://www.chem.utah.edu/directory/faculty/white.html. (33) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19. (34) Hockney, R. W.; Eastwood, J. W. In Computer Simulations Using Particles; Taylor and Francis Group: New York, 1988; pp 267−304. (35) Heyes, D. M.; Barber, M.; Clarke, J. H. R. J. Chem. Soc., Faraday Trans. 2 1977, 73, 1485−1496. (36) Yeh, I.-C.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155− 3162.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the Office of Naval Research through ONR Grant N-00014-09-1-1107 and in part by a grant of computer time from the Department of Defense (DOD) High Performance Computing Modernization Program at the Army Research Laboratory DOD Supercomputing Resource Center.

■

REFERENCES

(1) Long, J. W.; Dunn, B.; Rolison, D. R.; White, H. S. Chem. Rev. 2004, 104, 4463−4492. (2) White, R. J.; White, H. S. Langmuir 2008, 24, 2850−2855. (3) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471−2474. (4) Blumberger, J.; Bernasconi, L.; Tavernelli, I.; Vuilleumier, R.; Sprik, M. J. Am. Chem. Soc. 2004, 126, 3928−3938. (5) Blumberger, J.; Sprik, M. J. Phys. Chem. B 2004, 108, 6529−6535. (6) Blumberger, J.; Tateyama, Y.; Sprik, M. Comput. Phys. Commun. 2005, 169, 256−261. (7) Izvekov, S.; Mazzolo, A.; VanOpdorp, K.; Voth, G. A. J. Chem. Phys. 2001, 114, 3248−3257. (8) Izvekov, S.; Voth, G. A. J. Chem. Phys. 2001, 115, 7196−7206. (9) Lamperski, S.; Klos, J. J. Chem. Phys. 2008, 129, 164503. (10) Straus, J. B.; Calhoun, A.; Voth, G. A. J. Chem. Phys. 1995, 102, 529−539. (11) Calhoun, A.; Voth, G. A. J. Phys. Chem. 1996, 100, 10746− 10753. (12) Crozier, P. S.; Rowley, R. L.; Henderson, D. J. Chem. Phys. 2000, 113, 9202−9207. (13) Fedorov, M. V.; Kornyshev, A. A. J. Phys. Chem. B 2008, 112, 11868−11872. (14) Feng, G.; Zhang, J. S.; Qiao, R. J. Phys. Chem. C 2009, 113, 4549−4559. (15) Wernersson, E.; Kjellander, R. J. Chem. Phys. 2006, 125, 154702. (16) Philpott, M. R.; Glosli, J. N. J. Electrochem. Soc. 1995, 142, L25− L28. (17) Spohr, E. Electrochim. Acta 1999, 44, 1697−1705. (18) Glosli, J. N.; Philpott, M. R. Electrochim. Acta 1996, 41, 2145− 2158. (19) Siepmann, J. I.; Sprik, M. J. Chem. Phys. 1995, 102, 511−524. (20) Reed, S. K.; Lanning, O. J.; Madden, P. A. J. Chem. Phys. 2007, 126, 084704. (21) Willard, A. P.; Reed, S. K.; Madden, P. A.; Chandler, D. Faraday Discuss. 2009, 141, 423−441. (22) Reed, S. K.; Madden, P. A.; Papadopoulos, A. J. Chem. Phys. 2008, 128, 124701. (23) Griffiths, D. J. Introduction to Electrodynamics, 3rd ed.; Benjamin Cummings: Upper Saddle River, NJ, 1998. (24) Spohr, E. J. Chem. Phys. 1997, 107, 6342−6348. (25) Marcus, R. A. J. Chem. Phys. 1956, 24, 966−978. (26) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta, Bioenerg. 1985, 811, 265−322. (27) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. J. Comput. Phys. 1975, 17, 10−18. (28) Voter, A. F. In Radiation Effects in Solids; Sickafus, K. E., Kotomin, E. A., Uberuaga, B. P., Eds.; NATO Science Series; Springer: Netherlands, Dordrecht, 2007; Vol. 235, Chapter 1, pp 1−23. (29) Smalley, J. F.; Feldberg, S. W.; Chidsey, C. E. D.; Linford, M. R.; Newton, M. D.; Liu, Y.-P. J. Phys. Chem. 1995, 99, 13141−13149. (30) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269−6271. 4912


A Computationally Efficient Treatment of Polarizable Electrochemical

Recommend Documents