CatMAP: A Software Package for Descriptor-Based ... - Springer Link

Catal Lett (2015) 145:794–807 DOI 10.1007/s10562-015-1495-6

CatMAP: A Software Package for Descriptor-Based Microkinetic Mapping of Catalytic Trends Andrew J. Medford • Chuan Shi • Max J. Hoffmann Adam C. Lausche • Sean R. Fitzgibbon • Thomas Bligaard • Jens K. Nørskov

•

Received: 14 January 2015 / Accepted: 6 February 2015 / Published online: 17 February 2015 Springer Science+Business Media New York 2015

Abstract Descriptor-based analysis is a powerful tool for understanding the trends across various catalysts. In general, the rate of a reaction over a given catalyst is a function of many parameters—reaction energies, activation barriers, thermodynamic conditions, etc. The high dimensionality of this problem makes it very difficult and expensive to solve completely, and even a full solution would not give much insight into the rational design of new catalysts. The descriptor-based approach seeks to determine a few ‘‘descriptors’’ upon which the other parameters are dependent. By doing this it is possible to reduce the dimensionality of the problem—preferably to 1 or 2 descriptors—thus greatly reducing computational efforts and simultaneously increasing the understanding of trends in catalysis. The ‘‘CatMAP’’ Python module seeks to standardize

and automate many of the mathematical routines necessary to move from ‘‘descriptor space’’ to reaction rates for heterogeneous (electro) catalysts. The module is designed to be both flexible and powerful, and is available for free online. A ‘‘reaction model’’ can be fully defined by a configuration file, thus no new programming is necessary to change the complexity or assumptions of a model. Furthermore, various steps in the process of moving from descriptors to reaction rates have been abstracted into separate Python classes, making it easy to change the methods used or add new functionality. This work discusses the structure of the code and presents the underlying algorithms and mathematical expressions both generally and via an example for the CO oxidation reaction. Graphical Abstract

A. J. Medford (&) C. Shi M. J. Hoffmann A. C. Lausche S. R. Fitzgibbon T. Bligaard J. K. Nørskov Department of Chemical Engineering, Stanford University, Stanford, USA e-mail: [email protected]

Keywords Heterogeneous catalysis Surface reaction kinetics Kinetic modeling Ab initio calculations

A. J. Medford C. Shi M. J. Hoffmann A. C. Lausche T. Bligaard J. K. Nørskov SLAC National Accelerator Laboratory, SUNCAT Center for Interface Science and Catalysis, Menlo Park, USA

123

1 Introduction Microkinetic modeling (MKM) is a widely used technique in connecting an atomic or elementary-step level

CatMAP

description of a reaction pathway to the measureable kinetic behavior observed in experiments. There are different levels of kinetic modeling, varying in complexity from mean-field models with a single rate-determining step to stochastic simulations such as kinetic Monte Carlo [1–6]. At the high temperatures typical of most heterogeneously catalyzed reactions on metal surfaces, mean-field kinetic models are often found to describe the kinetics semiquantitatively, and there are a number of cases where it has been shown to describe trends in reactivity from one catalyst to the next very well [7–15]. Mean-field models are by far the simplest to use in the screening of catalytic processes, and will form the basis for the discussion in the present paper. In recent years the discovery of correlations between the energies of various intermediate states and transition states on the surfaces of heterogeneous catalysts [16–21], has led to the reduction in the dimensionality of the parameter space necessary to define the energetics of a reaction path or network [11, 22]. The combination of these techniques leads to the ability to predict trends in kinetic properties of heterogeneous catalysts across a range of surfaces using a ‘‘descriptor-based’’ analysis [10, 11]. This technique has been applied to a large number of chemical and electrochemical reactions, and often results in the prediction of novel catalysts with experimentally verifiable activity [8, 23–27]. Despite the success of this general approach, the implementation of such analyses is typically done on a case-bycase basis. While effective, the strategy of re-implementation is inefficient and in many cases prevents a proper analysis. In contrast, the prospect of automation suggests several major advantages: consolidation of knowledge, easy evaluation of sensitivity to assumptions, and a low barrier to test new ideas. Herein we present a Python module named the ‘‘Catalysis Microkinetic Analysis Package’’ (CatMAP) that provides a flexible framework for automatically constructing descriptor-based microkinetic analyses. CatMAP is designed to act as a platform for the development of new methods for mapping catalytic activity through descriptorbased analyses, and will eventually be integrated to the larger context of CatApp [28]. The paper discusses the interface and structure of CatMAP. In addition, each element of CatMAP is presented from a general perspective along with a simple CO oxidation example.

795

user can run the model and obtain a variety of outputs including the coverages of intermediates, reaction rates, turnover frequencies, selectivities, and ‘‘degree of rate control’’ [29, 30]. The structure of CatMAP, shown schematically in Fig. 1, utilizes the object orientation of Python to provide a modular and flexible interface, allowing advanced users to create their own classes that can extend the methods and assumptions outlined above. The class structure provided is: •

•

•

•

• •

ReactionModel: Messenger class that holds all attributes and common methods of the descriptor-based kinetic model. Parser: Converts the raw data into a structure compatible with the model. Could be extended to interface with databases, a file system, or other data structures. Scaler: Constructs the ‘‘scaling relations’’ which map descriptor space onto the parameter space necessary for solving the microkinetic model. Thermodynamics: Includes the thermal enthalpy, entropy, and any other free energy contributions that are not accounted for by the energies in the raw data, including corrections needed for electrochemical environments. Solver: Solves the microkinetic model, providing a map between the parameter space and a reaction rate. Mapper: Obtains initial guesses using Gibbs distributions and by using information from other points in descriptor space. Creates a ‘‘map’’ of the output as a function of descriptor space.

Implementation of new types of these classes requires minimal interaction with the other parts of CatMAP, allowing users to utilize some parts of the code without being tied to the assumptions in others. The code is opensource and freely available to facilitate communication and consolidation of knowledge between various researchers and groups. A link to the most recent version of the code can be found at http://suncat.slac.stanford.edu/#/ facility/software/catmap, and more thorough documentation, tutorials, and examples are provided there. In the next sections we discuss the methods used in each part of the code, along with examples for the CO oxidation reaction. The specific code to generate the CO oxidation example is provided along with the source code and several other examples online.

2 CatMAP Implementation 3 Generating Energy Inputs The CatMAP module relies on two input files from the user: a raw data table and a model configuration file. The module reads these files and provides an object-oriented API (application programming interface) from which the

The energetics of products, reactants, and intermediates are the central inputs to a MKM. In most instances these energetics will be used to construct rate and equilibrium

123

796

A. J. Medford et al.

Fig. 1 Schematic of the CatMAP code structure. Dark gray boxes correspond to classes in the CatMAP code, and white boxes represent different types of intermediate data along with some representative examples. Arrows indicate flow of data through the code from the raw input to final output

constants; however in the interest of preserving generality and enforcing thermodynamic consistency it is preferable to work with energies of individual states rather than the rate/equilibrium constants. In principle these energies may be obtained by experimental or theoretical means. In practice the energetics of adsorbed states are cumbersome to measure [31], and therefore CatMAP is geared towards the use of large data sets generated from electronic structure calculations such as density functional theory (DFT). This fact leads to the need to have an input format that is general enough to allow the user to input energies that are determined in various ways.

to accurately calculate (e.g. gas-phase O2 is very poorly described in DFT). In order to circumvent this issue it is useful to use a somewhat more general definition of formation energies. Let ~ e be a vector of ‘‘atomic reference energies’’, where there is one entry for each atomic species in the system of interest. In practice it is impossible to calculate absolute energies, so any measured or calculated energy is inherently relative to some reference energies ei. In this sense, computing formation energy is always equivalent to changing the atomic reference energies: X Ej ð~ eÞ ¼ Ej ð~ e Þ ni e i ð1Þ i

3.1 Formation Energy Approach The framework of thermodynamics requires that all energies be specified as relative quantities, and thus it is not possible to obtain energies of individual species that are completely independent of others. However, by taking the energies of all states relative to some common set of references it is possible to obtain ‘‘formation energies’’ for each species. Formation energies are, in general, computed by using a ‘‘reference state’’ for each element present in the state of interest. By IUPAC convention, the ‘‘standard state’’ for each element is the most stable pure allotrope of that element at standard temperature and pressure. In principle, the assignment of reference states is arbitrary, and it is only important to ensure that they are consistent within data sets. This becomes important when using theoretical methods to compute formation energies, since it is not uncommon for the ‘‘standard state’’ to be very difficult

123

where ~ e are defined as 0 relative to ~. e In order to use formation energies in a practical scheme it is required that all compatible energies use the same reference set so that the ~ e are consistent in computing relative energies between states with equivalent compositions. One approach for achieving this consistency is relating the atomic reference energies to the energies of specific molecular compounds. The IUPAC convention uses atomic reference energies equivalent to the most stable pure allotrope of a substance; however this is not required. Any set of formation energies will, in principle, be compatible provided their atomic references are related to the same molecular compounds. The example below demonstrates this. The use of formation energies as inputs to a model is also preferable to relative quantities such as reaction/activation energies because it ensures thermodynamic consistency. For example, if two elementary steps share an intermediate species in their final state, and both quantities

CatMAP

797

have a different error associated with them, then the formation energy of the shared species becomes ambiguous. Such ambiguities can often be resolved, but the formation energy approach avoids them and is hence preferred for the automated implementation.

eÞ ¼ Eslabþads ð~ e Þ Eslab ð~ e Þ Eads ð~

ð2Þ

~Þ is the ‘‘formation energy’’ of the adsorbed where Eads ðe ~ Þ is the energy of the adsorbate and the species, Eslabþads ðe ~ Þ is the energy of the bare slab, and ni surface slab, Eslab ðe is the number of molecular species i in the adsorbate. This is technically equivalent to Eq. (1), assuming that the reference energies, ~, e of the species in the surface are taken relative to the bare slab. It is also possible to compute ‘‘formation energies’’ using reaction energies of elementary steps. For this example, we will assume a very simple mechanism for CO oxidation over a catalyst surface:

Consider the CO oxidation reaction: 2COðgÞ þ O2ðgÞ ! 2CO2ðgÞ The energies of all relevant species can be calculated using electronic structure methods, or the experimental energies can be found in databases online [32]. Table 1 shows values calculated with the GPAW code [33] using the RPBE functional [34] (relative to the IUPAC standard states) and values from the CCCBDB [32] (corrected for zero-point energy and thermal enthalpy). It is clear that some species are in better agreement than others. In particular, the reference states agree perfectly (by definition), while there are relatively large discrepancies for the H2O(g) and CH4(g) molecules. It is well known that DFT performs poorly at calculating the triplet state of O2(g) [35, 36], and the RPBE functional does not include van der Waals forces necessary to describe graphite [34], so these are not ideal choices for reference states. Instead, we can choose to express the energies relative to carbon from CO(g), oxygen from H2O(g), and hydrogen from H2(g). Applying Eq. (1) gives the energies of all species in the new reference set as shown in Table 2. For this small set, the accuracy appears to have decreased; however, the discrepancies in the numbers in Table 2 are due to the shortcomings of DFT in calculating gas-phase O2 [35], CO2 [37], and graphite. Note that any reaction energies calculated from either table will remain constant, demonstrating that the choice of reference is arbitrary (assuming no corrections are applied). A similar approach can be obtained to calculate ‘‘formation energies’’ of adsorbed species:

Energies are given relative to graphite, O2(g), and H2(g)

ni e i

i

3.2 Example: Formation Energies for CO Oxidation

Table 1 Electronic potential energies of gas-phase species in CO oxidation reaction calculated using the RPBE functional and taken from experiment (from [32], corrected for zero-point and thermal enthalpy contributions)

X

COðgÞ þ ! CO O2ðgÞ þ 2 ! 2O CO þO ! CO2ðgÞ þ 2 This reaction has been studied previously with more complex models [38, 39], but we focus on this simple model here to illustrate the mechanics of CatMAP. In order to obtain the energetics for this reaction we will use the CatApp database [28]. The reaction energies for the elementary steps can be written in terms of the ‘‘formation energies’’ of each constituent species, leading to a relationship between the relative energies and the formation energies: DE1 ¼ ECO ECOðgÞ E ) ECO ¼ DE1

ð3Þ

DE2 ¼ 2EO EO2 ðgÞ 2E ) EO ¼ 0:5ðDE2 þ EO2 ðgÞ Þ ð4Þ where we have used the fact that E* is defined as 0 in the CatApp database, and ECO(g) is defined as zero when CO(g) is a reference state (Table 2). It is also necessary to determine the ‘‘formation energies’’ of transition-states in order to calculate the rate constants. This can be achieved in a similar way:

Species

Calculated energy (eV)

Experimental formation energy (eV)

CO(g)

-1.06

-1.10

CO2(g)

-4.14

-4.13

O2(g) H2O(g)

0 -2.37

0 -2.71

H2(g)

0

0

C(graphite)

0

0

CH4(g)

-1.45

-1.16

123

798 Table 2 Electronic potential energies of gas-phase species in CO oxidation reaction calculated using the RPBE functional and taken from experiment

Energies are given relative to CO(g), H2O(g), and H2(g)


Species

Calculated energy (eV)

CO(g)

0

0

CO2(g)

-0.71

-0.31

O2(g)

4.75

5.43

H2O(g)

0

0

H2(g)

0

0

C(graphite)

-1.31

-1.61

CH4(g)

-2.76

-2.77

DEa;2 ¼ EOO EO2 2E ) EOO ¼ DEa;2 þ EO2 ðgÞ ð5Þ DEa;3 ¼ EOCO ECO EO ) EOCO ¼ DEa;3 þ DE1 þ 0:5ðDE2 þ EO2 ðgÞ Þ

ð6Þ

where CO adsorption is assumed to have no barrier, and the relationships from Eqs. (3) and (4) have been used. For this example we will examine the reaction over FCC(111) transition-metal surfaces. The relevant reaction/activation barriers from CatApp and the corresponding computed ‘‘formation energies’’ are tabulated below. In order to properly describe the energetics of CO adsorption it is necessary to correct for the known over-binding [40, 41]. For simplicity we have taken the over-binding correction to be 0.25 eV for all metals, which is representative of the corrections for noble transition metals calculated with RPBE [41]. The ‘‘formation energies’’, Ei, in Table 2 (gasphase species) and Table 3 (adsorbed species) can be used as inputs to CatMAP, and will form the basis of subsequent analyses in this paper. The experimental gas-phase formation energies are used based on their improved accuracy and the fact that they present a convenient way to apply previously reported corrections for gas-phase CO2 [37] and O2 [35]. The specific format of the input to CatMAP will depend on the ‘‘Parser’’ class that is selected. The simplest possibility is a text table, and a parser for such tables is currently implemented. The required structure of the file is discussed in the CatMAP documentation. It is also possible to create custom ‘‘Parser’’ classes, allowing integration with databases or raw

Table 3 Reaction, activation, and formation energies of relevant adsorbed species and transition-states on FCC(111) metal surfaces

Formation energies are computed using Eqs. 3–6 and calculated gas formation energies from Table 2

123

Experimental formation energy (eV)

Surface

DE1

DE2

Pt(111) Pd(111)

-1.24 -1.39

-1.49 -1.63

data; more information is available in the online documentation. 4 Including Thermal Enthalpy, Entropy, and Electrochemical Corrections In most cases the energetics provided will be obtained using electronic structure theory, and hence correspond to 0 K. This leads to the need to include the thermal enthalpy ( r cp dT) and entropy associated with the finite temperatures at reaction conditions. This is a commonly encountered problem in chemistry, and a number of strategies exist for gas-phase species with well-known thermochemical properties [32]. For example, the contributions can be obtained using a theoretical approach such as the ideal gas equation, or an empirical approach such as the Shomate equation [42]. The case of adsorbates is less well defined, and it is often necessary to resort to a relatively crude approximation such as neglecting entropy and integrated heat capacity [43]. In cases where more information is available, it is possible to assume that all degrees of freedom can be approximated by vibrational modes; the vibrational partition function of statistical mechanics can then be applied to approximate the thermal enthalpy and entropy [44]. More advanced corrections are also possible [45]. The important thing for an automated implementation is to recognize that the method used may vary between models, and hence should be properly abstracted. The ‘‘Thermodynamics’’ class in CatMAP stores all information pertaining to finite temperature corrections.

DEa,2

DEa,3

ECO*

0.62 0.61

1.68 1.84

-0.99 -1.14

EO* 1.625 1.555

EO–O* 5.36 5.35 3.8

EO–CO* 1.355 1.515

Rh(111)

-1.60

-3.63

-0.94

0.74

-1.35

0.555

Ru(111)

-1.64

-4.87

-1.39

0.17

-1.39

-0.065

3.35

-0.155

0.405

Cu(111)

-0.36

-2.59

0.01

1.82

-0.11

1.075

4.75

1.495

Ag(111)

0.05

-0.63

1.25

2.69

0.3

2.055

5.99

2.365

Au(111)

0.10

0.49

2.49

3.38

0.35

2.615

7.23

3.055

CatMAP

Different correction schemes can be specified for gas-phase species and adsorbates. Currently gas-phase free energies may be estimated using Shomate parameters [42], the ideal gas approximation [44], or various simple schemes that assume fixed values for entropy/enthalpy. The free energy of adsorbates can be approximated using the harmonic approximation that treats all degrees of freedom as vibrational modes [44], or schemes that assume fixed values for entropy/enthalpy. The ‘‘Thermodynamics’’ class also handles the treatment of lateral adsorbate–adsorbate interactions, although current support is only experimental. Another key feature of the ‘‘Thermodynamics’’ class is the inclusion of corrections necessary for modeling electrochemical environments. Electrochemical reactions that involve the transfer of a proton–electron pair can be described in the framework of the Computational Hydrogen Electrode (CHE) model [35]. In this treatment, the definition of the reversible hydrogen electrode (RHE) 1 Hþ þ e $ H2ðgÞ 2 at 25 C, 1 bar H2(g), and all values of pH is used to define zero voltage. Applied bias vs. RHE, U, and an electrochemical transfer coefficient, b, can be entered into CatMAP as reaction parameters. In this framework, the free energy of the proton-electron pair is simply the free energy of H2(g) - eU, where e is the elementary positive charge. The free energy of any transition state containing a protonelectron pair is similarly corrected by -ebU. A simple scheme to account for the influence of solvation in electrochemical systems is also implemented as user-customizable static solvation corrections to adsorbates (e.g. CO* is stabilized by 0.1 eV due to weak hydrogen bonding and OH* is stabilized by 0.5 eV due to a strong interaction with the solvent). 4.1 Example: Free Energies for CO Oxidation For the example of CO oxidation considered here we make the simplest approximations. Zero-point energy and integrated heat capacity are neglected for gas and adsorbed species, entropy is neglected for adsorbed species, and all gas-phase species are assumed to have entropy of 0.002 eV/K: Ggas ¼ Egas ð0:002ÞT Gads ¼ Eads where Gi is the free energy of species i (gas = gas-phase, ads = adsorbate/transition-state) in eV, Ei is the ‘‘formation energy’’ of species i in eV, and T is the temperature in Kelvin. Despite the primitive nature of these approximations, we will show that it is still possible to accurately describe trends. This is primarily due to the fact that these

799

contributions are relatively minor, and cancellation of error occurs in relative quantities.

5 Constructing Scaling Relations The advent of descriptor-based analyses in heterogeneous catalysis was made possible by the discovery of correlations between adsorption energies of atomic species, intermediate species, and transition-states. These correlations, called ‘‘scaling relations’’ take many forms, including linear, piecewise linear, or injective [16, 20, 46]. CatMAP currently includes only linear scaling, but it is possible to include more complex correlations. The key feature of a ‘‘scaling relation’’ is that it provides a functional map between the descriptors and the energetics of the reaction scheme. Previous work on heterogeneous catalyst surfaces has shown two general classes of linear scaling relations: adsorption energy scaling relations, and transition-state scaling relations. Adsorption energy scaling relations express the energies of adsorbed intermediates as a linear function of the atomic species that they bind through [16]. This mapping also ensures linear correlations between the energies of all intermediates. Transition-state scaling relations express the energy of a transition-state species as a linear function of the initial or final state associated with it, or the activation energy of an elementary reaction as a linear function of its reaction energy (Brønsted–Evans–Polanyi relations) [20, 21]. The adsorption energy scaling relations provide a link between the ‘‘descriptors’’ (atomic adsorption energies or sometimes the adsorption energies of other intermediates) and the thermochemistry of the pathway, while the transition-state scaling relations connect the thermochemistry and kinetics of the pathway. It is clear that the scaling relations provide a link between the descriptors and the energies of adsorbed/transition states, although the procedure for obtaining all necessary energies requires several steps and is not ideal for automation. This can be improved by reformulation of the adsorption/transition-state energy scaling relations into a unified linear formation energy scaling relationship. This is achieved by recognizing that the transition-state energies can be given in terms of formation energies in an equivalent manner to the adsorbed states. Thus, all energies can be expressed as a vector, p, of length N where N is the number of adsorbed states plus the number of transition-states. The M descriptors are represented by a vector, d, of length M. In this scheme it is easy to see that there exists a linear map between the descriptor vector and the energy vector. This linear map is equivalent to an N 9 M matrix C such that: p ¼ Cd

ð7Þ

123

800


It is noted that in this formalism the descriptor vector will need to include one constant component to account for the intercept in the linear scaling relationships. For convenience this constant is normally 1, so that the coefficients of this constant term correspond directly to the intercepts of the linear equations. Furthermore, the energy vector p can be more generally interpreted as a ‘‘parameter vector’’ which could contain other input parameters to the kinetic model that might vary as a function of the descriptor space (interaction parameters, etc.). The final task necessary in this framework is determination of the ‘‘coefficient matrix’’ C. This can also be achieved using the mechanics of linear algebra. We can define a matrix of raw energies/parameters, R, which will be of size N 9 S where S is the number of total surfaces for which raw data exists. We further define a ‘‘descriptor matrix’’, D of size M 9 S, which consists of the raw values of the descriptors for each of the surfaces of interest. This leads to the following relation: R ¼ CD

ð8Þ

which represents the least-squares fitting normally used to obtain the scaling parameters. The coefficient matrix C which provides the optimal least-squares mapping between R and D is given by: C ¼ RDþ

ð9Þ

where D? is the Moore–Penrose pseudo-inverse of D [47, 48].

[50]

123

While this formulation is a cleaner solution from the standpoint of automation, it is not always ideal. In many cases the number of surfaces available to construct the fit is relatively low compared to the number of parameters. In this case the procedure outlined above may become numerically unstable, leading to ‘‘optimal’’ fits that are unphysical. For this reason, it is often necessary to include some constraints in the process of computing C. For example, intuition may dictate that a certain energy value depends only on one of the descriptors. In this case it is desirable to constrain the coefficients of the other descriptors to be zero. It is also intuitive that a surface that is more reactive towards a ‘‘descriptor’’ adsorbate should be more reactive towards other similar adsorbates, so that it is necessary to constrain the coefficient to be a positive value. Numerous other constraints might be expected from the physics of the problem [16, 49], which calls for a method of solving Eq. (8) subject to constraints on the values of C. One solution to the problem is to use a constrained minimization algorithm [50]. The algorithm is an iterative method that allows specification of an upper/lower bound for any value in the matrix C. The method is relatively simple, is fast for matrices of the size investigated in these applications (N * 10–100), and has been implemented in CatMAP. By treating scaling relations in this ‘‘unified linear formation energy scaling’’ framework it is possible to obtain the coefficient matrix through a well-defined algorithm, and use this coefficient matrix to quickly and easily move from ‘‘descriptor space’’ to ‘‘parameter space’’ via Eq. (8).

CatMAP

801

Another possibility is to recognize that Eq. (8) encapsulates the ‘‘adsorption energy scaling relations’’ [16] and the ‘‘transition-state scaling relations’’ [20]. The adsorption energy scaling relations relate atomic adsorption energies to the adsorption energies of other intermediates as: Eads ¼ aEd þ b

ð10Þ

where a is a slope obtained by linear regression or bondcounting arguments [16], Ed is the atomic adsorption energy of the atom through which species i is bound, and b is the intercept which relates to the site geometry. Meanwhile, transition-state scaling is defined as: ETS ¼ cEdiss þ n

ð11Þ

where Ediss is the energy of the dissociated state, and c and n are typically obtained by regression or from tabulated values [20, 21]. Since Ediss will always be a sum of adsorbed states, we can combine Eqs. (10) and (11) to give:

ETS ¼

X

c ð a i Ed þ b i Þ þ n i X X ¼ ca þ cb þ n E i d i i i

ð12Þ

This is analogous to Eq. (10), and is easily generalizable to the case of multiple descriptors (e.g. multiple atomic adsorption energies, Ed). Equations (10) and (12) can be rewritten in matrix/vector form to give Eq. (8). Hence, the coefficient matrix C can be determined using established linear-regression approaches or in some cases even extracted from tables in the literature (although it is noted that the constant term will depend on the atomic reference set used so care should be taken when using tabulated values). 5.1 Example: Scaling Coefficients for CO Oxidation For the example of CO oxidation we will use the CO* and O* binding energies as descriptors. In this case the adsorbate

123

802


scaling parameters are trivial since the only adsorbates, CO* and O* are descriptors of the system. This implies that aCO = aO = 1 and bCO = bO = 0 (in the notation of Eq. 10). Thus the only remaining task is to relate the transition-state energies of O-CO* and O–O* to the energies of the descriptors. First we determine these coefficients using linear ‘‘transition-state scaling’’, and show how to use this information to construct a ‘‘coefficient matrix’’, C (Eq. 8). The transition-state scaling coefficients can be obtained using a linear regression between the dissociated state energy and the transition-state energy (e.g. EOCO ¼ cOCO ðEO þ ECO Þ þ nOCO and EOO ¼ cOO ð2EO Þ þ nOO ) so that: 2 3 2 3 2 3 1 ECO 0 0 ECO 6 EO 7 6 0 7 0 74 1 6 7 6 5 ð12Þ 4 EOCO 5 ¼ 4 cOCO cOCO nOCO 5 EO 1 EOO 0 2cOO nOO The coefficients are determined using the values in Table 3 and linear regression, and the numerical values of the transition-state scaling coefficient matrix (CTS) are given in Eq. (13). For comparison, we also compute the coefficient matrix using multi-dimensional linear regression (CLR) as discussed in the previous section. The details are omitted for brevity, but the coefficient matrices and resulting scaling lines are shown in Eq. (13) and Fig. 2. 2 3 2 3 1 0 0 1 0 0 6 0 6 0 1 0 7 1 0 7 7 6 7 CTS ¼ 6 4 0:68 0:68 0:94 5; CLR ¼ 4 0:45 0:87 0:55 5 0 1:43 3:20 0:13 1:35 3:39 ð13Þ Using these scaling relations we can compute the relevant energetics of CO oxidation as a continuous function of the CO* and O* descriptor space. We note that since the gas-phase formation energies are not computed via scaling relations (i.e. they are constants) the gas-phase reaction energies are exactly constant regardless of the descriptor values. Free energy diagrams obtained from scaling relations at several points in descriptor space (with the transition-state scaling coefficient matrix, CTS, from Eq. 13) are compared to the raw data in Fig. 3. These free energy diagrams and scaling relations form the basis of the subsequent kinetic modeling.

the CatMAP implementation and discuss its implications and the method of solution. The reaction rates of the network of elementary steps are obtained by solving a mean-field microkinetic model to steady state. The differential equations governing the system are of the form: Y Y Y Y ri ¼ k þ hij pij k hij pil ð14Þ i i j

ohi ¼ ot

X

j

sij rj

l

l

ð15Þ

j

where ri are the rates of each elementary step, k i are the forward/reverse rate constants, hij and hil are the surface concentrations of the adsorbed reactants and products for elementary step i, pij and pil are the unitless pressures of the gas-phase reactants and products for elementary step i, and sij are coefficients for the stoichiometry of species i in elementary step j. This set of coupled non-linear ordinary differential equations has a steady-state solution such that: ohi ¼0 ot

ð16Þ

which must be solved subject to the site conservation constraint: X hi ¼ hTOT ð17Þ i

where hTOT is the normalized surface area. For the application of catalysis, the steady-state solution is relevant since it governs the long-term behavior of the catalyst. More complex types of kinetic models and solutions are possible [1, 3, 6], and CatMAP could include more advanced approaches in future versions. Some extensions to the simple expressions given in Eqs. (14)–(17) allow for a more sophisticated description of the surface reaction without significantly affecting the complexity. The first extension is the possibility of multiple sites. There are many examples where intermediates preferentially adsorb at different sites, leading to non-competitive adsorption [51–53]. In these cases it is desirable to enforce site conservation differently for different adsorbates: X hi ¼ hTOT ð18Þ j i2Sj

6 Solving the Kinetic Model The construction and use of microkinetic models is a wellestablished field within chemical engineering [1–3], and hence will not be covered in detail. Instead we briefly outline the general type of kinetic model that is applied for

123

where Sj is the set of indices corresponding to site j, and hTOT is the maximum possible coverage of site j. j In this case the solution of Eqs. (14)–(16) is no more complex, but significantly more physics can be captured. Another example of a relatively simple extension of Eqs. (14)–(17) is to allow the rate constants, k, to be

CatMAP

803

Fig. 2 Comparison of scaling lines determined from transition-state scaling (red) and multi-dimensional linear regression (blue). The transition-state scaling lines are more robust, but the direct fit gives lower mean absolute errors

Fig. 3 CO oxidation free energy diagrams computed for the (111) facet of Au (blue), Pt (red), and Ru (green). Solid lines correspond to energetics obtained using scaling relations, and dashed lines represent the data in Table 3

dependent on coverages. This allows at least semi-quantitative inclusion of the effect of adsorbate–adsorbate interactions, while retaining the general form of the mean-field equations. Allowing rate constants to be coverage dependent requires choosing a model for coverage dependence of intermediate energies and parameterization of the chosen model, which is not a trivial task and has been addressed elsewhere [51, 52]. These extensions will not be discussed in further detail here, but the CatMAP implementation is compatible with both multi-site models and models with coverage-dependent rate constants, though this will increase the computational overhead somewhat (at least a factor of N). The key assumptions associated with the mean-field model are: homogenous distribution of adsorbates, fast diffusion, and static active sites. In some situations it may be necessary to move to a more complex model in order to capture the physics at the surface. However, the mean-field

model and steady-state solutions have been found to be sufficiently accurate to describe trends in catalysis [10], and can easily and quickly be solved in an automated fashion. In an automated implementation it is fairly easy to construct expressions such as Eqs. (14) and (15) given some elementary steps as input. Yet, there are some considerations to take into account. The rate constants in Eq. (14) will, in general, vary by many ([20) orders of magnitude due to the exponential factor in the Arrheniustype equations. This leads to extreme numerical instability in double precision floating point arithmetic. It is often possible to eliminate or reduce this ‘‘stiffness’’ by intuition and mathematics; however, in an automated scheme this becomes considerably more difficult. The more general solution is to move to multiple-precision floating-point arithmetic. The CatMAP implementation uses this solution, relying on the Python mpmath library [54] in order to represent rate constants, coverages, rates, and intermediate quantities to an arbitrary precision. The required decimal precision varies from system to system, but in general a precision of 10-50 is sufficient. The computational cost associated with this more accurate representation of numbers is relatively low when the gmpy Python library is used to conduct the calculations in C. The steady-state solution to Eqs. (14) and (15) can be obtained using a numerical root-finding procedure. The simplest and most widely used method for multi-dimensional root-finding is Newton’s method. In the current application it has been found that Newton’s method exhibits good convergence behavior, although some minor modifications improve the speed and stability. One necessary modification is constraining the solutions for the coverages to conform to site conservation. This is achieved by defining a function that ensures that no coverages are negative, and that the sum of coverages is not greater than

123

804


hTOT. The function is applied each iteration of Newton’s method in order to ensure that the solutions conform to this requirement. It was also found that, unsurprisingly, convergence was faster and more stable when an analytical Jacobian of Eq. (15) was used. The analytical Jacobian for a single-site model with coverage-independent rate constants is given by: oh_ i X or i ¼ sij ohn ohn j

ð19Þ

Y X ohim Y Y X ohiq Y ori ¼ kþ pij hij k pil hil i i ohn ohn j6¼m ohn l6¼q m q j l ð20Þ ohim ¼ ohn

hTOT ifm ¼ dmn otherwise

ð21Þ

where dmn is the Kronecker delta. Extending Eqs. (19)– (21) to the case of a multi-site model is fairly simple, while the coverage-dependent case is considerably more complex and requires knowledge of the specific coverage-dependence model. Finally, it is worth noting that in the implementation Eqs. (14)–(21) must be automatically generated from the reaction expressions. Naturally this takes some computational time, and it is extremely advantageous to construct these functions one time and store/compile them, as opposed to constructing them at each iteration of Newton’s method. 6.1 Example: Rate Map for CO Oxidation The CatMAP package was employed to calculate the steady-state rates for CO oxidation as a function of the CO* and O* descriptor space based on the discussed approximations and assumptions. The resulting rate and coverage maps are shown in Fig. 4. The analysis correctly predicts that Pt, Pd, and Rh are the best catalysts, in agreement with experiment and previous analyses [38, 39].

7 Obtaining Initial Guesses The numerical solution of Eqs. (14)–(16) requires a sufficiently accurate initial guess of the coverage vector in order for Newton’s method to converge. Determining this initial guess by trial and error is highly impractical for most problems, and hence an automated solver should employ some systematic strategies. The first strategy that has proven successful is to use the Gibbs distribution in order to obtain an initial guess:

123

i exp G kB T hi ¼ P Gj j exp k T

ð22Þ

B

where Gi is the free energy of formation for species i, kB is Boltzmann’s constant, and T is temperature in Kelvin. It is intuitive that this strategy should provide a good guess since the steady-state coverages should approach the thermodynamic equilibrium coverages when kinetics effects do not play a key role. There is still some arbitrariness in the Boltzmann coverages, however, since in order to compute them it is necessary to assume that each atomic species is in thermal equilibrium with some reservoir [55] (i.e. the free energies must be calculated using a set of ‘‘atomic reference energies’’, see previous discussion and Eq. 1). It is relatively easy to have the automated implementation assume all possible cases and attempt to use them as an initial guess, and typically one is successful. The second strategy relies on using information from other similar systems. In a descriptor-based analysis this corresponds to using information from a nearby point in ‘‘descriptor space’’. The simplest way of using this information is to take the converged solutions from nearby points in descriptor space and attempt to use them as an initial guess. It was found that this is usually successful, but that in some cases it was necessary to use a ‘‘descriptorspace bisection’’ algorithm to obtain a successful initial guess using the solution from a nearby point. This simple algorithm attempts to use the solution from a successful point at a distance halfway between the successful point and the point of interest. This is repeated until a solution is obtained for a point closer to the point of interest, and then the bisection algorithm is re-applied from the new point. This approach can be somewhat time consuming, but it is typically only necessary for a limited number of points, and nearly always results in a valid solution. As more complex systems are studied, it may be necessary to implement more advanced strategies such as homotopy continuation in order improve the robustness and efficiency of mapping the descriptor space [56]. There are also other strategies that can be imagined for determining initial guesses (Langmuir isotherms, experimental knowledge, etc.) but these are more difficult to implement in an automated way. Nevertheless, new strategies could be developed, and even the strategies outlined above provide a multitude of possible initial guesses. This leaves the problem of how to determine the order in which to try initial guesses obtained by various strategies. The order will not affect the overall outcome of whether or not a valid solution can be obtained, but it can significantly affect the amount of time needed. One

CatMAP

805

(a)

(b)

(c)

Fig. 4 Turnover frequencies for CO oxidation on a range of metal (111) surfaces (a) and corresponding steady-state coverages of O* (b) and CO* (c) as a function of descriptor space. Rates and coverages

are calculated using CatMAP with the assumptions outlined in previous sections (T = 473 K, PCO = 1 bar, PO2 = 1/3 bar, PCO2 = 1/3 bar)

intuitive strategy is to compute the residual of each initial guess and begin by using the guess with the lowest residual. It is found that this is indeed successful in reducing the amount of time needed, although there is not a strict correlation between the initial residual and whether or not the model converges (i.e. it is possible that the initial guess with the highest initial residual is the only one which leads to convergence). Another issue associated with obtaining initial guesses, and the root-finding solution method in general, is the possibility of multiple steady-state solutions. A signature of this would be a discontinuity in the coverages as a function of descriptor space. This is a phenomenon that is rarely observed in systems studied using CatMAP. One notable exception is the case where a trivial solution exists (e.g. hi = 1, hi=j = 0). In such cases the solver will often converge to a solution that is obviously erroneous since the resulting rates are exactly zero. This is typically the result of an incomplete reaction network, and can be circumvented by adding elementary steps that remove the trivial solution. As new numerical methods for solving for kinetic steady states (e.g. homotopy continuation [56], advanced root-finding [57], numerical integration [58], or utilization of independent solvers [59]) are implemented it may become possible to circumvent these heuristics. The CatMAP package is designed to accommodate the implementation of advances via the ‘‘Mapper’’ and ‘‘Solver’’ classes, but even the simple strategies currently implemented in CatMAP are sufficiently fast and robust to enable the solution of relatively complex kinetic systems [53]. The specific time required for solving a CatMAP model is a complex function that depends on a variety of factors including the complexity of the reaction network, inclusion of coveragedependent energetics, the resolution and range over which the descriptor space is sampled, the convergence properties of the kinetic equations, and computational architecture. However, as a rough estimate, simple systems (N \ 10) can be solved with CatMAP on the order of minutes, while more complex systems (N [ 50) may require hours. In all

cases, the solution of mean-field models will require much less time than generating the energetic inputs from electronic structure theory or simulating kinetics with stochastic methods such as kinetic Monte Carlo.

8 Conclusion A general approach to creating descriptor-based microkinetic analyses is presented. The approach consolidates knowledge from many other works into a comprehensive overview of the methodology, showing that there is a clear path to relate ‘‘descriptors’’ of materials to their catalytic properties. An automated implementation of the approach, CatMAP, is introduced and discussed. These developments will lead researchers in the field of catalysis towards a more standard and easily understood methodology for producing these descriptor-based microkinetic models that have proved fruitful in the analysis of catalytic trends and the discovery of novel catalysts. Access to automated tools such as CatMAP will allow for increased creativity in the scientific process of designing new catalyst materials. Acknowledgments We acknowledge David Fenning and Maureen Tang for their help in naming CatMAP. Support from the DOE Office of Basic Energy Science to the SUNCAT Center for Interface Science and Catalysis is gratefully acknowledged. AJM is grateful for support by the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program.

References 1. Dumesic JA, Rudd DF, Aparicio LM et al (1994) The microkinetics of heterogeneous catalysis. AIChE J 40:1085–1087 2. Davis ME, Davis RJ (2003) Fundamentals of chemical reaction engineering. Chapter 7. McGraw-Hill, New York 3. Deutschmann O (2011) Modeling and simulation of heterogeneous catalytic reactions: from the molecular process to the technical system. Wiley, Weinheim 4. Stoltze P (2000) Microkinetic simulation of catalytic reactions. Prog Surf Sci 65:65–150. doi:10.1016/S0079-6816(00)00019-8

123

806 5. Reuter K, Scheffler M (2006) First-principles kinetic Monte Carlo simulations for heterogeneous catalysis: application to the CO oxidation at RuO2(110). Phys Rev B 73:045433. doi:10.1103/ PhysRevB.73.045433 6. Hoffmann MJ, Matera S, Reuter K (2014) kmos: a lattice kinetic Monte Carlo framework. Comput Phys Commun 185:2138–2150. doi:10.1016/j.cpc.2014.04.003 7. Honkala K, Hellman A, Remediakis IN et al (2005) Ammonia synthesis from first-principles calculations. Science 307:555–558 8. Andersson MP, Bligaard T, Kustov A et al (2006) Toward computational screening in heterogeneous catalysis: pareto-optimal methanation catalysts. J Catal 239:501–506. doi:10.1016/j. jcat.2006.02.016 9. Cheng J, Hu P (2008) Utilization of the three-dimensional volcano surface to understand the chemistry of multiphase systems in heterogeneous catalysis. J Am Chem Soc 130:10868–10869. doi:10.1021/ja803555g 10. Nørskov JK, Bligaard T, Rossmeisl J, Christensen CH (2009) Towards the computational design of solid catalysts. Nat Chem 1:37–46. doi:10.1038/nchem.121 11. Nørskov JK, Abild-Pedersen F, Studt F, Bligaard T (2011) Density functional theory in surface chemistry and catalysis. Proc Natl Acad Sci USA 108:937–943. doi:10.1073/pnas.1006652108 12. Grabow LC, Mavrikakis M (2011) Mechanism of methanol synthesis on Cu through CO2 and CO hydrogenation. ACS Catal 1:365–384. doi:10.1021/cs200055d 13. Rankin RB, Greeley J (2012) Trends in selective hydrogen peroxide production on transition metal surfaces from first principles. ACS Catal 2:2664–2672. doi:10.1021/cs3003337 14. Wang H, Schneider WF (2012) Comparative chemistries of CO and NO oxidation over RuO2 (110): insights from first-principles thermodynamics and kinetics. Mol Simul 38:615–630. doi:10. 1080/08927022.2012.671521 15. Vendelbo SB, Elkjær CF, Falsig H et al (2014) Visualization of oscillatory behaviour of Pt nanoparticles catalysing CO oxidation. Nat Mater 13:884–890. doi:10.1038/nmat4033 16. Abild-Pedersen F, Greeley J, Studt F et al (2007) Scaling properties of adsorption energies for hydrogen-containing molecules on transition-metal surfaces. Phys Rev Lett 99:016105. doi:10. 1103/PhysRevLett.99.016105 17. Logadottir A, Rod TH, Norskov JK et al (2001) The Brønsted– Evans–Polanyi Relation and the volcano plot for ammonia synthesis over transition metal catalysts. J Catal 197:229–231. doi:10.1006/jcat.2000.3087 18. Michaelides A, Liu Z-P, Zhang CJ et al (2003) Identification of general linear relationships between activation energies and enthalpy changes for dissociation reactions at surfaces. J Am Chem Soc 125:3704–3705. doi:10.1021/ja027366r 19. van Santen RA, Neurock M, Shetty SG (2010) Reactivity theory of transition-metal surfaces: a Brønsted-Evans-Polanyi linear activation energy-free-energy analysis. Chem Rev 110:2005–2048. doi:10.1021/cr9001808 20. Wang S, Temel B, Shen J et al (2010) Universal Brønsted-EvansPolanyi relations for C-C, C–O, C–N, N–O, N–N, and O–O dissociation reactions. Catal Lett 141:370–373. doi:10.1007/ s10562-010-0477-y 21. Wang S, Petzold V, Tripkovic V et al (2011) Universal transition state scaling relations for (de)hydrogenation over transition metals. Phys Chem Chem Phys 13:20760–20765. doi:10.1039/ c1cp20547a 22. Medford AJ, Vojvodic A, Hummelshøj JS, et al (2015) From the Sabatier principle to a predictive theory of transition-metal heterogeneous catalysis. J Catal (Accepted) 23. Jacobsen CJH, Dahl S, Clausen BS et al (2001) Catalyst design by interpolation in the periodic table: bimetallic ammonia

123


24.

25. 26.

27.

28.

29. 30.

31.

32.

33.

34.

35.

36.

37.

38.

39.

40.

41.

42.

synthesis catalysts. J Am Chem Soc 123:8404–8405. doi:10. 1021/ja010963d Studt F, Abild-Pedersen F, Wu Q et al (2012) CO hydrogenation to methanol on Cu–Ni catalysts: theory and experiment. J Catal 293:51–60. doi:10.1016/j.jcat.2012.06.004 Greeley J, Mavrikakis M (2004) Alloy catalysts designed from first principles. Nat Mater 3:810–815. doi:10.1038/nmat1223 Greeley J, Jaramillo TF, Bonde J et al (2006) Computational high-throughput screening of electrocatalytic materials for hydrogen evolution. Nat Mater 5:909–913. doi:10.1038/nmat1752 Hansgen DA, Vlachos DG, Chen JG (2010) Using first principles to predict bimetallic catalysts for the ammonia decomposition reaction. Nat Chem 2:484–489. doi:10.1038/nchem.626 Hummelshøj JS, Abild-Pedersen F, Studt F et al (2012) CatApp: a web application for surface chemistry and heterogeneous catalysis. Angew Chem Int Ed Engl 51:272–274. doi:10.1002/anie. 201107947 Campbell CT (1994) Micro- and macro-kinetics : their relationship in heterogeneous catalysis. Top Catal 1:353–366 Stegelmann C, Andreasen A, Campbell CT (2009) Degree of rate control: how much the energies of intermediates and transition states control rates. J Am Chem Soc 131:8077–8082. doi:10. 1021/ja9000097 Brown WA, Kose R, King DA (1998) Femtomole adsorption calorimetry on single-crystal surfaces. Chem Rev 98:797–832. doi:10.1021/cr9700890 Johnson RDI (2011) NIST Computational Chemistry Comparison and Benchmark Database. NIST Stand Ref, Database Number 101 Enkovaara J, Rostgaard C, Mortensen JJ et al (2010) Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. J Phys 22:253202. doi:10.1088/0953-8984/22/25/253202 Hammer B, Hansen L, Nørskov J (1999) Improved adsorption energetics within density-functional theory using revised PerdewBurke-Ernzerhof functionals. Phys Rev B 59:7413–7421. doi:10. 1103/PhysRevB.59.7413 Nørskov JK, Rossmeisl J, Logadottir A et al (2004) Origin of the overpotential for oxygen reduction at a fuel-cell cathode. J Phys Chem B 108:17886–17892. doi:10.1021/jp047349j Klu¨pfel S, Klu¨pfel P, Jo´nsson H (2012) The effect of the PerdewZunger self-interaction correction to density functionals on the energetics of small molecules. J Chem Phys 137:124102. doi:10. 1063/1.4752229 Peterson AA, Abild-Pedersen F, Studt F et al (2010) How copper catalyzes the electroreduction of carbon dioxide into hydrocarbon fuels. Energy Environ Sci 3:1311. doi:10.1039/c0ee00071j Jiang T, Mowbray DJ, Dobrin S et al (2009) Trends in CO oxidation rates for metal nanoparticles and close-packed, stepped, and kinked surfaces. J Phys Chem C 113:10548–10553. doi:10. 1021/jp811185g Grabow LC, Hvolbæk B, Nørskov JK (2010) Understanding Trends in catalytic activity: the effect of adsorbate-adsorbate interactions for CO oxidation over transition metals. Top Catal 53:298–310. doi:10.1007/s11244-010-9455-2 Mason S, Grinberg I, Rappe A (2004) First-principles extrapolation method for accurate CO adsorption energies on metal surfaces. Phys Rev B 69:161401. doi:10.1103/PhysRevB.69.161401 Abild-Pedersen F, Andersson MP (2007) CO adsorption energies on metals with correction for high coordination adsorption sites— a density functional study. Surf Sci 601:1747–1753. doi:10.1016/ j.susc.2007.01.052 Shomate CH (1954) A method for evaluating and correlating thermodynamic data. J Phys Chem 58:368–372. doi:10.1021/ j150514a018

CatMAP 43. Nørskov JK, Studt F, Abild-Pedersen F, Bligaard T (2014) Fundamental concepts in heterogeneous catalysis. Wiley, New York 44. McQuarrie D (2000) Statistical mechanics. University Science Books, Sausalito 45. Karplus M, Kushick JN (1981) Method for estimating the configurational entropy of macromolecules. Macromolecules 14:325–332. doi:10.1021/ma50003a019 46. Vojvodic A, Calle-Vallejo F, Guo W et al (2011) On the behavior of Brønsted-Evans-Polanyi relations for transition metal oxides. J Chem Phys 134:244509. doi:10.1063/1.3602323 47. Moore EH (1920) On the reciprocal of the general algebraic matrix. Bull Am Math Soc 26:385–397 48. Penrose R, Todd JA (1955) A generalized inverse for matrices. Math Proc Camb Philos Soc 51:406. doi:10.1017/S0305004 100030401 49. Hammer B, Norskov JK (1995) Why gold is the noblest of all the metals. Nature 376:238–240. doi:10.1038/376238a0 50. Axelsson O (1996) Iterative solution methods. Cambridge University Press, Cambridge 51. Xu Y, Lausche AC, Wang S et al (2013) In silico search for novel methane steam reforming catalysts. New J Phys 15:125021. doi:10.1088/1367-2630/15/12/125021 52. Lausche AC, Medford AJ, Khan TS et al (2013) On the effect of coverage-dependent adsorbate–adsorbate interactions for CO methanation on transition metal surfaces. J Catal 307:275–282. doi:10.1016/j.jcat.2013.08.002

807 53. Medford AJ, Lausche AC, Abild-Pedersen F et al (2014) Activity and selectivity trends in synthesis gas conversion to higher alcohols. Top Catal 57:135–142. doi:10.1007/s11244-013-0169-0 54. Johansson F (2010) mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.14) 55. Reuter K, Scheffler M (2003) First-principles atomistic thermodynamics for oxidation catalysis: surface phase diagrams and catalytically interesting regions. Phys Rev Lett 90:046103. doi:10.1103/PhysRevLett.90.046103 56. Li T-Y (1997) Numerical solution of multivariate polynomial systems. Acta Numer 6:399–436. doi:10.1017/S096249290000 2749 57. Gusma˜o GS, Christopher P (2014) A general and robust approach for defining and solving microkinetic catalytic systems. AIChE J. doi:10.1002/aic.14627 58. Chandler JP, Hill DE, Spivey HO (1972) A program for efficient integration of rate equations and least-squares fitting of chemical reaction data. Comput Biomed Res 5:515–534. doi:10.1016/ 0010-4809(72)90058-4 59. Coltrin ME, Kee RJ, Rupley FM (1991) Surface chemkin: a general formalism and software for analyzing heterogeneous chemical kinetics at a gas-surface interface. Int J Chem Kinet 23:1111–1128. doi:10.1002/kin.550231205

123