Framework for Scalable Adsorbate–Adsorbate Interaction Models

Jun 2, 2016 - We present a framework for physically motivated models of adsorbate–adsorbate interaction between small molecules on transition and ...
0 downloads 0 Views 10MB Size
Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Article

A Framework for Scalable Adsorbate-adsorbate Interaction Models Max J Hoffmann, Andrew J Medford, and Thomas Bligaard J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b03375 • Publication Date (Web): 02 Jun 2016 Downloaded from http://pubs.acs.org on June 13, 2016

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

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

Page 1 of 24

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

The Journal of Physical Chemistry

A Framework for Scalable Adsorbate-adsorbate Interaction Models Max J Hoffmann,∗,†,‡ Andrew J Medford,†,‡ and Thomas Bligaard†,‡ †Department of Chemical Engineering, Stanford University, Stanford, CA 94305, USA ‡SUNCAT Center for Interface Science and Catalysis, SLAC, National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA E-mail: [email protected] Phone: +1(650) 926-2233. Fax: +1(650) 926 4100

Abstract We present a framework for physically motivated models of adsorbate-adsorbate interaction between small molecules on transition and coinage metals based on modifications to the substrate electronic structure due to adsorption. We use this framework to develop one model for transition and one for coinage metal surfaces. The models for transition metals are based on the d-band center position and the models for coinage metals are based on partial charges. The models require no empirical parameters, only two first-principles calculations per adsorbate as input, and therefore scale linearly with the number of reaction intermediates. By theory to theory comparison with explicit density functional theory calculations over a wide range of adsorbates and surfaces we show that the root-mean squared error for differential adsorption energies is less than 0.2 eV for up to 1 ML coverage.

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Introduction For many surface chemical reactions adsorbate-adsorbate interactions can have significant influence on the rate of every reaction step.1–4 Adsorbate-adsorbate interactions can be expected to become particularly influential under operating conditions when high coverage of adsorbates implies low average separation between adsorbates. For an accurate prediction of the intrinsic kinetics and concomitantly the reactivity and selectivity of a catalyst it is therefore desirable to account for adsorbate-adsorbate interactions. Previously, in the context of first-principles derived microkinetic models significant advances have been made in the form of physically guided bond-order conservation models5 , a surface electronic structure mediated model for transition metals based on the d-band width by İnoğlu and Kitchin6 , linear and higher-order parametrizations as a function of mean-field coverages7 or the very general lattice gas Hamiltonian approach by deploying an expansion into adsorbates clusters of pairs, triplets, quadruplets, etc.8–10 The parametrization of lattice gas Hamiltonians via adsorbate cluster expansion can be applied irrespective of the underlying physical cause of interaction and cross-validation benchmarks of the training set allow for a systematic improvement of its accuracy. In particular the contributions from Schneider and coworkers have raised the bar to describe the interaction between one or two surface species with an impressive uncertainty of only 5 meV/site when compared to the corresponding electronic structure calculation.10 The generality and accuracy of the cluster expansion often makes it a natural first choice. However, when aiming for an atomistic microkinetic model of more complex, yet highly relevant, catalytic reactions such as CO methanation, partial methane oxidation, or the Fischer-Tropsch reaction with tens or even hundreds of different surface species interacting with one another the cluster expansion approach is fundamentally limited by the exuberant number of required calculations. Three or more adsorbate species is nowadays a non-trivial problem size for modern supercomputing capacities available for electronic structure calculations.10–12 2

ACS Paragon Plus Environment

Page 2 of 24

Page 3 of 24

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

The Journal of Physical Chemistry

The problem is exacerbated if the catalyst surface contains more than one adsorption site. The cluster expansion approach becomes impractical with increasing complexity because the number of possible clusters that have to be calculated grows in leading order at least as p2 , where p is the number of adsorbate species.13 While maintaining a large degree of generality the cluster-expansion approach typically relies only on one piece of information, the total energy, from the multitude of information that a converged DFT calculation of a particular arrangement of adsorbates on the surface offers. On the other end of the spectrum are more empirical approaches such as bond-order conservation models15 . These require a significant amount of intuition or experimental input about the interaction strength of a particular adsorbate. Furthermore the expressions for the adsorption energy are typically derived for specific adsorption geometries (e.g. hollow, top) instead of first-principles observables and there exists no straightforward way to transfer an expression from one adsorption site to another. The approach pursued here aims to better balance the generality and accuracy of a cluster expansion with the computational cost of more empirical approaches. The approach builds on the presupposition that nowadays in the context of first-principles derived kinetic model calculating adsorption energies at two well-defined coverages is of negligible cost compared to the development of a kinetic reaction network and identification of transition states. Cluster expansion approaches, on the other hand, quickly outgrow this cost if applied to three or more different species. One way of avoiding this prohibitive growth in the number of calculations is to separate the adsorbate-adsorbate interactions into two steps by explicitly accounting for auxiliary physical properties: what effect has an adsorbate on the substrate and what effect does such a modified substrate have on other adsorbates in the vicinity. This is the approach pursued here for the subclass of problems of first- and second row adsorbates on late transition and coinage metal16 surfaces. One would expect such an approach to work well if the interaction between adsorbates is mediated primarily by the surface eletronic structure. Using this general framework we develop physically motivated models in this paper that 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

allow quantitative predictions of adsorbate-adsorbate interactions from just a fixed number of first-principles calculations per species and site, i.e. one low coverage and one high coverage geometry. The total number of electronic structure simulations necessary to parametrize the model thus scales linearly with the number of adsorbates, making the computational cost of utilizing this approach in theory-based microkinetic studies a small constant factor more expensive than a similar study without inclusion of adsorbate-adsorbate interaction. The remainder of the paper is structured as follows: first we define a general linearized framework that allows for the aforementioned separation of adsorbate-adsorbate interactions. Then, by studying trends across adsorbates we develop suitable properties to be calculated for each surface atom that can act as faithful physical proxies in such an interaction model. Finally we validate the outcome of the interaction model by comparing explicit DFT calculations of mixed adsorbate-adsorbate interactions with the corresponding model.

Theoretical Background Defining the Problem The adsorbate-adsorbate interaction energy can be defined as the total energy of one or more adsorbates (a at site i, b at site j, c at site k) in close vicinity on a surface relative to the total energy of the same adsorbates but infinitely far away from each other. The central property this paper aims to predict is the interaction energy between two or more adsorbates on a periodic metal surface

Eint (ai, bj, ck, . . . ) = Eai,bj,ck − (Eai + Ebj + Eck + . . . ). where Eai,bj,ck is the adsorption energy of adsorbates in some close arrangement and Eai , Ebj , Eck are the adsorption energies of the same individual adsorbates but in the long distance limit from any other adsorbate. This adsorbate-adsorbate interaction is by construction 4

ACS Paragon Plus Environment

Page 4 of 24

Page 5 of 24

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

The Journal of Physical Chemistry

independent of any gas phase reference state. While predicting interactions between adsorbates for arbitrary co-adsorbate configurations seems like a daunting task, several assumptions can be made in the context of microkinetic models that allow the significant reduction of the requirements of such a model. First of all, adsorbates on single crystal surfaces adsorb on high symmetry sites. Thus an interaction model does not have to yield interaction energies for arbitrary positions of the adsorbate above a surface but only a discrete set of adsorption sites. Furthermore, for most surface species the highest expected coverage is 1 monolayer (ML) (by which we shall here mean one adsorbate per surface atom). This implies that atomic adsorbates are typically separated by more than one unit cell width which could help reduce the complexity of the electronic adsorbate-adsorbate interaction. Finally, in order to reduce the geometrical complexity of the problem all adsorbates are held in an upright orientation and adsorbate atom(s) can only relax along the z-direction while the two upper layers of the substrate can relax in all directions. In particular for adsorbate-adsorbate interactions involving OH a more detailed description that is beyond the scope of this paper maybe appropriate. A control sample geometry of OH and O in one case and N and CO in another on a (2x1) unit cell of Rh(111) show that full relaxation of adsorbates lowers the energy by 0.56 eV and 0.06 eV respectively due to tilting.

Different Causes of Adsorbate-adsorbate Interaction This section shall briefly summarize the different aspects and identify the most promising route. The aim is to predict adsorbate-adsorbate interaction from basic physical properties. The exact interaction of adsorbates in a given arrangement on a surface in some proximity is only captured by the exact solution of the full quantum mechanical equations. Different causes may govern the interaction for different combinations of adsorbates and substrates. It may prove cumbersome or even impossible to exactly decompose contributions from each

5

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

aspect since they are all in some way related within the self-consistent ground state solution. Hammer et al.17 divide the interaction between adsorbates on transition metals into four significant causes: direct interactions due to overlap of adsorbate wave functions, indirect interactions mediated by the electronic structure of the substrate, elastic interactions via displacements of surface atoms, and electrostatic interactions. The direct overlap of orbitals is most significant at very short adsorbate-adsorbate separations. At separations typical for nearest neighbor sites on fcc(111) metal surface (2.5-2.8Å) the direct overlap between adsorbate orbitals is expected to have only a minor influence except perhaps for more delocalized molecular orbitals as for example in sulfur or alkali metals18 . Elastic interactions are generally difficult to treat in a lattice picture since they typically involve directional preference that breaks from the lattice geometry. Before capturing elastic distortion in great detail, methodological developments towards how to include such displacements in microkinetic models would be a necessary prerequisite and are out of the scope of this paper. Electrostatic interactions describe the notion that an adsorbate-surface system can have significant electrostatic dipoles or higher-order multipoles which can interact with other adsorbate induced dipoles. Another notion is that an adsorbate creates an electric field normal to the surface much like a capacitor which is felt by other adsorbates. In practice it proves to be rather difficult to decompose such an interaction in a transferable way. In a simple way one can view two approaching adsorbates as two dipoles interacting with each other. However, since adsorbates in typical nearest-neighbor geometries share substrate atoms a clear separation between non-local electrostatic interaction and surface electronic structure mediated interaction is not guaranteed. If the creation of adsorbate induced dipoles is due to charge depletion from surface atoms the effect of overlapping charge depletion may be the dominating effect which connects the indirect electrostatic mechanism to a surface mediated electronic structure problem. The d-band model builds upon the Newns-Anderson

6

ACS Paragon Plus Environment

Page 6 of 24

Page 7 of 24

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

The Journal of Physical Chemistry

model, which does not take variations in the electrostatics into account. The d-band model itself, therefore also does not directly include electrostatic effects. The lowest-order version of the d-band model suggest a linear variation of adsorption strength with the position of the d-band with respect to the Fermi level. Even though this model in principle does not include electrostatic effects, the electrostatic contributions might in many cases also scale reasonably linearly with the Newns-Anderson hybridization strength, and monitoring the the total adsorption strength as a function of d-band position will therefore tend to include contributions from both electrostatics and the one-electron contributions. This mixing of effects might be a central reason for the success of the low-order versions of the d-band model. In the present study, however, as we parametrize the one-electron contribution to the adsorption strength variations as the variation of the adsorption strength with the variation in d-band center, the d-band center variation already contain a significant electrostatic contribution, which makes it impossible in the model to fully decompose the one-electron and the electrostatic contributions to the adsorbate-adsorbate interactions. Furthermore, for adsorbates consisting of more than one atom one faces the problem of how internal degrees of freedom should be handled in an adsorbate-adsorbate interaction model. In a closely-packed co-adsorption situation it is to be expected that relaxation of internal degrees of freedom of adsorbates systematically leads to a decrease of effective adsorbate-adsorbate interaction since relaxation of adsorbate internal degrees of freedom systematically leads for a better accommodation of neighboring adsorbates. In light of the already non-trivial level of theory proposed so far and the fact that these internal degrees of freedom will be quite cumbersome to account for in the next level of description such as microkinetic modeling we opt to keep all adsorbate atoms fixed at their low-coverage position with respect to the xand y-direction and only relax those in the z-direction. In particular for adsorbate-adsorbate interactions involving OH a more detailed description that is beyond the scope of this paper maybe appropriate. A control sample geometry of OH and O in one case and N and CO in another on a (2x1) unit cell of Rh(111) show that full relaxation of adsorbates lowers the 7

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

energy by 0.56 eV and 0.06 eV respectively due to tilting. We expect that the largest single contribution originates from surface electronic structure mediated interaction and moreover that the overall interaction is more or less proportional to this contribution. Hence focusing on electronic modifications is the most promising starting point.

Linear Surface Mediated Interaction Model We shall now motivate the principles and formalism of the surface mediated adsorbateadsorbate interaction framework. The idea behind a surface mediated adsorbate-adsorbate interaction framework is to separate the interaction into two steps: in the first step we describe the effect of an adsorbate on the electronic structure of the surfaces. In the second step we describe how this change of the electronic structure affects the adsorption energy of other adsorbates in the vicinity. The challenge of such an approach is that suitable auxiliary properties of the surface electronic structure need to be identified which can be superimposed. However, if a simple auxiliary property can be identified and assigned to each surface atom the identity of the adsorbate causing the change is no longer required to predict the adsorption energy of the affected adsorbate. In the following, generic adsorbates will be denoted by letters a, b; generic adsorption sites will be denoted by letters i, j, k; and generic surface atoms by m, n. For better readability adsorbates will be superscripts, while surface atoms and sites will be subscripts. Consider Fig. 1 displaying 3 oxygen atoms labeled A, B, C on an fcc(111) surface. Surface atoms are labeled using numbers 1, 2, 3 . . . and adsorption fcc/hcp sites are labeled 1’, 2’, 3’, . . . . Suppose we would like to predict the change in adsorption energy of oxygen atom A at site 21’ due to co-adsorbates B and C at sites 22′ and 11′ , respectively. We assume that for every adsorbate species a and every adsorption site i we have a measure for how strongly a the adsorbate interacts with every surface atom m denoted by V˜mi . The next section will

8

ACS Paragon Plus Environment

Page 8 of 24

Page 9 of 24

The Journal of Physical Chemistry

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

The Journal of Physical Chemistry

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

Page 10 of 24

as (1) Vmi =

X b

b b,max j6=i δj Vmj P ˜b k6=i Vmk

P

b V˜mj

and (2) δjb =

    1

   0

if adsorbate b occupies site j else.

In the given example the sum over b would only include b = O since oxygen is the only species in our system. The sum over sites would run over all sites except 21′ since we are not including the interaction of an adsorbate with itself. Without stating explicit values we would expect among all Vm,21′ the largest modifications at V12,21′ since atom 12 is interacting with both B and C. This will be followed by V7,21′ and V13,21′ since they may still be under the influence of both adsorbates nearby. This will be followed by V6,21′ and V19,21′ and all other surface atoms located further from B and C. Having a model for the modification of the surface electronic structure we shall next translate given modifications into a change in adsorption energy. To this end we need to fix an energy scale or conversion factor that converts a certain electronic structure modification into a change in adsorption energy. A canonical choice here is to use the energy difference in adsorption energy between a low coverage configuration (e.g. one adsorbate on a surface slab of (4 × 4) unit cells) denoted by Eia,ads−lo and the adsorption energy Eia,ads−hi at the highest attainable coverage (e.g. one adsorbate in a surface slab of (1 × 1) unit cells). We define this conversion factor as (3) 



∆Eia = Eia,ads−hi − Eia,ads−lo . 10

ACS Paragon Plus Environment

Page 11 of 24

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

The Journal of Physical Chemistry

Finally we assume that the influence of changes in the surface electronic structure on the a adsorption energy can be weighted by the same measure V˜mi , the interaction strength between

an adsorbate and the electronic structure of the surface atom. In the given example we A therefore sum over all non-zero elements of V˜m,21 ′ . The weighted average can then be written

as (4) Eia,model

=

∆Eia

Vmi ˜ a m V a,max Vmi mi

P

P

m

a V˜mi

.

Using this general framework the key missing ingredient is then to identify suitable electronic structure properties that can be extracted from reference first-principles calculations for specific combinations of substrates and adsorbates to populate the input parameters V˜ and V max .

Technical Settings The density functional theory (DFT) calculations in the following have been made using slab geometries of fcc(111) from Rh, Ir, Pd, Pt, Cu, and Ag consisting of 4 layers with the bottom 2 layers fixed at their bulk positions using the BEEF-vdW as xc functional as implemented in Quantum Espresso.19,20,21,22,23 The plane-wave cut-off energy was set to 500 eV, the density integration cut-off was set to 5000 eV and the non-centered k-point grid was set to commensurate sizes corresponding to (12 × 12 × 1) for the (1 × 1) slab. That is a (6 × 6 × 1) grid was used for the (2 × 2) slab, a (4 × 4 × 1) grid for the (3 × 3) slab, and a (3 × 3 × 1) grid for the (4 × 4) slab. For the self-consistent field convergence an energy threshold of 10−7 eV and for the geometry optimization a force threshold of 0.03 eV/A was chosen. The density of states projections were performed using a non self-consistent calculation with twice the number of k-points each axis. The density of states was sampled over an energy window of ±500 eV throughout this study with a step width of 0.01 eV. A 11

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 12 of 24

Page 13 of 24

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

The Journal of Physical Chemistry

energy is only approximately 0.2 eV while it is on the order of 1.0 eV for all other adsorbates. The conclusion is that for rhodium the d-band modifications can act as a reliable proxy for changes in the adsorption energy. Accordingly we identify V˜ with the shift of the d-band center of each surface atom in the (4 × 4) surface adsorption on transition metal surfaces. Furthermore we identify V max with the maximum shift of the d-band center of the surface atom at 1 ML coverage.

Charge Mediated Interactions

Figure 3: Density differences of an oxygen atom adsorbed on Rh(111) (left) and Cu(111) (right). The isosurface is drawn at a density of 3 × 10−5 e|Å|−3 . Blue indicates charge accumulation while red indicates charge depletion. The top row shows a top view of the geometry while the bottom row shows a side view of the same isosurface. The Rh atoms coordinated by oxygen display a charge reorganisation but the excess and depletion volume appear to be of similar size. On the other hand the Cu atom directly coordinated by oxygen feature a net depletion of charge. In the case of coinage metals the d-band is not a faithful descriptor of interaction trends as the d-band is not pinned to the Fermi level when the d-band is full. Instead we turn to qualitative differences in the rearrangement of charge upon adsorption between coinage and 13

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

0.01 0.05

0.04

0.02

0.02

0.05

0.04

0.01

0.04

0.04

0.00

0.07

-0.02

-0.08

0.02

-0.07

-0.06

0.02

-0.01

-0.04

0.05

-0.04

-0.01

0.04

-0.02

-0.08

0.02

0.01

0.01

0.04

0.01

0.01

0.02

0.02

0.01

0.04 0.02

-0.04

0.00

0.02 0.02

0.07

0.04

Page 14 of 24

-0.07 -0.21

-0.06

0.02

0.04

0.02

0.00

0.03

0.00

0.01

0.01

-0.21

0.01

0.02

-0.21

-0.21

0.01

0.02

0.02

0.02

0.04

0.01

0.00

0.01

0.03

0.00

0.01

Figure 4: An OH molecule adsorbed on the fcc site of Rh(111) (left) and Cu(111) right. The number on each surface atom indicates the Bader charge assigned to the respective atom in number of electrons. The Bader analysis confirms the intuition developed through fig. 3. The charge reorganisation almost cancels out in the case of Rh while on Cu a net charge transfer occurs from those atoms. transition metals as given by density difference plot as exemplified in fig. 3. The panels on the left show the density difference upon OH adsorbed in the fcc site on Rh(111) while the two right panels show the density difference of OH on Cu(111). Both density difference plots shows significant charge reorganization at the surface atoms directly coordinated by the OH adsorbate. By visual inspection however the excess and depletion volumes seems rather similar in the case of Rh atoms while the depletion volumes dominate at the Cu atoms. By using a Bader charge analysis24–26 for each surface atom as displayed in fig. 4 we can also quantify trends in charge depletion for different surfaces. While the charge attributed to each surface atom remains very close to neutral for Rh, there is a significant charge depletion in the case of Cu. The more general trend is depicted in fig. 5. Each point represents one adsorbate and site on Cu(111). The abscissa represents the amount of charge depleted from surface atoms according to the Bader analysis. The ordinate represents the difference in adsorption energy between the (4 × 4) (1/16 ML coverage) surface and the (1 × 1) (1 ML coverage) surface. The graph indicates a positive linear correlation between these two properties or in other words: the more charge that is transferred from the surface atom, the larger the change in adsorption energy. While for Cu the d-band center is not a reliable descriptor for adsorbate-adsorbate interactions 14

ACS Paragon Plus Environment

Page 15 of 24

we do observe a strong correlation between the charge transfer from surface atoms and adsorbate self-interaction. Accordingly we here identify V˜ with the amount of charge removed from surface atoms in the (4 × 4) surface adsorption. Furthermore we identify V max with the total charge removed from surface atoms. Note that unlike in the non-coinage case we only consider the removed charge (i.e. negative numbers). Cu(111)

2.0

O@fcc O@hcp

S@fcc S@hcp 1.5 ∆Eads [eV]

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

The Journal of Physical Chemistry

N@fcc N@hcp

HO@hcp HO@fcc C@hcp C@fcc

CO@hcp CO@fcc

1.0

CH@fcc CH@hcp 0.5

0.0 0.0

H@hcp H@fcc 0.1 0.2

NO@hcp NO@fcc 0.3

0.4 0.5 −∆q [e]

0.6

0.7

0.8

Figure 5: Relation between adsorbate self-interaction and charge transferred from the surface on Cu(111). The abscissa represents that the amount of charge depleted from surface atoms in the low coverage (1/16 ML) adsorption. The ordinate represents the change in adsorption energy between 1/16 ML coverage and 1 ML coverage.

Linearized Interaction Model

Figure 6: Two examples for the geometries used to validate the interaction model. Both adsorbates are on either the fcc or hcp site in a (2 × 1) (left) and a (2 × 2) (right) unit cell. The slab consists of four layers with the lowest two layers held fixed at the bulk positions. The adsorbate atoms are held fixed in the cartesian coordinates parallel to the surface (xand y-axis). In order to explicitly test the predicted interaction energies we compared the adsorbateadsorbate interaction energies for several geometries as obtained from the presented model 15

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 16 of 24

Page 17 of 24

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

The Journal of Physical Chemistry

to the corresponding value of a corresponding DFT calculation. The test set consists of all combinations of pair interaction of the adsorbates C, N, H, O, S, HO, C, CH, and CO. For each combination both adsorbates are either placed on the fcc or the hcp site of the fcc(111) surface and either at a nearest neighbor distance in a (1 × 2) or at a next nearest neighbor distance in a (2 × 2) surface unit cell (see fig. 6). In order to exclude the influence of directional bonding each adsorbate atom was only allowed to relax in the z-direction orthogonal to the surface plane. For the benchmark case of two adsorbates A and B on a surface slab S, the adsorbateadsorbate interaction energy can be explicitly written as E int (A + B) = E(A+B)@S − (EA@S + EB@S ) + ES where ES is the DFT total energy of the clean slab, EA@S the DFT total energy of A adsorbated on the slab, EB@S the DFT total energy of B adsorbated on the slab, and E(A+B)@S the DFT total energy of both A and B adsorbed on the surface. This property E int is compared in fig. S1 in the Supplementary Information. On the other hand the quantity that directly enters into a kinetic model and accounts for adsorbate-adsorbate interaction is the differential adsorption energy or the adsorption energy of a given adsorbate and a given coverage. We write the differential adsorption energy of adsorbate A after adsorbate B was already adsorbed as diff EA|B@S = E(A+B)@S − (EAgas + EB@S ).

By inserting E(A+B)@S from the first equation into the second we obtain diff ads EA|B@S = E int (A + B) + EA@S

17

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 18 of 24

or by replacing A and B we obtain diff ads EB|A@S = E int (A + B) + EB@S

Fig. 7 shows the comparison of all binding differential adsorption energies E diff between the surface mediated model and the corresponding DFT calculation. Geometries leading to positive differential adsorption energies (i.e. non-binding for at least one adsorbate) where discarded altogether. It should be noted however that on Ag NO and CO have negligible or even positive differential adsorption energies. In this regime a surface mediated model cannot be expected to be quantitatively correct and we therefore exclude these adsorbates from the comparison. A detailed comparison of all differential adsorption energies can by found in the Supporting Information. Overall the deviation of the model from the explicit DFT calculation is 0.18 eV (root mean squared error) or 0.12 eV (mean absolute error). From the group of transition metals there are noticeably more outliers for Pd and Pt than for Rh and Ir. This may in part be explained by the stronger interaction between adsorbates and Ir and Rh compared to Pd and Pt. The accuracy of prediction for coinage metals is similar as for transition metals; however Cu still has noticeably more outliers than Ag. In conclusion the overall agreement shows that a surface mediated interaction model reproduces the largest contribution to adsorbate-adsorbate interaction qualitatively and quantitatively. While there are infinitely many other configurations that could be tested this test aims to strike a balance between a non-biased evaluation for each parametrized adsorbate. Furthermore, this particular test focuses on the regime of a total coverage between 0.5 ML and 1 ML in which the changes in adsorption energy due to adsorbate-adsorbate interaction can be expected to be the largest. The presented test can therefore be expected to provide an upper bar for the absolute prediction errors. Furthermore, having extracted intermediate electronic structure descriptors which allow one 18

ACS Paragon Plus Environment

Page 19 of 24

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

The Journal of Physical Chemistry

to reliably and quickly estimate the magnitude of adsorbate-adsorbate interactions may lead to an improved rationalization of the influence of co-adsorbates on a given reaction network. This may play an important role for tailoring optimal reaction conditions (and thereby coverages) to achieve higher reactivities or selectivities.

Summary We have presented a framework which allows the prediction of adsorbate-adsorbate interactions on transition and coinage metal surfaces from only two DFT calculations per adsorbate. The model predicts the magnitude of interaction between adsorbates in an explicit lattice picture via an auxiliary electronic structure property. Decomposing the interaction between adsorbates in such a way allows the avoidance of explicit calculations of many co-adsorption geometries from expensive electronic structure calculations. For non-coinage transition metals we propose a model where the auxiliary property is the relative position of the d-band center. For coinage metals we have presented a model where the auxiliary property is the charge depletion as predicted by a Bader analysis. By construction these models recover the exact DFT result at 0 ML coverage and 1 ML coverage for any particular adsorbate interacting with itself and they provide an accurate picture of cross-interactions between different adsorbates. Explicit comparison with first-principles calculations of adsorbate-adsorbate interaction between different species at nearest neighbor and next nearest neighbor distances shows that the model yields accurate results with a RMSE error across the tested adsorbates of less than 0.2 eV. Generally speaking, adsorbate-adsorbate interaction between nearest neighbor adsorbates can be on the order of 1 eV (cf. Fig. S1 in the Supplementary Material). By using the presented framework with its very moderate computational effort one can therefore mitigate the large energetic deviation from the physical situation if not adsorbate-adsorbate interaction is considered. 19

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 20 of 24

Potential next steps are an adaptation of the model to higher index surfaces including the important (211) facet, application to alloys which would require a conceptual consolidation of the two physical properties, or even to transition-metal compounds such as oxides and carbides. However, the latter will require a different set of electronic structure descriptors to achieve predictions as concise and reliable as now possible for pure metals. The key advantage of the presented framework is that it allows the prediction of adsorbateadsorbate interactions for explicit adsorption arrangements while the number of required first-principles calculation only grows linearly with the number of adsorbate species. Firstprinciples kinetic models which faithfully account for adsorbate-adsorbate interaction involving 3 or significantly more adsorbate species therefore become feasible. Such models are necessary to facilitate an understanding of complex surface chemical reactions under conditions where the coverage is high.

Acknowledgement We thank the US DOE-BES for support through the Early Career program to SUNCAT Center for Interface Science and Catalysis.

Supporting Information Lattice constants, clean surface d-band centers, adsorption energies of individiual adsorbates, tabulated adsorbate-adsorbate interactions, detailed interaction elements for adsorbateadsorbate interactions.

20

ACS Paragon Plus Environment

Page 21 of 24

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

The Journal of Physical Chemistry

References (1) Feibelman, P. J. Theory of Adsorbate Interactions. Annual Review of Physical Chemistry 1989, 40 (1), 261–290 DOI: 10.1146/annurev.pc.40.100189.001401. (2) Nørskov, J. K. Chapter 1 - Adsorbate-adsorbate interactions on metal surfaces. In The Chemical Physics of Solid Surfaces; Woodruff, D. A. K. and D. P., Ed.; Coadsorption, Promoters and Poisons; Elsevier, 1993; Vol. 6, pp 1–27. (3) Einstein, T. L. Interactions between adsorbate particles. Handbook of Surface Science 1996, 1, 577–650. (4) Hellman, A.; Honkala, K. Including lateral interactions into microkinetic models of catalytic reactions. The Journal of Chemical Physics 2007, 127 (19), 194704 DOI: 10.1063/1.2790885. (5) Hansen, E.; Neurock, M. Predicting lateral surface interactions through density functional theory: Application to oxygen on Rh(100). Surface Science 1999, 441 (2–3), 410–424 DOI: 10.1016/S0039-6028(99)00873-0. (6) İnoğlu, N.; Kitchin, J. R. Simple model explaining and predicting coverage-dependent atomic adsorption energies on transition metal surfaces. Physical Review B 2010, 82 (4), 045414 DOI: 10.1103/PhysRevB.82.045414. (7) Grabow, L. C.; Hvolbæk, B.; Nørskov, J. K. Understanding Trends in Catalytic Activity: The Effect of Adsorbate–Adsorbate Interactions for CO Oxidation Over Transition Metals. Topics in Catalysis 2010, 53 (5-6), 298–310 DOI: 10.1007/s11244-010-9455-2. (8) Stampfl, C.; Kreuzer, H. J.; Payne, S. H.; Pfnür, H.; Scheffler, M. First-Principles Theory of Surface Thermodynamics and Kinetics. Physical Review Letters 1999, 83 (15), 2993–2996 DOI: 10.1103/PhysRevLett.83.2993. (9) Drautz, R.; Singer, R.; Fähnle, M. Cluster expansion technique: An efficient tool to search for ground-state configurations of adatoms on plane surfaces. Physical Review B 2003, 67 (3), 035418 DOI: 10.1103/PhysRevB.67.035418. 21

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 22 of 24

(10) Wu, C.; Schmidt, D. J.; Wolverton, C.; Schneider, W. F. Accurate coverage-dependence incorporated into first-principles kinetic models: Catalytic NO oxidation on Pt (1 1 1). Journal of Catalysis 2012, 286, 88–94 DOI: 10.1016/j.jcat.2011.10.020. (11) Frey, K.; Schmidt, D. J.; Wolverton, C.; Schneider, W. F. Implications of coveragedependent O adsorption for catalytic NO oxidation on the late transition metals. Catalysis Science & Technology 2014, 4 (12), 4356–4365 DOI: 10.1039/C4CY00763H. (12) Herder, L. M.; Bray, J. M.; Schneider, W. F. Comparison of cluster expansion fitting algorithms for interactions at surfaces. Surface Science 2015 DOI: 10.1016/j.susc.2015.02.017. (13) Stamatakis, M.; Vlachos, D. G. Unraveling the Complexity of Catalytic Reactions via Kinetic Monte Carlo Simulation: Current Status and Frontiers. ACS Catalysis 2012, 2 (12), 2648–2663 DOI: 10.1021/cs3005709. (14) Vlachos, D. G. Multiscale modeling for emergent behavior, complexity, and combinatorial explosion. AIChE Journal 2012, 58 (5), 1314–1325 DOI: 10.1002/aic.13803. (15) Shustorovich, E. Chemisorption phenomena: Analytic modeling based on perturbation theory and bond-order conservation. Surface Science Reports 1986, 6 (1), 1–63 DOI: 10.1016/0167-5729(86)90003-8. (16) McNaught, A. D.; Wilkinson, A. Compendium of Chemical Terminology: Iupac Recommendations: Gold Book, 2 Sub.; IUPAC International Union of Pure; Applied Chem: Oxford England ; Malden, MA, USA, 1997. (17) Hammer, B.; Nørskov, J. K. Theoretical surface science and catalysis—calculations and concepts. In Advances in Catalysis; Bruce C. Gates, H. K., Ed.; Impact of Surface Science on Catalysis; Academic Press, 2000; Vol. 45, pp 71–129. (18) Yxklinten, U.; Hartford, J.; Holmquist, T. Atoms embedded in an electron gas: The generalized gradient approximation. Physica Scripta 1997, 55 (4), 499 DOI: 10.1088/00318949/55/4/022.

22

ACS Paragon Plus Environment

Page 23 of 24

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

The Journal of Physical Chemistry

(19) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter 2009, 21 (39), 395502 DOI: 10.1088/0953-8984/21/39/395502. (20) Wellendorff, J.; Lundgaard, K.; Møgelhøj, A.; Petzold, V.; Landis, D.; Nørskov, J.; Bligaard, T.; Jacobsen, K. Density functionals for surface science: Exchange-correlation model development with Bayesian error estimation. Physical Review B 2012, 85 (23), 235149 DOI: 10.1103/PhysRevB.85.235149. (21) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D.; Lundqvist, B. Van der Waals Density Functional for General Geometries. Physical Review Letters 2004, 92 (24), 246401 DOI: 10.1103/PhysRevLett.92.246401. (22) Román-Pérez, G.; Soler, J. M. Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes. Physical Review Letters 2009, 103 (9), 096102 DOI: 10.1103/PhysRevLett.103.096102. (23) Sabatini, R.; Küçükbenli, E.; Kolb, B.; Thonhauser, T.; Gironcoli, S. de. Structural evolution of amino acid crystals under stress from a non-empirical density functional. Journal of Physics: Condensed Matter 2012, 24 (42), 424209 DOI: 10.1088/0953-8984/24/42/424209. (24) Tang, W.; Sanville, E.; Henkelman, G. A grid-based Bader analysis algorithm without lattice bias. Journal of Physics: Condensed Matter 2009, 21 (8), 084204 DOI: 10.1088/09538984/21/8/084204. (25) Sanville, E.; Kenny, S. D.; Smith, R.; Henkelman, G. Improved grid-based algorithm for Bader charge allocation. Journal of Computational Chemistry 2007, 28 (5), 899–908 DOI: 10.1002/jcc.20575. (26) Henkelman, G.; Arnaldsson, A.; Jónsson, H. A fast and robust algorithm for Bader decomposition of charge density. Computational Materials Science 2006, 36 (3), 354–360 DOI: 10.1016/j.commatsci.2005.04.010. 23

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 24 of 24

0

1 2 3 4 5 6 7 8

-2 -4 -6 ACS Paragon Plus Environment-8

-8

-6

-4

-2

0