A Lattice Kinetic Monte Carlo Solver for First-principles Microkinetic

Subscriber access provided by READING UNIV

Article

A Lattice Kinetic Monte Carlo Solver for First-principles Microkinetic Trend Studies Max J Hoffmann, and Thomas Bligaard J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00683 • Publication Date (Web): 22 Jan 2018 Downloaded from http://pubs.acs.org on January 23, 2018

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.

Journal of Chemical Theory and Computation 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 34 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

Journal of Chemical Theory and Computation

A Lattice Kinetic Monte Carlo Solver for First-principles Microkinetic Trend Studies Max J. Hoffmann∗,†,‡ and Thomas Bligaard†,‡ †Department of Chemical Engineering, Stanford University, Stanford, California 94305, United States ‡SUNCAT Center for Interface Science and Catalysis, SLAC, National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, California 94025, United States E-mail: [email protected]

Abstract Mean-field microkinetic models in combination with Brønsted-Evans-Polanyi like scaling relations have proven highly successful in identifying catalyst materials with good or promising reactivity and selectivity. Analysis of the microkinetic model by means of lattice kinetic Monte Carlo promises a faithful description of a range of atomistic features involving short range ordering of species in the vicinity of an active site. In this paper, we use the "fruit fly" example reaction of CO oxidation on fcc(111) transition and coinage metals to motivate and develop a lattice kinetic Monte Carlo solver suitable for the numerically challenging case of vastly disparate rate constants. As a result, we show that for the case of infinitely fast diffusion and absence of adsorbateadsorbate interaction it is, in fact, possible to match the prediction of the mean field theory method and the lattice kinetic Monte Carlo method. As a corollary, we conclude that lattice kinetic Monte Carlo simulations of surface chemical reactions are most likely to provide additional insight over mean-field simulations if diffusion limitations

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

or adsorbate-adsorbate interactions have a significant influence on the mixing of the adsorbates.

1

Introduction

Energy descriptor and energy scaling relation based microkinetic models differ fundamentally from models based on thermodynamic descriptors because the aim is to establish trends across materials or classes of materials rather than rationalizing one particular experiment on one catalyst surface. 1,2 Herein lies the unique ability to be able to recognize and predict entirely new catalytic materials. 3 To this date mean-field analysis (MFA) is the state-of-the art technique. Lattice kinetic Monte Carlo (kMC) has been successfully used to study the microkinetics of systems where intricate features such as adsorbate-adsorbate interactions, 4–7 slow diffusion 8–10 of surface species or complex multi-site catalyst surfaces, 11,12 require a description beyond mean-field microkinetics. Micro-kinetic models based on adsorption energies as descriptors invariably imply vastly disparate timescales. This poses a significant challenge for applying kMC throughout a typical phase space. In this paper, this challenge is directly tackled by the presented solver and demonstrated using the prototypical reaction of CO oxidation on an fcc(111) surface. Important contributions that aim at the acceleration of sampling lattice kinetic Monte Carlo models exist in the literature. Most approaches suggest a systematic bottom-up approach focusing on identifying connected low-energy states known as superbasins 13 or on judiciously increasing energy barriers of frequently visited transitions that prevent sufficient sampling of other critical reaction steps. 14,15 Towards the more algorithmic side are contributions focusing on optimal data-structures 16–18 and parallelization of parts of the execution model. 19 A comprehensive overview of different approaches is given in Ref. 20. The solver presented here was thus developed from a particular viewpoint: can we establish direct correspondence between MFA evaluation and kMC evaluation for a subset of

2


Page 2 of 34

Page 3 of 34 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


simplified energy descriptor based microkinetic models? Which acceleration techniques are needed in the context of energy descriptor based microkinetic models? A direct agreement between the two pictures can serve as a basis for further method development. Specifically, it allows to directly quantify the effect of the mean-field approximation. This may assist to systematically develop more detailed protocols for the purpose of screening heterogeneous catalyst materials. The remainder of this paper has the following structure: first we introduce a protocol that translates a given descriptor based mean-field microkinetic model into a lattice kMC model (Sec. 2.3) and introduce steps necessary to adequately sample the lattice kMC model for different regimes in descriptor space. More specifically we address the problem to automatically determine a sufficient number of kMC sampling steps for representing steady-state conditions (Sec. 2.4) and to decelerate the rate constants of fast processes (Sec. 2.5). Moreover, we present new techniques sampling regimes with total coverages either very close to one (Sec. 2.7, 2.8) or to zero (Sec. 2.9). We present the effect of each of those steps on the predicted turn-over frequencies (Sec. 3) resulting in quantitative agreement between the two approaches.

2 2.1

Theoretical Background Description of surface chemical reactions

We start by briefly summarizing a few of the key notions of heterogeneous catalysis. Detailed mechanistic understanding of the systems and their dynamics serves as a basis for the development of most microkinetic models. A common assumption is that the relevant catalyst surface can be represented on the atomistic level by a single-crystal surface facet. Moreover, the chemical reaction proceeds via a set of simple elementary reactions. The elementary reactions transform adsorbates that are otherwise bound to periodic sites defined by the potential energy surface created by the electronic structure of the underlying catalyst 3



Page 4 of 34

surface. The time spent during the transition is generally orders of magnitude shorter than in the intermediate state which is also referred to as rare-event dynamics. Furthermore, all changes due to these elementary reactions are in themselves local, that is their effect is contained to within a few unit cells. Recent studies 21,22 address the influence of finite size effects on nano-particles. The underlying interaction is not within the focus of this study. The notions introduced above motivate the concept of a list of elementary processes that describe the transition from one set of adsorbates to another or the transition from a detached molecule (in the gas phase or liquid phase) to an adsorbed particle(s) or vice versa. We can denote the most general elementary process using the following notation

X j

k+

i + Si,A Aj * ) j −

ki

X

− Si,A Ak k

i = 1, 2, 3 . . .

(1)

k

where Aj denotes an adsorbate or vacancy A and j denotes site (type) or possibly continuum ± phase, Si,A is the stoichiometry factor for the ith elementary process and adsorbate Aj , and j

ki+ /ki− are the forward/reverse rate constants of the ith elementary process. At this point in practice, the microkinetic modeler faces a choice about the level of description of the surface state: one can either maintain the explicit lattice representation or assume a homogeneous distribution of adsorbates across all identical sites which yields the mean-field approximation. Continuing from this choice onwards follow two different techniques. Treating the system as a set of coverages calls for the use of rate equations. Treating the system state as an explicit adsorbate configuration calls for the use of a Master equation and its solution via (lattice) kinetic Monte Carlo. The derivation of the corresponding equations are well-established within microkinetic modeling literature. 23 We will here briefly highlight the resulting equation for the purpose of establishing the correspondence between the two approaches. Here and in the following we assume a fixed list of elementary processes. Protocols that discover and generate such a list on-the-fly exist and are exemplified by adaptive kinetic Monte Carlo 24 , however, this is out of the scope of this study. When choosing to apply the MFA it is then expedient to use a coverage vector Θ for 4


Page 5 of 34 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


representing the surface state (e.g. Θ = (ΘCO , ΘO , Θ∗ ), which leads to a set of coupled differential equations

ri = ki+

Y

+ Si,A

ΘA j

j

− ki−

Y

j

˙A = Θ j

X

− Si,A

ΘAk

j

i = 1, 2, 3, . . .

(2)

k

− + Si,A − Si,A ri j j

i = 1, 2, 3, . . .

i

Note that the indices j, k in Eq. (2) run over those same terms that appear in the corresponding line of Eq. (1). These equations have to be solved for obtaining steady-state coverages and rates. When choosing the explicit lattice representation of the system it is expedient to denote the system state in terms of explicit lattice configurations u, v, . . . , where each configuration defines the occupation of each site of the surface section under consideration. To illustrate this for the example reaction under consideration one would first fix a certain lattice size (e.g. (100 × 100) unit cells) and one state could be the empty lattice, another an empty lattice with one CO adsorbate in the lower left corner, another an empty lattice with an CO adsorbate in the unit cell right of the lower left corner etc.. Since the transition between lattice configurations via elementary processes occur on timescales much shorter than the residence time in any particular lattice configuration, we assume that the next transition can be predicted solely based on the rate constants of all available elementary process of the current configuration and the past history of the system can be truncated. This assumption is known as the Markov approximation and allows us to describe the propagation of the probability density of the system using the Master equation

ρ˙ u (t) =

X

wuv ρv (t) − wvu ρu (t)

(3)

v

where ρu (t) is the probability to be in state u at time t, and wvu is sometimes referred to as transition matrix.

5



Page 6 of 34

The transition matrix wvu can be decomposed into a sum over elementary processes

wvu =

X

i wvu

(4)

i

i = where wvu

    χi ki+    

−

χi ki       0

(u → v) ∈ i (v → u) ∈ i else

where χi is a geometry factor compensating for the possibly multiple ways the same transition can be realized on a lattice and ki± are those same forward and reverse rate constants of the ith elementary process.

2.2

Lattice Kinetic Monte Carlo Background

Solving the Master equation for any but the smallest realizations of a surface lattice is not feasible since the transition matrix w grows quickly with the number of unit cells and species. To this end, lattice kinetic Monte Carlo provides a solution using stochastic sampling. There are various names for this (dynamic Monte Carlo, 25 n-fold way 26 ), and thorough introductions exist in literature. 27–29 Here we will only motivate the most high-level concepts. These however are important to illustrate the challenges a sampling protocol has to overcome for sampling the highly disparate rate constants and concomitant time scales posed by an activity trend study based on energy descriptors as is considered here. There are at least three different but equivalent formulations of lattice kinetic Monte Carlo. 30 We here focus on the so-called variable step size method (VSSM) that is a rejection free algorithm and is used throughout. 26 Starting from an initial state a list of all elementary processes is created. From this list, one process and site is randomly selected with a probability weighted by the rate constant of the corresponding elementary process. The process is executed by updating the lattice

6


Page 7 of 34 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


configuration accordingly and the simulated time is incremented by

∆tkMC =

− ln(r) . ktot

(5)

where r ∈]0, 1] is a uniform random number. Afterward, one kMC simulation cycle is complete and the algorithm proceeds by building or updating the list of available processes. After repeating this basic kMC loop possibly multiple millions of times, observables such as rates and coverages can be extracted from the resulting trajectory. The sampled averages are time-averages but they correspond to ensemble averages due to the ergodic theorem. 27 Despite the simplicity of the algorithm, in order to achieve high numerical efficiency for sampling any particular simulated system and list of elementary processes, customized code often has to be written in a fast low-level programming language. In this study, we leverage kmos, a general lattice kinetic Monte Carlo framework, that has been designed in order to perform the heavy-lifting in implementing a lattice kinetic Monte Carlo model. 31 At the end of this high-level introduction we will discuss what practical hurdles kMC sampling could face in practice, based on the fact that all input parameters such as rate constants and mechanisms are dictated by first-principles calculations and that the key property to sample in the context of heterogeneous catalysis is the turn-over frequency or production rate of a catalytic cycle. The main question here is which type of process has a high probability to be executed most frequently since the computational cost per kMC step varies only little across different elementary processes. As stated above in each kMC step one process is chosen among all available processes with a probability weighted with its respective rate constant. In regimes and corresponding configurations relevant to high catalytic activity, surface diffusion via hopping is one such possible process since reactive surface configurations often feature a surface covered with well mixed adsorbates and a significant number of vacant sites. Furthermore, the surface is not completely covered since an efficient catalyst allows for fast

7



removal of intermediates via a low barrier reaction step and empty surface sites are necessary to allow replenishment via new adsorbates from the gas or liquid phase. Thus, under all circumstances, one or several adsorbates will be allowed to diffuse to empty nearby sites. Since surface diffusion of adsorbates on transition metals usually has a diffusion barrier of 0.5 eV or less it is usually orders of magnitude faster than other crucial reaction steps. Thus, typically many diffusion steps will be executed before any crucial reaction step is executed and therewith sampled. Another bottleneck can be fast consecutive adsorption and desorption steps. Adsorption rate constants as predicted by kinetic gas theory are typically on the order of 106 −108 s−1 . In contrast, corresponding desorption rate constants can be similarly large. For weakly bound adsorbates it can take many thousands if not million of pairs of adsorption-desorption steps before neighboring adsorbates are selected to react. An important limit is the maximum number of kMC steps that can be sampled on a single core performance. Despite recent progress in parallelizing the evaluation of complex reaction pathways and possible co-adsorbate dependent expression for rate constants, 19 the core algorithm of kMC remains elusive to parallelization. Henceforth, the maximum number of kMC steps that can be executed for the most simple kMC models of around 104 -106 kMC steps per CPU and second remains a firm upper limit. As such, this limit is difficult to improve more then incrementally due to necessary frequent random memory accesses based on today’s hardware. A total number of approximately 1010 kMC steps per simulation (∼ 104 CPUsecond) therefore remains a firm practical constraint. Together this poses a serious challenge for mapping the reactivity via kMC sampling of relevant catalytic materials, which we will address using a combination of acceleration techniques, as explained in the following.

8


Page 8 of 34

Page 9 of 34 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


2.3

Mean-field to lattice model translation

We now turn to describe the protocol by which we translate mean-field oriented micro-kinetic model into a lattice oriented micro-kinetic model. The main input for the mean-field model is a formal representation of a list of reaction steps as established in eq. (1). For the purpose of evaluating the mean-field micro-kinetic model, we use the open-source Python package CatMAP that was previously developed by Medford et al. 32 . The purpose of CatMAP is to automate and simplify descriptor-based analysis across catalyst materials. CatMAP standardizes the setup of a micro-kinetic model including the reaction steps, adsorption energies, reaction energies, vibrational frequencies, thermodynamic corrections via an input file. Furthermore, from this input, CatMAP can construct scaling relations between adsorption energies and transition state energies and finally generate and solve the resulting kinetic equations. The translation protocol is straightforward. Starting from the list of elementary processes as defined in the CatMAP 32 input file, we iterate over each process and generate corresponding kMC elementary processes. The essential additional information not already contained in the CatMAP (*.mkm) input file is the shape of the unit cell and the positions of the adsorption sites. In the case of fcc(111) surfaces under study here, we use the (111) surface slab as the unit cell and use a single site (i.e. the fcc site) as the only active site. Having defined this additional input the translation routine inspects each elementary processes, reuses the same rate constants, and constructs corresponding kMC elementary processes multiplied by a possibly non-unity geometry factor. However, in the kMC picture sites refer to lattice sites instead of abstract active sites. In the case of two-site processes (i.e. dissociative oxygen adsorption), all pairs of nearest-neighbor sites are selected based on the Euclidean distance and corresponding kMC processes are defined. As a result of this translation procedure we obtain an input file suitable for evaluation and analysis using the kmos lattice kinetic Monte Carlo package.

9



2.4

Steady-state detection 1.2 Start steady-state sampling at step 54 or 5.4%. 1.0 sampled rate

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.8 0.6 signal EWMA controlled fraction final mean LCL/UCL

0.4 0.2 0.0

0

200

400 600 kmc steps

800

1000

Figure 1: Steady-state detection scheme or steady-state after Rossetti. 33 The original signal (blue) is shown together with its exponentially weighted moving average (EWMA, red). The controlled fraction (thick black line) indicates which part of the EWMA stays within the lower and upper control limit (black dashed line). Sampling of the kMC signal continues until less than a given fraction (e.g. 0.05) of the full trajectory is outside the control limit. A quantity of central importance to a microkinetic model is the production rate or turnover frequency at steady-state conditions. The idea is that the catalyst surface is in contact with reservoirs of both reactants and products and while never reaching equilibrium with these reservoirs, after some transient regime the catalyst surface will reach a state at which all time-averaged coverages and turn-over frequencies remain constant. Due to the stochastic nature of kMC, those quantities do not converge monotonously but instead fluctuate around a converging function which can render the automatic detection of said steady-state challenging. The task is to identify the part of the signal as transient, which has to be deleted from sampling, and identify the remainder of the signal as steady-state that also has to be comprehensive enough to provide for accurate averages. Here we adapt a protocol put forward by Rossetti et al., 33 which is based on exponentially weighted moving average (EWMA) charts. 34,35 The EWMA x¯(ti ) is defined for any function x(ti ) of discretized time ti as

10


Page 10 of 34

Page 11 of 34 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


x¯(ti ) = λx(ti ) + (1 − λ)¯ x(ti−1 ).

(6)

where λ is a decay parameter that we choose as λ = 0.05. The EWMA protocol further defines lower and upper control limits (LCL/UCL) c± as r ci± = x¯(t∞ ) ± Lσ

λ [1 − (1 − λ)2i ] 2−λ

(7)

here σ is the standard-deviation of the signal and L is a parameter representing a tolerance range around the mean value that we choose as L = 6. Given a series of previously sampled rates and coverages from ntot kMC steps the idea is then to identify the part of the EWMA x¯(ti ) that stays within the lower and upper control limit that if we delete the first i kMC steps, which we coin the controlled fraction pi . As long as L is chosen not too small, the function pi will always tend toward 1 so we can define a first kMC step n1st for which pi = 1. We can then define a convergence criterion pc demanding that n1st /ntot ≤ pc . We choose pc = 0.05, throughout. Once a chosen pc is obtained for all pairs elementary processes and coverages the initial fraction is discarded and the remaining fraction corresponding to 1−pc of the kMC trajectory is used to calculate averages. As illustrated in Fig. 1 the controlled fraction pi first assumes 1 at step 54 out of 1000 steps or a fraction of 0.054 that would be deleted. Since we chose pc = 0.05 the simulation would sample a few more kMC step until this threshold is reached.

2.5

Adaptation of Equilibrated Process Pairs

Since activity maps cover a wide range of input descriptor values (typically on the order of 2 or 3 eV) compared to a characteristic energy scale of kB T it is inevitable that rate constants between different elementary process are widely disparate and may differ by ten or more orders of magnitude. This is a severe challenge for sampling possibly limiting turn11



Page 12 of 34

over frequencies, which are thereby the key feature of an activity map, by means of lattice kinetic Monte Carlo. Dybeck et al. have recently recognized that an analysis in terms of pairs of forward and reverse elementary processes allow to efficiently divide the set of elementary processes pairs into quasi-equilibrated and non-equilibrated process pairs. 15 They then proceed by decreasing the rate constants of quasi-equilibrated process pairs in order to increase the likelihood of executing non-equilibrated processes. Once a non-equilibrated elementary process has been executed the rate constants of the quasi-equilibrated process are unscaled again to their original values since the system is now assumed to be in a different superbasin. The criterion used to determine whether a pair of elementary processes is quasiequilibrated is |ni,forward − ni,reverse | ≤δ ni,forward + ni,reverse

(8)

where ni,forward and ni,reverse are the number of times the forward and reverse elementary process of the ith elementary process pair have been executed and δ is a convergence parameter. They then proceed to repeatedly reduce the rate constants of quasi-equilibrated pairs of elementary processes while maintaining detailed-balance until a non-equilibrated elementary process is executed. Once a non-equilibrated process is executed all rate constants are unscaled to their original setting and the procedure is repeated. In this paper, we use the same idea of partitioning the reaction network into quasiequilibrated and non-equilibrated pairs of elementary processes. However, before reducing the frequency of quasi-equilibrated processes we establish that we have obtained steady-state averages by means of the EWMA technique introduced in Sec. 1. Once we obtain steady-state averages for this particular set of rate constants, all fast processes are decelerated and the procedure is repeated until every pair of elementary processes is sampled a minimum number of times since the last adaptation of rate constants. Note however that the rate constants of fast processes are not unscaled after each iteration. Thus, in the absence of adsorbateadsorbate interaction and with infinitely fast diffusion we assume that the following holds:

12


Page 13 of 34 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


• All configurations connected by diffusion steps (identical coverage) belong to the same superbasin. • Under steady-state conditions, non-equilibrated reaction channels only create small fluctuations in coverage around steady-state coverages and the escape times of all visited superbasins are similar. The threshold parameter to establish quasi-equilibrium δ = 0.05, quasi-equilibrated reaction pairs are decelerated by dividing the corresponding rate constants by a factor 5, the minimum number of events per pair elementary reaction pair is 50. The algorithm is graphically summarized in Fig.2. Begin Simulation

Populate Event Table

Sample Steady-State Averages using EWMA Approach

Have we obtained suﬃcient statistics from each pair of elementary processes Yes

No Rescale rate-constants of equilibrated frequent processes

Report Rates & Coverages

Figure 2: The basic flow chart of accelerating rate constants. After initializing the simulation, steady-state averages according to EWMA technique are sampled. If a minimum number of events are sampled for every pair of forward-reverse elementary reactions is sampled the simulation exits. Otherwise, all elementary reaction pairs that are quasi-equilibrated are decelerated and the next iteration cycle begins.

2.6

The Microkinetic Model

The microkinetic model under consideration here is CO oxidation on fcc(111) transition and coinage metal surfaces using Langmuir-Hinshelwood type of kinetics for the reaction step. The reaction steps are

13



Page 14 of 34

CO(g) +∗ * ) CO∗

(9)

O2(g) + 2∗s → 2Os

(10)

COs + Os → CO2(g)

(11)

The transition state energies are obtained using scaling relations as implemented in the CatMAP package 32 . The original formation energies have been published elsewhere. 36 Using the machinery implemented in CatMAP the transition energies that are used in the model are

EO−CO = 0.674 EO + 0.674 ECO + 1.833 eV

(12)

EO−O = 1.432 EO + 3.19 eV.

(13)

Rate constants are obtained using harmonic transition state theory

k=

kB T exp (−∆G/(kB T )) . h

(14)

The sticking coefficient is assumed as unity for CO and O2 throughout the simulated phase space. The difference in entropy for adsorption and desorption processes is calculated from tabulated parameters using the Shomate equation and using the frozen adsorbate approximation for the adsorbed state. The model is evaluated at a temperature of 600 K and partial pressures of pO2 = 0.33 bar and pCO = 1.0 bar. We note that these conditions may differ from typical experimental conditions however for the purpose of comparing two simulation techniques within the scope

14


Page 15 of 34 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


of micro-kinetic trend studies they are representative.

2.7

Diffusion Processes

In the case of mean-field evaluation of a microkinetic model, there is no need to explicitly account for diffusion between identical sites since the assumption of infinitely fast diffusion implies immediate equilibration of surface intermediates between identical sites. However, in the case of kMC diffusion needs to be explicitly accounted for by adding elementary processes since the absence would be equivalent to having a diffusion rate constant of k = 0. Thus, the microkinetic model is amended by adding the following elementary reactions

COs + ∗s → ∗s + COs

(15)

Os + ∗s → ∗s + Os .

(16)

Since this study focus on the limit of infinitely fast diffusion, all diffusion rate constants are set initially three orders of magnitudes faster than the fastest non-diffusion rate constants and afterward decelerated in the same way as outlined in Sec. 2.5. Furthermore, an additional exchange diffusion step is added

COs + Os → Os + COs

(17)

to better sample some high-coverage regimes as shown in Sec. 3.5.

2.8

Embedding Mean-field Boundary Conditions

For those regimes where the surface is covered with strongly binding species (poisoned regime), the kMC algorithm generates large time-steps if the desorption rate constant of the majority species is slow. This can present an additional challenge for sufficient sampling

15



Page 16 of 34

Figure 3: The active lattice surrounded by one row of pseudo MFT sites. These sites mimic a reservoir of species available for diffusion processes if insertion of species via regular processes would inhibit sufficient sampling of all elementary processes. non-equilibrium elementary processes even after decelerating all fast processes. To illustrate this using the CO oxidation kinetics, consider a situation where the desorption is slow (e.g. 10−26 s−1 ). Consider furthermore that the turn-over frequency of the rate limiting process was 10−16 s−1 . Thus for every desorption event on the order of 1010 CO oxidation events would have to be executed before assuming a quantitatively correct result. Besides being on the edge of computational feasibility, it is difficult to a priori impose the correct requirements without referring to the mean-field result. Thus in a practical scenario, where one assumes some fixed number of events per elementary-process pair as sufficient (100-1000), the overall rate would be grossly underestimated. To this end, we modify the boundary conditions of the active kMC lattice. Instead of the usual periodic boundary conditions the edge of the lattice is replaced by one row of reservoir sites. Those reservoir sites are able to execute all diffusion elementary processes modified by the mean-field coverages. Thus, if the regular diffusion process had the following update rules and rate constant kidiff were

COs + ∗s → ∗s + COs ;

k = kidiff

(18)

then the amended elementary process would be

MFTs + ∗s → MFTs + COs ;

k˜ = ΘCO kidiff

(19)

Corresponding processes are added for the reverse process and for every active surface 16


Page 17 of 34 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


species and site. The rate constants are initialized just like the corresponding regular diffusion rate constants and afterward iteratively adapted via the adaptation algorithm outlined in Sec 2.5. Careful testing using a kMC lattice of (20 × 20) sites

2.9

Sampling ultra-low coverage regimes (a) 1 adsorbate desorption diﬀusion adsorption

(b) 1 adsorbate

...

...

kMC steps, simulated time

> 1 adsorbates

> 1 adsorbates

Figure 4: One particle acceleration technique: in regimes with low coverages it is critical to sample configuration with more than one particle on the surface for sampling the complete catalytic cycle. Subfigure (a) indicates a typical trajectory of kMC steps involving many futile steps, whereas (b) shows a modified trajectory that skips those same futile steps for configurations where only one adsorbate is present on the surface and jumps directly to the next adsorption event. Another regime that is particularly difficult to sample using standard lattice kinetic Monte Carlo is that with dilute coverages. The challenge here is that even though one particle will be on the surface on average in every other kMC sampling step it is unlikely to have more than one particle on the surface. Consider the case where the average coverage of oxygen and CO is 10−9 ML, respectively. In this case, the probability of sampling both particles next to each other is therefore approximately 10−18 and therefore prohibitively low for sampling in 1010 kMC steps. However, having at least two particles on the surface is typically a necessary condition for sampling reaction typical between two reactants in a heterogeneous catalysis reaction. To this end, we avoid the sampling issue by blocking futile intermediate events for those configurations where there is only one particle on the surface. As such diffusion of this single particle among equivalent sites is disabled as it only generates more identical states due to the lattice symmetry of the system. Furthermore, desorption is disabled in order to 17



avoid exuberant desorption-adsorption sequence of one particle, which in itself cannot sample reactions. The only other steps possible and allowed in this situation is adsorption of another P particle. By virtue of the additivity of rates (ktot = i ki ), the waiting time for this second adsorption event is unchanged by the removal futile intermediate events. The basic idea is illustrated in Fig. 4. One key complication stems from the fact that simply freezing the single adsorbate on the surface artificially increases the probability that if a second particle adsorbs on the surface it will always find a first particle waiting. To avoid this error, we split each adsorption event in those one-particle configurations into two channels: one for the case that the first particle is still on the surface and one for the case that the next adsorbate impinges on an empty surface. The probability of executing each is weighted by the corresponding probability of finding each configuration based on mean-field coverages.

3

Results

Fig. 5 shows an overview of the resulting turn-over frequency and coverages from an evaluation of the mean-field calculation of the microkinetic model. This result will serve as the reference point for comparison of the kinetic solver after applying a number of different acceleration techniques. For high oxygen adsorption energies the surface is predominantly covered by oxygen, for CO adsorption energies the surface is predominantly covered by CO. Noticeably, in Fig. 5(b/c) there is a significant part of the activity map characterized by lower overall coverage. As exemplified by the regime at and around Ag there is still appreciable catalytic activity. While overall, most of the activity is seen near the transition regime between one majority species to another, we manifestly have to plot out a descriptor space of 2-4 eV for each descriptor to get a comprehensive picture of the activity space.

18


Page 18 of 34

Page 19 of 34 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


(a)

(b)

(c)

Figure 5: The coverage and turn-over-frequency resulting from evaluation of the mean-field version of the microkinetic model: the subfigures show (a) the turn-over frequency, (b) the occupation with CO species, and (c) the occupation with oxygen species.

19



Table 1: Geometry factors χi for four different lattice geometries fcc(111) with one site, fcc(111) with two sites (e.g. fcc and hcp), fcc(100), and fcc(110). All kMC simulations are run with all rate constants set to k = 1. Quantity (111/1) (111/2) (100) (110) ΘCO 0.33 0.33 0.33 0.33 0.33 0.33 0.33 0.33 ΘO CO Rads/des 0.33 0.66 0.33 0.33 −1 1 2 1 1 χ O2 Rads/des 0.66 0.66 0.44 0.22 −1 6 6 4 2 χ CO ox. 0.66 0.66 0.44 0.22 6 6 4 2 χ−1

3.1

Geometry factors

During the translation step from MFT to kMC, one elementary process in the MFT picture may translate into more than one process in the kMC picture due to the symmetry of the chosen surface facet. To illustrate this we consider the CO oxidation step on an fcc(111) facet assuming one active site per unit cell. The mean-field description of this process requires a CO adsorbate present on one site and an oxygen adsorbate on a neighboring site, i.e. implicitly assuming that every site is connected with exactly one neighboring site. However, when picturing this description on an fcc(111) we note that there are now equivalent nearest-neighbor sites resulting in six equivalent elementary processes in the kMC evaluation corresponding to one process in the MFT evaluation. For the purpose of elaborating which techniques have to be applied to achieve parity between the MFT evaluation and kMC evaluation, we have introduced geometry factors χi in the translation formalism indicated in Eq. 4. To determine the geometry factors translated the CO oxidation model using the automated procedure described in Sec. 2.3 with 4 different unit cell geometries: fcc(111) with one site per unit cell, fcc(111) with two sites per unit cell (fcc and hcp site positions), fcc(100) with one site, and fcc(110 with one site. For each of the resulting 4 kMC models, we perform independent kMC simulations setting all rate constants ki = 1. This way we ensure that every elementary process can be adequately sampled and we can read off the

20


Page 20 of 34

Page 21 of 34 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


corresponding geometry factors. For each geometry the resulting coverages where ΘCO = ΘO = 1/3 identical and as expected from the mean-field evaluation. However, the forward and reverse rates are geometry dependent and proportional to the geometrical multiplicity of the corresponding elementary process. The sampled coverages, rates, and geometry factors are summarized in Tab. 1. Consequently, geometry factors of 1, 1/6, 1/6 are used for CO adsorption/desorption, O2 adsorption/desorption, CO oxidation/reduction for the kMC elementary process for the remainder of this paper as is appropriate for the fcc(111) geometry.

3.2

Adaptation of Equilibrated Process Pairs

Fig. 6(a) shows the obtained kMC result after applying the steady-state detection technique as well as the quasi-equilibrated process acceleration technique. Fig 6(b) shows the ratio between the MFT and the kMC result (RMFT /RkMC ) on a logarithmic scale. For further analysis of these differences we have divided the diagram into regime I, II, and III. Regime I is characterized by high oxygen coverage and generally large deviations between the MFT and kMC sampling of up to 1024 . In the case of the largest difference, the kMC result underestimates the MFT result. Regime II is characterized by generally low coverage. For most parts, the kMC result underestimates the MFT by up to 108 in the low coverage case. Regime III is characterized by not too strongly poisoned surface coverages. Here the agreement between the kMC result and the MFT result is best and is within the same order of magnitude or better. As a result, we note, that strong poisoning by one species poses a challenge for kMC sampling even after successively decelerating quasi-equilibrated process until all elementary reactions can be sampled. The reason is that the long time intervals created by the kMC algorithm in the fully oxygen covered state can not be compensated by a sufficient number of CO oxidation steps since the total number of kMC steps is limited to about 1010 steps per simulation as discussed in Sec. 2.2. This situation is unique to the coverage representation of 21



(a)

II

I

III (b)

II

I

III Figure 6: The turn-over frequency after applying the steady-state detection technique and the quasi-equilibrated process acceleration technique: (a) the raw result from kMC sampling, (b) the logarithm of the ratio between the MFT and raw kMC rate.

22


Page 22 of 34

Page 23 of 34 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


lattice kMC. Considering that the vacancy coverage can be as low as 10−15 ML in regime I, the fully covered state could theoretically be avoided by using an active lattice of 1015 sites. However, this would be beyond the available memory of available hardware. Instead, we turn to the Mean-Field Embedding technique.

3.3

Embedding Mean-Field Boundary Conditions

(a)

V IV (b)

V IV

Figure 7: The kMC turn-over frequency without one-particle technique (top) and the differences with respect to the mean-field result (bottom). Fig. 7 shows the result of kMC sampling after applying the mean-field embedding technique described in sec. 2.8(a) and the ratio between the MFT sampling result and kMC sampling result on a logarithmic scale. At this point, most data-points in the sampled space 23



from kMC are in good agreement with the MFT results. As indicated by the regime IV the averaged turn-over frequency is at least within the same order of magnitude. A large deviation exists in the regime in the upper right corner of fig. 7 labeled as V. This regime is generally characterized by low coverages of both CO and O as low as 10−9 ML as shown in Fig. 5. Having such low coverage causes a probability of sampling CO and on neighboring sites on the order of 10−18 , which cannot be reliably sampled with the given number of 1010 kMC steps. Consequently, we apply the next acceleration techniques which avoid sampling of futile transitions in the low coverage regime altogether as described in Sec. 2.9.

3.4

Sampling ultra-low coverage regimes

Fig. 8(a) shows the kMC sampling result after blocking futile elementary reaction in the ultra-low coverage regime. The lower Fig. 8(b) show the linear ratio the MFT result and the kMC sampling result. The quantity describing the difference shown here is

∆=

    RMFT − 1 RkMC

RMFT ≥ RkMC

  −( RkMC − 1) else RMFT to ensure meaningful color bars. Most parts of this activity map now show good agreement for most a indicated as label VI. The only regime that still shows significant aberration from the MFT result is indicated under label VII. This part of the diagram is characterized by high CO coverage.

3.5

Exchange diffusion

Careful inspection of configuration snapshots under these conditions shows that the difficulty in sampling this regime stems from a fundamental asymmetry between the O poisoned and the CO poisoned regime. In the oxygen poisoned regime vacancies can enter the kMC lattice via the boundary conditions and diffuse readily through the other oxygen covered surface.

24


Page 24 of 34

Page 25 of 34 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


(a)

VI

VII (b)

VI

VII

Figure 8: The kMC sampled result after applying the low coverage sampling technique on top of deceleration of quasi-equilibrated processes and the mean-field embedding technique (top) and the difference between the

25



CO can adsorb in this single vacancy and complete the catalytic cycle by reacting with neighboring oxygen atoms. On the other hand in the CO poisoned regime, a single vacancy is not sufficient to adsorb O2 since two neighboring vacancies (divacancy) are needed. Thus, the only channel for the minority species oxygen to enter the CO poisoned surface exists if a vacancy happens to occur right next to the lattice boundary. Once oxygen enters at this position oxygen can only leave immediately through the boundary or react off with a CO molecule creating CO2 in the gas phase. To this end, we introduce an additional diffusion step that allows exchange between CO and oxygen adsorbates on neighboring sites (cf. Sec. 2.7). Fig. 9(a) shows the final kMC sampling result after applying all acceleration techniques including the exchange diffusion between CO and oxygen adsorbates. As can be seen in the subfigure 9(b) the quantitative agreement between the two approaches of evaluating microkinetic models is very good throughout the phase space. With few exceptions the turn-over frequencies match to within a factor of 2. Few data-points result in an agreement of only within a factor 4. Deviations on this order are still difficult to avoid even with carefully chosen sampling criteria due to the stochastic nature of the kMC sampling algorithm. Both simulations were carried out on identical standard cluster hardware using Intel Xeon CPUs with 2.66 GHz. Sampling the MFT result took 24.7 seconds, whereas sampling the final kMC result took 27.8 hours.

4

Discussion

We have presented both a translation algorithm and a solver that allows evaluating a standard microkinetic model either using the mean-field method or the lattice kMC method. We establish that for the conceptually important limit of infinitely fast diffusion and absence of adsorbate-adsorbate interaction (where adsorbate mixing should be ensured) both evaluations result in identical turn-over frequencies.

26


Page 26 of 34

Page 27 of 34 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


(a)

(b)

Figure 9: Subfigure (a) shows the kMC sampling result after applying the complete suite of techniques: automatic steady-state detection, deceleration of quasi-equilibrated elementary reactions, normalized geometry factors, ultra low coverage sampling, and exchange diffusion. Subfigure (b) shows the linear ratio between the kMC sampling and the MFT result. The white parts of the diagram indicate that the kMC result agrees with the MFT counterparts to within a factor of 2. The remaining outliers indicated slightly larger fluctuations which are difficult to avoid even with carefully chosen numeric settings.

27



The lattice kMC solver of the microkinetic model was motivated in the context of energy descriptor phase space using the fruit-fly reaction of CO oxidation on fcc(111) metal surfaces. Using the combination of the presented acceleration techniques yields quantitatively identical results when using the lattice kMC method. Naturally, the question arises how much of this will work in a similar way when turning to a different and possibly more complex microkinetic model. That is: how much of the deployed techniques is transferable to other reaction systems? In this regard, it is clear that the automatic detection should always be useful since it frees the modeler from manually testing the number of kMC steps that is sufficient, which may also differ vastly for different parts in phase space. The automatic deceleration of quasi-equilibrated elementary processes is critical to bridge the disparate time scales of different elementary processes. This often occurs in microkinetic models, and especially seems inevitable for energy descriptor based microkinetic models. Thus, this acceleration technique should be useful for almost any model but may require from fine-tuning of the particular convergence parameters. Furthermore, the existence of strongly poisoned regimes, as well as those characterized by low coverages (

A Lattice Kinetic Monte Carlo Solver for First-principles Microkinetic

Recommend Documents