Constant Potential, Electrochemically Active Boundary Conditions for

Subscriber access provided by Nottingham Trent University

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Constant Potential, Electrochemically Active Boundary Conditions for Electrochemical Simulation Kaitlyn Anne Dwelle, and Adam P. Willard J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.9b06635 • Publication Date (Web): 05 Sep 2019 Downloaded from pubs.acs.org on September 5, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Constant Potential, Electrochemically Active Boundary Conditions for Electrochemical Simulation Kaitlyn A. Dwelle and Adam P. Willard∗ Department of Chemistry, Massachusetts Institute of Technology, Cambridge, MA E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract In this manuscript we present a model for simulating active electrochemical systems using a classical molecular dynamics framework. We describe a computationally efficient method of enforcing the electrostatic properties of constant potential boundary conditions and demonstrate how this method can be adapted to support stochastic interfacial charge transfer processes. We highlight the utility of this model by simulating the nonequilibrium dynamics of a model battery system. We demonstrate the ability of this model to support the formation of a stable double structure, consistent with expectations from macroscopic equilibrium. We also illustrate how this model can be used to provide microscopic physical insight into the results of standard potential jump experiments.

1

Introduction

Chemical reactions at electrode interfaces are generally sensitive to the conditions of their local molecular environment. Heterogeneity in these conditions – for example due to thermal fluctuations in local electrolyte composition – cannot be resolved by most electrochemical measurements, yet can play a central role in determining the overall kinetics and thermodynamics of a given reaction. 1,2 This heterogeneity can be modeled using molecular dynamics (MD) simulations, however, incorporating the influence of chemical reactivity into MD simulations presents several challenges. In this paper, we describe an efficient model based on classical MD that is capable of simulating the dynamics of interfacial charge transfer in a fluctuating electrolyte solution confined between constant potential electrodes. We demonstrate the utility of this model by presenting results of nonequilibrium simulations of a generic reactive electrochemical system subject to an instantaneous change in the applied electrode potential. As we highlight, these simulations can be extended to compute direct experimental observables and relate them to the underlying molecular scale electrochemical dynamics. 2


Page 2 of 31

Page 3 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The reactions that drive most electrochemical systems take place at the electrode-electrolyte interface. The unique properties of this interface are primarily due to the physical characteristics of the electrode and its influence on the adjacent electrolyte. The electrode serves as a physical barrier that can support a tunable buildup of surface charge. This barrier breaks the translational symmetry of the electrolyte, promoting the emergence of anisotropic and possibly highly correlated interfacial molecular structure. The presence of surface charge provides an electrochemical potential gradient, and an associated potential drop, that ultimately leads to the formation of the electrical double-layer. In reactive systems, the electrode provides a source or sink of electrons, thereby facilitating (and sometimes catalyzing) the redox reactions that drive the flow of charge (and sometimes mass) in driven electrochemical applications such as batteries. Capturing these interfacial properties in a single model framework is challenging because their interactions span a wide range of characteristic time and length scales. A common solution to this challenge is to model the electrolyte as continuum that interacts empirically with the electrode. 3–5 Continuum modeling approaches have been widely used in the analysis and interpretation of electrochemical measurements because they are highly efficient and are easily extended to experimentally relevant time and length scales. However, because these models are highly parameterized and contain very few specific molecular details, they are not reliable as a predictive framework and therefore of limited use as a basis for molecular insight and design. Another common approach to modeling electrode-electrolyte interfaces is to utilize simulation methods based on first-principles electronic structure calculations. Ab initio molecular dynamics (AIMD), usually based on density functional theory (DFT), provides the ability to explicitly describe the electronic rearrangements involved in electrochemical reactions. 6 First-principles electronic structure and AIMD have been extensively used to compute reaction mechanisms 7,8 and energy barriers, 9,10 as well as to describe molecular structure and dynamics at the electrode-solvent interface. 11 Despite having generated an enormous amount

3



of valuable physical insight, these approaches are generally limited in scope to very small systems (typically 100s of atoms over 10s of ps) due to the inherent computational expense of electronic structure calculation. Connecting the results of these first-principles studies to experimentally relevant systems therefore often requires assumptions about the role of fluctuations, disorder, and molecular correlations in the extended system. Classical molecular dynamics (MD) simulation can effectively bridge the system size gap between continuum and first-principles modeling approaches. However, traditional force fields lack the functionality to model constant potential electrodes or to simulate chemical reactivity. Over the past 20 years, numerous methodological advances have targeted this lack of functionality. This includes the development of methods for simulating constant potential electrodes, 12,13 reactive force fields for simulating bond making/breaking, 14,15 and stochastic approaches to modeling interfacial electron transfer events. 16 Despite these separate advances, there have been no MD-based models that combine both tunable constant potential electrodes and the capability for interfacial electron transfer which is sensitive to fluctuations in the electrolyte system. In this paper, we present a model for elecrochemically active and electrostatically consistent electrodes held under constant potential conditions that is fully compatible with standard classical molecular dynamics. Specifically, this model (i) treats the constant potential boundaries in a computationally efficient manner, (ii) allows for outer-sphere oxidation and reduction at both electrodes, and (iii) can include the effect of ion intercalation into the electrode, so that both the total charge and number of redox-active particles in the electrolyte can fluctuate throughout the course of the simulation. We demonstrate the application of this model by simulating a nanoscale Li-ion-type battery system, such as illustrated in Fig. 1. The physical processes that determine the performance of this type of battery system, namely interfacial electron transfer and the transport of ions across the interface and through the electrolyte, cannot be simulated with standard classical MD. The functionality provided by our model enables a computationally efficient simulation of the microscopic dynamics of the

4


Page 4 of 31

Page 5 of 31

electrolyte-electrolyte interface under operating conditions.

—

e-

—

(1a) e-

+

+

(2)

+

+

—

(1b) μcathode

μanode

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∆V V0

V1

Figure 1: A diagram of a model battery system. Two electrodes are held at a constant potential difference, ∆V . V0 and V1 set the chemical potential of electrons in the electrodes, µanode and µcathode , respectively. Current is driven by the preference for oxidative formation of cations at the anode (1a at the left-hand electrode) and reductive elimination of cations at the cathode (1b at the right-hand electrode), and mediated by the diffusion of cations across the electrolyte (2). Thus, the overall performance is dominated by ion migration through the electrolyte and redox reactions at the electrodes. One of the primary challenges in modeling an electrochemical cell is efficiently handling the electrostatic interactions between the electrodes and the electrolyte. To capture the correct electrostatic behavior, proper treatment of the simulation cell boundary conditions and electrode polarizability are essential. Because of the long-range nature of Coulomb interactions, distance-based cutoff methods give the incorrect asymptotic behavior 17 and the periodicity of the system needs to be considered carefully. Here we focus specifically on simulation cells with parallel planar electrodes positioned at the boundaries of the z directions and periodically replicated in the other x and y directions. Figure 2 shows the difference in long range symmetry between a bulk system periodically replicated in all directions and a two-dimensional slab system appropriate for modeling an interfacial system. Summation techniques for the calculation of electrostatic forces and energies in a two-dimensional slab system exist, but are more computationally expensive and less numerically robust than a standard three-dimensional calculation. 18 Another challenge in modeling electrochemical systems is describing the surface charge 5



I.

II. +

+

+

x y

z

Page 6 of 31

+

+

+

+

+

+

+

+

+

Figure 2: A schematic comparison of I. a fully three-dimensional periodic system and II. a slab system with finite width in the z direction. Both systems are illustrated in twodimensions for convenience. of the electrode. A straightforward approach is to treat the electrode as a constant and uniformly charged surface. However, despite being easy to implement and computationally efficient, this approach lacks important physical effects, such as the ability to polarize in response to local charge fluctuations, that can qualitatively alter the static and dynamic properties of the electrode-electrolyte interface. 19 Under constant potential conditions, surface charge distribution polarizes in response to nearby charge density. A variational procedure, introduced by Siepmann and Sprik 20 and applied to electrochemical cells by Madden and coworkers, 12,21 enforces the constant potential condition by allowing the charge on individual electrode atoms to fluctuate based on the surrounding environment. While computationally expensive, this method has the distinct advantage of being able to enforce a constant potential condition in systems with non-planar electrodes. In systems with planar electrodes held at zero potential, it has long been understood that the electrostatic potential in the system can be obtained through the use of image charges. 22–24 Perram and Ratner have shown that the periodic repetition in the image charge solution can be included in an MD simulation using full periodic boundary conditions and a double cell. 25 More recently, Voth and coworkers have developed a computationally efficient procedure for enforcing constant potential boundary conditions based on use of explicit first order image charges and a uniform correction term. 13 Here, we build upon these efforts by reformulating the image charge approach for describing constant potential electrodes so that the forces and energies can be

6


Page 7 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


computed using the more efficient, full three dimensional Ewald summation. Methods for including redox reactions in classical molecular dynamics necessarily neglect most electronic detail. In this paper, we limit our scope to simple Marcus-type outer sphere electron transfer, where the active species is not electronically coupled to the electrode. The standard method for calculating electron transfer rates using classical MD simulation are based on the generation of free energy surfaces for the reactant and product states by sampling the equilibrium fluctuations in the vertical energy gap (i.e., the canonical Marcus theory reaction coordinate). With this method, the average kinetics of the electron transfer can be inferred from Marcus free energy curves and the average Marcus rate can be applied to treat the interconversion of products and reactants at an electrode boundary in nonequilibrium simulations. 13 To account for the effects of spatial and temporal variability in electron transfer rate, Subotnik and coworkers have developed a surface hopping method where the probability of electron transfer for a given redox active species is related to the instantaneous value of its vertical energy gap. This method inherently encodes the concentration dependence of the rate and allows for the study of correlations between electron transfer events. We build upon this method to include the mass transfer associated with the exchange of electrochemically active ions between the electrolyte and the electrode material. The structure of the paper is as follows: First, in section 2 we describe a computationally efficient method for enforcing the constant potential boundary conditions and compare the results and complexity to existing methods. We then describe the model for electron and mass transfer. The details of the implementation of the two methods are briefly described in section 2.4. Finally, in section 3 we present results from a nonequilibrium simulation of a model battery system.

7



2

Methods

2.1

Enforcing Constant Potential Boundary Conditions

The distribution of charge within a constant potential electrode can fluctuate in response to the dynamics of charged species within the electrolyte. The ability to properly describe these fluctuations is essential for understanding the equilibrium and nonequilibrium properties of the electrode-electrolyte interface. Here, we consider model systems that are bounded by two parallel, planar, polarizable electrodes, each with an independently assigned voltage. Each electrode imposes a boundary condition where the value of the electrostatic potential is equal to the assigned electrode voltage at all points along the electrode surface. We define our coordinate system so that the electrode boundaries lie along the z-axis, with one electrode at position z = 0 and the other at position z = d, as illustrated in Figs. 2 and 3. We decompose the constant potential conditions in to two components. One component that describes the polarization fluctuations that are native to neutral constant potential electrodes and one component that describes the uniform electric fields that result from charging the electrodes. This decomposition is possible due to the additive nature of electrostatic interactions. We enforce the neutral component of the constant potential condition at each boundary using the method of image charges reflected across the two planes at z = 0 and z = d. Reflecting all charges, including images generated by previous reflections, across both electrode planes leads to a periodic system consisting of the active system of interest, a reflected system, and infinite periodic replicas as illustrated in Fig. 3. The active and reflected system together compose an electrostatically neutral repeat unit, extending from z = −d to z = d, which is periodic in all three dimensions. Perram and Ratner showed using a Green’s function formalism that the total electrostatic energy in such a system is exactly twice that of a slab system sandwiched between two ideal metal boundaries. 25 This configuration ensures that the potential at z = 0 and z = d are zero, while also transforming the system to a

8


Page 8 of 31

Page 9 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fully three-dimensional periodic system (see Fig. 2I) where electrostatic interactions can be calculated using the more computationally efficient and robust family of fully three dimensional k-space summation techniques (e.g. pppm 26 ). Further computational speedup can be gained by exploiting the symmetry of the unit cell. 25 To enforce the electrostatic effect of a potential drop ∆V across the cell, we introduce a uniform field of magnitude ∆V /d between z = 0 and z = d along the z-axis which preserves a constant potential along each image plane. Notably, in an isolated system the absolute value of the reference potential is irrelevant and so we pin it to the z = 0 electrode. As such, this system setup corresponds to an electrostatic potential of Ψ = 0 at z = 0 and Ψ = ∆V at z = d, as shown in more detail in Appendix A. In theory, the field generated by the potential drop has magnitude −∆V /d in the region z = −d to z = 0 giving a repeat unit with a net field of zero. In practice, only the forces in the active system need to be calculated to propagate molecular dynamics since the configuration of the reflected system is completely determined by the configuration of the active system. -d

x y

z=0

z

d

repeat unit

+

+

+

-

-

active system reflected system

Figure 3: A schematic representation of the periodic system. The two parallel image charge boundaries at z = 0 and z = d lead to repeating symmetry due to a recursive introduction of image charges. The periodicity in the y dimension is directly analogous to the x dimension. The blue lines show the electrostatic potential drop due to the uniform field.

2.1.1

Active System Energy

The method of image charges ensures that the electrostatic potential is correct within the active system (i.e. the volume confined between the electrodes), but interactions between 9



Page 10 of 31

particles in the reflected system are not physically meaningful. The energy of the active system therefore includes all active-active and active-image coulomb interactions, but excludes all contributions from image-image interactions. The forces on the active particle are derived from this energy. Since force calculation can be restricted to atoms within the active system, propagating the dynamics of the active system using molecular dynamics is straightforward. Including short range, nonbonded energy terms and the electrostatic potential energy described above, the total active system potential energy Uact is given by,

Uact ({r}, ∆V ) = USR ({r}) + UCoul ({r}) + UF ({r}, ∆V ),

(1)

where {r} denotes the nuclear configuration of all members of the active system, USR is the potential energy of the system due to short range (e.g., Lennard-Jones-type interactions) within the active system (and its periodic replicas), UCoul is the electrostatic energy of the active system (including its interactions with image charges), and UF is energetic contribution of the potential drop on the active system. The short-range contribution is typically given by, USR ({r}) =

Nact X

wSR (ri , rj ) +

Nact X

wwall (ri ),

(2)

i

hiji

where the summation denotes a sum over all unique pairs of particles, wSR is the pair-wise shortwave interaction potential between particles, and uwall describes the interaction between particles and the electrode surface. The Coulomb contribution can be expressed as,

UCoul ({r}) =

Nact X

qi Φ(ri ),

(3)

i

where qi is the charge on species i and Φ(r) is the electrostatic potential (due to active species, inactive species, and their periodic images) evaluated at position r. the contribution

10


Page 11 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from the potential drop is given by,

UF ({r}, ∆V ) =

Nact X

∆V qi zi /d,

(4)

i

where zi is the position of species i along the z-axis. From a practical standpoint, we want to implement the system in such a way that the total active system energy can be calculated using existing standard molecular dynamics (MD) routines. For the Coulomb interactions, the active system energy is exactly half the total energy of the repeat unit which can be calculated using an Ewald summation or faster particle-mesh techniques. 25,27 For short-range interactions, treating the reflected system atoms as point charges with no short-range interactions means that the calculated total short-range energy can be used as-is. For the contribution to the potential due interactions with the constant field, only charges in the active system should be included.

2.1.2

Comparison to existing methods

As discussed in the introduction, there are several existing methods for including electrode boundary conditions in MD simulations of electrochemical systems. Table 1 shows a summary of the highlights and limitations of each method. The first two columns indicate whether a method correctly reproduces the phenomena of conducting electrodes at zero applied voltage and when the boundaries are held at an applied potential. The third column indicates whether the method can be applied to systems without planar electrodes. The last column uses the number of charged particles in a system as a proxy for the computational cost of the method. N is the number of charged particles in the active system and this column gives the total number of charged particles required in the total system. Since calculating electrostatic interactions scales roughly as M 3 , where M is the number of charged particles, it is often the most computationally expensive step of MD simulation. Although conceptually simple and computationally efficient, constant distribution of

11



charge on the electrode surfaces fails to reproduces the results of experimentally relevant conducting, or constant potential, boundary conditions. By adding an explicit term for the applied potential to the method developed by Perram and Ratner, 25 we are able to leverage the computational benefits of fewer charged particles and a fully three dimensional periodicity in a system with a potential bias. Table 1: Comparison of methods for including electrode boundary conditions Const. Potential at ∆V = 0 constant charge x Ratner 25 X 12 Madden X Voth 13 X This paper X

Const. Potential Non-planar System System at ∆V 6= 0 Electrodes Size Symmetry x X N 2D-slab – x 2N 3D X X –a 2D-slab X x 3N 2D-slab X x 2N 3D

a

This method adjusts the charges on the electrode atoms in the simulation through an iterative procedure. Each step of this iteration requires a full electrostatic calculation so the computational cost also includes a factor of the number of iterations required.

In order to verify the equivalency of this method to previous methods, we compare the energy contributions of a test charge in a simulation held at constant potential. We have recreated a version of the validation system, System I, used by Voth and coworkers where the image planes reside at z = 0 and z = 51.1 ˚ A. 13 The cell is 50 ˚ A in the x and y dimensions. A single charge of q = −e was placed at the center of the x and y directions and energy was measured as a function of z position. For comparison to the energies obtained by Voth and coworkers, an arbitrary energy constant was added to the results from these simulations since there are no explicit electrode atoms in this simulation and an addition of a constant energy term does not affect the magnitude of the calculated forces used to propagate dynamics in section 3. Figure 4 shows that the total system energy for System I as a function of the test charge position. In agreement with previous methods, the total energy of a charged particle between biased electrodes includes a nonlinear term due to image charge effects and a linear term whose slope is determined by the applied potential. 12


Page 12 of 31

Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


50

energy due to images

0

-50

e total en

-100

-150 0

rgy

ference ntial dif te o p ue to nergy d

e 10

20

30

40

50

Figure 4: Energy of a test charge as it is moved across an empty cell with constant potential boundary conditions. A single charge of q=e was placed at the center of the cell and total energy was computed as a function of z position. Orange (open) shapes are the results from the iterative method by Siepmann and Sprik,. 20 Black (closed) shapes are the energies obtained by Petersen et al. 13 and Blue lines are the result of the method described in this paper. Squares and the solid line represent the total energy of the active system. This energy is decomposed into two terms, long-dashed and circles mark the contribution to the total energy due to all image charge interactions and the short-dashed lines and triangles mark the contribution due to the applied potential.

2.2

Electrochemical boundary conditions

We model charge transfer as a stochastic process with a rate that fluctuates in response to the dynamics of the electrolyte. In this manuscript, we consider the most simplified case of coupled ion-electron transfer, analogous to the intercalation or deintercolation processes that drive the performance of Li-ion batteries. We model this process in the limit where adiabatic effects are fast compared to the simulation time step and thus are only described empirically. Specifically, we express the probability per unit time for creating a new ion or removing an existing ion (i.e., deintercalation or intercalation, respectively) as a function of the energy change to create a neutral particle at the interface and then pass a charge to it (or vice versa). The procedure we describe can be easily simplified to model the case where only charge (and not mass) is transferred. As an example of our procedure, consider the deintercalation of a cation, as illustrated in Fig. 5. We model this concerted process as the product of two sequential steps. We first insert a neutral dummy particle into the electrolyte at the electrode interface and then we 13



Page 14 of 31

pass an electron from the particle to the electrode to create a cation. We express the total probability for this event as a product of probabilities to perform each step independently. More formally, we express the probability for adding a new cation to the system as,

(1) (2) Pox = Pox Pox ,

(5)

(1)

where Pox is the probability to insert a neutral particle at a specific position in the system (2)

and Pox is the probability for outer-sphere electron transfer from the neutral particle to the electrode.

ΔUox

oxidation

(1)

0

ΔUred(2)

reduction

ΔUox(2) -

+

ΔUred(1)

Figure 5: The leftmost panel shows the system where a site (shown as a dotted circle) might be the site of a new cation. The middle panel shows the change in nuclear positions associated with the insertion of a neutral test particle. The final panel shows the electronic (1) rearrangement necessary to result in a positively charged ion in the simulation box. ∆Uox (2) and ∆red involve a change in nuclear configuration and thus are compared against a Boltz(2) (1) mann distribution to determine a probability of success. ∆Uox and ∆red involve a change in electronic configuration and are compared to the density of electronic states in the electrode. The probability for the first step is given by a standard grand canonical Monte Carlo acceptance criteria, (1) (1) Pox = min 1, exp −β∆Uox ,

(6) (1)

where 1/β = kB T – the Boltzmann constant times temperature – and ∆Uox is the energy change associated with the addition of the neutral particle. 28 This expression does not include a term for the chemical potential or fractional volume element of the neutral test particle since the neutral particle is replaced with a positively charged ion if the oxidation is successful. 14


Page 15 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The probability for the second step is taken from the method of Subotnik and coworkers, 16 where the energy to remove an electron from the simulation is compared to the distribution of unoccupied electronic states in a metal. Specifically,

(2) Pox =

Γ(z)

(2)

1 + exp Velectrode − β∆Uox

,

(7)

where Γ(z) is the distance dependent electronic coupling between the particle and electrode, (2)

Velectrode specifies the Fermi level of the electrode, and ∆Uox is the energy change associated with the charge transfer, which includes contributions from all induced image charges. We use Γ(z) = Γ0 exp(−∆z/λ), where Γ0 is the normalized maximum coupling strength, ∆z is the closest distance between the particle and the electrode surface, and λ is a decay parameter. We define the probability for the reverse process, i.e., cation reduction and intercallation, analgously to Eq. 5. Specifically, we define, 

 Pred =

(1)

(1) (2) Pred Pred

=

Γ(z)

(1)

1 + exp β∆Ured − Velectrode

h i (2)  min 1, exp −β∆Ured ,

(8)

(2)

where ∆Ured and ∆Ured correspond to the change in system energy associated with electron transfer from the electrode to the cation and the removal of the remaining neutral dummy particles, respectively. This method describes a process where the change in charge state is accompanied by an intercalation of the species into the electrode, removing it from the electrolyte. In a system where an electrochemically active species acts as a redox couple (i.e. both charge states are solvated species) the particle insertion/removal steps are no longer required and this method reduces to the model presented by Subotnik and coworkers. 16 Notably, this idealized model neglects several important details. In particular, the solidelectrolyte interphase (SEI), which can be hundreds of nanometers thick, is absent. Since

15



this phase can have its own energy barrier for intercalation 29,30 and even support part of the double layer formation, 31 its inclusion will be important in future models.

2.3

Simulation Details

All simulations were run using the open-source molecular dynamics software LAMMPS. 32 Simulations included two oppositely charged ionic species with the same mass and radius. The simulation box was 70 ˚ A by 70 ˚ A in the x and y directions with a 140 ˚ A electrode separation giving a total periodic z distance of 280 ˚ A including the explicit image charges. The initial concentration of each species was 1.2 mol/L . Solvent interactions were treated implicitly using a Langevin thermostat with a damping time of 30 fs and temperature of 302 K. Ion-ion interactions were treated using a Lenard Jones 12-6 potential and Coulomb interactions were calculated using the pppm method. 26 The ion-electrode short range interaction was included using a wall interaction derived by integrating over a three-dimensional half-lattice of Lennard Jones 12-6 particles. 33 Simulations were run with a timestep of 5 fs. Figure 6 shows a rendering of the active system (i.e. not including reflected system) for scale.

electrode

electrode

14 nm

7 nm

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

Figure 6: A snapshot of the active system. Blue particles are negative ions and orange particles are positive ions. The system is 7 nm by 7 nm by 14 nm.

2.4

Implementation Details

Separate fixes were developed to handle the image charge placement and the electrochemical boundaries. For the image charges, the fix ensures that every charged atom in the active 16


Page 17 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


system has a corresponding oppositely charged particle reflected across the image plane. The fix also keeps track when charged particles are added or removed from the system and adjusts the image system accordingly which is necessary for the correct computation of the (2)

(1)

change in energy due to the addition of a charge (e.g. ∆Uox and ∆Ured ). For the electrochemical reactions, the routine is summarized below: 1. Attempt frequency: an appropriate frequency must be chosen to control how often a redox event is attempted. Since each trial requires a full system energy calculation, we choose a frequency that allows for a well-sampled system while avoiding computational waste due to unnecessary attempts. This frequency can be adjusted to model the dynamics of specific chemical systems. 2. Type selection: it is first decided whether a redox attempt will be an oxidation or a reduction through a random selection. The probability of a reduction attempt is freduction and oxidation is 1 − freduction . 3. Site selection: The trial redox event is attempted at a randomly selected site. For computational efficiency, the statistics of Γ(z) are used to weight the site selection such that it is exponentially more likely that a site closer to the electrode is selected for evaluation. For oxidation, the selected site becomes the center of mass for the insertion of the test particle. For reduction of an existing particle, the selected site serves as the center of a sphere with radius rsearch , within which any active cation may be reduced. The exact relationship between freduction , rsearch , and the volume concentration of active and inactive species is not considered in this work. Simulations in this work use freduction = 0.9 and rsearch = 2.8 ˚ A 4. Acceptance criteria: equations 5 and 8 are used to calculate the probability of acceptance based on the energy differences. Note that the Γ(z) term is included in the site selection statistics and is thus not included explicitly in the evaluation of equations 5 and 8. This probability is then compared to a pseudo-random number and on a 17



successful attempt the number of particles in the active system changes by one. On a failed attempt, the system is reset to the previous configuration.

Results and Discussion Emergent Double Layer Formation Normalized Density

3.1

Negative Ion Density

3 2

2500

1 0

0

14

28

42

56

70

84

98

112

136

140

Normalized Density

Distance (Å)

2000

Positive Ion Density

3

1500

2 1 0

0

14

28

42

56

70

84

Distance (Å)

98

112

136

140

1000

Applied Field (MV/cm)

3

Normalized Density

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 31

Net Charge Density 2

500

0 -2

0

14

28

42

56

70

84

98

112

136

140

0

Distance (Å)

Figure 7: The equilibrium distribution of ions as a function of distance from the left electrode. Upper and middle panels show the density of positive an negative ions normalized by the average density with no applied field respectively. Field strength increases as the color changes from teal to dark blue. The dashed black lines show the linearized PoissonBoltzmann solution for the highest applied field in the region near the left electrode. Lower: The net charge density obtained by subtracting the negative ion density from the positive ion density. One important benchmark is that the model is efficient enough to reach equilibrium conditions under an applied field. For an electrochemical cell, this field arises from the differential charge buildup on the electrode surfaces when held at a potential. Migration of the ions to screen the field through the formation of a double layer structure can take times ranging from microseconds to minutes (or longer) depending on the total concentration of 18


Page 19 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mobile charge carriers, viscosity of the electrolyte, and distance between electrodes since ions may need to diffuse nearly the entire length of the cell. Gouy-Chapman theory proposes that the distribution of these screening ions can be understood by taking into account the electrostatic interactions between charged ions and the surface as well as the entropic contribution of ion position. 34 The resulting interfacial potential can be solved using the Poisson-Boltzmann equation, and results in a diffuse layer over which the interfacial field is screened. More sophisticated models include the effects of finite-volume charge carriers and short-range interactions on the double layer structure. 35,36 For simplicity, we will compare the results of our simulations to the linearized PoissonBoltzmann solution even though we expect the lack of excluded volume interactions to lead to a more compact predicted double layer. Figure 7 shows the average equilibrium distribution of ions for simulations with an applied field ranging from 300 to 2700 MV/cm. Densities are normalized to the bulk concentration at equilibrium and the net charge density is computed by subtracting the negative charge density from the positive charge density. As expected, the Gouy-Chapman solution, shown for the highest field as a dashed black line, predicts a more compact double layer structure since charge can build up at the interface to an arbitrary density whereas real systems, and this model, are limited by the finite volume of charge-carrying ions.

3.2

Current-Voltage Response

One strength of this model in comparison to those that use a fixed reaction rate is the ability to better understand the non equilibrium behavior before a steady state is reached. In order to verify that this microscopically defined model correctly reproduces the expected macroscopic response, we computed an ensemble of trajectories where the system was first equilibrated for 140 ns with ∆V = 0mV and at time t = 0 the potential was instantaneously increased to ∆V = 350mV. The current density is computed by keeping track of the number of successful oxidations 19



∆V (mV)

400 300 200 100 0

Current Density (A/cm2)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

0

200

400

0

200

400

600

800

1000

1200

1400

600

800

1000

1200

1400

Time (ps)

300 250 200 150 100 50 0 -50

Time (ps)

Figure 8: The response of the model to an applied potential step. Upper: the difference in potential between the two electrodes as a function of time. Lower: the solid black line plots the current density through the system measured as a function of time. The red dashed line is a fit to a single exponential decay plus an offset. and reductions at both the left and right boundary. The total average flux through the − 1 + #e system is then 12 abs #e− in,left out,right = 2 abs ([nred,left − nox,left ] + [nox,right − nred,right ]) per time interval. nox and nred are the number of successful reductions or oxidations at the given boundary. It should be noted that while this count of electron transfer events gives the Faradaic current through the system, this definition does not include the capacitive current generated by the movement of charge within the cell. At steady-state, however, the total current in an electrochemical cell has no contribution from nonfaradaic, or capacitive, current. A future model which includes explicit enforcement of a voltage differential as discussed in the previous section should address these issues. Figure 8 shows the current density as a function of time averaged over 1000 independent trajectories. At t < 0, fluctuations around zero are an indication that even at equilibrium, there is still electrochemical activity with deviations from zero driven by thermal fluctuations. At t = 0, the spike initial current in driven by the large difference between the initial electrochemically active ion concentration at the interface and the concentration at steady state. The decay to steady-state is fit well by an exponential decay plus an offset equal to 20


Page 21 of 31

the steady-state current density, shown in the figure by the orange dashed line. Typically, current decay due to a diffusion-limited depletion of active concentration at the interface is characterized by a t−1/2 dependence, 34 however we observe an exponential decay which is characteristic of bulk depletion of a reactant due to a first order reaction. This is most likely an indication that the chosen set of parameters favors a charge transfer limited current instead of a diffusion limited current at early times. The current density for the simulation is notably much higher than existing electrochemical cells. Factors contributing to this include the inherently smoother dynamics of the implicit solvent and the large surface area to volume ratio of the nanoscale cell. Incorporation of more chemical detail into the solvent dynamics will lead to slower reaction rates. The rate can also be tuned by varying the electronic coupling decay parameter λ or the trial frequency for charge transfer events. Additionally, the current model describes a system with zero overpotential, but its inclusion would also

Cation Density

2

1400

1.5

1200

1

0.5 0

0

1000 28

56

84

Distance (Å)

112

140

Anion Density

2

600 400

1.5 1

200

0.5 0

800

time (ps)

Normalized Density

provide another physically relevant parameter for obtaining more realistic current densities.

Normalized Density

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

28

56

84

Distance (Å)

112

140

0

Figure 9: The early-time concentration density as a function of Upper: the normalized concentration profiles of the electrochemically active cations. At t = 0 (lightest) the voltage in stepped from zero to 350mV. The concentration then decays at the left electrode and increases at the right electrode, driving the change in current density seen in figure 8. The darkest line is the density averaged over the t = 1.1ns to t = 1.4ns after the application of the potential difference. Lower: The normalized density profiles of the electrochemically inactive counterion (anion). 21



The response of the system to the potential jump can also be understood through the response of the ion concentration profiles. Figure 9 shows the ion densities as a function of time over the first 1400 ps after the application of the voltage difference. The cation density at the left interface decays and the density at the right interface increases, setting up a concentration gradient across the cell. The initial spike and decay in observed current is due to the large charge transfer rate at the beginning when the chemical potential of the system at the interface is far away from the chemical potential enforced by the electrode boundary conditions. Additionally, at the left interface, the larger initial concentration of positive ions at the interface leads to an increase in the probability of any one ion being reduced. As the concentration at the interfaces approaches the steady state concentration, the reaction rate at the interface decays to the steady state rate. The anion density changes much less dramatically than the cation density since any changes are second order effects due to the Coulombic interactions with the density profile. The small increase in ion concentration at the interface is due to image charge stabilization, the same interactions which lead to the decay in free energy near the interface for an isolated ion shown in figure 4.

4

Conclusions

In this paper, we have implemented an electrochemically active model of a complete electrochemical cell with molecular detail. By adding a term to account for the voltage bias, we were able to leverage the efficiency of an existing method 25 to compute accurate electrostatic interactions. The electrochemical reactions were included as stochastic events using both the vertical energy gap due to a charge transfer and the change in energy due to a particle insertion/removal to determine the reaction probability. These methods, applied to a simple test system reproduce the expected equilibrium and nonequilibrium response and provide a toolset for understanding more complicated systems where electrostatic interactions, solvent structure, and electrochemical driving play a role in complicating the ‘textbook picture’ of

22


Page 22 of 31

Page 23 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the electrochemical interface.

23



5

Page 24 of 31

Appendix A: Electrostatic potential at the boundary: details of result used in section 2.1

For simplicity, and without loss of generality, we can consider the metal boundaries to be at z = 0 and z = d. For each particle in the active system at position ract = (xi , yi , zi ) and charge qi , the reflected particle is at position rref = (xi , yi , −zi ) and has charge −qi . The electrostatic potential Ψ in the active system in this system is given by eq 9.

Ψ(x, y, z) = zF +

∞ ∞ X X

∞ X N X

qi p

A + (z − (zi + 2dn))2

l=−∞ m=−∞ n=−∞ i=0

(9)

−qi

+p

A + (z − (−zi + 2dn))2

Note that the definition of A is given by equation 10 and involves the x and y contributions to the distance which are the same for active charges and image charges.

A = (x − (xi + lx0 ))2 + (y − (yi + my0 ))2

(10)

The sums over l, m and n account for periodic images in the x, y, and z dimensions respectively and N is the total number of atoms in the active system. x0 and y0 are the length of the periodic cell in the x and y directions respectively. The first term accounts for the applied potential while the first term comes from including the reflection of each particle. By breaking the summation over n into a sum from 1 to infinity, we can rewrite the expression as shown in equation 11, where the n = 0 term is accounted

24


Page 25 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for by ψn=0 =

∞ P

∞ N P P

l=−∞ m=−∞ i=0

√

qi A+(z−zi )2

Ψ(x, y, z) = zF + ψn=0 +

−qi

+√

A+(z−(−zi ))2

. Note that ψn=0 = 0 when z = 0

∞ ∞ ∞ X N X X X

qi p

A + (z − (zi + 2dn))2

l=−∞ m=−∞ n=1 i=0

qi

+p

A + (z − (zi − 2dn))2 −qi

(11)

+p

A + (z − (−zi + 2dn))2 −qi

+p

A + (z − (−zi − 2dn))2

At the z = 0 plane, the first and fourth terms inside the the summation will cancel as will the second and third terms. Thus Ψ = 0 at z = 0 as desired. Since the system is infinitely periodic, we can rewrite equation 9 in terms of a repeat unit spanning z = 0 to z = 2d. In this case, a repeat unit will include a particle at position zi and reflection −zi + 2d. Again splitting the sum over n into two parts gives equation 12. With the shifted frame, the n = 0 ∞ ∞ N P P P −qi √ qi +√ . Note that here ψn=0 = 0 term is now ψn=0 = 2 2 l=−∞ m=−∞ i=0

A+(z−zi )

A+(z−(2d−zi ))

when z = d.

Ψ(x, y, z) = zF + ψn=0 +

∞ ∞ ∞ X N X X X

qi p

l=−∞ m=−∞ n=1 i=0

+p +p +p

A + (z − (zi + 2dn))2 qi A + (z − (zi − 2dn))2 −qi

(12)

A + (z − (−zi + 2d + 2dn))2 −qi A + (z − (−zi + 2d − 2dn))2

Evaluated at z = d, the first and fourth terms inside the summation add to zero as do the middle two terms. This gives a total potential of Ψ = dF = ∆V everywhere along the z = d plane as desired.

25



Acknowledgement This work was supported by the US Department of Energy through the grant number DESC0018094. KD was partially supported through the Toyota Research Institute and the National Science Foundation through the Graduate Research Fellowship Program. AW acknowledges the Research Corporation for Scientific Advancement (Cottrell Scholar).

Supporting Information Available Description of the procedure for selecting the absolute value of Velectrode in Eqs 7 and 8 and details on obtaining the code used to modify LAMMPS.

References (1) P´eraud, J.-P.; Nonaka, A.; Bell, J. B.; Donev, A.; Garcia, A. L. Fluctuation-enhanced electric conductivity in electrolyte solutions. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, 10829–10833. (2) Bentley, C. L.; Kang, M.; Unwin, P. R. Nanoscale surface structure-activity in electrochemistry and electrocatalysis. J. Am. Chem. Soc. 2018, 141, 2179–2193. (3) Bazant, M. Z.; Chu, K. T.; Bayly, B. J. Current-voltage relations for electrochemical thin films. SIAM J. Appl. Math 2005, 65, 1463–1484. (4) Brown, M. A.; Bossa, G. V.; May, S. Emergence of a stern sayer from the sncorporation of hydration interactions into the gouy-chapman model of the electrical double layer. Langmuir 2015, 31, 11477–11483. (5) Pilon, L.; Wang, H.; D’Entremont, A. Recent advances in continuum modeling of interfacial and transport phenomena in electric double layer capacitors. J. Electrochem. Soc. 2015, 162, A5158–A5178. 26


Page 26 of 31

Page 27 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(6) Car, R.; Parrinello, M. Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett. 1985, 55, 2471–2474. (7) Ayala, R.; Sprik, M. Ligand field effects on the aqueous Ru(III)/Ru(II) redox couple from an all-atom density functional theory perspective. J. Chem. Theory Comput. 2006, 2, 1403–1415. (8) Cheng, T.; Xiao, H.; Goddard III, W. A. Full atomistic reaction mechanism with kinetics for CO reduction on Cu(100) from ab initio molecular dynamics free-energy calculations at 298 K. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, 1795–1800. (9) Oberhofer, H.; Blumberger, J. Charge constrained density functional molecular dynamics for simulation of condensed phase electron transfer reactions. J. Chem. Phys. 2009, 131, 064101. (10) Adriaanse, C.; Cheng, J.; Chau, V.; Sulpizi, M.; Vandevondele, J.; Sprik, M. Aqueous redox chemistry and the electronic band structure of liquid water. J. Phys. Chem. Lett. 2012, 3, 3411–3415. (11) Leung, K. Electronic structure modeling of electrochemical reactions at electrode/electrolyte interfaces in lithium ion batteries. J. Phys. Chem. C 2013, 117, 1539–1547. (12) Reed, S. K.; Lanning, O. J.; Madden, P. A. Electrochemical interface between an ionic liquid and a model metallic electrode. J. Chem. Phys. 2007, 126, 084704. (13) Petersen, M. K.; Kumar, R.; White, H. S.; Voth, G. A. A computationally efficient treatment of polarizable electrochemical cells held at a constant potential. J. Phys. Chem. C 2012, 116, 4903–4912. (14) Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M. et al. The ReaxFF reactive 27



force-field: Development, applications and future directions. npj Comput. Mater. 2016, 2, 15011. (15) Ostadhossein, A.; Kim, S.-Y.; Cubuk, E. D.; Qi, Y.; van Duin, A. C. Atomic insight into the lithium storage and diffusion mechanism of SiO2/Al2/O3 electrodes of Li-ion batteries: ReaxFF reactive force field modeling. J. Phys. Chem. A 2016, 120, 21142127. (16) Ouyang, W.; Saven, J. G.; Subotnik, J. E. A surface hopping view of electrochemistry: non-equilibrium electronic transport through an ionic solution with a classical master equation. J. Phys. Chem. C 2015, 119, 20833–20844. (17) Spohr, E. Effect of electrostatic boundary conditions and system size on the interfacial properties of water and aqueous solutions. J. Chem. Phys. 1997, 107, 6342–6348. (18) Widmann, A. H.; Adolf, D. B. A comparison of Ewald summation techniques for planar surfaces. Comput. Phys. Commun. 1997, 107, 167–186. (19) Merlet, C.; P´ean, C.; Rotenberg, B.; Madden, P. A.; Simon, P.; Salanne, M. Simulating supercapacitors: can we model electrodes as constant charge surfaces? J. Phys. Chem. Lett. 2013, 4, 264–268. (20) Siepmann, J. I.; Sprik, M. Influence of surface topology and electrostatic potential on water/electrode systems. J. Chem. Phys. 1995, 102, 511–524. (21) Reed, S. K.; Madden, P. A.; Papadopoulos, A. Electrochemical charge transfer at a metallic electrode: A simulation study. J. Chem. Phys. 2008, 128, 124701. (22) Torrie, G. M.; Valleau, J. P.; Patey, G. N. Electrical double layers. II. Monte Carlo and HNC studies of image effects. J. Chem. Phys. 1982, 76, 4615–4622. (23) Lorrain, P.; Corson, D. Introduction to Electromagnetic Fields and Waves; Freeman: San Francisco, 1962; Chapter 4.

28


Page 28 of 31

Page 29 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(24) Griffiths, D. J. Introduction to Electrodynamics, 3rd ed.; Prentice Hall: Upper Saddle River, New Jersey, 1999. (25) Perram, J. W.; Ratner, M. A. Simulations at conducting interfaces: boundary conditions for electrodes and electrolytes. J. Chem. Phys. 1996, 104, 5174. (26) Plimpton, S. J.; Pollock, R.; Stevens, M. Particle-mesh Ewald and rRESPA for parallel molecular dynamics simulations. Proc of the Eighth SIAM Conference on Parallel Processing for Scientific Computing. Minneapolis, MN, 1997. (27) Toukmaji, A. Y.; Board, J. A. J. Ewald summation techniques in perspective: a survey. Comput. Phys. Commun. 2003, 95, 73–92. (28) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from Algorithms to Applications, 2nd ed.; Academic Press: San Diego, 2002; Chapter 5. (29) Abe, T.; Fukuda, H.; Iriyama, Y.; Ogumi, Z. Solvated Li-ion transfer at interface between graphite and electrolyte. J. Electrochem. Soc. 2004, 151, 1120–1123. (30) Borodin, O.; Bedrov, D. Interfacial structure and dynamics of the lithium alkyl dicarbonate SEI components in contact with the lithium battery electrolyte. J. Phys. Chem. C 2014, 118, 18362–18371. (31) Leung, K.; Leenheer, A. How voltage drops are manifested by lithium ion configurations at interfaces and in thin films on battery electrodes. J. Phys. Chem. C 2015, 119, 10234–10246. (32) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117, 1–19. (33) Magda, J. J.; Tirrell, M.; Davis, H. T. Molecular dynamics of narrow, liquid-filled pores. J. Chem. Phys. 1985, 83, 1888–1901.

29



(34) Bard, A. J.; Faulkner, L. R. Electrochemical Methods: Fundamentals and Applications, 2nd ed.; John Wiley & Sons, Inc.: New York, 2001; Chapter 3. (35) Gupta, A.; Stone, H. A. Electrical double layers: effects of asymmetry in electrolyte valence on steric effects, dielectric decrement, and ion-ion correlations. Langmuir 2018, 34, 11971–11985. (36) Goodwin, Z. A.; Feng, G.; Kornyshev, A. A. Mean-field theory of electrical double layer in ionic liquids with account of short-range correlations. Electrochim. Acta 2017, 225, 190–197.

30


Page 30 of 31

Page 31 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry —

oxidation

— +

+

+ -

+

—

∆V

reduction

31


+

Constant Potential, Electrochemically Active Boundary Conditions for

Recommend Documents