Surface Energy as a Descriptor of Catalytic Activity - The Journal of

Sep 27, 2016 - For these particular materials, the surface energies in vacuum are nearly unchanged upon .... JDFTx: Software for joint density-functio...
0 downloads 0 Views 912KB Size
Article pubs.acs.org/JPCC

Surface Energy as a Descriptor of Catalytic Activity Published as part of The Journal of Physical Chemistry virtual special issue “Mark S. Gordon Festschrift”. Houlong Zhuang,† Alexander J. Tkalych,‡ and Emily A. Carter*,§ †

Department of Mechanical and Aerospace Engineering, ‡Department of Chemistry, and §School of Engineering and Applied Science, Princeton University, Princeton, New Jersey 08544, United States S Supporting Information *

ABSTRACT: Computational searches for catalysts of the hydrogen evolution reaction commonly use the hydrogen binding energy (HBE) as a predictor of catalytic activity. Accurate evaluation of the HBE, however, can involve large periodic supercell slab models that render high-throughput screening relatively expensive. In contrast, calculations of other relevant surface properties, such as the surface energy, work function, and potential of zero charge (PZC), require only small surface unit cells and are hence less expensive to compute. Correlations between catalytic activity and these surface properties warrant exploration because of this reduced computational cost. Here, we use density functional theory in conjunction with three different exchange-correlation functionalsthe local density approximation (LDA), the Perdew−Burke− Ernzerhof (PBE) generalized gradient approximation, and the PBEsol functional (a reparameterization of the PBE functional)to calculate the lattice constants, surface energy, cohesive energy, and work function of six common catalysts: three metals (Au, Pd, and Pt) and three transitionmetal carbides (TMCs; WC, W2C, and Mo2C). The three exchange-correlation functionals produce identical trends, and PBEsol yields results between those calculated using LDA and PBE and most often closer to experiment. We therefore use PBEsol to obtain the surface energy, work function, and PZC of nine novel hybrid catalysts, each containing a metal monolayer on a TMC substrate. Importantly, a volcano-shaped correlation between the experimental exchange current density and the theoretical surface energies emerges. We also investigate solvation effects on the surface energy and work function using a polarizable continuum model within the framework of joint density functional theory. For these particular materials, the surface energies in vacuum are nearly unchanged upon exposure to an aqueous solution, prior to any reaction with water. The volcano-shaped correlation observed between the exchange current densities and the surface energies is not observed for the work function or PZC. Our work thus reveals opportunities for more rapid computational screening of reduced Pt-loading catalysts using the surface energy as a computationally efficient catalytic descriptor.



INTRODUCTION More than half of industrially produced hydrogen originates from the water gas shift reaction.1 This hydrogen is used in the manufacturing of numerous important chemicals, like ammonia and methanol, as well as in processes like catalytic reforming of petroleum.1−3 Despite the clear importance of these uses of hydrogen, the additional possibility of utilizing hydrogen as a clean energy carrier drives much of the current research. The high energy costs and carbon dioxide emissions of the water gas shift reaction motivate research into improving efficiency and lowering costs of alternative hydrogen production processes. The electrolysis of water consists of two half reactions: the hydrogen evolution reaction (HER) and the oxygen evolution reaction. Nicholson and Carlisle first observed the electrolysis of water in 1800 when they inserted platinum wires from each end of a voltaic pile into a small tube filled with water.4 Over two hundred years later, platinum and expensive metal oxides such as IrO2 and RuO2 remain among the most effective electrocatalysts for water-splitting.5 However, the scarcity and associated costs of these materials highlight the necessity of © XXXX American Chemical Society

dramatically reducing, if not eschewing entirely, the use of these metals. There has been a tremendous amount of work investigating the possibility of replacing Pt as the material of choice for the HER. To this end, a number of materials have been studied and shown to have potential. Ni-based binary transition-metal alloys,6,7 transition-metal chalcogenides such as MoS2,8,9 transition-metal nitrides,10 borides,11 and phosphides12 all exhibit significant HER activity. One of the most promising classes of materials for this reaction are transition-metal carbides (TMCs).13,14 A recent study of various TMCs at elevated temperatures (200−400 °C) found that the onset reduction potential for the HER decreases (becomes more negative, i.e., harder to reduce) in the order WC > Pt = Mo2C > NbC > TaC.15 However, the activity of WC has been known since the 1970s13 to be much lower than Pt in solution because Received: September 24, 2016

A

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

Article

The Journal of Physical Chemistry C

therefore is of great interest. The present work aims to identify potential correlations between catalytic activity and different surface properties such as the surface energy, work function, and potential of zero charge (PZC), all of which can be calculated using only surface unit cells. These properties were chosen because of their possible relevance to surface catalysis. The surface energy is influenced by the degree of coordinative unsaturation of the surface atoms; a higher surface energy is generally indicative of a more reactive surface. The work function is related to the ease of electron transfer, which is important for the reduction reactions of interest here. Lastly, the PZC can be regarded as the work function as modified by the presence of the liquid−solid interface. We consider here 15 catalysts, including six conventional catalysts (Au, Pd, Pt, WC, W2C, and Mo2C) and nine novel hybrid catalysts (the nine possible monolayer combinations of the three metals on the three TMC substrates). The exchange current densities of almost all of these catalysts (except Pt/ Mo2C, Pd/W2C, and Au/W2C) have been measured in recent experiments.18 Variations of the experimental exchange current densities with our computed surface energies reveal a volcanoshaped correlation, which shows that (1) Pt exhibits an optimal surface energy, consistent with its superior catalytic activity among these 15 catalysts, and (2) the catalytic activities of the other 14 catalysts are strongly dependent on the deviation of their surface energies from the optimal surface energy. On the basis of our predicted surface energies of these novel hybrid catalysts, we predict that Pt/Mo2C should exhibit catalytic activity similar to Pt and Pt/WC. Notably, no such volcanoshaped correlation is found for either the work function or the PZC, although an almost constant shift exists between these two properties for all of the catalysts considered.

of the former’s predicted stronger binding of hydrogen atoms.16,17 Nevertheless, the TMCs are still viewed as auspicious for reducing Pt loading in electrocatalytic systems. This is because many TMCs, including WC, W2C, Mo2C, TaC, VC, and TiC, have electronic properties similar to those of Pt-group metals and are stable under electrochemical conditions.18,19 These properties make the TMCs particularly well-suited to act as substrates for Pt-group metal overlayers. Extensive work by Chen and co-workers has demonstrated that Pt overlayer coverages as low as 1 monolayer (ML) on a WC substrate leads to a catalyst as active as bulk Pt.18,20−22 Further research from the same group has shown that a coverage of 1.5 ML of Pt on a Mo2C substrate recovers most of the HER activity of bulk Pt.18 Another study found that ML Pd is also active toward the HER on both Mo2C and WC substrates.23 Given these findings, and the many possible combinations of other metal overlayers and TMC supports, it would be of great interest to develop a technique that could quickly screen for the suitability of such hybrid material HER catalysts. Calculation of the hydrogen binding energy (HBE) is one of the most effective computational tools for predicting the efficacy of a material surface toward the HER. Sabatier’s principle states that the activity of a catalyst is related to the binding interactions between the surface and the reactive intermediates: surfaces having too weak an interaction with an intermediate will fail to retain the intermediate long enough for a reaction, whereas surfaces having too strong an interaction with an intermediate will fail to release the intermediate for further reaction.24 This qualitative principle has been used with impressive success for predicting the relative reactivities of various catalytic surfaces for a variety of reactions. In 1958, Parson found that the maximum exchange current density for the HER was attained when the hydrogen adsorption free energy is close to zero.25,26 Numerous groups have illustrated the truth of Sabatier’s principle using so-called “volcano plots.” These plots demonstrate that those catalysts for which the Gibbs free energy of adsorbed atomic hydrogen is closest to zero are generally the most active. Indeed, the most active catalyst for the HER, Pt, has a Gibbs free energy of hydrogen atom binding that is very nearly zero. Computational studies have repeatedly shown that the volcano plot relationship holds for various pure metal catalysts and have used calculated HBE values to successfully predict new HER catalysts, e.g., BiPt and MoS2.27−29 Rigorous, thorough calculations of the HBE are timeconsuming. Such calculations involve consideration of a variety of possible hydrogen adsorption sites on the surface. Vibrational frequencies of the optimized structures must then be calculated to ensure that the optimized structures are true minima and to calculate thermal corrections. Although this approach has proven successful, it is not the only possible predictor of catalytic performance. There are many other factors besides the HBE that determine the activity of an electrocatalyst, such as conductivity, crystallinity, the local environment of the electrolyte at the surface, and surface roughness. For example, work by Schmickler and co-workers suggests that, in addition to following Sabatier’s principle, it is important that an HER catalyst have a d-band that spans the Fermi level, and that there is a strong long-range interaction between the hydrogen 1s orbital and the catalyst’s d-orbitals.30 An accurate descriptor retaining the general applicability of the HBE without the associated computational expense



METHODS Density functional theory (DFT) calculations were performed within the Vienna Ab initio Simulation Package (VASP, version 5.3.5).31 The projector-augmented-wave (PAW) method32,33 within VASP was employed for all calculations. The 2s2p states of C; the 4d5s states of Mo and Pd; and the 5d6s states of W, Pt, and Au were treated as valence electrons, using the default PAW potentials in VASP to describe the nuclei and all (frozen) core electrons. We employed three different exchangecorrelation (XC) functionals: the local-density approximation (LDA),34,35 the generalized gradient approximation within the Perdew−Burke−Ernzerhof (PBE) form,36 and PBE with a reparametrization for solids (PBEsol).37 We used periodic slab models for calculating surface properties. Each slab had two symmetric surfaces. A vacuum spacing of 30 Å was used to isolate the periodic images from each other. We examined the dependence of the surface energies on slab thickness and found that nine layers of Pd, Pt, and Au and 11 layers of WC, W2C, and Mo2C are sufficient to converge the surface energies to an accuracy of 6 mJ/m2. We adopted these surface slabs with converged surface energies for subsequent work function and PZC calculations. We used Γpoint-centered Monkhorst−Pack38 k-point grids for all calculations. The k-meshes for the bulk crystals were 18 × 18 × 18 for conventional (four-atom) cubic unit cells of facecentered cubic Au, Pd, and Pt; 30 × 30 × 30 for the primitive (two-atom) unit cell of hexagonal-close-packed (hcp) WC; 27 × 27 × 18 for the primitive (three-atom) unit cell of hexagonal W2C; and 15 × 12 × 15 for the orthorhombic (12-atom) unit cell of Mo2C. We used a k-point grid of 33 × 33 × 1 for the Au, B

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

Article

The Journal of Physical Chemistry C

Table 1. Lattice Constants a0, b0, and c0 (Å); Surface Energy, Esurf (J/m2); Bulk Cohesive Energy, Ecoh (eV/atom); and Work Function, Φ (eV), of Au, Pd, Pt, WC, W2C, and Mo2C Obtained with the LDA, PBE, and PBEsol XC Functionalsa a0 Au

Pd

Pt

WC

W2C

Mo2C

b0 b

4.052 4.157c 4.081d 4.06552 3.841b 3.944c 3.874d 3.85952 3.897b 3.968c 3.916d 3.91252 2.885b 2.919c 2.898d 2.90258 3.019b 3.057c 3.033d 3.00159 4.682b 4.740c 4.704d 4.75258

c0















2.814b 2.844c 2.827d 2.84958 4.619b 4.670c 4.637d 4.73659 5.168b 5.232c 5.190d 5.21058



5.986b 6.060c 6.011d 6.01158

Esurf

Ecoh b

b

b

Φ

PZC

1.125 0.700c 0.976d 1.5053 1.850b 1.331c 1.675d 2.0553 1.935b 1.464c 1.806d 2.4853 3.804b 3.353c 3.686d

4.311 3.036c 3.722d 3.8154 5.062b 3.747c 4.476d 3.8954 7.107b 5.516c 6.402d 5.8454 9.927b 8.369c 9.047d

5.57 5.14c 5.29d 5.26 ± 0.0455 5.67b 5.24c 5.39d 5.3−5.656 6.09b 5.71c 5.83d 6.10 ± 0.0657 5.19b 4.97c 5.06d

4.99d

3.290b 2.906c 3.194d

9.929b 8.261c 8.973d

5.28b 4.99c 5.10d

4.68d

3.441b 2.989c 3.310d

8.696b 7.010c 7.737d

4.98b 4.64c 4.78d

4.25d

5.01d

5.50d

4.59d

a PZCs (V) for each of these materials, calculated using PBEsol, are also listed. Only the lowest surface energy is shown for WC, W2C, and Mo2C surfaces. Available experimental data (italicized) are shown for comparison. bLDA. cPBE. dPBEsol.

solvent) is divided into three contributions that arise from the electronic subsystem (solute), the liquid subsystem (solvent), and the coupling interactions between the two subsystems. Each of these contributions can be dealt with at different levels of approximations, for example Schwartz et al. employed a more accurate but also expensive quantum Monte Carlo method to treat the electronic system.45 In our work, the electronic system is treated at the DFT level. We then choose PCMs to approximate the liquid as a dielectric response. Correspondingly, there exists a dielectric cavity at the interface between the electronic and liquid subsystems. Two empirical parameters nc and σ enter the cavity shape function, s, that determines the shape of the cavity:

Pd, Pt, WC, and W2C surface calculations. The extremely dense k-point grid used in these surface calculations was used solely to demonstrate efficiency; a coarser k-point grid (20 × 20 × 1) in fact is sufficient to reach an accuracy of 1.0 meV/atom for the total energy. For the Mo2C surfaces with much larger in-plane lattice constants, we used a k-point grid of 12 × 15 × 1. We truncated the plane-wave basis at a kinetic energy of 700 eV, which, together with the above k-point grids, converged all total energies to within 1.0 meV/atom. We fully optimized the atomic coordinates until the interatomic forces were below 0.01 eV/Å. The Brillouin zone was integrated using the Methfessel− Paxton method,39 with a smearing width of 0.05 eV during total energy evaluation and geometry optimization. To calculate the bulk cohesive energy Ecoh, each of the constituent atoms is placed in the center of a rectangular vacuum box with asymmetric side lengths of 20, 21, and 22 Å. This asymmetric box is necessary to obtain nonspherically symmetric orbital occupations corresponding to the electronic ground states of these isolated atoms. We determined the DFT ground states of C, Mo, Pd, W, Pt, and Au atoms to be 3P, 7S, 1 7 3 S, S, D, and 2S, respectively. All of these DFT ground states are consistent with experimental atomic spectroscopic data,40 except for W, for which DFT produces the wrong ground state (7S instead of 5D)17 caused by errors in the XC functionals.17,41 We used the 7S state in the ensuing cohesive energy calculations because this state is 0.62 eV lower than the 5D state when using the PBE functional.17 Because they are less commonly used, we briefly describe joint density functional theory (JDFT) and the polarizable continuum model (PCM) employed here, referring the reader to the original references42−44 for further details. In the JDFT framework, the free energy of a solvated system (solute and

s=

⎧ log(n/nc) ⎫ 1 ⎬ erfc⎨ ⎩ σ 2 ⎭ 2

(1)

where nc is the critical electron density at which the dielectric cavity forms, separating the solvent and the solute, and smoothly switches the dielectric constant from the solvent value to the vacuum value (ε = 1); σ is the width controlling the transition from liquid to solute. We selected the recent parametrizations of nc and σ for our PCM calculations.46 These parameters were obtained by fitting calculated solvation energies of a number of molecules to experimental data. The JDFT-PCM method has been implemented in the VASPsol47 and JDFTx48,49 software packages. In this work, we mainly utilized the former package for solvation calculations with water as the solvent (i.e., dielectric constant ε = 78.0). We used JDFTx to crosscheck several surface energies in both vacuum and liquid environments and compared the surface energies with those calculated with VASP. To reduce computational cost within JDFTx, the optimized geometries C

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

Article

The Journal of Physical Chemistry C

where Eslab is the total energy of a surface periodic slab, Ebulk refers to the energy per atom of the bulk metal, N denotes the number of atoms in the surface slab, and A0 is the crosssectional area of the surface slab unit cell. Table 1 lists the predicted Esurf of the three metals. The LDA functional yields surface energies closest to experimental values for Au, Pd, and Pt, consistent with numerous previous theoretical studies (e.g., ref 60). Jellium calculations for metal surfaces have shown that the LDA functional benefits from error cancellation (overestimation of exchange is compensated for by underestimation of correlation) while the PBE functional underestimates both contributions.61,62 Our results suggest that this error cancellation also occurs for real metal surfaces such Au, Pd, and Pt. As with the lattice constants, the PBEsol surface energies lie between the LDA and PBE values, and all functionals predict similar trends in surface energies: Pt > Pd > Au. Corresponding trends are seen in the bulk cohesive energies, Ecoh, calculated as the energy difference between a bulk crystal unit cell and that of its constituent (isolated) atoms. The cohesive energy is the energy cost to break all bonds in the bulk solid. Because the same number of bonds are broken to create each of the (111) surfaces of Au, Pd, and Pt, all three functionals produce the same trend in Ecoh as for Esurf. The cohesive energy therefore serves as a useful indicator for expected trends in energies of analogous (same Miller index) metal surfaces. We next evaluated the surface energies of the three TMCs. We calculated the surface energy of W- or C-terminated WC using the same methodology adopted in our previous work on WC.17 The method for computing the surface energy of the other two TMCs M2C (M = Mo and W) is similar:63

from VASP were directly used as input configurations for the JDFTx calculations. Furthermore, we needed to adopt coarser k-point grids in the JDFTx calculations, i.e., 12 × 12 × 1 for all surface slabs other than the bare Mo2C surface and its hybrid structures with metal monolayers, for which 3 × 5 × 1 was used. These settings represented the most accurate calculations we could afford. The JDFTx calculations employed the PBEsol functional and ultrasoft pseudopotentials from the Garrity− Bennett−Rabe−Vanderbilt pseudopotential library.50 The kinetic energy cutoff of the plane wave basis was set to 30 hartree. The PBEsol surface energies calculated with VASP and JDFTx are very similar (see Supporting Information, Tables S1 and S2), even though the parameters used in the two programs are different. This good agreement validates the JDFTx parameters used in the ensuing solvation calculations.



RESULTS To benchmark our choice of physical models, we first compare the lattice constants (a0, b0, and c0) of bulk Au, Pd, Pt, WC, W2C, and Mo2C calculated using the LDA, PBE, and PBEsol XC functionals. Experimental lattice constants for each of these materials are also listed in Table 1. Figure 1 illustrates the

Esurf =

1 (Eslab − NMμM ‐NCμC ) 2A 0

(3)

where Eslab is the energy of a M2C surface slab with M- or Ctermination and NM and NC are the number of M and C atoms in the surface slab, respectively. The chemical potentials μM and μC then satisfy the following condition: μ M C (bulk) = 2μM + μC

Figure 1. Top and side views of primitive unit cells of (a) WC, (b) W2C, and (c) Mo2C. The two inequivalent W sites in W2C are denoted by I and II, respectively.

(4)

2

where μM2C (bulk) refers to the energy per formula unit of bulk M2C. The surface energy can be written as a function of the chemical potential μC:

structures of hcp WC, hexagonal W2C, and orthorhombic Mo2C, with space groups P6̅m2, P3̅1m, and Pbcn, respectively. Table 1 lists the theoretical lattice constants obtained with the three XC functionals. We observe typical underestimation and overestimation of most of the experimental lattice constants by the LDA and PBE XC functionals, respectively.51 The PBEsol functional frequently remedies the shortcomings of the LDA and PBE XC functionals, more often leading to the best agreement with experimental values for the six bulk materials. We then computed the surface energy of the three metals and three TMCs using the LDA, PBE, and PBEsol functionals. We focused on the (111) surface of the metals, the (0001) surfaces of WC and W2C, and the (100) surface of Mo2C. These surfaces were selected in order to eventually compare with their HER activities, as determined in ref 18. For the three pure metals, the surface energy, Esurf, is calculated as 1 Esurf = (Eslab − NE bulk ) 2A 0 (2)

Esurf =

⎛N ⎞ ⎞ 1 ⎛ 1 ⎜Eslab − NMμ M C (bulk) + ⎜ M − NC⎟μC ⎟ 2 ⎝ 2 ⎠ ⎠ 2A 0 ⎝ 2 (5)

In the C-rich limit, μC refers to the free energy of a C atom within graphite.We used Grimme’s D3 correction64 when computing μC to treat the interlayer van der Waals interactions in graphite. The surface energy can also be expressed as a function of the chemical potential μM: Esurf =

1 (Eslab − NCμ M C (bulk) + (2NC − NM)μM ) 2 2A 0 (6)

In the M-rich limit, μM is the free energy of an M atom within ground-state, body-centered-cubic W or Mo. Table S3 lists all of the calculated energies of the surfaces with different terminations and in the M- or C-rich limits. We accounted for the two types of terminations for the W2C (0001) surface, i.e., the WI and WII terminations as illustrated in Figure 1b. The D

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

Article

The Journal of Physical Chemistry C

We next calculated the surface energies of the nine novel catalysts using the optimized heterostructures as input geometries. As a point of reference for our calculations, we assume that a heterostructure with a metal monolayer on a TMC substrate could be realized by cleaving another heterostructure with N (e.g., N = 8 in this work) metal layers that have already been deposited on the TMC substrate. Figure S1 illustrates the heterostructure and surface slabs of hybrid catalyst and pure metal formed upon cleavage. The surface energy of each hybrid catalyst then can be defined as the energy cost needed to reduce the number of metal layers to a single layer, that is

surface energy of the WI-terminated surface is independent of the W- or C-rich condition because the surface slab is stoichiometric. The W2C surface with the WI-termination has an energy that is much lower (∼1.5 J/m2 lower) than that with the WII termination. Given the geometry depicted in Figure 1b, this finding is not surprising: cleavage at the weakest, longest bonds, the metal−metal bonds, produces this surface. The lowest energies of each surface computed with the three XC functionals were extracted from Table S3 and are displayed in Table 1. Although experimental surface energies are not yet available for these TMCs, our theoretical results agree well with prior DFT values for WC17,65 and Mo2C.66 The TMC surface energies are predicted to be significantly higher than those of the three elemental metals. This comparison is again consistent with the much higher cohesive energies of the TMCs (Table 1). The metal-terminated surfaces of the three TMCs, including the W-terminated WC, W I-terminated W 2C, and Moterminated Mo2C surfaces, are all more energetically favorable than their C-terminated counterparts. We henceforth refer to these surfaces simply as the WC, W2C, and Mo2C surfaces and base our subsequent simulations of hybrid metal monolayer− TMC catalysts on these surfaces as substrates. We next computed the work function, Φ, of the three elemental metals and the three TMCs. Table 1 lists the calculated work function using the three XC functionals. Our results are in excellent agreement with previous calculations, especially for the three metals calculated using the LDA and PBE functionals.60,67 The PBEsol functional most accurately reproduces the work function for Au and Pd, among the three functionals considered. Indeed, the calculated Φs of these two materials are nearly identical to the experimental values measured from photoemission. By contrast, the LDA functional best reproduces the measured work function of Pt, although the difference in pure metal work function values between all three XC functionals is not terribly large ( Pd > Au. Note that this ordering is exactly the same as the trends in surface and cohesive energies already presented. Table 1 also reveals that the TMC work functions are consistently predicted to be smaller than those of the elemental metals. The physical origin of this trend may be that surface dipoles present on the metal-terminated TMC surfaces (metal atoms slightly positively charged) decrease the energy cost to extracting the electron from the material into the vacuum. The above benchmark tests of the LDA, PBE, and PBEsol XC functionals on the bulk and surface properties of the six conventional catalysts show that the PBEsol functional consistently produces results intermediate between those of the LDA and PBE functionals. We therefore expect that the PBEsol functional would result in similarly reliable, intermediate surface properties for the nine novel hybrid catalysts. We thus performed all of our subsequent calculations using the PBEsol functional. Consideration of different stacking patterns of a metal monolayer on top of TMC substrates revealed that the most stable heterostructure always exhibits close-packed stacking. Table S2 lists the layer-stacking sequences and optimized interlayer distances between the metal monolayer and the TMC substrates.

Esurf =

1 (E2L/TMC + 2E7L − E16L/TMC) 4A 0

(7)

where E2L/TMC and E16L/TMC are the energies of the surface slabs with one and eight metal layers deposited on each surface of a TMC substrate slab, respectively. These thicknesses are not arbitrary but are those that lead to convergence of the surface energy in our systems here, similar to how one should evaluate the surface energy of a pure material (i.e., converging it with respect to slab thickness). A factor of 1/2 is implicitly included in these two terms because of the two symmetric interfaces in the supercells. E7L is the energy of a seven-layer surface slab resulting from the cleavage on one side of the hybrid catalyst slab. Equation 7 contains a factor of 1/4 rather than of 1/2 for the surface area term (A0) because four additional surfaces were formed as a result of cleavage on both sides of the hybrid catalyst slab. The surface energy of a heterostructure based on eq 7 is the same as the cleavage or interface energy defined for the interface between two dissimilar systems: one a pure metal and the other a monolayer metal on a TMC substrate. The Supporting Information provides a derivation showing that this definition reduces to the conventional definition of surface energy for a pure material, provided the slab model is sufficiently thick. Numerical tests given in the supplement verify that our models meet this requirement, thus yielding a unique definition of surface energy for our heterostructures. Table 2 reports the surface energies calculated using this definition. Table 2. Surface Energy, Esurf (J/m2); Work Function (eV); and PZC (V) of Hybrid Catalysts Calculated with the PBEsol Functional Au/WC Pd/WC Pt/WC Au/W2C Pd/W2C Pt/W2C Au/Mo2C Pd/Mo2C Pt/Mo2C

Esurf

Φ

PZC

1.241 1.765 2.052 1.208 1.483 1.885 1.143 1.513 1.809

5.43 4.14 5.04 4.82 4.44 4.77 4.59 4.49 4.83

5.09 3.71 4.69 4.49 4.09 4.47 4.28 4.15 4.54

The surface energies of these novel catalysts are significantly smaller than those of pure TMC surfaces, which is to be expected because the addition of a metal monolayer leads to passivation of the dangling bonds at the TMC surfaces. Exposed metal atoms in a monolayer film have a surface energy that is much lower than that of exposed carbon atoms with localized dangling bonds. E

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

Article

The Journal of Physical Chemistry C

compute the free energy of adsorption. Furthermore, expensive supercell calculations are often necessary to account for coverage variations. In contrast, surface energy calculations involve only three calculations on surface unit cells: calculations of the energies of the hybrid system and the two pure materials. The surface energy could therefore provide a more efficient descriptor than the HBE for high-throughput computational searches. This hypothesis was verified by calculating the surface energies of six hybrid electrocatalysts with known current densities, as shown in Figure 2. The experimental exchange current densities for Pt/Mo2C, Pd/W2C, and Au/W2C hybrid catalysts have not yet been measured. The “volcano” relation established in Figure 2 allows us to predict their exchange current densities by relating their computed surface energies to the “volcano” trend line. This trend line suggests that Pt/Mo2C will be a better catalyst than Pd/ W2C or Au/W2C for the HER. Indeed, the surface energy of Pt/ Mo2C indicates that it could be a hybrid catalyst that exhibits activity nearly identical to that of pure Pt. We also explored whether any similar correlation exists between catalytic activity and other surface properties such as the work function and the PZC. This investigation was inspired by the finding several decades ago that the HER’s catalytic activity correlated with the work function of metallic electrodes.69 We therefore computed the work functions of the nine hybrid monolayer metal/TMC catalysts (Table 2). The variation of the experimental exchange current densities with work functions of all 15 materials is displayed in Figure 3. We

Having calculated these 15 surface energies, we next explored the correlation between the measured exchange current densities, i0, and the computed surface energies. Figure 2

Figure 2. Experimental exchange current densities (adopted from ref 18) as a function of calculated surface energies at the DFT PBEsol level. Surface energies for Pt/Mo2C, Pd/W2C, and Au/W2C are shown as dashed lines. The crossing points of these three lines with the trend line of the volcano plot are the corresponding anticipated exchange current densities.

illustrates the variation of i0 (on a logarithmic scale)18 with the computed surface energies. A volcano-shaped correlation between log(i0) and Esurf is clearly evident, similar to the relation between log(i0) and the HBE shown in ref 18. The peak of this volcano plot in this work corresponds to the material with the highest experimental current density (in this case, pure Pt). This peak is then referenced to the surface energy of Pt. The surface energies of the other materials studied in this work are then referenced to this Pt surface energy. The resulting theoretical volcano plot is therefore uniquely specified. Two key features emerge from Figure 2: (i) the catalytic activity peaks around the surface energy of Pt and (ii) surface energies that are either too large (e.g., the TMCs) or too small (e.g., Au) correspond to poor catalytic activities. Higher surface energies generally correspond to higher surface chemical reactivities.68 This is in keeping with the notion of the surface energy being a measure of the degree of coordinative unsaturation of a material’s surface atoms. Those systems with a large number of coordinatively unsaturated atoms will be more reactive than those with fewer such atoms. Coordinatively unsaturated atoms tend to bind intermediates too strongly, inhibiting their subsequent transformation and/or release from the surface. This effect is an illustration of Sabatier’s principle: those surfaces having either too weak or too strong an interaction with reactive intermediates will demonstrate diminished catalytic activity. The extent of this interaction is certainly affected by the degree of coordinative unsaturation of the surface atoms. Surface energy calculations require dramatically reduced computational time compared to HBE calculations because of the former’s reduced complexity. Rigorous HBE calculations necessitate identification of the most stable binding site for H on a surface. Once this site has been determined, the binding energy is then calculated. Frequency calculations must also be performed to ensure that a stable minimum was found and to

Figure 3. Experimental exchange current densities (adopted from ref 18) as a function of calculated work function Φ and PZC. Filled and empty symbols denote Φ and PZC, respectively. Square symbols refer to Φ and PZC of the three catalysts (Pt/Mo2C, Pd/W2C, and Au/ W2C) with their predicted exchange current densities based on Figure 2 intersection points.

used the predicted exchange current densities that were extracted from Figure 2 for Pt/Mo2C, Pd/W2C, and Au/ W2C. The scatter in the data demonstrates no strong correlation exists between current density and work function. Finally, we studied how solvation affects the surface energy and the work function. Tables S1 and S2 list the predicted surface energies of all 15 catalysts in aqueous solution. Because these surface energies are almost the same as they are in vacuum, the liquid environment appears to only trivially affect the surface energies. As our solvation model is just a polarizable F

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

Article

The Journal of Physical Chemistry C continuum, our conclusion applies only to surfaces in which the electrode−electrolyte interface is an “innocent” one, in which water coexists but does not react or dissociate on the electrode surface. For such interfaces, our volcano-shaped relation should remain valid. This is in keeping with the experimental finding that these hybrid systems are highly robust under electrochemical cycling. This inertness suggests that these surfaces are negligibly affected by interaction with the bulk electrolyte. Tables S1 and S2 also once again show that the solvated surface energies calculated with VASP are similar to the results calculated by JDFTx. We also examined the potential correlation of the catalytic activity with another important surface property of catalysts in solution, namely, the PZC, which can be regarded as the “work function” of an electrode immersed in a liquid. This “work function” in solution is equal to the work function in vacuum corrected by a potential drop at the electrode/liquid interface. In principle, we can use a periodically infinite surface slab with implicit solvent to calculate the PZC. However, convergence of the calculation is slow. In practice, implicit cations/anions are added to the simulation supercells by setting their ionic concentrations, as implemented in both JDFTx and VASP. These ions induce Debye screening in the simulation supercell, which guarantees that the electrostatic potential rapidly decays to zero.70 In this work, we set the total concentration of monovalent ions to 1 M (0.5 M H+ and 0.5 M HSO4−) to model the 0.5 M H2SO4 solution used in the experiment.18 To demonstrate the effect of ions on the electrostatic potential, Figure S2 displays the electrostatic potential profiles for the WC(0001) surface in the absence and presence of ions in the simulations. We thus verify that the electrostatic potential is indeed zero outside the electrode if the ions are present in the simulation. This zero potential far away represents an absolute reference for the electrostatic potential. The added ions also lead to a determinate electrostatic potential, and the PZC is the negative of the calculated Fermi level of a solvated system as a result. The calculated PZC of the metals and TMCs and of the hybrid catalysts are reported in Tables 1 and 2, respectively. Figure 3 also displays the variation of the experimental exchange current density with the calculated PZC. We observe a nearly constant shift of the PZC with reference to the work function, i.e., the PZC can simply be written as71 PZC = Φ/e − c0. We determined the average shift c0 to be 0.36 V with a standard deviation of 0.07 V. Similar to the work function, the PZC exhibits no volcano-shaped correlation with the exchange current densities; therefore, neither electronic propertythe work function or the PZCis as good a descriptor of the catalytic activity as is the surface energy.

hybrid catalysts, each of which consists of one metal monolayer on a TMC substrate. We explored the correlation between experimental catalytic activity and three surface properties of 15 traditional and novel catalysts: the surface energy, work function, and PZC. These properties were selected for study because of their close relationship with many properties thought to influence catalytic activity. We discovered a strong volcano-shaped correlation between the catalytic activity and the surface energy, indicating that the surface energy can be used as a novel, efficient descriptor of catalytic activity. This simple descriptor is successful because it reflects the degree of coordinative unsaturation of the surface atoms and thus the reactivity of the surface. Further validation of this new descriptor is underway; HBE calculations to be reported elsewhere on, e.g., Pt/MoC2, confirm its placement near the peak of the volcano and its potential high activity for HER. We further expect this descriptor could be useful beyond HER, to other catalytic reactions such as the oxygen reduction reaction, but verification of this speculation is the subject of future work. We also performed electrochemical simulations for the 15 materials based on JDFT. We found that the surface energy in vacuum of these 15 materials is essentially unchanged by the presence of an aqueous solution. This suggests that the vacuum surface energy is a sufficiently accurate metric, provided the water at the interface remains undissociated. Finally, we predicted the PZCs for most of the 15 catalysts, all of which await experimental verification. We also confirmed the linear relationship between the PZC and work function. However, unlike the surface energy, there exists no obvious correlation between the catalytic activity and either the work function or the PZC. Our work thus identified an instructive indicator of catalytic activity of the HER. We therefore propose that the surface energy can be used as a descriptor to search for Ptefficient catalysts exhibiting surface energies similar to that of the pure Pt surface.

CONCLUSIONS We employed DFT to study the bulk and surface properties of six common catalysts, including three elemental metals (Au, Pd, and Pt) and three TMCs (WC, W2C, and Mo2C) using the LDA, PBE, and PBEsol XC functionals. We showed that the PBEsol functional consistently produces properties between those of the LDA and PBE functionals and generally closest to those of experiment. We therefore expect the PBEsol functional to partially compensate for the shortcomings of the LDA and PBE functionals. Indeed, the PBEsol functional gives rise to lattice constants and work functions in excellent agreement with experiment. We therefore employed only the PBEsol functional to further investigate surface properties of nine novel





ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b09687. Surface energies and PZCs of six common catalysts computed with both VASP and JDFTx; stacking sequences, interlayer spacings, surface eneriges, and PZCs of nine novel catalysts; surface energies of TMCs with different surface terminations and in metal-rich or carbon-rich limits; derivation of equivalences between surface energy expressions and numerical validation; electrostatic profiles of WC(0001) surface with and without the presence of implicit ions (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We are grateful to the Air Force Office of Scientific Research for funding (Grant FA9550-14-1-0254). We acknowledge use of the TIGRESS high-performance computer center at Princeton University. We also acknowledge use of the G

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

Article

The Journal of Physical Chemistry C

(19) Kelly, T. G.; Chen, J. G. Metal overlayer on metal carbide substrate: unique bimetallic properties for catalysis and electrocatalysis. Chem. Soc. Rev. 2012, 41, 8021. (20) Weidman, M. C.; Esposito, D. V.; Hsu, Y.-C.; Chen, J. G. Comparison of electrochemical stability of transition metal carbides (WC, W2C, Mo2C) over a wide pH range. J. Power Sources 2012, 202, 11. (21) Kimmel, Y. C.; Esposito, D. V.; Birkmire, R. W.; Chen, J. G. Effect of surface carbon on the hydrogen evolution reactivity of tungsten carbide (WC) and Pt-modified WC electrocatalysts. Int. J. Hydrogen Energy 2012, 37, 3019. (22) Esposito, D. V.; Chen, J. G. Monolayer platinum supported on tungsten carbides as low-cost electrocatalysts: opportunities and limitations. Energy Environ. Sci. 2011, 4, 3900. (23) Kelly, T. G.; Hunt, S. T.; Esposito, D. V.; Chen, J. G. Monolayer palladium supported on molybdenum and tungsten carbide substrates as low-cost hydrogen evolution reaction (HER) electrocatalysts. Int. J. Hydrogen Energy 2013, 38, 5638. (24) Sabatier, P. La Catalyse en Chemie Organique; Béranger: Paris et Liége, 1920. (25) Parsons, R. The rate of electrolytic hydrogen evolution and the heat of adsorption of hydrogen. Trans. Faraday Soc. 1958, 54, 1053. (26) Parsons, R. Hydrogen evolution on platinum electrodes. The heats of activation for the component reactions. Trans. Faraday Soc. 1960, 56, 1340. (27) Nørskov, J. K.; Bligaard, T.; Logadottir, A.; Kitchin, J. R.; Chen, J. G.; Pandelov, S.; Stimming, U. Trends in the Exchange Current for Hydrogen Evolution. J. Electrochem. Soc. 2005, 152, J23. (28) Greeley, J.; Jaramillo, T. F.; Bonde, J.; Chorkendorff, I. B.; Norskov, J. K. Computational high-throughput screening of electrocatalytic materials for hydrogen evolution. Nat. Mater. 2006, 5, 909. (29) Jaramillo, T. F.; Jørgensen, K. P.; Bonde, J.; Nielsen, J. H.; Horch, S.; Chorkendorff, I. Identification of Active Edge Sites for Electrochemical H2 Evolution from MoS2 Nanocatalysts. Science 2007, 317, 100. (30) Quaino, P.; Juarez, F.; Santos, E.; Schmickler, W. Volcano plots in hydrogen electrocatalysis − uses and abuses. Beilstein J. Nanotechnol. 2014, 5, 846. (31) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169. (32) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953. (33) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758. (34) Ceperley, D. M.; Alder, B. J. Ground State of the Electron Gas by a Stochastic Method. Phys. Rev. Lett. 1980, 45, 566. (35) Perdew, J. P.; Zunger, A. Self-interaction correction to densityfunctional approximations for many-electron systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 5048. (36) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (37) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett. 2008, 100, 136406. (38) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188. (39) Methfessel, M.; Paxton, A. T. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 40, 3616. (40) Sansonetti, J. E.; Martin, W. C. Handbook of Basic Atomic Spectroscopic Data. J. Phys. Chem. Ref. Data 2005, 34, 1559. (41) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Insights into Current Limitations of Density Functional Theory. Science 2008, 321, 792. (42) Fattebert, J.-L.; Gygi, F. First-principles molecular dynamics simulations in a continuum solvent. Int. J. Quantum Chem. 2003, 93, 139.

COPPER supercomputer at the Air Force Office of Scientific Research High Performance Computing Center. We thank Dr. Kiran Matthew, Dr. Ravishankar Sundararaman, and Prof. Richard Hennig for helpful discussions on the VASPsol and JDFTx calculations. We also thank Ms. Nari Baughman and Dr. Johannes M. Dieterich for critical reading of this manuscript.



REFERENCES

(1) Sharma, S.; Ghoshal, S. K. Hydrogen the future transportation fuel: From production to applications. Renewable Sustainable Energy Rev. 2015, 43, 1151. (2) Balat, M. Potential importance of hydrogen as a future solution to environmental and transportation problems. Int. J. Hydrogen Energy 2008, 33, 4013. (3) Bakenne, A.; Nuttall, W.; Kazantzis, N. Sankey-Diagram-based insights into the hydrogen economy of today. Int. J. Hydrogen Energy 2016, 41, 7744. (4) Nicholson, W.; Carlisle, A. Philos. Mag. Ser. 1 1800, 7, 337. (5) Trasatti, S. Electrocatalysis by oxides  Attempt at a unifying approach. J. Electroanal. Chem. Interfacial Electrochem. 1980, 111, 125. (6) Raj, I. A. Nickel-based, binary-composite electrocatalysts for the cathodes in the energy-efficient industrial production of hydrogen from alkaline-water electrolytic cells. J. Mater. Sci. 1993, 28, 4375. (7) Fosdick, S. E.; Berglund, S. P.; Mullins, C. B.; Crooks, R. M. Evaluating Electrocatalysts for the Hydrogen Evolution Reaction Using Bipolar Electrode Arrays: Bi- and Trimetallic Combinations of Co, Fe, Ni, Mo, and W. ACS Catal. 2014, 4, 1332. (8) Hinnemann, B.; Moses, P. G.; Bonde, J.; Jørgensen, K. P.; Nielsen, J. H.; Horch, S.; Chorkendorff, I.; Nørskov, J. K. Biomimetic Hydrogen Evolution: MoS2 Nanoparticles as Catalyst for Hydrogen Evolution. J. Am. Chem. Soc. 2005, 127, 5308. (9) Li, Y.; Wang, H.; Xie, L.; Liang, Y.; Hong, G.; Dai, H. MoS2 Nanoparticles Grown on Graphene: An Advanced Catalyst for the Hydrogen Evolution Reaction. J. Am. Chem. Soc. 2011, 133, 7296. (10) Fernández, E. M.; Moses, P. G.; Toftelund, A.; Hansen, H. A.; Martínez, J. I.; Abild-Pedersen, F.; Kleis, J.; Hinnemann, B.; Rossmeisl, J.; Bligaard, T.; Nørskov, J. K. Scaling Relationships for Adsorption Energies on Transition Metal Oxide, Sulfide, and Nitride Surfaces. Angew. Chem., Int. Ed. 2008, 47, 4683. (11) Vrubel, H.; Hu, X. Molybdenum Boride and Carbide Catalyze Hydrogen Evolution in both Acidic and Basic Solutions. Angew. Chem. 2012, 124, 12875. (12) Popczun, E. J.; McKone, J. R.; Read, C. G.; Biacchi, A. J.; Wiltrout, A. M.; Lewis, N. S.; Schaak, R. E. Nanostructured Nickel Phosphide as an Electrocatalyst for the Hydrogen Evolution Reaction. J. Am. Chem. Soc. 2013, 135, 9267. (13) Armstrong, R. D.; Bell, M. F. Tungsten carbide catalysts for hydrogen evolution. Electrochim. Acta 1978, 23, 1111. (14) Hwu, H. H.; Chen, J. G. Surface Chemistry of Transition Metal Carbides. Chem. Rev. 2005, 105, 185. (15) Meyer, S.; Nikiforov, A. V.; Petrushina, I. M.; Köhler, K.; Christensen, E.; Jensen, J. O.; Bjerrum, N. J. Transition metal carbides (WC, Mo2C, TaC, NbC) as potential electrocatalysts for the hydrogen evolution reaction (HER) at medium temperatures. Int. J. Hydrogen Energy 2015, 40, 2905. (16) Esposito, D. V.; Hunt, S. T.; Stottlemyer, A. L.; Dobson, K. D.; McCandless, B. E.; Birkmire, R. W.; Chen, J. G. Low-Cost HydrogenEvolution Catalysts Based on Monolayer Platinum on Tungsten Monocarbide Substrates. Angew. Chem. 2010, 122, 10055. (17) Zhuang, H.; Tkalych, A. J.; Carter, E. A. Understanding and Tuning the Hydrogen Evolution Reaction on Pt-Covered Tungsten Carbide Cathodes. J. Electrochem. Soc. 2016, 163, F629. (18) Esposito, D. V.; Hunt, S. T.; Kimmel, Y. C.; Chen, J. G. A New Class of Electrocatalysts for Hydrogen Production from Water Electrolysis: Metal Monolayers Supported on Low-Cost Transition Metal Carbides. J. Am. Chem. Soc. 2012, 134, 3025. H

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

Article

The Journal of Physical Chemistry C (43) Petrosyan, S. A.; Rigos, A. A.; Arias, T. A. Joint DensityFunctional Theory: Ab Initio Study of Cr2O3 Surface Chemistry in Solution. J. Phys. Chem. B 2005, 109, 15436. (44) Andreussi, O.; Dabo, I.; Marzari, N. Revised self-consistent continuum solvation in electronic-structure calculations. J. Chem. Phys. 2012, 136, 064102. (45) Schwarz, K. A.; Sundararaman, R.; Letchworth-Weaver, K.; Arias, T. A.; Hennig, R. G. Framework for solvation in quantum Monte Carlo. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 201102. (46) Gunceler, D.; Letchworth-Weaver, K.; Sundararaman, R.; Schwarz, K. A.; Arias, T. A. The importance of nonlinear fluid response in joint density-functional theory studies of battery systems. Modell. Simul. Mater. Sci. Eng. 2013, 21, 074005. (47) Mathew, K.; Sundararaman, R.; Letchworth-Weaver, K.; Arias, T. A.; Hennig, R. G. Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways. J. Chem. Phys. 2014, 140, 084106. (48) Ismail-Beigi, S.; Arias, T. A. New algebraic formulation of density functional calculation. Comput. Phys. Commun. 2000, 128, 1. (49) Sundararaman, R. JDTFx. http://jdftx.sourceforge.net. (50) Garrity, K. F.; Bennett, J. W.; Rabe, K. M.; Vanderbilt, D. Pseudopotentials for high-throughput DFT calculations. Comput. Mater. Sci. 2014, 81, 446. (51) Jones, R. O.; Gunnarsson, O. The density functional formalism, its applications and prospects. Rev. Mod. Phys. 1989, 61, 689. (52) Davey, W. P. Precision Measurements of the Lattice Constants of Twelve Common Metals. Phys. Rev. 1925, 25, 753. (53) de Boer, F. R.; Boom, R.; Mattens, W. C. M.; Miedema, A. R.; Niessen, A. K. Cohesion in Metals; North Holland: Amsterdam, 1988. (54) Kittel, C. Introduction to Solid State Physics; Wiley: Hoboken, NJ, 2004. (55) Hansson, G. V.; Flodström, S. A. Photoemission study of the bulk and surface electronic structure of single crystals of gold. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 18, 1572. (56) Murata, Y.; Starodub, E.; Kappes, B. B.; Ciobanu, C. V.; Bartelt, N. C.; McCarty, K. F.; Kodambaka, S. Orientation-dependent work function of graphene on Pd(111). Appl. Phys. Lett. 2010, 97, 143114. (57) Derry, G. N.; Ji-Zhong, Z. Work function of Pt(111). Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 39, 1940. (58) Page, K.; Li, J.; Savinelli, R.; Szumila, H. N.; Zhang, J.; Stalick, J. K.; Proffen, T.; Scott, S. L.; Seshadri, R. Reciprocal-space and realspace neutron investigation of nanostructured Mo2C and WC. Solid State Sci. 2008, 10, 1499. (59) Pierson, H. O. 6 - Carbides of Group VI: Chromium, Molybdenum, and Tungsten Carbides. In Handbook of Refractory Carbides and Nitrides; William Andrew Publishing: Westwood, NJ, 1996; p 100. (60) Singh-Miller, N. E.; Marzari, N. Surface energies, work functions, and surface relaxations of low-index metallic surfaces from first principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 235407. (61) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Tests of a ladder of density functionals for bulk solids and surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 075102. (62) Michaelides, A.; Scheffler, M. An Introduction to the Theory of Crystalline Elemental Solids and their Surfaces. In Surface and Interface Science; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2014; p 13. (63) Reuter, K.; Scheffler, M. Composition, structure, and stability of RuO2(110) as a function of oxygen pressure. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 65, 035406. (64) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (65) Siegel, D. J.; Hector, L. G., Jr; Adams, J. B. Adhesion, stability, and bonding at metal/metal-carbide interfaces: Al/WC. Surf. Sci. 2002, 498, 321.

(66) Politi, J. R. d. S.; Vines, F.; Rodriguez, J. A.; Illas, F. Atomic and electronic structure of molybdenum carbide phases: bulk and low Miller-index surfaces. Phys. Chem. Chem. Phys. 2013, 15, 12617. (67) Fall, C. J.; Binggeli, N.; Baldereschi, A. Work-function anisotropy in noble metals: Contributions from d states and effects of the surface atomic structure. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 61, 8489. (68) Totten, G. E.; Linag, H. Surface Modification and Mechanisms: Friction, Stress, and Reaction Engineering, 1st ed.; CRC Press: Boca Raton, FL, 2004. (69) Conway, B. E.; Bockris, J. O. M. The d-Band Character of Metals and the Rate and Mechanism of the Electrolytic Hydrogen Evolution Reaction. Nature 1956, 178, 488. (70) Letchworth-Weaver, K.; Arias, T. A. Joint density functional theory of the electrode-electrolyte interface: Application to fixed electrode potentials, interfacial capacitances, and potentials of zero charge. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 075140. (71) Trasatti, S. Work function, electronegativity, and electrochemical behaviour of metals. J. Electroanal. Chem. Interfacial Electrochem. 1971, 33, 351.

I

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