Controlled-Potential Simulation of Elementary Electrochemical

Article Cite This: J. Phys. Chem. C XXXX, XXX, XXX−XXX

pubs.acs.org/JPCC

Controlled-Potential Simulation of Elementary Electrochemical Reactions: Proton Discharge on Metal Surfaces Georg Kastlunger, Per Lindgren, and Andrew A. Peterson* School of Engineering, Brown University, Providence, Rhode Island 02912, United States S Supporting Information *

ABSTRACT: The simulation of electrochemical reaction dynamics from first principles remains challenging, since over the course of an elementary step, an electron is either consumed or produced by the electrode. For example, the hydrogen evolution reaction begins with a simple proton discharge to a metal surface, but with conventional electronic structure methods, the simulated potential, which is manifested as the metal’s workfunction, varies over the course of the simulation as the electron is consumed in the new metal−hydrogen bond. Here, we present a simple approach to allow the direct control of the electrochemical potential via charging of the electrode surface. This is achieved by changing the total number of electrons in the self-consistent cycle, while enforcing charge neutrality through the introduction of a jellium counter charge dispersed in an implicit solvent region above the slab. We observe that the excess electrons localize selectively at the metal’s reactive surface and that the metal workfunction responds nearly linearly to the variation in electronic count. This linear response allows for control of the potential in simulations with a minimal computational penalty compared to standard electronic structure calculations. This scheme can be straightforwardly implemented with common electronic structure calculators (density functional theory in the present work), and we find this method to be compatible with the commonly used computational hydrogen electrode model, which we expect will make it useful in the construction of potential-dependent free-energy diagrams in electrochemistry. We apply this approach to the proton-deposition (Volmer) step on both Au(111) and Pt(111) surfaces and show that we can reliably control the simulated electrode potential and thus assess the potential dependence of the initial, transition, and final states. Our method allows us to directly assess the location along the reaction pathway with the greatest amount of charge transfer, which we find to correspond well with the reaction barrier, indicating this reaction is a concerted proton−electron transfer. Interestingly, we show that the Pt electrode has not only a more favorable equilibrium energy with adsorbed hydrogen but also a lower intrinsic barrier under thermoneutral conditions.

■

mechanisms correctly in first-principles simulations, the ability to control the electrochemical potential is crucial. As a consequence, electrochemical reactions that involve charge transfer at the interface between the electrode and molecular species in the electrolyte are particularly difficult to simulate. A fundamental limitation lies in the size of the simulation unit cells; one must always consider the trade-off between capturing the macroscopic behavior of the system and creating a computationally feasible system. Since charge is transferred to/from the electrode in the course of such a reaction, the electrochemical potential, which can be interpreted as the Fermi level, varies in a canonical simulation. This leads to an artificial situation, where the initial, transition, and final states are all described at different applied potentials. A simple but still theoretically challenging reaction is the electrochemical proton discharge to a surface (H+ + e− + * → H*, where * represents a surface site or adsorbed species),

INTRODUCTION

The development of atomistic simulation tools and the growth in computational power in the past decades have made it possible to study reaction dynamics at solid surfaces from an ab initio perspective.1,2 In this context, density functional theory (DFT)3,4 has proven extraordinarily useful, allowing for calculations of systems consisting of a reasonably high number of atoms. A wide variety of surface/adsorption reactions have been studied and understood at the atomistic level,1,2,5 and in silico discovery of relevant materials and processes is already in practice. Dealing with electrochemical environments, however, is still challenging within conventional electronic structure theory methods that are bound to the canonical ensemble, that is, a fixed number of electrons in the periodic unit cell. In experimental measurements, the potential of the working electrode is controlled by the external circuitry and is subject to factors such as the alignment of the electrochemical potentials of the involved electrodes6 and the structure of the electrode− electrolyte interface.7,8 To capture reaction energetics and © XXXX American Chemical Society

Received: March 13, 2018 Revised: May 16, 2018

A

DOI: 10.1021/acs.jpcc.8b02465 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

it is subject to increasing computational demands imposed by the number of electrons in the system, where the scaling is O(N2e ) − O(N3e ),30 where Ne is the number of electrons (or basis functions) in the system. On the other hand, explicit charging schemes can directly prevent the change in workfunction but demand a strategy for dealing with a net charge in periodic unit cells. Proposed approaches for this include applying a homogeneous background charge and tackling the loss of a field-free region by relating the electrode’s Fermi level to the potential inside a water bulk structure,22 introducing a charged plane (counter electrode) of opposite charge in the unit cell,20,21 introducing a polarizable continuum with infinite dielectric constant that mimics a perfect conductor with a dissolved counter charge,23 and utilizing implicit solvent approaches, where the distribution of the counter charge is defined on the basis of the Debye screening lengths of the electrochemical double layer.25,26,28,29 In the current work, we develop and employ a relatively simple method to explicitly control the reactive workfunction by employing a jellium slab, whose constant-charge background is tempered by an implicit solvent; we will refer to this as the “solvated jellium” method. While similar in spirit to predecessor methods, the current method offers a straightforward approach that can in principle be implemented in any generic DFT code that is able to handle an implicit solvent in the unit cell while keeping a physically reasonable behavior. In the following, we will present the background on this approach as well as its application to describing the proton-deposition reaction in a potential-controlled manner.

which is known as the Volmer reaction in the context of the hydrogen evolution reaction (HER).9,10 Theoretical investigations into this reaction are commonly based on a periodic metallic “slab” that represents the electrode; this slab is infinite (periodic) in the x and y Cartesian directions and has a surface exposed to vacuum between layers in the z direction. Additionally, at least one icelike water layer is added to serve as the proton donor. A common technique is to include a surplus of at least one hydrogen atom in this water layer. In a ground-state electronic structure calculation, this leads to a charge reorganization compared to the separated slab and water systems, resulting in a positively charged acidic water structure and a net negative charge (surplus of electrons) localized on the electrode, as one would expect from chemical intuition. Interestingly, the amount of charge separation is typically in the range of 0.6−0.7,11−13 where it is still under debate if the noninteger value is a consequence of the well-known selfinteraction error of DFT or due to the hybridization of the acidic water network with the metallic surface.12 Since the electrode potential is generally defined as its Fermi energy μ with respect to a reasonable reference potential (i.e., the workfunction), one specific electrochemical potential is treated in such a setup. Over the course of the Volmer reaction, the surplus of electrons localizes on the transferring proton forming an adsorbed hydrogen atom and leading to two uncharged subsystems: the water layer and the slab with adsorbed hydrogen. The charge distribution in the system thus changes drastically, and this reorganization causes the workfunction, that is, the simulated electrode potential, to vary in the course of the reaction. Various methods to simulate electrochemical reactions have been proposed. Perhaps, the most prominent and widely used approach is the computational hydrogen electrode (CHE), in which the potential dependence of the reaction thermodynamics is dealt with in postprocessing steps only.14−17 The method accounts for the chemical potential of a p r o t o n / e l e c t r o n p a i r b a s e d o n t h e eq u i l i b r i u m 1 μ[H+] + μ[e−] = 2 μ[H 2(g)] at 0 VRHE (where RHE is the reversible hydrogen electrode) and corrects the driving force with the deviation of the applied potential from the equilibrium situation. This method has proven to be reliable for a large variety of applications; however, it is only designed to capture differences between stable elementary intermediates and does not offer a means to directly capture reaction dynamics. That is, when charge transfer directly takes place in the dynamics of a simulation (such as in the calculation of a barrier), it is necessary to explicitly account for the electrode potential. Nonetheless, the CHE approach provides a thermodynamic foundation that should still be valid in any approach developed to account for the potential-dependent dynamic systems; we will compare our methods to the CHE model later in this study. Methods directly addressing the change in electrode potential in the course of an electrochemical reaction have been published in the last decade and include extrapolating the change in electrode potential to an infinitely large slab,7,18 applying classical capacitor models,12,13 referencing electrochemical reactions to comparable nonelectrochemical processes,19 and varying the number of electrons in the electrode in the course of the reaction.20−29 All approaches have advantages and drawbacks; for example, the cell-size extrapolation scheme directly addresses the root cause of the workfunction change, being the small unit cell sizes. However,

■

THEORY “Solvated Jellium” Method. The solvated jellium method presented here is based on the introduction of an overall charge-neutral “jellium slab” to the unit cell. Our approach is named after the origin of the term jellium,31 which referred to a “uniform distribution of positive chargewe may call it a ’positive jelly’and a compensating number of electrons”. In the same spirit, we add (spatially unconstrained) electrons to our system that are compensated with a uniform distribution of positive charge in the vacuum region above the slab. This simple addition allows for the number of electrons to be varied continuously in a straightforward manner; of course, this approach also allows the opposite strategy of removing electrons from the system and compensating with a negatively charged jellium region. By itself, such a setup would create an artificial potential profile between the electrode and the jellium counter charge. We circumvent this problem by filling the region above the explicit atomic system with a dielectric continuum, which screens the large electric field and creates a reasonable potential reference inside the fluid.25−29 Applying this method, we consistently observe that the applied charge is located dominantly at the slab surface that faces the explicit water layer, that is, the electrode surface as would be expected experimentally. Therefore, a realistic situation resembling a direct charging of the surface by an applied potential is created, where the back side of the electrode, being an artifact of the simulation, is not affected and only the occupation of the metal surface states in the region of interest is varied. It is worth briefly discussing limitations in the atomistic modeling of electrochemical processes involving charge transfer. Hypothetically, if one had enough computational power, a second explicit slab could be introduced in the unit cell, directly simulating the “counter electrode”; such an approach would B


Article


Figure 1. Electrostatic potential resulting from three different approaches for the description of a charged periodic system, namely, addition of a homogeneous background charge (black), a capacitor-type approach (red), and our approach (blue). In all cases, 0.083 electrons per surface atom have been added to the system. The red block on the right-hand side represents the defined counter electrode analogous to the method of Lozovoi and Alavi.20 The light blue background schematically displays the region, which is filled by the implicit solvent and also contains the jellium background charge, which is depicted by the blue pattern. The implicit solvent and the broad counter charge are the characteristic features of SJM compared to the other two methods.

overlap between the jellium slab and electrons. However, the large dipole between the slab and the jellium background would create a large electric field, which would, in principle, lead to a polarization of the electron density. As a consequence, an adsorbed water layer would experience the highest electric field, which would lead to an artificial charging of the bilayer. A realistic potential profile can, however, be created by immersing the jellium background in an implicit solvent, changing the role from a counter electrode in vacuum to mimicking the electrochemical double layer. The charging scheme in the solvated jellium method (SJM) is based on the same foundation as jellium, being a uniform constant background charge region and electrons/holes, which distribute based on the system’s electrostatics. Since, however, the unit cell also contains an atomic slab/surface system, the added electrons/holes do not create a homogeneous electron gas but occupy the Kohn−Sham states resulting from the SCF cycle. Finally, immersing the counter charge region in a dielectric continuum characterized by the permittivity of the solvent is based on its role as a simplified electrochemical double layer, which consists of ions drawn to the surface inside the electrolyte. In the current work, we based the implicit solvent on the effective potential cavity method as proposed by Held and Walter.33 In principle, any implicit solvent model based on explicitly solving the electrostatic problem over real space is applicable, however. The basis of such approaches is the (generalized) Poisson equation

require a large increase in computational efforts due to the vast increase in the number of electrons in the system and the scaling of electronic structure codes. If all reactions at our working electrode were exactly balanced with a reaction at the counter electrode, we would maintain charge neutrality in the system, but it would not be straightforward to separate the energetics of the elementary processes, that is, to correctly assign a reaction barrier energy to the process at the working electrode. Schemes would still be necessary to control the potential in such a construction, for example, by introducing an artificial dipole in the artificial vacuum region between the back sides of the two slabs. In other words, there would still be a large amount of art and arbitrariness in the setup of such simulations, and in practice, it is much simpler and more computationally efficient to focus on the half-reactions taking place at a single electrode. Our solvated jellium approach represents a viable trade-off in this line of thought. Here, we give details on how the two components of this method are applied in practice. A more detailed description of the algorithm and a workflow diagram, as well as a performance benchmark example, can be found in the Supporting Information (SI). The variation in electrons (and therefore the electrical potential) is achieved by introducing a jellium slab in the vacuum region, which exists only on the reactive side. A dipole correction32 is located between the jellium counter charge and the reverse of the periodic slab to isolate this counter charge (and other field effects) on the reactive side. While the localized static background charge region of this jellium is fixed in the unit cell with respect to its charge magnitude and size, the charge with opposite sign is added to or subtracted from the electron reservoir of the atomic system and its localization is optimized in the course of the electronic structure calculation (e.g., the self-consistent field (SCF) cycle in DFT). Jellium itself is a charge-neutral entity, and its inclusion in a periodic unit cell does not cause divergence. Thus, after convergence, the integral of the localized background charge is of equal magnitude but opposite sign to the excess charge of the explicit system. In the current implementation, the lower limit of the jellium has been chosen to be at least 3 Å (two van der Waals radii of oxygen) above the explicit system, where the electron density has fallen off to negligible values, that is, there is no spatial

∇(ϵ(r)∇Φ(r)) = −4πρ(r)

(1)

where ρ(r) represents the charge density of the system, Φ(r) is the resulting electrostatic potential, and ϵ(r) is a smooth permittivity function describing the solute’s cavity and accounting for the screening due to a solvent, which is represented on the same grid. In SJM, 1 is extended to include the charge density introduced into the unit cell by the jellium slab. Hence, ρ(r) is changed into ρ(r) = ρexplicit (r) + ρjellium (r)

(2)

with ρexplicit(r) containing the charge densities from the explicit system (with both electrons and ions, giving a net charge of C


Article

The Journal of Physical Chemistry C zero) and ρjellium(r) being the charge density introduced by the added jellium (containing added electrons/holes and the “background slab”, again with a net zero charge). In the slab geometries common in catalyst models, the simulated electrode has a back side that is also exposed to vacuum, but is meant to mimic an infinitely deep bulk system. Since we are addressing the charging of the electrode− electrolyte interface, it is essential that the structure and electrostatics at the back side of the slab do not change when the chemical potential is varied. Therefore, we exclude the solvent from the vacuum on the back side of the slab; combined with the dipole correction mentioned earlier, we will show that this results in a constant potential profile in this region. In Figure 1, we show comparisons of the electrostatic potentials that result from two previously applied means of adding charge to the periodic system, namely, a capacitor-type approach (red) and a homogeneous background charge (black), and the solvated jellium approach of this work (blue). The capacitor approach consists of defining a charge layer in the unit cell, which neutralizes the net charge created by the addition of electrons or holes in the SCF cycle, as first described by Lozovoi and Alavi.20 The region between the counter slab and the explicit atomic system is treated as vacuum (ϵrel = 1). Since the narrow counter charge slab has to incorporate the whole neutralizing charge, a high charge density is created therein. This charge density combined with the absence of a screening medium in the region between the atomic system and the neutralizing slab creates an artificially high electric field. As a consequence, introduced electron/holes in the SCF cycle localize on top of the explicit water bilayer instead of the metallic electrode, where the latter is the desired behavior when simulating an electrochemical half-cell system. A homogeneous background charge scheme does not suffer from this limitation. However, the uniformly distributed background charge does not provide a region where the electric field is screened entirely. Hence, this approach lacks a reasonable reference for the chemical potential. Taylor et al. addressed this issue and developed a methodology known as the doublereference method.22 Although this approach leads to a reasonable reference, two limitations are still apparent: first, it demands rather large simulations because the interior must be filled with water, and second, there will always be a spurious interaction of the electron density inside and behind the metal slab with the homogeneous background charge. In the solvated jellium method, the implicit solvent reduces the electric field on top of the explicit water, thus avoiding artificial charging of the explicit water. Further, the immersed jellium background completely screens the electric field and provides a reference for the chemical potential. In Figure 2, we show the electrostatic potential difference between charged and neutral unit cells with the solvated jellium approach. The electrostatic potential profile is in good qualitative agreement with the Gouy−Chapman−Stern model;34 the innermost region of the double layer has a large, constant electric field (Stern layer), whereas the electric field in the slipping plane (diffuse layer) decays exponentially to a field-free region. This field-free region corresponds to the electrolyte bulk and serves as a reference for the chemical potential. Therefore, this relatively simple charge neutralization scheme produces realistic double-layer potential profiles. In all of the following cases, the absolute electrode potential Φe with respect to the solvent bulk has been determined as

Figure 2. Schematic view of the unit cell setup in an SJM calculation (top) and electrostatic potential with respect to the potential where no extra charge has been added to the unit cell (bottom).

Φe = Φw −

μ e

(3)

where μ is the metal Fermi level, e is the elementary charge, and Φw represents the electrostatic potential deep in the implicit solvent, where the net charge has been screened completely. Both quantities Φw and μ are subject to the same arbitrary internal reference chosen for the electrostatic potential, making Φe independent of this choice and dependent on factors such as metallic species, surface charge, thickness of the electrochemical double layer, and the dielectric constant of the solvent. That is, eΦe is the workfunction measured on the top side of the slab.

■

RESULTS AND DISCUSSION

Controlling the Potential via Selective Charging. By varying the amount of charge and observing the change in potential by 3, it is straightforward to extract a double-layer capacitance within this model. These double-layer capacitances, calculated without explicit water molecules, were calculated as 18 and 22 μF/cm2 for Au(111) and Pt(111), respectively (for further details, see the Supporting Information). This is in reasonable agreement with the experimental values that show a wide range between 10 and 40 μF/cm2 for Au(111)35−37 and between 9 and 38 μF/cm237−40 for Pt(111) and the computational result based on the application of JDFT,25,29 where a value of 19 μF/cm2 for Pt(111) has been reported. It should be emphasized, however, that the double-layer capacitance is highly dependent on the electrolyte in the experimental setup, where the value can differ by a factor of 2.36 When a water bilayer is adsorbed on the metal electrodes, an implicit surface dipole is introduced, influencing the workfunction and therefore the resulting electrode potential Φe.20 As shown in Table 1, the H-down bilayer structure increases the potential at zero applied charge on the electrode, indicated as Φe,0. This is primarily a consequence of the polar O−H bonds laying perpendicular to the electrode surface. The implicit dipole of this bond leads to a partial negatively charged oxygen pointing toward the vacuum (implicit solvent) above the electrode, thus increasing Φw relative to the bare surface situation arising from the electrostatic penalty for an electron to be placed there. In the H-up geometry, the effect is opposite, D


Article


potential with respect to the known potential of a reference electrode. Translating this into an atomistic description requires that a change in potential should be caused by a change in the electrode’s Fermi level with respect to a known reference potential, which by the present methodology is manifested via a variation in the electrode’s surface charge (this can be considered as an interchange of action and reaction compared to experimental conditions). However, the distribution of charge within the system is set by the ground-state electronic configuration, that is, in the case of DFT, by the self-consistent field cycle, and any added charge does not necessarily localize at the electrode surface. In the present study, we keep the geometry of the water bilayer on top of the electrode constant (at the relaxed geometry of the uncharged situation) as we vary the charge to avoid the need to account for the change in surface dipole due to a reorientation of the water bilayer; this constraint will be lifted when we examine reactive systems later in this work. The amount of charge that localizes on the surface of the slab versus within the water layer can be assigned (somewhat crudely) via a Bader analysis.42 This analysis is shown in the lower left panel of Figure 3. With this analysis, we can clearly see a linear, charge-selective region for charges of low and intermediate magnitudes, highlighting the physical region of allowed charging. Next, we relate the number of electrons in the system to the electrode potential, which is related to the electronic chemical potential (μ) by 3. The top panel of Figure 3 shows the electrostatic potential referenced to μ for the H-down structure,

Table 1. Calculated Potentials at Zero Surface Charge (Φe,0) and Work Functions (Γ) on the Back Side of Au(111) and Pt(111) for the Bare Slabs in an Implicit Solvent and Including the Water Bilayers system

interface

Φe,0, V

Γ, eV

Au(111)

bare H-down H-up bare H-down H-up

4.69 5.99 4.54 5.04 6.37 4.43

5.20 5.24 5.23 5.74 5.74 5.73

Pt(111)

since the positively charged H stabilizes an electron in near-field vacuum, hence reducing Φw and by extension Φe,0. As mentioned above, the back side of the electrode, which represents the transition into the semi-infinite electrode crystal, is not altered in terms of occupation or energetics of its surface states. This region is solvent-free, and thus, its potential Φvac serves as a measure of the vacuum workfunction Γ defined as

Γ = e Φvac − μ

(4)

As seen in Table 1, these calculated work functions are essentially unaffected by the charging scheme and agree well with previous computational studies of uncharged systems.41 Next, we study where the applied charge localizes within this system, i.e., on the electrode surface, in the electrode bulk, or on the water layer. In an experimental setup, the applied voltage is established via the circuitry which sets the electrode’s

Figure 3. Change in electrode potential (Φe) and localization of the excess charge via the present method. Top: electrostatic potential referenced to the system’s Fermi level (μ) of the H-down structure as it was used for determining Φe in the lower right panel. Steps of 0.033 e/atom (0.2 e/cell) are shown in both charging directions up to a charge of 0.167 e/atom (1 e/cell). Bottom: resulting metal electrode surface charge (left) and electrode potential (right) at varying applied charge. The H-down and H-up bilayer situations are depicted in black and red, respectively. E


Article


Figure 4. Density of states projected onto the pz states of the water bilayer adsorbed on Au(111) for the H-down water structure.

Figure 5. Left panel: electron density change with respect to the uncharged electrode surface when 0.2 (red) and 1.0 (blue) electrons have been added (top panel) or extracted (bottom panel) to/from an electrode−H2O down system consisting of six surface metal atoms, corresponding to 0.03 and 0.17 electrons per surface atom. The dotted black curves in the top and bottom panels correspond to the electron density created by the water bilayer LUMO and HOMO, as determined by subdiagonalization of the water subspace in a calculation containing the full slab + water system with an LCAO dzp basis set implemented in GPAW.46 Right panel: integration over the electron density difference (in electrons per surface atom) shown in the left panel.

as they were used for determining Φe. As can clearly be seen, the workfunction on the back side of the electrode stays constant at a value of approximately 5.2 eV as it has also been determined for all of the systems shown in Table 1, while the effect of the applied charge manifests itself in the potential in the metal−water interface region. In fact, this asymmetry in charging of the metal slab is not a consequence of adding an explicit water layer on top of the metal slab, but is rather due to the created electric field on the reactive side, as shown in more detail in the Supporting Information. Additionally, the electrostatic potentials related to all of the curves shown in Figure 3 as well as in Figure 6 and a more thorough analysis addressing the asymmetry in the charge localization on the metal slab on the basis of the density of states are provided in the Supporting Information. In the lower right panel of Figure 3, we show the response of this potential, Φe, to the applied charge on the Au(111) surface. An approximately linear relationship between the two is apparent over a wide range of several volts, and the standard hydrogen electrode potential (0 VSHE), located at 4.44−4.85 V,43,44 falls nicely into this linear region. The power of this linear response should be emphasized; in practice, this means that potential control by adjusting the added charge is straightforward, as adjustments to the system potential will typically take only one extra electronic structure calculation if an estimate of the slope of this curve is known. Additionally, using the Kohn−Sham states from a previous run as a starting guess for the self-consistent field cycle reduces the numerical

effort significantly. This minimizes the computational cost of this method, which in a worst case would carry a penalty of 3× a normal electronic structure calculation, but in practice, it is much closer to 1×, as we describe later. Before moving to a description of proton deposition with this model, we first turn to the issue of the nonlinear regions encountered in Figure 3. For applied charges beyond the positive charge limit, the extra charge is shared between the electrode’s surface and the water subsystem. This is a physically expected behavior, since the Fermi edge of the electrode is varied within the highest occupied molecular orbital (HOMO)−lowest unoccupied molecular orbital (LUMO) gap of the adsorbed bilayer. Once a resonance between the Fermi level and the frontier orbitals of the adsorbed water is created, any additional noninteger charge will be added to both the metal and the pinned water orbitals.45 The shift in the range of the linear region between the two H2O orientations indicates that the alignment of the H2O frontier orbitals changes due to the orientation of its implicit dipole, where the H-down structure frontier orbitals are shifted to higher eigenenergies with respect to μ than their H-up counterparts. To further investigate this, we examined the position of the molecular eigenstates relative to μ, at various levels of applied charge. Figure 4 shows the density of states projected onto the atomic pz eigenstates of the bilayer oxygens (PDOS) relative to μ above the Au(111) surface. (A projection of this kind is reasonable for studying the resonance of the water HOMO with μ, since the molecular HOMO has been determined as a F


Article


Figure 6. Top panels: electrode surface charge with respect to the applied charge. Bottom panels: electrochemical potential with respect to the applied charge. The left panel relates to the Au(111) surface, whereas the right panel shows the results for Pt(111) electrode. In both panels, the red full circles represent H+ + * and the black open circles depict H* geometries. The gray-shaded areas represent the charges (potentials), where the results deviate from a selective charging of the metal electrode. The horizontal dashed line represents the absolute potential of the standard hydrogen electrode at pH 0.

π-MO consisting solely of these atomic eigenstates.) In the situation where no extra electrons have been included, the electrode’s chemical potential lies inside the HOMO−LUMO gap of the water bilayer, which explains why addition or subtraction of electrons leads to a selective filling or emptying of the electrode’s surface states, respectively. Reducing the number of electrons in the unit cell leads to a positive charge on the surface and a reduction of the metal Fermi energy. As a consequence, the Fermi energy approaches the molecular HOMO until they are in resonance. From this point on, reducing the electrons in the unit cell always leads to an extraction of electrons, that is, positive charging, of both the electrode and the adsorbed bilayer due to the pinned molecular HOMO. This is in agreement with the discussion of Björketun and colleagues,45 who correctly suggested that only a certain range of charges can be selectively applied in electrochemical reaction simulations. Figure 5 shows the change in electron density (Δnz) with respect to the uncharged slab and the traces of an integration over it, where the two cases of extracting 0.03 and 0.12 e− per electrode surface atom are depicted in the bottom panel. In agreement with the findings from the PDOS analysis in Figure 4, the electron density difference at low positive charges exhibits partial extraction of electrons from the metallic surface states. Additionally, the electron density in the water bilayer is polarized due to the newly created positive charge. Hence, the electron density increases below the oxygen center of the water layer and decreases above it. When higher positive charges are applied, the situation changes, since the resonance region described in the last paragraph is reached. In this case, the change in electron density in the water bilayer is not only due to polarization, but also because of direct charging. This phenomenon manifests itself as the negative double peak that has the same shape as the molecular HOMO. When electrons are added to the surface system, Δnz is comparable to the positive charging within a certain range, as shown in the top panel of Figure 5. However, once a certain threshold is surpassed, the excess electrons cannot be held by the atomic cores anymore and are, therefore, pulled out into the solvent + background region of the unit cell. This behavior, obviously, does not result from a pinning of the water layer LUMO, but is rather related to a failing of the charging scheme due to the high applied charge and low slab workfunction.

As it will be shown later, this, however, does not impede the applicability of the model for the study of the HER Volmer reaction, since the electrode potential in this case is far off the potentials applied in an experimental situation. Finally, it should be noted that the limits of the region of selective charging is subject to the well-known gap problem of DFT,47,48 when limiting the description to conventional exchange−correlation functionals. (Björketun et al.45 showed that this can be partially alleviated by employing functionals with a Hubbard-like correction.) Therefore, these results should not be taken as providing a reliable range, where the electrode can be charged selectively, but rather as a tool where the applied charge selectively varies the electrode’s Fermi level unless the frontier orbitals of the water bilayer are pinned to it. This makes the charging scheme a viable tool for selfconsistent, potential-dependent descriptions of HER, which are commonly performed at negative potentials near 0 VSHE, where the Fermi level is far from the frontier orbitals. Modeling Proton Deposition. The electrochemical proton-deposition (Volmer) reaction, H+ + e− + * → H*, is arguably the simplest and most fundamental reaction in electrocatalysis, but has proven to be difficult to describe via electronic structure calculations. Here, we employ our solvated jellium approach to describe this reaction on Au(111) and Pt(111) surfaces, two commonly used metals in the study of HER. Since HER shows a significantly higher rate under acidic reaction conditions,49−52 in the following, we will study the system where the initial state contains an additional hydrogen in the water bilayer so that the proton donor is H3O+, denoted H+ + *, and the final state, where the hydrogen is adsorbed on the metal surface (H*). We first describe the physical response in the system without added charge (but with the added implicit solvent). In the H+ + * state, 0.64 (Au) and 0.63 (Pt) excess electrons are localized on the metal surfaces, as determined by a Bader charge analysis; see our earlier discussion of the noninteger nature of this charge separation. In both cases, the positive counter charge localizes on the H5O2+ dimer in the water layer. As a consequence, the H-up configuration spontaneously rotates to evade the electrostatic penalty arising from the interaction between the electrode’s surface charge and the lone pair on the oxygens. The H-down structure, on the other hand, is stable, where the H5O2+ dimer is created from one water molecule G


Article


Figure 7. Grand potential energy relative to the initial state for the Volmer reaction under acidic conditions. The voltage values given on the righthand side are shown relative to 4.4 V, which we approximate as 0 VSHE. IS, TS, and FS are the initial, transition, and final states, respectively.

HOMO is located comparatively near μ, and the addition of the H-atom and along with it the positive charge stabilize it and shift it further away from μ. As a consequence, the positive charge, which has to be applied to establish a resonance, is increased in the latter case, where μ is lower and accordingly Φe is higher. The deviation of the negative applied charges on the other side is not related to the resonance of μ and the water LUMO, but is rather related to the low Φe in this case, as will be discussed in the following paragraph. Given this, it is straightforward to vary the potential between about 2 and 7 V for both Au and Pt, but at potentials beyond this range, the excess electrons tend to distribute between the HOMO of the water and the metal surface or localize in the implicit solvent above the slab, leading to a less pronounced change and an artificial rise in the electrochemical potential once a Φe value of ≈1 V is reached, respectively. As earlier, the standard hydrogen electrode potential of 4.44−4.85 V43,44 falls nicely into this linear region. Now, we turn toward simulating the kinetics of the reaction at constant potential by calculating the reaction barrier along the minimum energy pathway (MEP) between these two endstate geometries of the Volmer step. We performed climbingimage nudged elastic band (CI-NEB) calculations,53,54 where the initial and final states at each single defined potential have been connected with five interior images. We chose a final state consisting of an adsorbed proton in an ontop position located beneath the former H5O+2 dimer. This structure, although unstable in the canonical scheme, has been identified as the first stable site on the minimum energy pathway of the Volmer reaction and is stable if the electrode’s surface is negatively charged (corresponding to potentials below about 0 VSHE, where HER takes place). In the shown results, we focus on the calculation of the socalled “grand” potential energy, Ω, in which the DFT-calculated energy (EDFT) containing Ne extra electrons is corrected for the energy of the extra Ne electrons such that comparisons can be made between the energies of states with varying numbers of electrons along a reaction pathway

pointing down and one lying parallel to the electrode’s surface, with the hexagonal structure persisting. Compared to the adsorption of pure water, the slab−oxygen distances decrease by about 0.5 Å, as a consequence of the electrostatic attraction between the two subsystems and the limited size of the periodic cell. We are aware that unit cell effects can lead to artificial behavior in the simulation of electrode reactions on periodic surfaces. For the proton-deposition process, we expect this issue to be minor. However, for reactions involving larger adsorbates, such as CO2 or COOH, we expect this issue to be more severe and we plan to study this effect in future work. The same behavior for both the H-up and H-down structure was observed with the implicit solvent absent, and in this case, the charge separations (0.62 and 0.60 for Au and Pt, respectively) were in exact agreement with the literature.11 In contrast, the H* setup did not result in a notable charge transfer between the electrode and the water bilayer for all of the studied adsorption sites of H. The most stable adsorption sites in the uncharged situation have been determined as the face-centered cubic (fcc) and the hexagonal close-packed hollow positions for the H-down and H-up bilayers, respectively. In the adsorption in the ontop bonding site, where the hydrogen is located beneath a flat lying H2O in the H-down configuration, the hydrogen spontaneously desorbs at positive potentials with respect to 0 VSHE and the structure relaxes into the H+ + * configuration described in the last paragraph. This bonding site, however, cannot be discarded and will be studied further below, since it represents the most reasonable intermediate adsorption site in the course of HER. We can achieve potential control of this reaction dynamically, and we show the relationship between the applied charge and the resulting potential for the initial and final states in Figure 6. Similar to the case without the added proton (Figure 3), a region exists where the charging of the metal surface is selective and the potential response is linear. Matching physical intuition, in the H+ + * geometry, there is a surplus of surface electrons and a correspondingly lower potential at any applied charge within the linear-response region. That is, even at zero applied charge, the H+ + * geometry is at a significantly lower potential than the H* geometry and electrons must be extracted from the former to match the potentials of the two states. Interestingly, the deviation from the linear, selective region at positive applied charges happens at the same surface charge for both H+ + * and final H*. This behavior can be explained by the energetic localization of the water HOMO with respect to the slab Fermi edge in the two cases. While in the H* system, the molecular

Ω = E DFT + Nee Φe

(5)

We note that the DFT energy also contains solvation effects of the implicit solvent. In the current work, we focus only on the potential energy and neglect nuclear free-energy corrections to these terms. For simplicity in the current work, we take the absolute potential of the standard hydrogen electrode at pH 0 to be 4.4 V, although as noted before, there is experimental H


Article


Figure 8. Minimum energy pathways (black) and electrons added or subtracted along the bands (gray) for the Volmer step on Au(111) and Pt(111) at 0 VSHE.

Figure 9. Applied charge to the atomic systems at given potentials. The results of IS and FS relaxed at the respective potentials are shown as circles and squares, respectively. The charge on TS, resulting from constant potential CI-NEB calculations, are shown by triangles. The dashed lines in the cases of IS and FS correspond to a linear fit of the data, while for TS, a quadratic fit was applied.

barrier energies are higher on Au than on Pt. These energetic differences manifest themselves in a later transition state on Au(111), where the reaction is almost barrierless, i.e., the Butler−Volmer56 symmetry factor is close to unity. Interestingly, we find that the highest rate of electron transfer occurs at roughly the transition state in both cases, with little electron transfer near the end states. These results indicate that the model captures the electrochemical reaction, which can clearly be thought of as a concerted proton−electron transfer, and that any subsequent surface diffusion step exhibits low or moderate potential dependence. As a final benchmark, we now compare these results to the most prominent model of accounting for potential dependence in electrochemistry, namely, the CHE model. CHE predicts that for a reaction containing the transfer of an electron, a change in the potential of 1 V would lead to a change in the reaction free energy by 1 eV (with the free-energy change following NeeΔU, with ΔU being the potential change). In this work, we found that ΔΩ only changes by 0.6 eV for Au and 0.50 eV for Pt over a range of 1 V. Where does this apparent deviation come from? An explanation can be found in Figure 9. This shows the amount of charge required to achieve a certain potential in both the initial and final states. Interestingly, the difference in these curves gives the number of electrons that must be “injected” into the system to hold the potential constant over the course of the reaction, that is, the number of electrons transferred over the course of the reaction. This is not found to be 1 e−, but 0.59 and 0.49 e− for Au and Pt,

uncertainty in this choice and it is not guaranteed to be consistent with the electronic structure calculations. In principle, an internal reference for 0 VSHE could be defined,7 which we leave as a task for future studies. The resulting grand potential energy diagrams are shown in Figure 7. First, we note one large difference between Au and Pt: although the overall Volmer step is slightly downhill at 0 VSHE on Pt, it is still strongly uphill on Au at this potential, in line with understandings from earlier studies via a thermodynamic (CHE) model,5,55 that is, Pt binds the H* intermediate with a near-optimal strength, while Au is too noble. As the applied potential (Φe) becomes more negative, the reaction becomes more downhill in a linear fashion with respect to potential, due to the difference in the response to potential of the initial and final states. Similarly, as Φe becomes more negative, the reaction barrier decreases, but not as strongly as the reaction energetics, as would be expected for a reaction intermediate in which incomplete charge transfer has occurred, and not in a linear fashion. Interestingly, according to this analysis, Pt is a better catalyst from the standpoint of not only the H* intermediate but also the transition-state energy: even when the elementary step is thermoneutral, Pt exhibits a lower barrier (200 meV). Figure 8 illustrates the location of the reaction barrier and charge transfer across the reaction coordinate. These examples were chosen at the equilibrium potential (0 VSHE), but the minimum energy pathways of other electrode potentials are included in the Supporting Information. The reaction and I


Article


With this, we analyzed the Volmer (proton-deposition) reaction on Au and Pt fcc(111) surfaces and found the potential dependence of this elementary reaction. We were able to observe a systematic decrease in the reaction barrier with applied (negative) potential and saw that the susceptibility of the barrier to potential is less than that of the stable intermediate H*, as expected. Additionally, we observed that the thermoneutral barrier, the barrier at the potential when the elementary step is thermoneutral, is intrinsically larger on the Au surface than on the Pt one, further suggesting the reason for the superior activity of Pt in this reaction.

respectively, which is roughly consistent with the Bader analysis of the charge separation described earlier. This behavior and the partial charge transfer indicate that the chosen initial state is already a (stable) intermediate state of the Volmer reaction of HER and that the migration of the proton through a water network would have to be regarded to completely describe the reaction. On the basis of these empirical results and a more detailed mathematical comparison, which is given in the SI, it can be confidently concluded that the presented methodology is in agreement with the predictions of the CHE model, while also directly accounting for the potential-dependent charge transfer, which gives access to studies of reaction kinetics. Figure 9 also contains the respective charge at the transition state as determined from the MEPs resulting from CI-NEB calculations (shown in the Supporting Information). As can be clearly seen, the charge transfer at the transition state is only partial with respect to the total transferred charge. In general, the transition state is closer to FS for the Au slab compared to Pt, which comes from the fact that it is located later on the MEP. Additionally, the charge on TS at varying potentials does not follow a linear trend as IS and TS do. This is a consequence of the transition state moving along the MEP, where it approaches IS and FS for cathodic and anodic potentials, ultimately leading to activationless and barrierless reactions, respectively. This indicates a change in the Butler−Volmer symmetry factor at varying potentials, which will be studied further in subsequent work.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b02465. Computational details; comparison with the computational hydrogen electrode (CHE) model; workflow in the SJM model; computational performance evaluation; constant potential climbing-image nudged elastic band minimum energy pathways of the Volmer step on Au(111) and Pt(111); workfunction vs applied charge for different water geometries; density of states description of charge localization; charge localization without explicit water; implicit solvent model; and double-layer capacitance studies (PDF)

■

■

CONCLUSIONS In an electrochemical reaction, the effect of the applied potential is to tune the chemical potential of the electrode relative to that of the chemicals with which it reacts. It is essential that this behavior is reflected in an atomistic simulation of such a process. We have presented a simple methodology that introduces a jellium slab immersed in an implicit solvent, which allows us to apply a net charge on the explicit system and thus controls the potential. Due to the screening of the created electric field, the net charge localizes on the electrode surface facing the adsorbed icelike water layer. Therefore, the electrode’s Fermi level is selectively altered with respect to the eigenlevels of the adsorbed water, whereas the back side of the electrode, which mimics the transition into the semi-infinite electrode bulk, is not altered in terms of occupation. Importantly, we observe a linear response (within reasonable physical limits) of the potential to the applied charge. This linear response allows us to quickly adjust the potential to a desired value with minimal computational effort not to exceed a factor of 3, which in practice is much less due to the relatively constant slope of the response as well as the ability to use the previous Kohn−Sham states as an initial guess to the selfconsistent field cycle. We found that the applied charge selectively localizes on the electrode’s surface within a certain range, which has been identified as the potential region in which the metal Fermi level is located within the adsorbate’s HOMO−LUMO gap. At potentials beyond this region, any applied charge (or potential) will be distributed between the electrode and the H2O bilayer, due to resonance and pinning of the adsorbate’s frontier orbitals to μ. The equilibrium potential of the standard hydrogen electrode falls nicely into the linear-response region between added charge and potential.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Andrew A. Peterson: 0000-0003-2855-9482 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS G.K. and P.L. were funded by the Office of Naval research under award N00014-16-1-2355 during this work. The authors thank Tim Mueller, Yogesh Surendranath, Tim Lian, Leanne Chen, Ali Alavi, and Jens Nørskov for fruitful discussions, which finally led to the described model.

■

REFERENCES

(1) Gross, A.; Schnur, S. Catalysis in Electrochemistry; John Wiley & Sons, Inc., 2011; Chapter 5, pp 165−196. (2) Nørkov, J. K.; Abild-Pedersen, F.; Studt, F.; Bligaard, T. Proc. Natl. Acad. Sci. 2011, 108, 937−943. (3) Hohenberg, P.; Kohn, W. Phys. Rev. 1964, 136, B864−B871. (4) Kohn, W.; Sham, L. J. Phys. Rev. 1965, 140, A1133−A1138. (5) Nørskov, J. K.; Bligaard, T.; Rossmeisl, J.; Christensen, C. H. Nat. Chem. 2009, 1, 37−46. (6) Tripkovic, V.; Björketun, M. E.; Skúlason, E.; Rossmeisl, J. Phys. Rev. B 2011, 84, No. 115452. (7) Rossmeisl, J.; Skúlason, E.; Björketun, M. E.; Tripkovic, V.; Nørskov, J. K. Chem. Phys. Lett. 2008, 466, 68−71. (8) Rossmeisl, J.; Chan, K.; Ahmed, R.; Tripkovic, V.; Bjorketun, M. E. Phys. Chem. Chem. Phys. 2013, 15, 10321−10325. (9) Safizadeh, F.; Ghali, E.; Houlachi, G. Int. J. Hydrogen Energy 2015, 40, 256−274. (10) Zheng, Y.; Jiao, Y.; Jaroniec, M.; Qiao, S. Z. Angew. Chem., Int. Ed. 2015, 54, 52−65. (11) Björketun, M. E.; Tripkovic, V.; Skúlason, E.; Rossmeisl, J. Catal. Today 2013, 202, 168−174. J


Article

The Journal of Physical Chemistry C (12) Chan, K.; Nørskov, J. K. J. Phys. Chem. Lett. 2015, 6, 2663− 2668. (13) Chan, K.; Nørskov, J. K. J. Phys. Chem. Lett. 2016, 7, 1686− 1690. (14) Nørskov, J. K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J. R.; Bligaard, T.; Jónsson, H. J. Phys. Chem. B 2004, 108, 17886−17892. (15) Peterson, A. A.; Nørskov, J. K. J. Phys. Chem. Lett. 2012, 3, 251− 258. (16) Rossmeisl, J.; Qu, Z.-W.; Zhu, H.; Kroes, G.-J.; Nørskov, J. K. J. Electroanal. Chem. 2007, 607, 83−89. (17) Keith, J. A.; Jerkiewicz, G.; Jacob, T. ChemPhysChem 2010, 11, 2779−2794. (18) Skúlason, E.; Tripkovic, V.; Björketun, M. E.; Gudmundsdóttir, S.; Karlberg, G.; Rossmeisl, J.; Bligaard, T.; Jónsson, H.; Nørskov, J. K. J. Phys. Chem. C 2010, 114, 18182−18197. (19) Akhade, S. A.; Bernstein, N. J.; Esopi, M. R.; Regula, M. J.; Janik, M. J. Catal. Today 2017, 288, 63−73. (20) Lozovoi, A.; Alavi, A.; Kohanoff, J.; Lynden-Bell, R. M. J. Chem. Phys. 2001, 115, 1661−1669. (21) Lozovoi, A. Y.; Alavi, A. Phys. Rev. B 2003, 68, No. 245416. (22) Taylor, C. D.; Wasileski, S. A.; Filhol, J.-S.; Neurock, M. Phys. Rev. B 2006, 73, No. 165402. (23) Otani, M.; Sugino, O. Phys. Rev. B 2006, 73, No. 115407. (24) Jinnouchi, R.; Anderson, A. B. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, No. 245417. (25) Letchworth-Weaver, K.; Arias, T. A. Phys. Rev. B 2012, 86, 75140. (26) Mathew, K.; Sundararaman, R.; Letchworth-Weaver, K.; Arias, T. A.; Hennig, R. G. J. Chem. Phys. 2014, 140, No. 084106. (27) Sakong, S.; Naderian, M.; Mathew, K.; Hennig, R. G.; Groß, A. J. Chem. Phys. 2015, 142, No. 234107. (28) Goodpaster, J. D.; Bell, A. T.; Head-Gordon, M. J. Phys. Chem. Lett. 2016, 7, 1471−1477. (29) Sundararaman, R.; Schwarz, K. J. Chem. Phys. 2017, 146, No. 084111. (30) Goedecker, S. Rev. Mod. Phys. 1999, 71, 1085−1123. (31) Ewald, P. P.; Juretschke, H. In Structure and Properties of Solid Surfaces; Gomer, R., Smith, C. S., Eds.; University of Chicago Press, 1953; pp 82−120. (32) Bengtsson, L. Phys. Rev. B 1999, 59, 12301−12304. (33) Held, A.; Walter, M. J. Chem. Phys. 2014, 141, No. 174108. (34) Wang, H.; Pilon, L. J. Phys. Chem. C 2011, 115, 16711−16719. (35) Kolb, D. M.; Schneider, J. Electrochim. Acta 1986, 31, 929−936. (36) Vasiljevic, N.; Trimble, T.; Dimitrov, N.; Sieradzki, K. Langmuir 2004, 20, 6639−6643. (37) Bartlett, P. N.; Cook, D. A. J. Electroanal. Chem. 2015, 746, 18− 24. (38) Pajkossy, T.; Kolb, D. M. Electrochim. Acta 2001, 46, 3063− 3071. (39) Hou, Y.; Aoki, K. J.; Chen, J.; Nishiumi, T. Univers. J. Chem. 2013, 1, 162−169. (40) Hou, Y.; Aoki, K. J.; Chen, J.; Nishiumi, T. J. Phys. Chem. C 2014, 118, 10153−10158. (41) Singh-Miller, N. E.; Marzari, N. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, No. 235407. (42) Bader, R. F. W. Acc. Chem. Res. 1985, 18, 9−15. (43) Trasatti, S. Pure Appl. Chem. 1986, 58, 955−966. (44) Hansen, M. H.; Jin, C.; Thygesen, K. S.; Rossmeisl, J. J. Phys. Chem. C 2016, 120, 13485−13491. (45) Björketun, M. E.; Zeng, Z.; Ahmed, R.; Tripkovic, V.; Thygesen, K. S.; Rossmeisl, J. Chem. Phys. Lett. 2013, 555, 145−148. (46) Larsen, A. H.; Vanin, M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Phys. Rev. B 2009, 80, No. 195112. (47) Perdew, J. P.; Levy, M. Phys. Rev. Lett. 1983, 51, 1884−1887. (48) Perdew, J. P. Int. J. Quantum Chem. 1985, 28, 497−523. (49) Sheng, W.; Myint, M.; Chen, J. G.; Yan, Y. Energy Environ. Sci. 2013, 6, 1509.

(50) Shinagawa, T.; Takanabe, K. J. Phys. Chem. C 2015, 119, 20453−20458. (51) Haghighat, S.; Dawlaty, J. M. J. Phys. Chem. C 2016, 120, 28489−28496. (52) Rossmeisl, J.; Chan, K.; Skúlason, E.; Björketun, M. E.; Tripkovic, V. Catal. Today 2016, 262, 36−40. (53) Jónsson, H.; Mills, G.; Jacobsen, K. W. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds.; World Scientific, 1998; Vol. 385. (54) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. J. Chem. Phys. 2000, 113, 9901−9904. (55) Nørskov, J. K.; Bligaard, T.; Logadottir, A.; Kitchin, J. R.; Chen, J. G.; Pandelov, S.; Stimming, U. J. Electrochem. Soc. 2005, 152, J23− J26. (56) Holewinski, A.; Linic, S. J. Electrochem. Soc. 2012, 159, H864− H870.

K


Controlled-Potential Simulation of Elementary Electrochemical

Recommend Documents