pH in Grand Canonical Statistics of an ... - ACS Publications

Article pubs.acs.org/JPCC

pH in Grand Canonical Statistics of an Electrochemical Interface Martin Hangaard Hansen†,‡ and Jan Rossmeisl*,† †

Nano Science Center, Department of Chemistry, University of Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark Center for Atomic-scale Materials Design, Department of Physics, Technical University of Denmark, 2800 Kongens Lyngby, Denmark

‡

S Supporting Information *

ABSTRACT: We present an atomic-scale model of the electrochemical interface, which unfolds the effects of pH and electrode potential using a generalized computational hydrogen electrode. The liquid structure of the solvent is included with the use of ab initio molecular dynamics to sample thousands of microstates with varying numbers of protons and electrons. The grand canonical probability weight function at constant pH and electrode potential is calculated a posteriori. The only inputs to the model are the fundamental assumptions of an equilibrated solvent, charge neutrality of the interface, and the dimensions of the system. The structures are unbiased outputs, and several atomic-scale quantities are calculated for our model system, water/Au(111), as weighted averages. We present the potentials of zero charge, Gibbs isotherms, and differential capacities as a function of pH. The potential of maximum entropy is also calculated, which to our knowledge has not previously been done with any first-principles method. The model predicts a non-Nernstian pH behavior for the potential of maximum entropy.

■

INTRODUCTION Electrochemical energy conversion and storage will become an increasingly important part of the future energy system.1,2 The energy conversion itself takes place via the electrochemical reactions in the interface between the electrolytes and electrodes. Further development of such technologies therefore requires better atomic-scale insight into the structure electrochemical interface.3−5 Macroscopic parameters of the electrolyte such as pH and concentrations of ions govern the rates of some reactions through changes happening in the interface on the molecular scale,6−13 and it is therefore necessary to understand how changing the electrolyte as well the surfaces affects the processes at the interface. Despite the growing demand for computational solutions, present methods still face a tough challenge in solving the problem of relating the atomicscale electrochemical interface models to the surrounding macroscopic environment.14,15 Various computational methods have been proposed previously for modeling the electrochemical interface using density functional theory (DFT).16−27 The methodology presented in the present paper distinguishes itself by treating both the protons and electrons grand canonically and focuses on determining the structure of the interface. The method is based on DFT and a generalization of the computational hydrogen electrode (CHE) using the work function as an absolute electrode potential scale and varying the numbers of protons and electrons explicitly.28 An a posteriori scheme is presented, which investigates pH-dependent structures of the aqueous electrochemical interface, with water/Au(111) as the case study. First, the theory of the © XXXX American Chemical Society

generalized computational hydrogen electrode (GCHE) is introduced, second the method for sampling of the interface partition functions is presented, and, finally, observables calculated as weighted averages are evaluated and compared to experimental results from literature.

■

THEORY Observables or properties of an atomic-scale system can be calculated as averages over all the states in which the system can be. The probability weight function in this sampling is the Gibbs factor, if the system is grand canonical. For a quantity of interest, A, this is written ⟨A⟩ =

∑ Aipi

(1)

i

pi = ZG =

1 exp[−(Ei − ZG

∑ μs ni ,s)/kBT ] s

∑ exp[−(Ei − ∑ μs ni ,s)/kBT ] i

s

(2)

(3)

where ZG is the grand partition function, pi is the probability of state, i, and Ei is the energy of the state.29 The sum ∑sμsni,s denotes the energy term for adding species s to the system in Received: September 7, 2016 Revised: November 10, 2016 Published: December 5, 2016 A

DOI: 10.1021/acs.jpcc.6b09019 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C equilibrium with a reservoir with chemical potential μs of that species. If an ensemble containing all relevant states is known and if the energy of all the states is known as a function of pH and electrode potential surrounding the interface, then ⟨A⟩ can be calculated. The interface system can be modeled as an extended metal/water slab with variable numbers of protons and electrons;20 see Figure 1.

(4) The interface region is charge neutral and large enough to screen all charges. Macroscopically, it is clear that the interface region is charge neutral, since it is surrounded by and in equilibrium with charge conducting media. A limitation to assumption 4 is that the surrounding media should be conductive enough that the Debye screening is on the order of the system size or smaller. Then it is a practical question to use sufficiently large models in the computations, to allow sufficient room for charges and counter charges. Assumptions 1, 2, and 4 add a constraint on the interface that leads to eURHE = ϕe− − ϕSHE + 2.303kBT · pH

(6)

where ϕe− is the work function and ϕSHE is the work function of the electrode at pH = 0 and USHE = 0 V. The reader is referred to our recent papers for further details on this.28,33 Equation 6 reflects the connection between the work function and the electrode potential, which is well-known in the experimental literature34,35 as an absolute electrode potential. From simulations, the work function can be obtained from the planar-averaged electrostatic potential and the Fermi level using ϕe− = E H(zmax ) −E F Figure 1. Snapshot from a molecular dynamics simulation plotted with the Hartree potential (black line) and the Fermi level (green line). The work function is the difference between the vacuum level on the right side and the Fermi level. The boundary conditions are periodic in the plane and nonperiodic with a dipole correction30 out of the plane.

⟨E H(z)⟩ =

2

− nμ0 H /2 2

− n(ϕSHE − ϕe− − 2.303kBT · pH)

(4)

(9)

Equation 6 acts as a constraint on the macroscopic interface, so in the limit of large unit cells, it must hold. In the present work, effects of unit cell size were not evaluated, so the rest of the analysis assumes that the unit cell is large enough and that this relation holds for every microstate.

Eint(n , URHE) = E(n) − E(n = 0) 2

(8)

Eint(n , ϕe− , pH) = E(n , ϕe−) − E(n = 0)

and the integral energy, Eint, of each state is accordingly

− n(μH0 /2 − eURHE)

∫x ,y EH(x , y , z)dxdy

where EH(x,y,z) is the electrostatic potential, which can be calculated for each structure by DFT (see Figure 1). Inserting eq 6 into eq 5, the integral energy of a microstate is written

In the original computational hydrogen electrode (CHE) scheme,31 the system exchanges protons and electrons in equal numbers, n; thus, ∑sμsni,s is reduced to n(μe− + μH+), where the i index is removed to simplify notation. The chemical potentials are given by μe− + μH+ = μH0 /2 − eURHE

(7)

■

(5)

METHOD To create an ensemble of interface structures and to calculate the grand partition function, ZG, ab initio molecular dynamics (AIMD) was employed. The electrode surface was represented by a four atomic layer Au(111) slab in a 3 × 4 orthogonal super cell. The electrolyte layer included 24 water molecules ± the above-mentioned hydrogen atoms. The system was centered with 9 Å vacuum on either side. The vacuum interface was kept by a constrained layer of eight water molecules in a hexagonal structure with a dipole moment out of the plane close to zero. The position of the constrained outer layer was chosen to approximately conform to the number density of water in bulk at standard conditions. The volume was counted from the surface atoms to the plane of the constrained water, and the constrained water molecules were counted in. Trajectories containing from −6 to +6 extra hydrogen atoms were made. Since the unit cell had 12 surface atoms, we write the coverages or excesses as n/N, where N = 12. The hydrogens were included as specifically adsorbed H* and/or one solvated proton, while n < 0 resulted in chemisorbed OH and/or O.

where E(n) is the energy obtained from the density functional theory and μ0H2 is the chemical potential of H2(g) at standard conditions, obtained from a DFT calculation and thermodynamic tables.31 Zero point energy differences are currently omitted. Usually, the integral free energy, Gint, is calculated using the enthalpy from the energy, Eint, of the ground state only, and the entropy contributions are taken from thermodynamic tables and configurational entropy of adsorbates on the surface. Instead, one may also calculate Eint of all microstates and get the free energies from the grand partition function, ZG, as we aim to do. Three assumptions are made in the original CHE: (1) The bulk electrolyte is in equilibrium. (2) The interfaces are in equilibrium with the bulk electrodes and the electrolyte surrounding them. (3) The energy of each state is independent of the electrostatic field and boundary conditions.32 If Eint is dependent on electric fields, as it may be in the case of water structures, the GCHE must be used.28 This replaces assumption 3 with assumption 4: B


Article

The Journal of Physical Chemistry C

Figure 2. Calculated Eint vs USHE and URHE (bottom and top axis, respectively) for all states, sampled from molecular dynamics. Coverage (n/N) denotes the number of H atoms added to the neutral water system divided by the number of surface atoms, and this value is denoted for every state by color. The two plots show projections of Eint onto pH = 0 (left panel) and pH = 14 (right panel).

Figure 2 shows the integral energy, Eint, of all the individual snapshots sampled from molecular dynamics at pH = 0 in the left panel and pH = 14 in the right panel. Eint is calculated using eq 9, where E(n,ϕe−) is the total energy from the molecular dynamics sampling; i.e., both kinetic and potential energy is included in E. Note in Figure 2 how the points fall on lines with the slope coming from the last term in eq 9. Note also that the crossing points of these lines are mostly at the same point on the URHE scale, while a different pH changes the work function at the crossing points, since the work function scales with the USHE. The forces and the electronic structure in every time step were calculated by density functional theory using the GPAW code.36,37 To reduce the total energy oscillations, the convergence criteria for the density was tightened to 2·10−5 Å3 to improve the self-consistency in calculation of the forces, still without sacrificing too much efficiency. The exchangecorrelation contribution to the energy was represented using the RPBE functional.38 Single electron wave functions were expanded in localized atomic orbitals (LCAO) with the doubleζ-polarized basis set.39 The k-point sampling was 3 × 3 × 1. In the out-of-plane direction, the boundary conditions were nonperiodic and the Bengtsson dipole correction30 was used for the electrostatic potential, in order to obtain the work function, ϕe−, for the a posteriori analysis. The LCAO basis sets do affect the accuracy of the work functions, since the electron spill out from the interface is not captured as accurately. The calculations with LCAO consistently underestimate the calculated work functions, which was seen by comparing a few structures calculated with LCAO to calculations with the more complete grid-based basis sets, which revealed differences of 0.3−0.5 eV. Therefore, we argue that the structure corresponding to 0 V versus SHE would have a work function 0.3−0.5 eV smaller in our DFT calculation than the 4.44 eV (or 4.85 eV) that is reported in experiments.35,40 Thus, a lower value of 3.9 eV was tentatively used for ϕSHE. The time step for the dynamics was 0.5 fs, and a Berendsen thermostat41 was employed with characteristic time constant of 1000 fs, to equilibrate the temperature of all trajectories toward 300 K. A total number of 90 000 states were sampled, giving each trajectory a sample time of approximately 2.5 ps. Prior to sampling, the structures had been initiated as relaxed icelike layers and annealed to around 300 K with coarser parameters

for the force calculation and more aggressive thermostat settings for at least another 5 ps and equilibrated with sample parameters for 1 ps before sampling. In the following, it is shown how the pH-dependent grand canonical averages were obtained from all the states sampled in the molecular dynamics trajectories. Calculating any atomic-scale quantity requires an ensemble with the correct averages of surface coverage and proton excess, at any given electrode potential and pH. This is not already achieved by averages of a raw ensemble, consisting of all states from the trajectories. It is our choice how long simulation time is given to each trajectory and thereby how many states are tried at each coverage. To make the ensemble unbiased with respect to n, the probability weight function was calculated using a grand canonical Metropolis Monte Carlo algorithm that walks around in all the states that are sampled. The probability, πij, to walk from an accepted state j to a trial state i and add i to the new ensemble is a piecewise function of the difference in energy ΔEijint ≤ 0 ⇒ πij = 1

(10)

ΔEijint > 0 ⇒ πij = exp( −ΔEijint /kBT )

(11)

ΔEint ij

where = − By moving between trajectories, the trial step can add or remove 0, 1, 2, or 3 protons + electrons, all at equal attempt probabilities. Added hydrogens are equally often as solvated protons or adsorbed hydrogens. Such a set of trial steps yields a quasi-ergodic Markov chain, and therefore it gives unbiased relative probabilities of occurrence for the states in the annealed ensemble.42 The resulting weight function is equivalent to just sampling equally many states with each number of particles. Eint of each state is only calculated relative to states with the same work function. This is achieved by binning the states according to their work function, resolved in a number of small intervals or bins. In min int essence, ΔEint i (ϕi) = Ei (ϕi ∈ b) − E (ϕ,b), where b denotes min the interval and E (ϕ,b) is the energy minimum as a function of work function. Emin(ϕ,b) is approximated in each bin simply as min{Eint(ϕ ∈ b)}. This works well if the bins are small enough, but the smaller the bins are, the more sampling is needed to get the converged averages in all bins. For the present calculations, 128 bins were used. Using significantly more bins increases the noise in the results, and it is undesirable C

ΔEint i

ΔEint j .


Article


Figure 3. Calculated Eint vs USHE and URHE (bottom and top axis respectively) for all atomic structure states sampled from molecular dynamics and annealed to 300 K with a grand canonical Monte Carlo sampling. n/N denotes the number of H atoms added to the neutral water system, and this value is shown for every state by color as indicated by the color bar. The two plots show projections of Eint onto pH = 0 (left panel) and pH = 14 (right panel).

the solvent layer, but no solution has yet been implemented to do this. Static atomic-scale properties of the interface, such as structure, can now be calculated from the annealed ensembles. Since the ensemble is made from AIMD, the structures are liquid and relevant at 300 K. Dynamic properties, such as diffusion constants, that are directly dependent on the trajectories are not currently calculated. In addition, the present method is still biased with regard to pressure, since a constant value is assumed for the volume. To alleviate this bias, one could test trajectories with different volumes or implement constant pressure dynamics for the relevant volume, which is from the constrained top layer to the surface. The latter would probably be the more desirable, since each additional trajectory has to be equilibrated separately.

to use significantly fewer bins in order to resolve results such as the Gibbs isotherms on the potential axis. For each pH value, an ensemble was generated using the procedure described above. The energies in the annealed ensemble follows a Gibbs distribution as would be expected.29 An example of this distribution can be found in the Supporting Information. The energy, Eint, is plotted versus electrode potential for ensembles resulting from the annealing in Figure 3, (a) for pH = 0 and (b) for pH = 14. Comparing Figure 3 with Figure 2 gives an overview over which states remain after the Metropolis Monte Carlo annealing. As an example, consider that states with high number of hydrogens at a very anodic potential are nonexistent in the annealed ensemble, as would be expected in the real system. The Monte Carlo annealing effectively filters away all states, which are irrelevant or improbable, because they are unstable at that pH and potential. As explained in previous papers,15,28,33 the pH effect enters on the atomic scale through the work function in the relevant set of states at a given electrode potential. For example, states at 0 V versus RHE and pH = 0 have a calculated work function of ϕSHE, while states at 0 V versus RHE at pH = 14 have a calculated work function of ϕSHE + 2.3 kBT·14 = ϕSHE + 0.8 eV. If the energies of some structures depend on the work function, they will be more or less likely to appear in the annealed ensemble lower or higher pH’s. In addition, different frequencies of occurrence in the dynamics simulation carry over in the Monte Carlo annealing, thereby capturing entropy dependence on the interface dipole and the work function. Water structures and surface species that are field-dependent on the atomic scale are thereby pH-dependent on the macroscopic scale. It should be noted, however, that adding or removing adsorbates or ions does not change the number of configurations in the system in the present implementation of the model. Protons also have the same number of configurations as adsorbed hydrogens. In reality, the number of configurations may differ significantly, but if the protons are located in the single layer closest to the surface, the number of configurations is probably on the same order of magnitude as adsorbed hydrogen. It would also be possible to assign different attempt probabilities for protons or ions to account for any difference if the number of configurations could be estimated in

■

RESULTS Simulations of several experimentally observable quantities are presented in this section, and the trends are thereafter compared to experiments in the Discussion. The calculated quantities include the Gibbs isotherms, differential capacities, and the potentials of zero charge. They are evaluated using grand canonical averages based on eq 1. First, the Gibbs isotherms are presented. They are calculated as the average value of n/N from the annealed ensemble. This is equivalent to n/N (pH, ϕe−)= 1 ZG(pH, ϕe−)

∑ i

ni exp[(−Eiint(ni , pH, ϕe−)/kBT )] N

(12)

The surface coverages for hydrogen, θH, or oxide species, θOH, or the proton excess, ΓH+, may be inserted instead of n/N to calculate their respective isotherm. In practice, the protons were distinguished from adsorbed hydrogens using a threshold distance to the surface. The proton excess is plotted in Figure 4, and the oxide species coverage is plotted in Figure 5. An equivalent plot of the hydrogen surface coverage can be found in the Supporting Information, but it is virtually identical at all pH values from 0 to 14. The isotherms are not directly observable in experiments, but they can be inferred from cyclic voltammetry (CV) by D


Article


Integrating the experimental differential capacity from the point of zero total charge (PZTC) to U leads the total charge, σtotal(U). This is written U

σtotal(U ) =

∫PZTC C dU

(14)

The differential capacity, C, may be evaluated from the present calculations using

Q = e n/N Figure 4. Calculated proton excess isotherms at pH = 0 and 14. In this region, nothing is chemisorbed on the surface, while the proton excess lines predict the stability of a proton in the double layer. The blue and red solid lines show an isotherm expression (eq 17), fitted to the calculated data. The dashed black line shows the isotherm calculated using the original computational hydrogen electrode (CHE), thus neglecting the thermodynamic constraint on the work function.

C=

dQ dU

(15)

(16)

where e is the unit charge. C is plotted for three different pH values in Figure 6a. In Figure 6b, C is plotted for pH = 7 together with experimental data in neutral solution.

■

DISCUSSION Experimental cyclic voltammograms on Au(111)43−47 show an oxidation region, which is highly asymmetric in the potential axis,45−47 and a hydrogen evolution region. No underpotential deposition of hydrogen is observed on Au(111).46 The presently calculated Gibbs isotherms contain the distinct regions, which are due to either hydrogen species, no species, or oxide species. A wide potential range of zero surface coverage comes out as expected on Au(111). On the RHE scale, the surface adsorption isotherms of hydrogen are almost identical at pH = 0 and pH = 14 as expected from past computational studies.32 This can be seen from the dashed black line in Figure 4, which is the isotherm calculated with the original computational hydrogen electrode, using eq 2, only counting in the single state with lowest total energy for each coverage and inserting the total energy difference from eq 6. This leaves out the constraint on the work functions, thereby assuming the total energy of the proton in the double layer is independent of the electric field, ignoring pH effects. The pH dependence in the isotherms can also be seen in Figure 7, where the equilibrium potentials of proton insertion and hydrogen chemisorption have been plotted for various pH values. The equilibrium potentials are shown

Figure 5. Oxide species isotherm at pH = 0 and pH = 14. It is not shown here which species are on the surface; only the surface specific number of hydrogen atoms that have been removed − n/N is shown. The dashed black line shows the isotherm calculated using the original computational hydrogen electrode (CHE).

calculating the differential capacity from the current, j, and the scan speed, v: j C= (13) v

Figure 6. Specific capacities, calculated from numerical differentials of ⟨n/N⟩. (a) Simulated specific capacities co-plotted at various pH. Three main features appear. From left to right, the first sharp peak indicates equilibrium potential of adsorbed H, the small second peak is proton excess, and the three peaks between 1 V and 2 V vs. RHE indicate chemisorbed OH or O. (b) Black lines: Simulated specific capacity at pH = 7. Blue dashed lines: Experimental differential capacities from CV’s in 0.1M NaClO4 by Hamelin et al.43 in the oxidizing region, and in the HER region by Wang et al.44 measured in 0.5M H2SO4 , plotted on the RHE axis. The high potential regions, marked by red dotted lines are very large and are therefore divided by 10 in the plot. E


Article


It is also possible that the most experimentally observed pH effects are not due to pH only, but rather governed by the counterions. Increasing the pH from acidic to alkaline will increase the chemical potential of positive counterions in the bulk solution, which will result in a replacement of the protons for cations in the interface. Introduction of counterions could also allow neutral ion pairs to account for dissociated or associated states in the interface region. These factors are not accounted for at the present stage, so to better model the alkaline regime, it would be the logical next step to add cations to the ensemble and include it in the analysis in the future. Considering the oxidizing region, rather little chemisorption of OH has been reported to occur on Au(111) in the ORR region in alkaline.50 In agreement with this, the calculated Gibbs isotherms, shown in Figure 5, show that low coverages of OH start to adsorp at around +1.0 V versus RHE. The original CHE predicts a slightly later OH chemisorption, as shown by the dashed black line. The discrepancy is likely due to a lack of stable reference states with 0 coverage and relevant high work functions in the generalized CHE ensemble. The cyclic voltammograms reported by Climent et al.51 show similar features to the simulated CVs in Figure 6b. Any features from counterions would also be omitted in the present study. There are some important additional contributions to the total charge, which are not due to chemisorption or ion insertion, and these are discussed in the following. The total charge, σtotal, which can be quantified from experiments using eq 14, can be divided into several contributions

Figure 7. Equilibrium potentials of hydrogen adsorption and of proton insertion to the interface plotted versus pH together with the potential of maximum entropy (PME). The former was evaluated by fitting an isotherm espression (eq 17) to the data in Figure 4.

together with the calculated potential of maximum entropy, which will be elaborated upon at the end of the discussion. The equilibrium potentials in Figure 7 were extracted from the calculated n/N from the DFT calculations by fitting the following isotherm expression n/N =

n max eq N (1 + exp[(eURHE − eURHE )/kBT ])

(17)

where n was either protons (seen in Figure 4) or adsorbed hydrogen atoms. In contrast to the chemisorbed species isotherm, the calculated proton excess is not overlapping at various pH values, as seen in Figure 4. The equilibrium potential for insertion of protons shows a slightly sub-Nernstian shift. The formation energy of states with a proton in the interface is affected by the work function, as reported in past studies,22,48,49 and reproduced in the present study (see Supporting Information). Therefore, a non-Nernstian shift in the equilibrium potential for proton insertion of close to 59 mV per pH unit was expected. Here it appears, however, that the proton gains a stabilization of no more than 0.15 eV as the pH is increased from 0 to 14, which could mean that it does not feel a very large part of the potential despite an increase in work function of 0.82 eV. To compare with previous methods, the dashed black line in Figure 4 shows the isotherm at undefined pH, using the original CHE, thereby ignoring the constraint on work function and effects of the interface dipole. The original CHE thus takes the energy of the most stable proton state and the most stable reference state,22 whereas the generalized CHE ensemble is constrained to the proton and reference states at the work functions that are relevant around the given potential. The large discrepancy clearly shows that dipole and pH effects on proton or other ion insertion energies should not be neglected. It should be mentioned that protons in the double layer do not give away a full electron to the surface in the calculations. Bader charge analyses of a few sample states including a solvated proton indicated that only 0.3−0.4 electrons were transferred to the slab, and that was also the result if the structures were recalculated with a more complete grid based basis set. This may reduce the pH effect on the proton stability, although the solvated proton still carries a significant positive charge. Finally, it may be that most of the variations in the electric field are located further away from the surface than the most stable location of the protons, which would also reduce the pH effect on the proton stability.

σtotal = σF − FθH + FθOH

(18)

where Fθs is charge attributed to electrosorption of species s. The free charge is σF, which is observed in CV experiments as the capacitive current from metal electrode in aqueous media. The free charge contains a contribution, σdl, which is electronic surface charge, screening ions in the double layer. Thus, σdl = FΓOH− − FΓH+ in the model, due to the assumption of charge conservation. States containing solvated OH− were not put into the dynamics, so ΓOH− is 0. Let us say there is another contribution, σdip, to the free charge that is not accounted for in the model: σF = σdl + σdip (19) σdip may arise from relocation of image charges in the electrode that screen water dipoles in the double layer. It must be a smaller contribution, since the presence of ions give rise to much larger dipoles than individual water molecules.52 In the simulated cyclic voltammograms, the only contribution that is not included is σdip. Thus, it is not counted in the calculated differential capacities (shown in Figure 6a and 6b). The potential where σdip is zero may be captured, however, as will be explained in the following. The potential of zero free charge (PZFC) and the potential of zero total charge (PZTC) are both quantities of interest to interpret cyclic voltammetry experiments, and in the following it is shown how they are obtained. Experimentally, the PZTC and the PZFC sometimes show non-Nernstian behavior (a potential shift of ≈2.3 kBT/e per pH unit on the RHE potential scale).53 The PZFC has been correlated closely with the potential of maximum entropy (PME).51,53 The explanation for the non-Nernstian shift may be found in the distribution of work functions and interface dipoles sampled by the molecular dynamics. F


Article


non-Nernstian behavior, assuming constant coverage. The calculated work functions of maximum entropy is not independent of coverage but varies as shown in Figure 9 for subsets with various coverages.

The net dipole moment of the water layer is the result of the integrated dipole moment in the direction normal to the surface, which is set up by the orientation of water molecules, by electronic counter charge and positions of ions. The molecular dynamics does not include any external fields, and thus it accesses all of the configuration space of up and down orientation of the molecules. The work function varies on a rather short time scale, indicating reorientation happens often and with negligible energy barriers (see the Supporting Information for plots of the work function versus time). The variation in net dipole moments of the water behaves somewhat like a random walk.54 The distribution of dipole moments projected onto the axis normal to the surface is reflected in the distribution of work functions, shown in Figure 8. This distribution resembles a normal distribution as expected

Figure 9. Peak of the work function distribution as a function of coverages of oxide species (−n/N) and of hydrogen θH with and without one proton in the double layer.

The work function is related to the potential of maximum entropy by eUSHE = ϕe− − ϕSHE. If ϕSHE = 3.9, then the potential of maximum entropy is −0.6 V versus SHE according to the data in Figure 9. The PME evaluated self-consistently with the surface coverage at various pH values is plotted in Figure 7. The self-consistently evaluated PME stays in the region with high hydrogen coverage, and therefore the ideal non-Nernstian behavior comes out. The experiments by Climent et al. concludes that the PME is around 0.2 V versus SCE, measured in 0.1 M HClO4, while the calculated PME in the subset with one proton and a high hydrogen surface coverage is −0.6 V versus SHE, which means the discrepancy is around 1 V. The discrepancy may be due to the presence of counterions, which have not yet been included in the calculations. It is straightforward to do so within the current methodology, disregarding the demands on computations to expand the ensemble. Since adding protons shifts the PME toward lower potentials, we expect adding a negative ion will shift the PME toward higher potentials and that will improve agreement with the experiment.

Figure 8. Histogram of calculated work functions in the ΓH+ = 0, θH = 0 subset of structures obtained directly from the molecular dynamics sampling.

from a system performing a random walk. The shape is converged rather quickly, after sampling just a few picoseconds, indicating that the interface dipole quickly becomes uncorrelated from the starting structure. The peak of the distribution shows which macrostate of the net dipole moment has the highest multiplicity, if the order parameter is the work function or interface dipole. Therefore, the peak of the distribution corresponds to the work function of maximum entropy. This may explain the observation of a potential of maximum entropy (PME) in electrochemical experiments by Climent et al.51,55 In these experiments, a very small heating of the interface is induced by a laser pulse on the order of microseconds, resulting in a shift toward a higher entropy macrostate. This shift is observed as a small transient in the measured potential, as the surface cools and returns to the equilibrium structure. The shift is explained by the following relation ∂U ∂ΔS =− ∂T q ∂q T

■

CONCLUSION We presented a more unbiased way of calculating pH and potential-dependent properties for the electrochemical interface on the atomic scale that does not assume a Helmholtz model. Fewer assumptions are necessary in this methodology, due to the application of the generalization of the computational hydrogen electrode to provide the energies of interface states as functions of pH and potential. Ab initio molecular dynamics (AIMD) was used to sample states with various numbers of adsorbed O, OH, H, and (H+ + e−). This grand canonical ensemble was used to calculate the static properties of the water/Au(111) interface as grand canonical averages. So far, only intrinsic pH effects were captured, but it is straightforward to include effects of chemical potentials of counterions within the methodology. Currently, the main limitation to the method is the high demand on time and computational resources for generating sufficiently large and well-equilibrated ensembles for obtaining accurate and well-converged results. AIMD is limited in terms of the time scales and configurations it can access, and it may be too slow to generate sufficiently equilibrated and

(20)

where q is charge and U is the electrode potential. This is derived for an ideal polarizable electrode from the electrocapillary equation.56,57 Equation 20 is used to identify the potential of zero potential transient with the potential of maximum entropy by measuring ΔU. Since the macrostate with the highest entropy is found at a certain potential, an interface already in this macrostate will show no shift and thus no potential transient. In the present model, the PME is constant on the SHE potential scale because the work function is linked directly to the SHE scale by eq 6. Therefore, the calculated PME will have G


Article


(5) Calle-Vallejo, F.; Koper, M. T. First-Principles Computational Electrochemistry: Achievements and Challenges. Electrochim. Acta 2012, 84, 3−11. (6) Hori, Y.; Murata, A.; Takahashi, R. Formation of Hydrocarbons in the Electrochemical Reduction of Carbon Dioxide at a Copper Electrode in Aqueous Solution. J. Chem. Soc., Faraday Trans. 1 1989, 85, 2309−2326. (7) Gattrell, M.; Gupta, N.; Co, A. A Review of the Aqueous Electrochemical Reduction of CO2 to Hydrocarbons at Copper. J. Electroanal. Chem. 2006, 594, 1−19. (8) Sheng, W.; Gasteiger, H. A.; Shao-Horn, Y. Hydrogen Oxidation and Evolution Reaction Kinetics on Platinum: Acid vs Alkaline Electrolytes. J. Electrochem. Soc. 2010, 157, B1529−B1536. (9) Markovic, N. M.; Sarraf, S. T.; Gasteiger, H. A.; Ross, P. N. Hydrogen Electrochemistry On Platinum Low-Index Single-Crystal Surfaces in Alkaline Solution. J. Chem. Soc., Faraday Trans. 1996, 92, 3719−3725. (10) Durst, J.; Siebel, A.; Simon, C.; Hasche, F.; Herranz, J.; Gasteiger, H. A. New Insights into the Electrochemical Hydrogen Oxidation and Evolution Reaction Mechanism. Energy Environ. Sci. 2014, 7, 2255−2260. (11) Subbaraman, R.; Tripkovic, D.; Strmcnik, D.; Chang, K.-C.; Uchimura, M.; Paulikas, A. P.; Stamenkovic, V.; Markovic, N. M. Enhancing Hydrogen Evolution Activity in Water Splitting by Tailoring Li+−Ni(OH)2−Pt Interfaces. Science 2011, 334, 1256−1260. (12) Zheng, J.; Sheng, W.; Zhuang, Z.; Xu, B.; Yan, Y. Universal Dependence of Hydrogen Oxidation and Evolution Reaction Activity of Platinum-Group Metals on pH and Hydrogen Binding Energy. Sci. Adv. 2016, 2, e1501602. (13) Varela, A. S.; Ju, W.; Reier, T.; Strasser, P. Tuning the Catalytic Activity and Selectivity of Cu for CO2 Electroreduction in the Presence of Halides. ACS Catal. 2016, 6, 2136−2144. (14) Schnur, S.; Groß, A. Properties of Metal−water Interfaces Studied from First Principles. New J. Phys. 2009, 11, 125003. (15) Nielsen, M.; Björketun, M. E.; Hansen, M. H.; Rossmeisl, J. Towards First Principles Modeling of Electrochemical Electrode− Electrolyte Interfaces. Surf. Sci. 2015, 631, 2−7. (16) Filhol, J.-S.; Neurock, M. Elucidation of the Electrochemical Activation of Water over Pd by First Principles. Angew. Chem., Int. Ed. 2006, 45, 402−406. (17) Otani, M.; Sugino, O. First-Principles Calculations of Charged Surfaces and Interfaces: A Plane-Wave Nonrepeated Slab Approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 115407. (18) Taylor, C. D.; Wasileski, S. A.; Filhol, J.-S.; Neurock, M. First Principles Reaction Modeling of the Electrochemical Interface: Consideration and Calculation of a Tunable Surface Potential from Atomic and Electronic Structure. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 165402. (19) Taylor, C.; Kelly, R. G.; Neurock, M. Theoretical Analysis of the Nature of Hydrogen at the Electrochemical Interface between Water and a Ni(111) Single-Crystal Electrode. J. Electrochem. Soc. 2007, 154, F55−F64. (20) Skulason, E.; Karlberg, G. S.; Rossmeisl, J.; Bligaard, T.; Greeley, J.; Jonsson, H.; Nørskov, J. K. Density Functional Theory Calculations for the Hydrogen Evolution Reaction in an Electrochemical Double Layer on the Pt(111) Electrode. Phys. Chem. Chem. Phys. 2007, 9, 3241−3250. (21) Jinnouchi, R.; Anderson, A. B. Electronic Structure Calculations of Liquid−Solid Interfaces: Combination of Density Functional Theory and Modified Poisson−Boltzmann Theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 245417. (22) Rossmeisl, J.; Skúlason, E.; Björketun, M. E.; Tripkovic, V.; Nørskov, J. K. Modeling the Electrified Solid−Liquid Interface. Chem. Phys. Lett. 2008, 466, 68−71. (23) McCrum, I. T.; Janik, M. J. pH and Alkali Cation Effects on the Pt Cyclic Voltammogram Explained Using Density Functional Theory. J. Phys. Chem. C 2016, 120, 457−471.

uncorrelated ensembles, but the methodology may as well rely on other sampling algorithms, which could be more efficient. Other limitations to the current version of the method include the following: (1) the assumption of high electrolyte conductivity, (2) counter ions and neutral pairs should be accounted for in the ensembles, (3) the effect of pressure is not yet included. Limitation 2 is straightforward to solve but will come at additional computational expense because additional trajectories must be added. Alleviating limitation 3 could probably be done using a form of constant pressure dynamics, which is not currently included. Further developments in software and hardware may provide the possibility of faster and more uncorrelated sampling of the interface and improvements to accuracy. Improvements to accuracy should include van der Waals forces and also use larger unit cells to fit solvation shells of ions. Given sufficient time and resources, the method may have the potential for making reference databases effectively predicting interface Pourbaix diagrams with structures of single crystal interfaces, which would be very useful for development of electrochemical interfaces for fuel cells, electrolysis cells, batteries, and other electrochemical devices.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b09019. Histograms of work functions calculated in each AIMD step, plotted for trajectories with various numbers, n, of hydrogen atoms added, and of those how many are specifically adsorbed; histogram showing the total energy distribution of an ensemble after MC steps; adsorbtion isotherm for atomic hydrogen at pH = 0 and pH = 14; plot of work function versus time for a molecular dynamics trajectory (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Martin Hangaard Hansen: 0000-0003-0818-1515 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Axel Groß, Marc T. M. Koper, and Karsten W. Jacobsen for fruitful discussions. The work is funded by the Danish Council for Strategic Research via the NACORR Project No. 12-132695, by the Carlsberg Foundation via case No. CF15-0165, by the Technical University of Denmark, and by the University of Copenhagen.

■

REFERENCES

(1) Schlögl, R. Sustainable Energy Systems: The Strategic Role of Chemical Energy Conversion. Top. Catal. 2016, 59, 772−786. (2) Dunn, B.; Kamath, H.; Tarascon, J.-M. Electrical Energy Storage for the Grid: A Battery of Choices. Science 2011, 334, 928−935. (3) Gasteiger, H. A.; Kocha, S. S.; Sompalli, B.; Wagner, F. T. Activity Benchmarks and Requirements for Pt, Pt-Alloy, and non-Pt Oxygen Reduction Catalysts for PEMFCs. Appl. Catal., B 2005, 56, 9−35. (4) Koper, M.; Wieckowski, A. Fuel Cell Catalysis: A Surface Science Approach; The Wiley Series on Electrocatalysis and Electrochemistry; Wiley: New York; 2009. H


Article

The Journal of Physical Chemistry C (24) Bonnet, N.; Marzari, N. First-Principles Prediction of the Equilibrium Shape of Nanoparticles under Realistic Electrochemical Conditions. Phys. Rev. Lett. 2013, 110, 086104. (25) Todorova, M.; Neugebauer, J. Extending the Concept of Defect Chemistry from Semiconductor Physics to Electrochemistry. Phys. Rev. Appl. 2014, 1, 014001. (26) Zeng, Z.; Greeley, J. Characterization Of Oxygenated Species At Water/Pt(111) Interfaces From DFT Energetics And XPS Simulations. Nano Energy 2016, 29, 369. (27) Zhao, M.; Anderson, A. B. Predicting pH Dependencies Of Electrode Surface Reactions In Electrocatalysis. Electrochem. Commun. 2016, 69, 64−67. (28) Rossmeisl, J.; Chan, K.; Ahmed, R.; Tripkovic, V.; Björketun, M. E. pH in Atomic Scale Simulations of Electrochemical Interfaces. Phys. Chem. Chem. Phys. 2013, 15, 10321−10325. (29) Gould, H.; Tobochnik, J. Statistical and Thermal Physics: With Computer Applications; Princeton University Press: Princeton, NJ, 2010. (30) Bengtsson, L. Dipole Correction For Surface Supercell Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 12301−12304. (31) Nørskov, J. K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J. R.; Bligaard, T.; Jónsson, H. Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode. J. Phys. Chem. B 2004, 108, 17886−17892. (32) Rossmeisl, J.; Nørskov, J. K.; Taylor, C. D.; Janik, M. J.; Neurock, M. Calculated Phase Diagrams for the Electrochemical Oxidation and Reduction of Water over Pt(111). J. Phys. Chem. B 2006, 110, 21833−21839. (33) Hansen, M. H.; Jin, C.; Thygesen, K. S.; Rossmeisl, J. Finite Bias Calculations to Model Interface Dipoles in Electrochemical Cells at the Atomic Scale. J. Phys. Chem. C 2016, 120, 13485−13491. (34) Kolb, D. M. Electrochemical Surface Science. Angew. Chem., Int. Ed. 2001, 40, 1162−1181. (35) Trasatti, S. The Absolute Electrode Potential: An Explanatory Note (Recommendations 1986). Pure Appl. Chem. 1986, 58, 955. (36) Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-Space Grid Implementation of the Projector Augmented Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 035109. (37) Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dulak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A. Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. J. Phys.: Condens. Matter 2010, 22, 253202. (38) Hammer, B.; Hansen, L. B.; Nørskov, J. K. Improved adsorption energetics within density-functional theory using revised PerdewBurke-Ernzerhof functionals. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 7413−7421. (39) Larsen, A. H.; Vanin, M.; Mortensen, J. J.; Thygesen, K. S.; Jacobsen, K. W. Localized Atomic Basis Set in the Projector Augmented Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 195112. (40) Kötz, E.; Neff, H.; Müller, K. A UPS, XPS, and Work Function Study of Emersed Silver, Platinum, and Gold Electrodes. J. Electroanal. Chem. Interfacial Electrochem. 1986, 215, 331−344. (41) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (42) Valleau, J. P.; Cohen, L. K. Primitive Model Electrolytes. J. Chem. Phys. 1980, 72, 5935−5941. (43) Hamelin, A.; Sottomayor, M.; Silva, F.; Chang, S.-C.; Weaver, M. J. Cyclic Voltammetric Characterization of Oriented Monocrystalline Gold Surfaces in Aqueous Alkaline Solution. J. Electroanal. Chem. Interfacial Electrochem. 1990, 295, 291−300. (44) Wang, T.; Liu, L.; Zhu, Z.; Papakonstantinou, P.; Hu, J.; Liu, H.; Li, M. Enhanced Electrocatalytic Activity for Hydrogen Evolution Reaction from Self-Assembled Monodispersed Molybdenum Sulfide Nanoparticles on an Au Electrode. Energy Environ. Sci. 2013, 6, 625− 633.

(45) Rodriguez, P.; Feliu, J. M.; Koper, M. T. Unusual Adsorption State of Carbon Monoxide on Single-Crystalline Gold Electrodes in Alkaline Media. Electrochem. Commun. 2009, 11, 1105−1108. (46) Perez, J.; Gonzalez, E. R.; Villullas, H. M. Hydrogen Evolution Reaction on Gold Single-Crystal Electrodes in Acid Solutions. J. Phys. Chem. B 1998, 102, 10931−10935. (47) Chen, A.; Lipkowski, J. Electrochemical and Spectroscopic Studies of Hydroxide Adsorption at the Au(111) Electrode. J. Phys. Chem. B 1999, 103, 682−691. (48) Roudgar, A.; Groß, A. Water Bilayer on the Pd/Au(111) Overlayer System: Coadsorption and Electric Field Effects. Chem. Phys. Lett. 2005, 409, 157−162. (49) Schnur, S.; Groß, A. Challenges in the First-Principles Description of Reactions in Electrocatalysis. Catal. Today 2011, 165, 129−137. Theoretical Catalysis for Energy Production and Utilization: From First Principles Theory to Microkinetics. (50) Štrbac, S.; Adžić, R. The Influence of OH− Chemisorption on the Catalytic Properties of Gold Single Crystal Surfaces for Oxygen Reduction in Alkaline Solutions. J. Electroanal. Chem. 1996, 403, 169− 181. (51) Climent, V.; Coles, B. A.; Compton, R. G. Laser-Induced Potential Transients on a Au(111) Single-Crystal Electrode. Determination of the Potential of Maximum Entropy of DoubleLayer Formation. J. Phys. Chem. B 2002, 106, 5258−5265. (52) Schmickler, W. A Jellium-Dipole Model for the Double Layer. J. Electroanal. Chem. Interfacial Electrochem. 1983, 150, 19−24. (53) García-Aráez, N.; Climent, V.; Feliu, J. M. Potential-Dependent Water Orientation on Pt(111) Stepped Surfaces from Laser-Pulsed Experiments. Electrochim. Acta 2009, 54, 966−977. (54) Eminyan, M.; Rubin, K. Introduction à la simulation des systèmes physiques: Un manuel interactif pour comprendre la physique par l’informatique; InterEditions: Paris, 1994. (55) Climent, V.; Coles, B. A.; Compton, R. G. Laser Induced Current Transients Applied to a Au(111) Single Crystal Electrode. A General Method for the Measurement of Potentials of Zero Charge of Solid Electrodes. J. Phys. Chem. B 2001, 105, 10669−10673. (56) Harrison, J.; Randles, J.; Schiffrin, D. The entropy of formation of the mercury-aqueous solution interface and the structure of the inner layer. J. Electroanal. Chem. Interfacial Electrochem. 1973, 48, 359− 381. (57) Benderskii, V.; Velichko, G. Temperature jump in electric double-layer study. J. Electroanal. Chem. Interfacial Electrochem. 1982, 140, 1−22.

I


pH in Grand Canonical Statistics of an ... - ACS Publications

Recommend Documents