Following Molecules through Reactive Networks: Surface Catalyzed

Subscriber access provided by CLEMSON UNIV

Article

Following Molecules through Reactive Networks: Surface-Catalyzed Decomposition of Methanol on Pd(111), Pt(111), and Ni(111) Zeb C. Kramer, Xiang-Kui Gu, Dingyu D.Y. Zhou, Wei-Xue Li, and Rex T. Skodje J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/jp503056u • Publication Date (Web): 20 May 2014 Downloaded from http://pubs.acs.org on May 29, 2014

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

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

Page 1 of 58

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

The Journal of Physical Chemistry

Following Molecules through Reactive Networks: Surface-Catalyzed Decomposition of Methanol on Pd(111), Pt(111), and Ni(111) Zeb C. Kramer,a Xiang-Kui Gu,b Dingyu D. Y. Zhou,a Wei-Xue Li,*b and Rex T. Skodje*a

a) Department of Chemistry and Biochemistry, University of Colorado, Boulder, CO 80309 b) Key State Laboratory of Catalysis, Dalian Institute for Chemical Physics, Dalian 116023, China Abstract:

We present a model of the surface kinetics of the dehydrogenation reaction of

methanol on the Pd(111), Pt(111), and Ni(111) metal surfaces.

The mechanism consists of ten

reversible dehydrogenation reactions that lead to the final products of CO and H2.

The rate

coefficients for each step are calculated using ab initio transition state theory that employs a new approach to obtain the symmetry factors.

The potential energies and frequencies of the reagents

and transition states are computed using plane wave DFT with the PW91 exchange correlation functional.

The mechanism is investigated for low coverages using a global sensitivity analysis

that monitors the response of a target function of the kinetics to the value of the rate coefficients. On Pd(111) and Ni(111), the reaction COH → CO + H is found to be rate limiting and overall rates are highly dependent upon the decomposition time of the COH intermediate.

Reactions at

branches in the reaction network are also particularly important in the kinetics. A stochastic atom-following approach to pathway analysis is used to elucidate both the pathway probabilities in the kinetics and the dependence of the pathways on the values of the key rate coeﬃcients of the mechanisms.

On Pd(111) and Ni(111) there exists significant competition between the

pathway containing the slow step and faster pathways that bypass the slow step. A discussion is given of the dependence of the model target’s probability density function on the chemical pathways. RTS: WXL:

[email protected] [email protected]

1 ACS Paragon Plus Environment


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

Page 2 of 58

Keywords: Kinetics, catalysis, pathway, sensitivity analysis, methanol, density functional analysis.

I.

Introduction The development of kinetic mechanisms that can accurately describe large networks

of chemical reactions is typically an iterative process that requires repeated revision before satisfactory agreement with experiment is achieved. Since many model parameters come with large uncertainties due to limitations in theory or experiment, it can be difficult to determine which parts of the mechanism are in the most immediate need for improvement.

Such obsticals

are particularly apparent for the treatment of surface catalyzed reactions where in situ measurements of species concentrations and computational determination of rate coefficients are quite challenging.

The possible participation or influence of multiple crystal faces, defect sites,

and lateral interactions between adsorbates can lead to added complexities and ambiguities in the kinetic modeling.

In this work, we introduce an approach to analyze surface reactions which

we have recently employed with some success in treating gas phase combustion reactions.1

2 3 4

In this theoretical approach, the reaction rate coefficients are obtained using quantum mechanical computation of the potential energy barriers and frequencies combined with transition state theory (TST).5 analysis.

The behavior of reaction network itself is understood using global sensitivity

In a manner described in detail in Sec. III, we determine the particular rate coefficient

that most strongly influences an observable “target” function of the kinetics from its “sensitivity indices.” The sensitivity index is a function of both the direct variation of the target with respect to the rate coefficient (or “factor”) as well as to the size of its uncertainty range.

Once

the key reaction step is identified as the most sensitive, it can be improved or further studied. While the sensitivity of the target function to a particular reaction step can be identified in a mathematically well defined manner, it is not always physically obvious why a given reaction affects a target function so strongly.

Hence, in this work we introduce a “stochastic pathway

analysis” in which individual chemical moieties are followed through the chemical network.

It

is shown how the pathway analysis can provide unique insight into the origin of the sensitivity of 2 ACS Paragon Plus Environment

Page 3 of 58

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


target functions. In the present work, we investigate the dehydrogenation reaction of methanol on several catalytic metal surfaces.

While the mechanism for this reaction is relatively simple compared

to that observed for more complicated alcohols, it provides us with a computationally tractable system to test our methodology.

Furthermore, the metal surface catalyzed decomposition of

methanol is of great interest in fuel cell research6 due to methanol’s high hydrogen to carbon ratio, its availability, and its ability to be transported in liquid form.

It is an important process

in both direct methanol fuel cells (DMFCs)7 and in the chemistry of methanol steam reforming (MSR) to hydrogen gas, CH3OH + H2O → 3H2+ CO2. Both DMFCs and the MSR process have been studied as a fuel source for low emissions vehicles8 9. The decomposition of methanol is also interesting because it shares many chemical characteristics with more complex carbohydrates used in fuel cell research10

11

and its surface decomposition mechanism may

include many of the same elementary reactions as the surface catalyzed mechanisms involving more complicated molecules. Therefore, surface reactions of methanol are of practical interest and serve as a prototype system for larger fuel cell components. There have been many studies, both experimental and theoretical, of the surface catalyzed decomposition of methanol on various metal surfaces, particularly Pt, Pd, and Cu, even on Pd nanocrystals14.

11 12 13

and

The mechanism of the decomposition reaction of methanol on

Pd(111), Pt(111), and Ni(111) has been characterized using density functional theory (DFT) in several studies15

11 8 16

.

The reaction is known to proceed through a sequence of reversible

hydrogen transfer reactions that ultimately lead to the formation of CO and H2 products which can desorb from the surface. Greeley and Mavrikakis11 have noted that there may be competition between the reaction pathways in methanol decomposition, an issue that we shall further explore in this study.

As we shall discuss later, the pathway initialized by the transfer of the methyl

hydrogen from methanol typically dominates over pathways that begin with hydrogen extraction from the hydroxyl hydrogen.

Experimentally, it was found that on Pt(111) the methanol

preferred to undergo complete dehydrogenation to adsorbed CO and H atoms , and no methoxy intermediate could be found on the surface.17

On the other hand, for Pd(111), the methoxy

intermediate was in fact found; the methoxy ultimately decomposed to adsorbed CO and H atoms.18 With increase of partial pressure up to millibar pressures, methanol decomposition via C-O bond scission occurred on Pd(111), although the rate was far slower than that of 3 ACS Paragon Plus Environment


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

Page 4 of 58

dehydrogenation to CO and hydrogen.19

The key element of a theoretical treatment of the surface chemistry of the methanol reaction is the determination of the energies and frequencies of the transition states and reagent states using quantum chemistry. These computations are typically carried out using primarily plane wave periodic density functional theory (DFT) methods. well.

This shall be our approach as

While density functional methods have not yet achieved the level of chemical accuracy

that can be provided for small gas phase reactions of small molecules using pure quantum methods, DFT calculations are rapidly developing and even now can be progressively improved through the use of more sophisticated functionals.

Indeed, a motivation of the present work is

to identify key reaction steps that may ultimately warrant higher level and hence more expensive DFT calculations.20 Once the DFT calculations are carried out for all the transitions states (TSs) and reagent state configurations relevant to the surface chemistry, the rate coefficients are determined using transition state theory (TST) which provides a statistical approximation for the rates.

Since

TST is commonly used for surface reactions, no detailed review of this method is necessary.21

22

We do, however, point out two important issues that will require special attention to implement the method.

First, since some vibrational modes for adsorbed molecules are extremely “loose”,

the usual normal mode harmonic oscillator approximation for the partition function can become quite inaccurate.

In such circumstances we shall treat such modes as rotors.

Second, since the

crystal faces possess symmetries that lead to multiple equivalent reaction paths, the usual symmetry number used for gas phase reactions must be modified.23

A general approach to this

issue will be discussed. After the kinetic model for the dehydrogenation reaction of methanol has been constructed, we shall then switch to the main focus of the work which is the sensitivity and pathway analysis of the reaction network.

Our purpose here is to present global sensitivity and

pathway analyses in the context of a surface reaction problem, and to discuss the implications and connections between the two types of analyses on a non-trivial, but still manageably small, reaction network. Our emphasis will be to illustrate the power of the theoretical algorithm to reveal the chemical pathways.

Both global sensitivity and pathway analyses provide important 4 ACS Paragon Plus Environment

Page 5 of 58

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


insights into the kinetic model and provide a connection between kinetic behavior and the topology of the reaction network. The remainder of this paper is organized as follows.

In Section II, we construct rate

coefficients for the elementary steps in the dehydrogenation mechanism of methanol on the Pd(111), Pt(111), and Ni(111) metal surfaces.

The ten forward (dehydrogenation) and ten

backward (recombination) reactions will be modeled using TST that makes use of the DFT calculations.

In Section III, we review the global sensitivity analysis implementation and then

apply the global sensitivity analysis methods to the decomposition of methanol on the metal surfaces in order to identify the key reactions in the decomposition process.

The methanol

decomposition mechanism constructed in this work has a relatively simple structure and is therefore a useful trial system for applying the pathway analysis method.

In Section IV, we use

stochastic pathway analysis (SPA) to explore the chemical pathway picture of the methanol decomposition reaction. Also in this section, we will compare the pathway analysis to the first and second order global sensitivity analyses to provide a clear interpretation of the physical origin of the sensitivity.

A key concept introduced here is the notion of pathway switching in

which altering the rate coefficients changes the path that chemical moieties follow through the reaction network.

II.

In Section V, we summarize the conclusions of this work.

Construction of the Kinetic Model The decomposition of methanol on Pt(111) was found by Greeley and Mavrikakis11

13

to proceed by a series of H-atom elimination reactions through the intermediates CH2OH, CH3O, CH2O, CHOH, CHO, and COH.

The reaction process is initiated by either the

CH3OH→CH2OH+H or CH3OH→CH3O+H.

It was found, based on the computation of

reaction barriers, a particular set of 10 elementary reactions were most important.

In our model,

we also find these 10 reaction steps are responsible for the dehydrogenation reaction for all the catalytic surfaces studied here, i.e. Pt(111), Pd(111), and Ni(111).

We compute their

temperature dependent rate coefficients using full TST (including the temperature dependent prefactor) which does present some quantitative differences (roughly one order of magnitude in either direction) with the generic 1013/sec prefactor.

Furthermore, we include the

recombination reactions for each step, e.g. CH3O+H→ CH3OH.

A schematic diagram



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

illustrating the reaction steps included in our model is shown in Fig. 1. Fig. 1, we identify each reaction with an index from 1 to 10.

Page 6 of 58

Using the notation of

The forward rate coefficients (for

dehydrogenation) are denoted by ki (1/sec) and the backward (recombination) rate coefficients are k-i (area/sec).

We omit the explicit inclusion of the adsorption and desorption steps of the

kinetics, which may under some conditions be important, since our focus is on the behavior of the surface reactions themselves.

A. Density Functional Calculations & Energetics Periodic plane wave DFT calculations were performed using the Vienna ab initio Simulation Package (VASP)24

25

. The exchange-correlation energy was calculated within the

generalized gradient approximation (GGA)26 using the Perdew-Wang 1991 functional (PW91)27. A plane-wave basis set with a kinetic energy cutoff of 400 eV was used.

A periodic 3 × 3 unit

cell with nine atoms in each layer, corresponding to an adsorbate surface coverage of

1 of a 9

monolayer, was chosen to represent the system.

The surface Brillouin zone was sampled with a

5 × 5 × 1 Monkhorst-Pack k -point grid28.

Four layer slabs were used to model the surface

with a 12 Å vacuum separating slabs.

The top two layers of the surface and the adsorbate

atoms were allowed to relax in optimization calculations, while the bottom two slabs were held in their bulk configuration.

Transition states were calculated using the climbing image nudged

elastic band method (CI-NEB)29 with 8 images separating the initial and final states. The transition state was taken as the highest image along the minimum energy path defined by the NEB calculation. Vibrational frequencies were calculated for the equilibrium and transition state structures by diagonalizing force constant matrix, including only the adsorbate degrees of freedom, obtained from finite differences of the potential performed by VASP.

For further

details on the methodology, we refer the reader to earlier work of Gu and Li8. The equilibrium geometries for the most stable configurations of the nine chemically relevant species are shown in Fig. 2.

The structures obtained from the geometry optimization

are found to be consistent with those of Greenley and Mavrikakis11 for the Pt(111) surface.

The

top row of structures are consistent with the optimized geometries of all species on the Ni(111) and Pd(111) surfaces and all but three molecules for Pt(111).

The structures for the methoxy

radical (CH3O), the formyl radical (HCO), and formaldehyde (CH2O) were found to be different on Pt(111) compared to Ni(111) and Pd(111).

As illustrated in Fig. 2, the methoxy radical is 6

ACS Paragon Plus Environment

Page 7 of 58

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


bound to the three-fold site for Ni(111) and Pd(111), while it is bound to the on-top site for Pt(111).

The carbon of the formal radical bound to the on top site for Pt(111) but formed a

bridge structure for Pd(111) and Ni(111).

The formaldehyde binding sites is seen to form two

distinct bridge structures for Pt(111) versus Pd(111) and Ni(111).

The absolute value of the

binding energies for each species on each surface is shown in Table 1. Using the force constant matrix obtained through a finite difference approximation from the analytical gradients, the normal mode vibrational frequencies corresponding to the adsorbate degrees of freedom were obtained for all the species on each of the three catalytic surfaces. frequencies obtained are reported in the Supporting Material.

The

The vibrational normal mode

frequency are found to be reasonably similar on the different surfaces.

The high frequency

modes (>1000cm-1) are typically the same to within 5% for all surfaces while somewhat larger differences develop for the lower frequencies.

The very low frequency modes ( /?

(2.10)

using a least squares proceedure for all the forward and reverse reactions over the temperature range 500-1000 K.

Since most of the Arrhenius plots were fairly linear, we have also fit the

results to the conventional two parameter Arrhenius form, 789 = :@@ < =>/? The results obtained are reported in Table. 3.

(2.11)

The zero point and thermal effects are seen to

have a significant influence on the energy of activation.

Indeed, comparing the classical barrier

heights reported in Table 1 to the extracted activation energies we see differences of up to a factor of 2.

The use of TST has an even more profound influence on Aeff.

While the

pre-exponential factor is very roughly in the realm of the conventional 1013/sec value we can see variations in Aeff in Table 3 of up to factors of 100! It is readily apparent that the model chemical mechanism of the reaction is very different for different metals, although the reaction network topology is the same for each.

However, on

the Pd and Pt metals, we clearly see that reaction 1 has a very small forward rate constant relative to reaction 5.

This is consistent with the findings of Greeley and Mavrikakis11 for

reaction on Pt. Both Pt(111) and Ni(111) have one large rate coefficient, for step 10 on Pt and step 7 on Ni, that results from very small reaction barriers, see Fig. 4a-c.

For step 7 on Ni, the

barrier is low enough that the zero point corrected barrier height is actually negative, suggesting that the HCOH adsorbate is unstable on Ni(111).

We included reaction 7 on Ni with its TST

rate coefficient within our global sensitivity and pathway analyses without much concern, since the order of magnitude of the rate coefficient for step 7 is so much greater than all the other reactions, even within our rate coefficient sample space, it can approximately be thought of as happening “instantly".

This is consistent with the HCOH being unstable.

Lastly, it is

important to compare the rate coefficients for reactions at branches in the reaction mechanism. Since reactions 1 and 5 have very similar rate coefficients on Ni, significant competition is expected at this junction. Similarly, reactions 6 and 9 on Pd are quite close together. As temperature increases, the all the rate coefficients become closer in value, and we expect 12 ACS Paragon Plus Environment

Page 13 of 58

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


pathway competition to increase as the temperature increases. The potential energy surface along various reaction paths is depicted in Fig. 5 for the three catalytic surfaces.

It is often conjectured that reaction pathway on surfaces should simply

follow the lowest barrier route through the reaction intermediates.

By this rule of thumb, path 1

should dominate Pd(111), path 3 should dominate Pt(111), and path 1 should be the most prominent for Ni(111).

As we shall see below, the actual reaction pathways followed are

distributed in more complicated and temperature dependent manner.

III.

Global Sensitivity Analysis of Methanol Decomposition A. Theoretical Method A useful way to understand the kinetics of the reaction network is to employ the methods

of global sensitivity analysis. literature,31

32 33

Since these techniques have been discussed extensively in the

including our own previous work,1

2 3 4

we present here only a brief review.

Further details are presented in the Appendix. The analysis employs a target function of the kinetic trajectories, τ, which represents some aspect of the chemistry that is of particular interest. We solve the kinetics equations from a fixed set of initial conditions, i.e. initial coverages

θ(t=0)=(θ1,θ2,..)t=0 to obtain the trajectories for each choice of values for the rate coefficients, i.e. "# "+

= A 7B9 C = DEE FGrndm>0 to compute the time of the reaction, trxn, without the need to converge the master equation with respect to ∆t.

The sum in the denominator of eq. (4.2) is the sum of all the

transition probabilities out of species Ai at time t.

A second random number is then used to

select which species the followed atom will migrate to upon reaction.

Since the present

problem is not computationally demanding, we employed the simple Markov chain method.

B. Pathway Analysis of the Sensitivity Indices The specific utility of the pathway data that is generated from this algorithm will naturally depend on the particular scientific objective of the study at hand.

Here, we are

interested in understanding the origin of the observed sensitivity of our target function to the value of the rate coefficients.

Since the target represents the speed of passage through the full

reaction network, we are interested only in retaining information about the statisitics of the pathways taken and associated passage times.

For the followed-atom we select the carbon atom

which passes from the MeOH reagent to the final CO product through various intermediate pathways.

Because of the simplicity of the present network, we specify the pathway using an

N-string of indices (J1, J2, J3,…,JN) where Ji labels the species containing the carbon atom, and Ji can range from 1-8.

There is only a single reaction that capable of the transformation Ji→Ji+1

so the reaction index is unnecessary in the present system. The final species is always CO so JN=8.

Each unique string of indices defines a path, and the probability of the mth pathway is

simply

Pm =

nm N tot

(4.3)

where nm is the number of Markov chain trajectories corresponding the mth pathway and Ntot is the total number of stochastic trajectories.

Associated with each pathway “m” is a distributions

of arrival times of the carbon atom to the final product, Ti(m) with i=1,…,nm.

Under each set of

initial conditions, such as temperature and coverage, there is a mean value of the arrival time .

Hence, the average value of the arrival times can be computed from the pathway

probabilities as

r ≈ 8tu = ∑v \v ∙ 〈87q9〉

(4.4)

Naturally, the 80% conversion time target, τ, is not exactly equal to Tave, but we expect that it is 19 ACS Paragon Plus Environment


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

Page 20 of 58

similar enough to the target that one can use it to understand the sensitivities.

In this

expression, it is only assumed that the reaction is run to completion and all contributing pathways are included.

Thus, in this picture, the sensitivity of the target is represented in terms

of a pathway representation.

If only a few paths are important, the summation (4.4) can be

accurately truncated to a small number of terms and the sensitivities may become easy to understand.

C. Pathway results For each initial condition studied, 15000-50000 stochastic trajectories were run and the statistics were extracted for the chemical pathways and the passage times. It was found that only a handful of pathways had a significant probability for the conditions studied.

In fact, most of

the results could be understood in terms of just four paths: CH3OH→CH2OH→HCOH→COH→CO

Path 1

CH3OH→CH2OH→HCHO→HCO→CO

Path 2

CH3OH→CH2OH→HCOH→HCO→CO

Path 3

CH3OH→CH3O→HCHO→HCO→CO

Path 4

These paths are depicted graphically in Fig. 12.

It should be noted that only paths built with

forward reactions showed significant probablities and thus influenced the target.

The reverse

reactions tended to be slow, and uncommon in the pathway analysis, because of the high barriers to recombination.

They become even less important at low coverages.

The probabilities for

the 4 dominant paths are shown in Fig. 13 at T=400-700 K for Pd(111), Ni(111), and Pt(111). The sum of probabilities for these four paths is within 2% of unity in all cases.

(The one

exception is for Ni(111) at the highest temperature where the reverse of reaction 5 contributes about 3% of the path 1 probability.) at the lower temperatures. path 1.

For the case of Pd(111), it is seen that path 1 is dominant

At the higher temperatures, paths 2 and 3 become competive with

For Pt(111), for all temperatures only path 3 contributes significantly and is seen to

possess a probability of near unity. probability at all temperatures.

For Ni(111), path 1 and path 4 have roughly equal

Clearly, the pathway depends strongly on the metal substrate

and, for Pd(111), on the temperature. The mean values of the “passage time”48 along each pathway have also been numerically 20 ACS Paragon Plus Environment

Page 21 of 58

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


extracted.

These are the times for the carbon atom to first appear in the CO molecule.

As we

see in Table 4, there are very large differences in the passage times between the paths.

In

Pd(111) at 400 K, e.g., we see the dominate path 1 is about 500 times slower than the next most probable pathway.

For Pt(111) on the other hand, despite the clear dominance of path 3, all the

paths exhibit roughly the same passage time.

For Ni(111), the paths 1 and 4 were roughly

equally probable, yet at 400 K path 1 was two order of magnitue slower than path 4.

The

relative passage times narrows to about a factor of 10 at 700 K, but yet the pathway probabilities are nearly constant functions of temperature. It is very interesting to observe the arrival time (i.e. passage time) probability distribution function (PDF) obtained from the stochastic simulation, which we express in terms of the distribution of log(t).

In Fig. 14, we show the PDF for the logarithm of the arrival time

distributions for Pd(111) and Ni(111) which shows a clear bimodal distribution.

The arrival

times can be separated into components for each pathway since the pathway is recorded for each trajectory; this decompostion is also shown in Fig. 14.

For Pd(111) it is seen that the peak at

long times is due path 1 while the faster arrival time peak is a combination of paths 2 and 3, which have nearly equal mean arrival times.

For Ni(111), the long arrival time peak is again

due to path 1 while the faster trajectories arrive via path 4.

We conjecture that the arrival time

PDF may provide a diagnostic for multiple paths in more general kinetic systems. In terms of the pathway analysis, the sensitivity of a target function to a given rate coefficient may result from either the variation of the target vs. ki along the important paths, or by the variation of the pathway probabilities vs. ki.

If one has, as in the present problem, a

reaction pathway that consists of a set of n-first order processes

: →?! w →? x →?y … \p{i|H_ then the mean passage time from A to Product is given by

〈r〉~ ∑

r

where τi is mean passage time for each step.

(4.5)

If a reaction step has a single product, i.e. X→Y,

then

r = ?

(4.6)



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

Page 22 of 58

If there are multiple products formed with rate coefficients ki,j where j=1,..m, then we have

r = ∑~

(4.7)

L! ?,L

When the chemistry of a given process is dominated by a single pathway then the sensitivities may be understood by examining the behavior of the target as the rate coefficients are varied along that path.

In particular, for cases where there is a clear rate limiting step, one term

dominates in the sum and the first order sensitivity is due almost entirely to that step.

This is in

keeping with our chemical intuition of the rate limiting step controlling the overall rate of reaction.

The peaks shown in the log-arrival time PDF, Fig. 14, can be easily related to the

values of the rate coefficients for the rate limiting steps. Pd(111) lies at 10-6 sec at 500 K. which has 1/k8=1.6×10-6 sec. ×10-8 sec.

For example, the longest time peak for

The slowest step along the dominant path 1 is reaction 8

At 700 K this peak shifts to ~10-8 sec in the PDF while 1/k8=1.7

The “faster” peak for the arrival time distribution of Pd(111) occurs at 10-9 sec for

T=500 K and 5×10-10 sec at T=700 K.

The SPA reveals that this feature was due to two

pathways, 2 and 3, which showed apparently nearly equivalent arrival time distributions.

The

slowest steps along path 2 are reactions 5 and 9 which have rate coefficients that are nearly equal.

The rate limiting step along path 3 is reaction 5 again.

agreement with the observed peak of the PDF. dominate.

The 1/k5 time is again in good

For Ni(111), pathways 1 and 4 were found to

The longer time peak of the arrival time PDF corresponds to pathway 1, and the

peak position is in good agreement with the 1/k8 for the rate limiting step, reaction 8. pathway 4, the rate limiting step is reaction 2.

Along

The predicted delay time 1/k2 is found to

correlate well with the faster peak positions at both 500 K and 700 K.

Finally, the

dehydrogenation reaction on Pt(111) occurs through pathway 3 at all temperatures.

The rate

limiting step along path 3 is reaction 5 and 1/k5 is found to correspond well with the peak positions observed for the arrival time PDF. It is interesting to note that arrival time heirarchy revealed in Fig. 14 using the stochastic trajectory simulation also clearly appears in the conventional kinetic simulations generated using eq. (3.1).

In Fig. 15 the production rate, dθCO/d(Logt)∝ t×dθCO/dt, is plotted

for Pd(111), Ni(111) and Pt(111) at T=500 K and 700 K.

As one sees by comparison of Fig. 15

to Fig. 14, log-arrival time PDF is nearly identical to the conventional production rate.


Page 23 of 58

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


In the present problem, we find that the first order sensitivities of the 80% CO conversion time target can often be understood in terms of the rate limiting step(s) along the dominant pathway.

Consider the first order sensitivity indices of Pd(111) at low temperatures, a problem

that is clearly dominated by pathway 1.

Reaction 8 is found to be several orders of magitude

slower than any other reactions along this pathway, i.e. 5, 6, and 7.

In keeping with the

previous discussion, the first order sensitivity index S8 is clearly dominant.

For Pt(111), Fig. 13

shows that the reaction proceeds almost exclusively along pathway 3.

The slowest step along

this path is reaction 5 with reaction 6 being the second slowest (by a factor of 2.4 at T=500 K). Consistent with this, the S5 is found to be the largest 1st-order index while S6 is the second largest.

The case of reaction on Ni(111) shows S8 to be largest first-order sensitivity with S5

and S1 being next most sensitive with values near 0.1.

The pathway analysis has revealed that

paths 1 and 4 dominate reaction on Ni(111) and reaction 8 is rate limiting on path 1. role of reactions 1 and 5 are not so obvious.

However,

Reaction 2 is actually somewhat slower than

reaction 1 along path 4, and reaction 5 is much faster than reaction 8 along path 1.

We believe

that the significance of reactions 1 and 5 is related to a phenomenon we term “pathway switching”.

That is, the value of the rate coefficients k1 and k5 determines the pathway

probabilities (Pm in eq. 4.3)

rather than the transit time along a path.

This pathway switching

is essential to understanding the origin of correlation that unlies the 2nd order sensitivity coefficients. The 2nd order sensitivity indices, Si,j, were found to be numerically important for the Ni(111) and Pd(111) cases, but unimportant for Pt(111).

We note first that Pt(111) is

dominated by one pathway (path 3) at all temperatures, while multiple pathways contribute for Ni(111) and Pd(111).

Consider first the Ni(111).

order indices were important for reactions 1, 5, and 8. and S5,8 are numerically large.

As with the 1st order sensitivity, the 2nd Specifically, we see in Fig. 11b that S1,8

As shown in Fig. 16, varying the values of the rate coefficients

k1 and k5 has a major influence on the pathway probabilities at all temperatures.

On the other

hand, reaction 8 has little effect on the pathway probabilities over its uncertain range. to understand the effect of k1 and k5.

The reagent CH3OH can dehydrogenate to CH2OH via

reaction 5 or to CH3O via reaction 1. respectively.

It is easy

These are the first steps of pathways 1 and 4,

If k1 is increased while k5 is held fixed at its nominal model value, the probability

of proceeding along path 4 grows like k1/(k1+k5).

When k5 is increased, the probability of 23



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

taking path 1 grows like k5/(k1+k5).

Page 24 of 58

While it is obvious why reactions 1 and 5 influence the

pathway probabilities, it is less clear why a correlation develops to reaction 8 via S1,8 and S5,8. Here, the specific choice of the τ80% target plays a role.

We see in the PDF shown in Fig. 14

that the trajectories moving along path 4 reach the CO target 50-100 times faster than those moving along path 1.

Hence, the trajectories along path 4 arrive very quickly and the 80%

target is finally achieved only when the slower path 1 trajectories arrive.

The time for which

those path 1 trajectories arrive is dominantly controlled by the k8 rate limiting step.

Hence, the

correlation develops because k1 and k5 determine the pathway probabilities and k8 determines the arrival time of the last portion of the product.

Clearly, if the target was set differently, say 10%

conversion to CO, then the 2nd order sensitivity to reaction 8 would be significantly lessened since 10% CO could be easily achieved by passage along path 4 alone. The 2nd order sensitivity indices on Pd(111) were found to be large for S6,8, S7,8, S8,9, and S8,10.

We recall that for Pd(111) path 1 was most probable at low temperatures while paths 2 All the large 2nd order indices

and 3 become increasing probable at higher temperatures. correlate to reaction 8.

Once again, we believe this is a consequence of reaction 8 being the rate

limiting step along path 1 which is the slowest path (see Fig. 14, upper panel). and 10 are important for the pathway switching mechanism. these rate coefficients are shown in Fig. 17.

Reaction 6, 7, 9,

The pathway probabilities versus

We see from the figure that reactions 6 and 9

control the pathway switching between paths 1 and 2.

This effect is more pronounced at T=700

K than it is at T=500 K in keeping with the fact that pathway probabilities are more commensurate at higher temperatures.

Likewise, reactions 9 and 10 appear to control the

pathway switching between paths 1 and 3; again the effect is stronger at high temperature. Fig. 18, we show the pathway probabilities versus k8. whatsoever of the pathway probabilities on k8.

In

As expected, there is seen to be no effect

The correlation develops between reaction 8 and

reaction 6, 7, 9, and 10 due to the combined effect of pathway switching and the rate limiting nature of reaction 8 on path 1. As a final illustration of the influence of pathway analysis in the kinetic simulations, we consider the distribution of target values, τ, obtained from the global sampling of the rate coefficients over the uncertainty hypervolume.

The distribution of target values is

quantitatively described by a target probability distribution function, or tPDF, P(τ).

In previous 24


Page 25 of 58

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


studies of combustion problems, we have found that the tPDF resembled a Gaussian or skewed Gaussian function.1 instead bimodal. P(log(τ)).

3 37

Here, for the Ni(111) and Pd(111) cases, we find that the tPDF is

To aid in clarity, we represent the distribution logarithmically, i.e. we plot

As shown in Fig.19, the log-tPDF for these two case shows a large peak at large τ

and a smaller peak at smaller τ.

The bimodal form of the distribution is apparently a

consequence of the multiple contributing pathways that describe the kinetics for the Ni(111) and Pd(111) surface reactions.

Reactions on the Pt(111), for which only one pathway is important,

gives a unimodal distribution of target values.

A rough approach linking the chemical

pathways to the target distribution is to compute the target distribution where all the pathways except a single one has been “disabled”.

This is accomplished by setting the rate coefficients

to zero except along the selected pathway, and then the nonzero rate coefficients along that pathway are sampled.

This process is repeated for each contributing pathway.

For Ni(111),

we see from Fig. 19 that the slow (large τ) peak is due exclusively to the trajectories passing along path 1. 4.

The faster (small τ) peak is due exclusively to the trajectories moving along path

For Pd(111), the slow peak is due to path 1 and the fast peak is combination of the

unresolvable contributions of paths 2 and 3.

While the distributions shown in Fig. 19 bear

resemblence to the arrival time PDF of Fig. 14, we should point out that there are important differences.

The target distribution reflects the variation of the 80% CO production target at the

rate coefficients are varied.

The arrival time distribution presented earlier, in contrast, shows

the spread of arrival times for trajectories at fixed values of the rate coefficients.

V.

Conclusion In this report we have discussed the decomposition of methanol on three metal surfaces,

Pd(111), Pt(111), and Ni(111).

In the first part of this work, we constructed a quantitative

kinetic model for the surface chemistry.

A set of DFT calculations were conducted to establish

the energies and frequencies of the reagents and transition states relevant for all the important surface reactions.

Using this information, rate coefficients for the forward (dehydrogenation)

and reverse (recombination) reactions were calculated using TST.

The relaxation of surface

atoms has a significant effect on the energetics and, hence, on the rate coefficients.

Special care

was devoted to the computation of the symmetry factors which account for interchange



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

Page 26 of 58

symmetry, chirality, and reaction path degeneracy. In the second portion of the work, the sensitivity of the surface chemistry to the values of the rate coefficients was investigated using a combination of global sensitivity analysis and stochastic pathway analysis.

It was found that, within the sampled uncertainty space of the rate

coefficients (a factor of 5 for each coefficient), the model was insensitive to changes in the initial methanol surface coverage.

The decomposition of methanol proceeds as a sequence of

(forward) dehydrogenation steps leading to the CO product with little significant reverse reaction noted.

On the surfaces Pd(111) and Ni(111), first order sensitivity analysis showed that the

80% CO conversion target time was very sensitive to the final reaction step COH → CO + H, since this is the rate limiting step within the mechanism on those surfaces.

Some first order

sensitivity is also observed for the reactions at pathway branching points in the mechanism. The importance of these reactions originates from their ability to redirect reaction flux along different pathways within the mechanism, i.e. pathway switching.

For Pt(111), it was seen that

the conversion of methanol to CH2OH the clear limiting step along the pathway within the rate constant uncertainty space used. Second-order global sensitivity of the target function was found to be numerically important for reactions on Pd(111) and Ni(111) but unimportant for reaction on Pt(111).

In

both Ni(111) and Pd(111), the main second order effects were correlations between the COH → CO slow step (reaction 8) and the reactions at pathway junctions.

In Pd(111), the results of the

sensitivity analysis were sensitive to the temperature, however, the Ni(111) sensitivities were approximately independent of temperature.

In the latter case, the independence was due to

weak temperature dependence in the pathway branching probabilities k1/(k1+k5) and k5/(k1+k5). An important part the present work involves of the use of the stochastic pathway analysis (SPA) to interpret the origin of the sensitivities and correlations.

A pathway is defined by

following a specific atom through a sequence of reactions until a stopping criterion for the trajectory is met.

The pathway is labeled by the species, reactions, and times at which and to

which the chemical transitions occur.

While in principle there can be a very large number of

contributing pathways, in the present problems in which we followed the C-atom, it was found that only 4 paths were important.

The target function was chosen for simplicity to be the time

needed for 80% of the methanol reagent to be converted to the CO product.

Thus, pathways

were followed until the C-atom appeared in the CO-product for the first time.

The pathway 26


Page 27 of 58

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


probabilities were extracted to high accuracy using a straightforward dynamical Monte Carlo algorithm.

It was found that target sensitivity to the various rate coefficients originated in one

of two ways.

First, the sensitive reaction could be a rate limiting step along an important

pathway to the target. important pathways. scenario.

Or second, the sensitive reaction could control the switching between The first order sensitivities tended to arise via the rate limiting step

The results for Ni(111) and Pd(111) suggest that the COH intermediate is quite stable

and much of the reactive flux is restricted at this step in the kinetics.

On the other hand, the

interaction (or correlation) observed in the 2nd order sensitivity tended to involve pathway switching working in concert with a rate limiting step. control pathway switching.

Often, there are a pair of reactions that

For example, the relative values of the rate coefficients for

reactions 1 and 5 determine whether path 1 or path 4 is preferred on Ni(111).

In future work on

other kinetic systems, additional sources of target sensititivies may well be discovered.

Indeed,

the SPA method appears to be a very robust method to interpret the behavior of complicated reaction networks. In addition to providing a means to understand the origin of target sensitivity, the SPA provides unique insight into other reaction observables.

For example, the distribution of arrival

times at the CO product was found to show a clear imprint of the pathway structure.

In

addition, the target distribution function exhibited a bimodal distribution for reactions that possessed pathway switching.

Acknowedgements We are grateful to the National Science Foundation of the United States and the Chinese Academy of Science for support of this work. D. Zhou was supported by the Division of Chemical Sciences, Geosciences, and Biosciences, the Office of Basic Energy Sciences, the U.S. Department of Energy, under contract number DE-AC02-06CH11357.

WXL was supported by

the NSFC (21173210, 21225315) and 973 Program (2013CB834603).

Supporting Information Available The detailed results of the DFT calculations are available as supporting information.

This 27



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

Page 28 of 58

information is available free of charge via the Internet at http://pubs.acs.org.

Appendix This section gives an overview of the global sensitivity analysis method, based on the work by Sobol49 and the HDMR method of Rabitz34

35 36 37

.

In our previous work, we have

used this methodology as a screening method to select rate coefficients in large combustion mechanisms that were most in need of “improvement” through high level computation.

In this

method, each rate coefficient parameter has an associated uncertainty, ∆ki , which may represent actual error bars or may be introduced artificially for purposes of analysis.

The goal of global

sensitivity analysis is to determine how the uncertainty in all the model parameters (i.e. all the rate coefficients) translates to the uncertainty in the model result as measured using the variance. Therefore, using Monte Carlo sampling, a large number of random sets of model parameters are selected from the uncertainty hypervolume of rate coefficients, {k1n , k 2n , L} , where n stands for the n th sample. to completion. obtain τn.

The kinetic trajectory for each member of the statistical sample is propagated The value of the target observable is evaluated for that set of rate coefficients to

The initial conditions of the kinetic model should be taken to be the same for each

member of the random sample, i.e the temperature and initial methanol coverage, so only the rate coefficient are being varied.

This generates a set of values of the target that define a probability

distribution of the target variable that exhibits a total variance σ2.

The main objective of the

method is to assign components of the total variance to specific rate coefficients (factors) or to the correlations between factors. We have previously described the advantages of global sensitivity analysis over the more familiar local sensitivity analysis in kinetics. Specifically, while global sensitivity analysis covers the full parameter uncertainty range, local sensitivity analysis only represents the response of the target at the nominal (i.e. zero uncertainty) value k0, e.g. through its local deriviatives, ∂τ/∂ki.

Since the target may be a highly nonlinear function of k, and the “true”

value of k may be far from the nominal model value, then local sensitivity analysis can prove very deceptive. According to the HDMR representation we can expand the target function as a function


Page 29 of 58

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 each of its parameters,

(

)

nk

nk nk

i =1

i > j j =1

τ k1 , k 2 , L , k nk = τ 0 + ∑Fi (k i ) + ∑∑Gij (k i , k j )

(A1)

where Fi (k i ) are functions that represent the dependence of the target on a single rate coefficient, and G ij (k i , k j ) the are functions that represent dependence with respect to the coupling of two rate coefficients. Higher order terms can be used, but it has been observed that many physical systems can accurately be represented by the truncation of Equation 8 at the second order term and so we shall not consider expansion beyond the 2 nd order. The term τ 0 is a constant term. The F and G component functions are expressed in terms of an orthogonal polynomial basis, so that they satisfy the relations, 〈 Fi 〉 = 0, 〈 Fi F j 〉 = δ ij 〈G n ,m 〉 = 〈 Fi G n ,m 〉 = 0

(A2)

〈G n ,m G j ,k 〉 = δ nj δ mk

M where 〈L〉 denotes an average (volume integral) over all the uncertainty hypervolume.

In this

study, orthonormal Legendre polynomials were used, L0 ( x ) = 1 L1 (x ) = 3 (2 x − 1) L2 ( x ) = 5 (6 x 2 − 6 x + 1)

(A3)

M The first and second order terms in eq. A1 are expressed as, ηi

Fi (k i ) = ∑α li Ll (κ i )

(A4)

l =1

ηi η j

Gij (ki , k j ) = ∑∑βlijδ Ll (κi )Lδ (κ j )

(A5)

l=1 δ =1

where κ i represents a scaled variable needed since the Legendre polynomials are defined in the



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

Page 30 of 58

space x ∈ [0,1] . The κ i factor is defined as,

κi =

k i − k iMIN k iMAX − k iMIN

(A6)

The terms kiMIN and k iMAX are the minimum and maximum values considered in the uncertainty sample space of rate coefficient ki . The sampling space of the rate coefficients is defined by a probability distribution function Pk (k ) , where k stands for the set of the rate coefficients in the mechanism. In the orthogonal polynomial representation, the functions in eq. (A1) can be determined knowing

Pk (k ) , by using the hyper-volume integrals

τ0 = ∫

Ω FULL

Pk (k )τ (k )dk1 dk 2 L dk N

(A7)

Fi (k i ) = ∫ Pk (k )τ (k )dk1dk 2 L dk i −1 dk i +1 L dk N − τ 0

(A8)

Ω −i

Gij (k i , k j ) =

∫

Ω −i − j

Pk (k )τ (k )dk1 dk 2 L dk j −1 dk j +1 L dk i −1 dk i +1 L dk N

− Fi (k i ) − F j (k j ) − τ 0

(A9)

where Ω is the full rate coefficient uncertainty space, and Ω −i signifies the reduced rate constant uncertainty space at a specific set value of ki . It follows that the partial variances can be defined by the low dimensional integrals, ηi

Pi ( k i ) Fi 2 (k i )dk i = ∑α li 2

σ i2 = ∫

All ki Space

σ =∫ 2 ij

∫

(A10)

l =1

Pi (k i ) Pj (k j )G

All k j All ki Space

2 i, j

(k , k )dk dk i

j

i

ηi η j

j

= ∑∑β lijδ 2 l

(A11)

δ

where we have made use of the orthonormality of the Legendre polynomials. The accuracy and convergence of the expansion can be judged based on the ability of eq. (A1) to reproduce the observed probability density function (PDF) of the target function when sampled over the rate coefficient uncertainty space. While it is possible to obtain the expressions for the polynomial expansion of the target using the integration formuale above, even for mechanisms of moderate size the multidimensional integrations can become computationally expensive or prohibitive.

The 30


Page 31 of 58

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


scheme we use to obtain the coefficients α and β in eqs. (A4) and (A5) is based on linear regression.1

3 35

The uncertainty space of k is randomly sampled according to its probability

distribution Pk (k ) a number N times.

This can be done using Monte Carlo sampling and

the random sampling HDMR (RS-HDMR) technique of Rabitz el al.35

The model is run for

each of the N set of samples of k , N target values are generated, and the total variance is determined according to, 2 σ TOT = 〈τ 2 〉 − 〈τ 〉 2

(A12)

Since there are now N samples of the target over k-space, the computation of the means in eq. (A12) is quite straightforward as a sum over a discrete set of data.

The coefficients in the

component functions are determined by a linear least squares fitting.

In this work, the

coefficients α for each individual first order function of the form eq. (A4)

were fit to the

residuals of the sampled target values, τ i Ri1 = τ i − τ 0

where τ0=.

(A13)

The second order coefficients, β , of eq. (A5) were determined by linear

regression of the residuals of the sampled target values from the mean and the first order fits, nk

R = τ i − τ 0 − ∑Fi (ki ) 2 i

(A14)

i =1

where the Fi(ki) functions have been already determined.

A more accurate method, which is

also more computationally demanding, is to fit all the coefficients simultaneously.

The first and

second order sensitivities can then be determined from,

σ i2 2 σ TOT

(A15)

σ ij2 S ij = 2 σ TOT

(A16)

Si =

using eqs. (A10) and (A11).

An example of the fitted component functions has been presented

in Fig. 8. The probability distribution function for each of the 10 forward rate coefficients was assumed to be uniform, pk (ki ) = i

k

MAX i

1 − k iMIN

(A17) 31



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

Page 32 of 58

We adopt a definition of the maximum and minimum values of ki according to, k iMAX = kiTST (| ∆logk i |)

(A18)

kiMIN = kiTST (| ∆logki |)

(A19)

−1

The ranges chosen here go from a factor of 5 larger to a factor of 5 smaller than the nominal value. Although the true uncertainty may be larger, this value was deemed sufficient as a preliminary test of the global sensitivity of the given mechanism.

The reverse reactions are

related to the forward rate coefficients through the equilibrium constant, kib =

ki f K eq

(27)

and throughout the sensitivity analysis process, the nominal values of the equilibrium constant are held constant.


Page 33 of 58

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


Table 1: Adsorbate absolute binding energies and reaction barrier heights obtained by DFT calculation. Absolute Value of Adsorbate Binding Energies (kcal mol-1) Adsobate Pd(111) Pt(111) Ni(111) CH3OH 5.8 4.6 5.8 CH3O 45.9 36.7 61.3 HCHO 14.8 10.1 18.9 HCO 57.9 52.3 54.2 CH2OH 43.3 48.0 38.0 HCOH 75.2 78.2 68.5 COH 107.9 108.6 103.8 CO 48.9 41.3 44.5 † H 15.2 11.1 13.4 Forward and Reverse Adsorbate Dehydrogenation Barriers (kcal mol-1) Reaction Pd(111) Pt(111) Ni(111) 24.7 24.2 18.7 6.0 18.0 31.1 CH3OH➙CH3O+H 13.8 23.3 4.6 14.5 19.6 15.9 CH3O➙HCHO+H 8.3 27.0 12.7 30.2 7.8 16.8 HCHO➙HCO+H 6.5 38.7 6.5 28.4 4.4 34.1 HCO➙HCO+H 12.7 18.2 14.3 21.4 17.1 15.4 CH3OH➙CH2OH+H 13.4 24.2 13.6 18.7 9.5 17.1 CH2OH➙HCOH+H 8.8 27.2 13.1 25.1 1.4 20.5 HCOH➙COH+H 19.1 44.3 21.0 33.4 22.1 45.2 COH➙CO+H 15.9 19.4 18.7 8.8 18.7 25.4 CH2OH➙HCHO+H 15.2 26.5 3.5 9.0 3.5 27.4 HCOH➙COH+H * All reported energies are the absolute energies and do not contain the zero point correction. Values are reported to roughly to a precision of 0.3 kcal mol-1. † Energy is with respect to H2 molecule and so corresponds to the reverse reaction of the dissociative adsorption of H2. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.



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

Page 34 of 58

Table 2. Symmetry prefactors to the transition state theory rate coefficients shown in Equations 2.5 and 2.7. Parentheses are put in when Pt(111) factors differ from those of Pd(111) and Ni(111). The ‘*’ indicates a free rotor partition function is used for in-site rotation of one of the reactants. RXN

Forward

Reverse

1*(*)

2

2

2*(*)

6

2/3 (1/3)

3

2 (1)

2/3 (2)

4(*)

2 (6)

2

5*(*)

3

1/3

6

1

2/3

7

1

1/3

8

1

1

9

1

2/3 (1/3)

10

2

2/3


Page 35 of 58

Table 3. The fit parameters of the forward and reverse TST rate coefficients to the form 79 = b 7− 9. The forward rate coefficients are expressed in units of s-1 and the reverse rate

Pt(111)

Pd(111)

coefficients in units of cm2 s-1. The ‘Ea’ parameter is expressed in kcal mol-1 to facilitate comparison with the reported reaction barrier heights and T is in K. The parameter ‘Aeff’ is the pre-exponential term obtained by fitting the rate coefficient to a standard two parameter Arrhenius expression. The units of Aeff and A are chosen so that rate coefficients have the proper cgs units.

Ni(111)

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


Rxn 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10

Aeff 6.05×1012 5.56×1013 4.85×1013 2.37×1013 1.25×1012 5.04×1012 8.47×1012 5.08×1012 4.57×1012 1.66×1014 1.36×1012 9.55×1012 2.11×1013 1.13×1013 3.61×1012 3.04×1012 2.38×1013 7.89×1012 6.46×1011 5.59×1013 3.59×1012 5.49×1013 3.31×1013 1.10×1013 2.75×1012 2.13×1013 3.52×1013 1.15×1013 6.85×1012 2.69×1013

Forward Reaction A δ 8 1.93×10 1.424 9 2.27×10 1.393 1.35×1011 0.810 2.36×1013 0.320 2.80×107 1.472 4.37×1010 0.656 1.51×1011 0.554 4.21×1010 0.657 3.71×109 0.980 2.96×1011 0.874 1.27×108 1.270 5.52×108 1.341 2.36×1011 0.621 3.26×109 1.120 2.95×108 1.295 1.55×107 0.726 2.02×1011 0.667 6.82×1010 0.650 2.43×109 0.764 3.93×1011 0.681 4.46×107 1.553 2.91×109 1.360 3.06×1011 0.646 2.83×1012 0.188 7.53×107 1.447 1.59×1011 0.676 3.34×1010 0.955 5.80×1010 0.728 6.38×109 0.960 2.75×1010 0.950

Ea 19.3 9.4 4.4 4.8 7.7 9.8 6.3 15.1 10.6 10.7 13.9 0.9 9.6 3.4 10.1 9.6 10.3 17.0 14.1 -0.3 12.0 15.3 4.6 3.3 12.1 5.8 -1.9 17.8 9.0 10.1

Aeff 1.33 0.0388 0.0566 0.271 0.0149 0.501 0.00323 0.0229 0.0579 0.249 0.0113 4.71×10-3 0.0140 0.0760 9.68×10-3 0.0612 5.92×10-3 0.0107 9.87×10-4 9.37×10-4 0.604 0.0209 9.67×10-3 0.0174 0.0260 0.0305 0.0513 0.0171 0.0463 0.037

Reverse Reaction A δ 0.0105 0.679 0.0163 0.111 0.0153 0.174 0.3818 -0.055 0.0143 0.014 0.0169 0.138 0.0118 -0.180 0.0197 0.330 0.01465 0.202 0.1314 0.100 -4 7.92×10 0.369 0.0257 -0.229 5.97×10-4 0.431 0.5677 0.275 0.0125 -0.040 0.2460 1.881 0.0158 -0.139 9.06×10-3 0.0226 3.17×10-3 -0.159 3.34×10-5 0.4572 4.06×10-5 1.307 0.0353 0.061 0.0197 -0.089 0.0319 -0.090 0.0227 0.029 0.0388 -0.021 0.0520 0.013 6.30×10-4 0.760 3.19×10-4 0.672 3.51×10-4 0.629

Ea 23.5 23.3 26.7 39.5 18.1 24.5 28.1 43.3 18.6 26.2 6.7 15.6 30.7 29.3 21.8 19.6 26.0 33.3 9.5 9.2 28.5 16.3 17.4 35.1 15.8 17.4 20.7 41.1 23.5 25.7



Pt(111)

Pd(111)

Table 4. The average methanol to carbon monoxide passage time along each pathway in the mechanism observed in the stochastic kinetics trajectories.

Ni(111)

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

Page 36 of 58

Temp. (K) 400 500 600 700 400 500 600 700 400 500 600 700

Average CO Conversion Time (s) Path 2 Path 3

Path 1 8.76 × 10-5 1.65 × 10-6 1.18 × 10-7 1.80 × 10-8

1.05 × 10-9 * 1.23 × 10-3 1.19 × 10-5 5.20 × 10-7 5.70 × 10-8

1.69 × 10-7 1.47 × 10-8 2.82 × 10-9 8.44 × 10-10 7.78 × 10-7 * 3.78 × 10-8 5.46 × 10-9 1.65 × 10-9 5.47 × 10-6 1.25 × 10-7 1.35 × 10-8 2.50 × 10-9

1.61 × 10-7 1.41 × 10-8 2.81 × 10-9 8.61 × 10-10 6.25 × 10-7 3.98 × 10-8 6.06 × 10-9 1.55 × 10-9

3.01 × 10-11 *

Path 4

2.94 × 10-9 * 1.94 × 10-9 * 6.74 × 10-10 5.41 × 10-7 3.15 × 10-8 4.15 × 10-9 1.04 × 10-9 2.55 × 10-5 4.93 × 10-7 3.59 × 10-8 5.46 × 10-9

The numbers with a ‘*’ correspond to rare event pathways (occur in less than 0.01% of the 15000 trajectory ensemble), and so the average conversion times may have significant error.


Page 37 of 58

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


Figure 1: The reaction network for the decomposition of methanol that occurs through 10 reversible dehydrogenation reactions. The H-atoms are eliminated in each successive forward reaction step.



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

Page 38 of 58

Figure 2: The reactant adsorbates in the most stable adsorption sites on Pd(111) and Ni(111) are shown on the first row. Some of the equilibrium structures of Pt(111) were visibly different from those on Pd(111) and Ni(111), and these structures are shown in the second row.

Figure 3: The optimized transition states configurations for each reaction on Pd(111). The structures for the other metals are quite similar.


Page 39 of 58

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


Figure 4: The classical barrier energies (no zero point correction) for the dehydrogenation reactions are shown superimposed on the reaction network for Pd(111), Pt(111), and Ni(111). The forward (green) and reverse (red) activation energies calculated at a temperature of 500 K are shown for Pd(111), Pt(111), and Ni(111), respectively. All energies shown are in kcal mol −1 .



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

Page 40 of 58

Figure 5. The potential energy surfaces of methanol dehydrogenation on Pd(111), Pt(111), and red line Ni(111). The blue line CH3OH→CH3O→HCHO→HCO→CO; CH3OH→CH2OH→HCOH→COH→CO; purple dash line CH3OH→CH2OH→ HCHO→HCO→CO; green dot line CH3OH→CH2OH→HCOH→HCO→CO. These correspond to paths 4, 1, 2, and 3 in Fig. 12 below.


Page 41 of 58

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


Figure 6. The rate coefficients for the 10 forward reactions are depicted in Arrhenius plots for the three catalytic surfaces. It is seen that reaction 8 is the slowest or next slowest step on all surfaces.

Figure 7. Schematic of the global sensitivity analysis process. The parameter sample space is randomly sampled. A probability density function is defined, and therefore a variance can be obtained. The variance is then decomposed into contributions from each parameter and each pair of parameters, when taken to 2nd order.



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

Page 42 of 58

Figure 8. An illustration of the regression fitting of the component functions for the most sensitive reactions. In panel (a) we show the function F8(k8) for the Pd(111) reaction along with a subset of the 50,000 random points used in the fitting. Panel (b) shows the function F5(k5) for the Pt(111) surface.


Page 43 of 58

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


Figure 9. The first order sensitivity coefficients of the 80% conversion time target computed for the three surfaces. The sensitivities are seen to be very different for the three surfaces.



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

Page 44 of 58

Figure 10. The sum of second-order sensitivity indices, ∑ UJ K ,J for all pairs of reactions versus temperature.


Page 45 of 58

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


Figure 11. The contribution of specific pairs of correlated reactions that most strongly affect the second order sensitivity. The Ni(111) has only two important correlated pairs, (1,8) and (5,8) while Pd(111) has four strongly interacting pairs, (6,8), (7,8), (8,9) and (8,10).

Figure 12. The four most important pathways followed by a carbon atom through the reaction network.



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

Page 46 of 58

Figure 13. The probabilities for the four most important pathways computed at temperatures 400, 500, 600, and 700 K.


Page 47 of 58

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


Figure 14. The probability distribution function (PDF) for log(arrival time) obtained from stochastic simulation of the 80% conversion time target. The black curve is the total PDF for all trajectories. The colored curves are the components of the total PDF due to the specific pathways indicated



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

Page 48 of 58


Page 49 of 58

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


Figure 15. The production rate of CO computed using conventional simulation of concentrations at two temperatures .



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

Page 50 of 58

Figure 16. The influence of k1 and k5 on the pathway probabilities for reaction on Ni(111). Here, kmin1 and kmax1 indicate the minimum and maximum values for k1 within its uncertainty range; kmin5 and kmax5 are defined similarly. The pathway probabilities are computed holding all the rate coefficients fixed at the nominal model values except for either k1 or k5.


Page 51 of 58

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


Figure 17. The influence of k6, k7 , k9, and k10 on the pathway probabilities for reaction on Pd(111). Here, kmin6 and kmax6 indicate the minimum and maximum values for k6 within its uncertainty range. The notation is similarly defined for the other reactions. The pathway probabilities are computed holding all the rate coefficients fixed at the nominal model values except the coefficient that is being explicitly varied.



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

Page 52 of 58

Figure 18 The influence of k8 on the pathway probabilities for reaction on Pd(111). Here, kmin8 and kmax8 indicate the minimum and maximum values for k8 within its uncertainty range. The pathway probabilities are computed holding all the rate coefficients fixed at the nominal model values except for the coefficient that is being explicitly varied.


Page 53 of 58

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


Figure 19. The log(target) probability distribution function (tPDF) for Pd(111) and Ni(111) at two temperatures. The normalized distribution of target values for the full simulation is shown in black. The normalized distributions obtained using single pathways are shown as colored curves. Since the single pathway results are approximate (see text) the peak positions do not exactly align.



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

Page 54 of 58

1 Skodje, R. T.; Tomlin, A. S.; Klippenstein, S. J.; Harding, L. B.; Davis, M. J. Theoretical Validation of Chemical Kinetic Mechanisms: Combustion of Methanol. J. Phys. Chem. A, 2010, 114, 8286-8301. 2 Klippenstein, S.J.; Harding, L.B.; Davis, M.J.; Tomlin, A.S.; Skodje, R.T. Uncertainty Driven Theoretical Kinetics Studies for CH3OH Ignition: OH2+CH3OH and O2+CH3OH. P. Comb. Inst. 2011, 33, 351-357. 3 Davis, M.J.; Skodje, R.T.; Tomlin, A.S. Global Sensitivity Analysis of Chemical-Kinetic Reaction Mechanisms: Construction and Deconstruction of the Probability Density Function. J. Phys. Chem. A 2011, 115, 1556-1578. 4 Zhou, D. D. Y.; Davis, M. J.; Skodje, R. T. Multitarget Global Sensitivity of n-Butanol Combustion. J. Phys. Chem. A 2013, 117, 3569-3584. 5 Truhlar, D. G.; Hase, W. L.; Hynes, J. T. Current Status of Transition State Theory. J. Phys. Chem. 1983, 87, 2664-2682. 6 Holladay, J. D.;Wang, Y.; Jones, E. Review of Developments in Portable Hydrogen Production Using Microreactor Technology. Chem. Rev. 2004, 104, 4767-4790. 7 Burstein, G. T.; Barnett, C. J.; Kucernak, A. R.; Williams, K. R. Aspects of the Anodic Oxidation of Methanol. Catal. Today, 1997, 38, 425-437. 8 Gu, X. K.; Li, W. X. First-Principles Study on Origin of the Different Selectivities for Methanol Steam Reforming on Cu(111) and Pd(111). J. Phys. Chem. C 2010, 114, 21539-21547. 9 Trimm, D. L.; Önsan, Z. I. Onboard Fuel Conversion for Hydrogen-Fuel-Cell-Driven Vehicles. Cat. Rev. – Sci. Eng. 2001, 43, 31-84. 10 Cortright, R. D.; Davda, R. R.; Dumesic, J. A. Hydrogen from Catalytic Reforming of Biomass-Derived Hydrocarbons in Liquid Water. Nature 2002, 418, 964-967. 11 Greeley J.; Mavrikakis, M. Competitive Paths for Methanol Decomposition on Pt(111). J. Am. Chem. Soc. 2004, 126, 3910-3919. 12 Hoffmann, J.; Schauermann, S.; Johànek, V.; Hartmann, J.; Libuda, J. The Kinetics of Methanol Oxidation on a Supported Pd Model Catalyst: Molecular Beam and TR-IRAS Experiments. J. Catal. 2003, 213, 176-190. 13 Greeley, J.; Mavrikakis, M. A First-Principles Study of Methanol Decomposition on Pt(111). J. Am. Chem. Soc. 2002, 124, 7193-7201. 14 Yudanov, I. V.; Matveev, A. V.;. Neyman, K. M.; Rösch, N. How the C-O Bond Breaks during Methanol Decomposition on Nanocrystallites of Palladium Catalysts. J. Am. Chem. Soc. 2008, 130, 9342-9352. 54 ACS Paragon Plus Environment

Page 55 of 58

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


15 Desai, S. K.; Neurock, M.; Kourtakis, K. A Periodic Density Functional Theory Study of the Dehydrogenation of Methanol over Pt(111). J. Phys. Chem. B 2002, 106, 2559-2568. 16 Jiang, R.; Guo, W. Y.; Li, M.; Fu, D. L.; Shan, H. H. Density Functional Investigation of Methanol Dehydrogenation on Pd(111). J. Phys. Chem. C 2009, 113, 4188-4197. 17 Sexton, B. A. "Methanol Decomposition on Platinum (111)", Surf. Sci. 1981, 102, 271-281. 18 Davis, J. L.; Barteau, M. A. , "Spectroscopic identification of alkoxide, aldehyde, and acyl intermediates in alcohol decomposition on Pd(111)", Surf. Sci. 1990, 235, 235-248. 19 Kaichev, V. V.; Miller, A. V.; Prosvirin, I. P.; Bukhtiyarov, V. I. ," In situ XPS and MS study of methanol decomposition and oxidation on Pd(111) under millibar pressure range", Surf. Sci. 2012, 606, 420-425. 20 Burke, K. Perspective on Density Functional Theory. J. Chem. Phys. 2012, 136, 150901. 21 Campbell, C. T.; Sun, Y. K.; Weinberg, W. H. Trends in Preexponential Factors and Activation Energies in Dehydrogenation and Dissociation of Absorbed Species. Chem. Phys. Lett. 1991, 179, 53-57. 22 Tautermann, C. S.; Sturdy, Y. K.; Clary, D. C. Reaction Rates of All Hydrogenation Steps in Ammonia Synthesis over a Ru(0001) Surface. J . Catal. 2006, 244, 199-207. 23 Pollak, E,; Pechukas, P. Symmetry Numbers, not Statistical Factors, should be used in Absolute Rate Theory and in Broensted Relations. J. Am. Chem. Soc. 1978, 100, 2984-2991. 24 Kresse G.; Furthmüller J. Efficiency of Ab Initio Total Energy Calculations for Metals and Semiconductors using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15-50. 25 Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total Energy Calculations using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter 1996, 54, 11169-11186. 26 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865-3868. 27 Perdew, J. P.; Wang, Y. Accurate and Simple Analytic Representation of the Electron-Gas Correlation Energy. Phys. Rev. B: Condens. Matter 1992, 45, 13244-13249. 28 Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrators. Phys. Rev. B: Condens. Matter 1976, 13, 5188-5192. 29 Henkelman, G.; Uberuaga, B. P.; Jonsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901-9904. 55 ACS Paragon Plus Environment


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

Page 56 of 58

30 Fernandez-Ramos, A.; Ellingson, B. A.; Neana-Paneda, R.; Marques, J. M. C.; Truhlar, D. G. Symmetry Numbers and Chemical Reaction Rates. Theor. Chem. Acc. 2007, 118, 813-826. 31 Saltelli, A.; Ratto, M.; Tarantola, S.; Campolongo, F. Update 1 of: Sensitivity Analysis for Chemical Models. Chem. Rev. 2012, 112, PR1-PR21. 32 Zador, J.; Zsely, I. G.; Zador, I.; Turanyi, T. Local and Global Uncertainty Analysis of Complex Chemical Kinetic Systems. Reliab. Eng. Syst. Safe. 2006, 91, 1232-1240. 33 Saltelli, A.; Chan, K.; Scott, E. M. Sensitivity Analysis; John Wiley & Sons: Chichester, West Sussex, England, 2000. 34 Li, G. Y.; Rosenthal, C.; Rabitz, H. High Dimensional Model Representations. J. Phys. Chem. A 2001, 105, 7765-7777. 35 Li, G. Y.; Wang, S. W.; Rabitz, H. Practical Approaches to Construct RS-HDMR Component Functions. J. Phys. Chem. A 2002, 106, 8721-8733. 36 Ziehn, T.; Tomlin, A. S. GUI-HDMR – A Software Tool for Global Sensitivity Analysis of Complex Models. Environ. Modell. Softw. 2009, 24, 775-785. 37 Ziehn, T.; Hughes, K. J.; Griffiths, J. F.; Porter, R.; Tomlin, A. S. A Global Sensitivity Study of Cyclohexane Oxidation under Low Temperature Fuel-Rich Conditions using HDMR Methods. Combust. Theor. Model. 2009, 13, 589-605. 38 Laidler, K. J. Chemical Kinetics; Harper Collins: New York, 1987. 39 Campbell, C. T. Finding the Rate-Determining Step in a Mechanism. J. Catal. 2001, 204, 520-524. 40 Radhakrishnan, K.; Hindmarsh, A. C. Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations. LLNL Report UCRL-ID-113855 1993. 41 McQuarrie, D. Stochastic Approach to Chemical Kinetics. J. Appl. Prob. 1967, 4, 413-478. 42 Arkin, A.; Ross, J. Statistical Construction of Chemical Reaction Mechanisms from Measured Time-Series. J. Phys. Chem. 1995, 99, 970-979. 43 Cao, Y.; Gillespie, D. T.; Petzold, L. Multiscale Stochastic Simulation Algorithm with Stochastic Partial Equilibrium Assumption for Chemically Reacting Systems. J. Comput. Phys. 2005, 206, 395-411. 44 Gillespie, D. T. Stochastic Simulation of Chemical Kinetics. Annu. Rev. Phys. Chem. 2007, 58, 35-55.


Page 57 of 58

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


45 Mevel, R.; Javoy, S.; Coudoro, K.; Dupre, G.; Paillard, C.-E. Assessment of H2-CH4-Air Mixtures Oxidation Kinetic Models used in Combustions. Int. J. Hydrogen Energy 2012, 37, 698-714. 46 Dumesic, J. A. Analyses of Reaction Schemes using De Donder Relations. J. Catal., 1999, 185, 496-505. 47 Gillespie, D. T.; Hellander, A.; Petzold, L. R. Perspective: Stochastic Algorithms for Chemical Kinetics. J. Chem. Phys. 2013, 138, 170901. 48 Van Kampen, N. G. Stochastic Processes in Physics and Chemistry; North Holland Personal Library, Elsevier: Oxford, UK, 1992. 49 Sobol, I. M. Global Sensitivity Indices for Nonlinear Mathematical Models and their Monte Carlo Estimates. Math. Comput. Simulat. 2001, 55, 271-280.

TOC graphic 57 ACS Paragon Plus Environment


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

Page 58 of 58


Following Molecules through Reactive Networks: Surface Catalyzed

Recommend Documents