- Email: [email protected]
Error propagation in first-principles kinetic Monte Carlo simulation
Accepted Manuscript Error propagation in first-principles kinetic Monte Carlo simulation Sandra Döpking, Sebastian Matera PII: DOI: Reference:
S0009-2614(17)30161-6 http://dx.doi.org/10.1016/j.cplett.2017.02.043 CPLETT 34556
To appear in:
Chemical Physics Letters
Please cite this article as: S. Döpking, S. Matera, Error propagation in first-principles kinetic Monte Carlo simulation, Chemical Physics Letters (2017), doi: http://dx.doi.org/10.1016/j.cplett.2017.02.043
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Error propagation in first-principles kinetic Monte Carlo simulation Sandra D¨ opking, Sebastian Matera∗ Institute for Mathematics, Freie Universit¨ at Berlin, Arnimallee 6, 14195 Berlin, Germany
Abstract First-principles kinetic Monte Carlo models allow for the modeling of catalytic surfaces with predictive quality. This comes at the price of non-negligible errors induced by the underlying approximate density functional calculation. On the example of CO oxidation on RuO2 (110), we demonstrate a novel, efficient approach to global sensitivity analysis, with which we address the error propagation in these multiscale models. We find, that we can still derive the most important atomistic factors for reactivity, albeit the errors in the simulation results are sizable. The presented approach might also be applied in the hierarchical model construction or computational catalyst screening. Keywords: stochastic simulation, sensitivity and uncertainty analysis, Quasi Monte Carlo, heterogeneous catalysis, first-principles kinetic Monte Carlo 2010 MSC: 65C05, 65C40, 82C99, 62-07 First-principles kinetic Monte Carlo (1p-kMC) have become popular for modeling the chemical kinetics in heterogeneous catalysis[1]. This popularity results from the predictive quality of the approach, enabling an understanding of the detailed surface chemistry without the possible masking by the fitting of parameters. In 1p-kMC, the mechanism and corresponding energy barriers are instead obtained from first-principles electronic structure methods and kinetic Monte ∗ Corresponding
author Email addresses: [email protected] (Sandra D¨ opking), [email protected] (Sebastian Matera)
Preprint submitted to Chemical Physics Letters
February 15, 2017
Carlo simulation then provides unbiased estimates of the solution of the resulting master equation. Albeit being predictive, practically applicable electronic structure methods, such as Density Functional Theory (DFT), provide only approximate results. Commonly reported uncertainties for DFT barriers are in the range of 0.2 − 0.3 eV. At low to medium temperatures, the corresponding rate constants carry then an uncertainty of some orders of magnitude. Hitherto, the propagation of these errors to the 1p-kMC simulations results has not been systematically addressed, only in the context of the much cheaper but also less accurate meanfield microkinetic modeling[2]. With the focus on the former, we employ Quasi Monte Carlo (QMC) based global uncertainty quantification (GUQ) to estimate the errors in the final 1p-kMC results. We further demonstrate a novel approach to global sensitivity analysis (GSA), which allows to quantify the error induced by each uncertain energy barrier using the same set of QMC points. This allows a global view on the main driving forces in the mechanism and generally provides a route to an error-aware analysis of such multiscale models. We exemplify this approach on the CO oxidation on the RuO2 (110) surface using the 1p-kMC model by Reuter[3], which is one of the most thoroughly investigated 1p-kMC models[4–9]. It is therefore a good test example, also because limitations of this model have been discussed in the literature[10–12]. For this model, we indeed find large possible errors, but also that many qualitative mechanistic conclusions seem to be not much affected by the DFT inaccuracies. Uncertainty Analysis. For an uncertainty analysis, the input parameters xi , the barriers or the rate constants in our case, must be considered as random variables, with a given joint probability density function (PDF). The output of a DFT code is deterministic, so this PDF must be interpreted in a information theoretical setting and reflects our knowledge about these parameters. We consider the CO oxidation on RuO2 (110) for CO and O2 partial pressures of 1 bar and a temperature of 600 K, in order to make contact with existing theoretical studies[4–6, 9], which investigate the behavior around these reaction
2
1
1
0.5
0
0
coverage
log10 TOF (s-1cell-1)
2
average nominal -1 TOF
CObr
COcus
Obr
Ocus
ebr
ecus
Figure 1: Turnover-frequencies (TOF) and coverages for CO occupied, O occupied and and empty bridge and cus sites for the CO oxidation on RuO2 . Grey bars represent the parameter averaged values and the error bars represent the standard deviation resulting from parameter uncertainty. For comparison, the values at the nominal parameter settings are displayed as green bars.
conditions. For these conditions, we also expect a rapid change of the reactivity when varying the energetic parameters within the DFT error; something which is particularly challenging for GUQ and GSA in high dimensions and therefore a good test for the methodology. We provide additional results for a CO partial pressure of one mbar in the supplement. The model considers two types of adsorption sites, bridge (br) and coordinatively unsaturated (cus) sites, which are arranged on a square lattice as alternating columns. The mechanism includes unimolecular adsorptions of CO on both types of sites, dissociative oxygen adsorption on nearest neighboring sites and the corresponding reverse elementary steps. Further, both adsorbates can make diffusive hops to neighboring empty sites and adsorbed CO and O on neighboring sites can react to gaseous CO2 . Altogether, this makes 22 uncertain parameters, for which rate constants (RC) and energy barriers are summarized in the supplementary material[13]. For simplicity, we assume that errors appear uncorrelated and the logarithm (to the base 10) of every RC is uniformly dis-
3
tributed between ±2 around its nominal value. This corresponds to uncertain DFT barriers with a uniformly distributed error in the bound [−0.25, 0.25] eV. The only exception from this rule are the adsorption steps, which have no barriers in the nominal setting. For these, we assume the nominal value to be the upper bound and allow variation only to smaller values. The lower bounds for these steps correspond then to a by a factor hundred smaller RCs. By the linear dependence of these RCs on the partial pressures, we therefore also address the dependence on these (at least to some extend). In order to access the impact of parameter uncertainty, we need to sample the parameter space and perform one kMC simulation for each sampling point. The results for these points are then used to calculate parameter-averaged key quantities, like the mean or the variance. The later then serves as a measure of the uncertainty of the simulation results[14]. We employ uniformly distributed Quasi Monte Carlo (QMC) points for this sampling of parameter space. This choice was motivated by the more even distribution of the points and a faster convergence rate in medium high dimensions compared to classical Monte Carlo sampling. In detail, we used 218 ≈260, 000 points from Sobol’s sequence for QMC[15, 16] and each kMC simulation was performed for a lattice of 20 × 20 units cells using the kmos kMC package[17]. Focusing on stationary operation, we employed 107 kMC steps for initial relaxation to reach the steady state and further 107 steps for sampling. This amounts to a total cost of ∼1000 CPU-hrs for the whole data set. The results from this parameter averaging are depicted in Figure 1. Averages are represented as gray bars and the error bars indicate the corresponding standard deviation (STD). For comparison, we have added the nominal values as green bars (i.e. those values which are obtained with the DFT barriers from Reuter[3]). The most important output from a 1p-kMC simulation is the turn-over frequency (TOF), i.e. the number of CO to CO2 reactions per second and surface unit cell. As for the RCs, we consider the logarithm of the TOF, i.e. what is shown is the average and the STD of log10 [TOF/(s−1 cell−1 )]. The nominal and averaged values for this quantity differ by 0.5 and the STD is 2.5, 4
1.5
CO des. cus CO ad. cus O2 des. cus O2 ad. cus Ocus/COcus Obr/COcus O ad. br/cus CO ad. br O diff. br CO des. br Ocus/CObr O2 des. br/cus O2 ad. br O2 des. br O diff. br/cus O diff. cus/br O diff. cus CO diff. br CO diff. br/cus CO diff. cus/br CO diff. cus Obr/CObr
1.25
Si
1
0.75
0.5
0.25
0 1000
100000
10000
sample size Figure 2: Convergence of the sensitivity measure Si with the number of Quasi-Monte-Carlo points. Shown are the sensitivities of the log10 [TOF/(s−1 cell−1 )] with respect to the logarithms of the rate constants of all elementary steps for the model for the CO oxidation on RuO2 (110).
which corresponds to a factor three respectively 300 on a linear scale. Additionally, Fig. 1 shows the results for the coverages on both site types, i.e. the relative number of empty (ecus/ebr), CO-covered (COcus/CObr) or O-covered (Ocus/Obr) cus and br sites. We find a similar picture for the CO and oxygen coverages on cus sites as for the TOF. For these the STD is not much smaller than one and the difference between nominal and average is of similar size. For the coverages on the br sites, the STD is a little smaller, but still sizable. Only for the coverages of empty sites, we find a good agreement between nominal and averaged values and also the STD is rather small. So, we can only safely predict that the surface will be not empty. Whether the surface is active or inactive, CO or oxygen covered or in a mixed phase is outside the accuracy of the model.
5
Sensitivity Analysis. Having large uncertainties, the question naturally arises, which microscopic errors are responsible for these. Derivative based local sensitivity analysis (LSA)[14] applies only to models which behave at least close to linear over the whole range of parameters. With RCs varying by orders of magnitude, we can not expect the model to behave linear and we need a global approach, which can properly deal with the existing non-linearities. We therefore employ a probability distribution based global sensitivity analysis (GSA), which allows to extract sensitivity information from the data for GUQ. For this, we consider the distribution based sensitivity index (DSI) Si = dxi f (xi )D(xi ) = dxi f (xi ) dy|F (y) − F (y|xi )|,
(1)
where f (xi ) is the marginal PDF of the input parameter xi . F (y) is the cumulative distribution function (CDF) of the output y and F (y|xi ) is the CDF, which arises when we fix xi . D(xi ) is the distance between F (y) and F (y|xi ) measured in the L1 -norm, which equals a Wasserstein metric in one dimension[18]. If y is independent of xi both CDFs would agree and Si = 0. If y depends sensitively on xi , we would see large differences, when we vary xi and Si would be large. To estimate Si from data, we decompose the data set along xi -axis into M subsets, such that we have the same number sampling points in each subset. We approximate Si by M 1 emp. (y)|, dy|F emp. (y) − Fm Si ≈ M m=1
(2)
where F emp. (y) is the linearly interpolated empirical CDF which results from emp. the whole data set and Fm (y) is the linearly interpolated empirical CDF[19]
resulting from the data in the mth subset. The above approach is similar to the approach proposed by Plischke et al.[20] with the difference of CDF instead of PDF based distance measures D(xi ) and the way to estimate it. The DSI Si can be interpreted as the uncertainty, which is caused by the uncertainty of the parameter xi . This interpretation can be motivated by the metric D(xi ), which measures the distance between two distribution in terms of an effective shift on the y-axis. Indeed, if both distributions have the same shape 6
but differ only in the mean, this distance is the absolute difference between both means. As F (y) is the average of F (y|xi ), the DSI Si measures the fluctuations of these shifts, i.e. the variability induced by the uncertainty of the parameter xi [13]. We want to investigate the logarithm to the base ten of the TOF, and the DSI must be interpreted in this setting. A DSI of one therefore corresponds to an induced uncertainty of one order of magnitude in the TOF and a DSI of two corresponds to two orders of magnitude, etc.. The DSI shares this interpretation with the so-called total sensitivity index[14]. However, in sampling approaches for the latter, the number of required data points scales linearly with the number of parameters[20]. Another GSA approach consist of averaging local derivative information[21], which has the drawback that derivatives are extremely cumbersome to estimate from 1p-kMC[9]. Polynomial Chaos based approaches might overcome both limitations[22], but we need to be able to represent the output as a low order polynomial of the input variables. The sensitivity indices for the TOF (precisely log10 [TOF]) with respect to all RCs for M = 16 are shown in Figure 2 for different sample sizes. Adsorptions are represented by solid lines and crosses, desorptions by dashed lines and open circles. Already at rather small sample sizes, the six most important elementary steps are properly identified with well enough converged DSIs for most purposes. The five most important steps are the ad- and desorptions of CO (blue) and O2 (red) on cus sites as well as the reaction of both adsorbates from adjacent cus sites (solid black). A little less important is the Obr + COcus reaction (Obr/COcus, solid gray). At larger sample sizes above 20, 000, we find that the O2 adsorption on two adjacent bridge and cus sites (O ad. br/cus, magenta), the CO adsorption on br sites (CO ad. br, turquoise) and the oxygen diffusion between two bridge sites (O diff. br, brown) have DSIs of 0.1 or higher. Also, that other reactions seem to have a DSI below 0.1, is already established at this relatively small sample size. Notably, those smallest Si seem to converge the slowest. However, with the variability interpretation of the DSI, these correspond to an induced uncertainty of at most 20% for the non-logarithmic 7
TOF. Considering the two order of magnitude total uncertainty, they seem to be sufficiently well converged. Not only the convergence with respect to the number of QMC points is important, but also the choice of M . In the supplementary information[13], we provide such convergence plots as Figure 2 for the choices M = 4 and M = 64. All three choices of M converge to the almost the same values for the DSIs when increasing the the sample size. We therefore concluded that M = 16 is a safe choice and introduces only a small bias. Also this allows us to judge to some extend about the accuracy of those DSIs, which converge slowly with increasing sample size. At smaller choices of M these converge significantly faster and they all converge to the same values as for M = 16. Thus our choice of M and of the number of QMC points provides a sufficient accuracy to estimate the DSIs. Another important aspect is the impact of the necessarily approximate kMC simulations from which we generated the data set. Here, we have two sources of errors: I) a systematic error due to an insufficient number of relaxation steps and II) a statistical error due to the estimation of the stationary TOF from a finite number of kMC steps. We have performed tests for 105 /105 , 106 /106 and 107 /107 relaxation/sampling steps, see the supplementary information[13]. We still see significant quantitative differences between 106 /106 and 107 /107 relaxation/sampling steps. We attribute this difference to insufficient relaxation to steady state, as a study with 107 relaxation but only 106 sampling steps show only very small quantitative deviations from the 107 /107 case. To be sure we recalculated the first 8000 QMC points with 108 /108 relaxation/sampling steps and the corresponding DSI estimates agree well with those for the afore mentioned cases with only 107 relaxation steps. We conclude that the approach seems to be rather robust against sampling errors. However, the systematic error due to insufficient relaxation has a significant impact. Discussion. The huge uncertainty of 1p-kMC results could have been expected. The TOF is a homogeneous function of grade one of the RCs, i.e. if we multiply all RCs with the same number the TOF changes by this number[9]. If we want
8
2
distribution based main effect local
sensitivity index
1.5
1
0.5
0
ff. br
us
s
cu
O
r
.b
ad
r/c
.b s
cu
O
s/
/C
di
br
O
O
O
C
ad O2
us
s
cu
.c
cu
O
C
s.
ad O2
us
s
cu
.c
ad
s.
de
de O2
O
C
O
C
Figure 3: Comparison of the distribution based sensitivity index (gray) with the main effect sensitivity index (green) and local sensitivity analysis (blue) for the CO oxidation on RuO2 (110). The two global sensitivity measures give almost identical results, while local sensitivity analysis assign zero sensitivity to the last three reactions.
9
the TOF to be accurate by a factor of two, the important energetic parameters need to be accurate to at least 40 meV. Also, our variation means that we change the relative binding energy of both adsorbates by up to 0.5 eV. We can not expect that the surface stays oxygen covered as at the nominal settings. Indeed, we expect that the surface gets CO covered when increasing the relative binding energy of CO; this happens at nominal parameter settings, when playing around with the partial pressures[3]. That the coverage of empty sites is relatively unaffected by the uncertainty, can be explained by the strong binding of both adsorbates. Regardless, how we change the binding energy, both like to bind to the surface and the coverage of empty sites stays small. The DSI indicates the impact of a rate constant’s variability and therefore connect to the notion of rate-determining steps, usually identified by a LSA[5, 9, 23]: When the local sensitivity is low for the whole parameter space, the corresponding DSI will be small. When the local sensitivity is large over a wide range of parameter settings the corresponding DSI will be large. In contrast to the classical approach, GSA therefore provides a more complete picture, identifying those steps which are most important over a large range of parameter settings. For the CO oxidation on RuO2 (110), we find a clear hierarchy in the importance of different elementary steps, with the strongest impact of reactions involving cus-sites. Among these the uncertainty of the COcus desorption has the highest impact with a respective effective factor of ∼20 of induced uncertainty for the TOF. Oxygen desorption and CO adsorption on cus are responsible for a factor ∼6, each, and Ocus/COcus and O ad. cus for ∼4 and ∼3, respectively. The Obr/COcus is the only reaction with an effective factor > 2, which involves br sites. For a more detailed discussion, figure 3 shows a comparison of the DSI (gray bars) with the variance based main effect (ME, unnormalized, green bars))[14] and local sensitivity analysis (blue) for the nine most important RCs[13]. The ME can be estimated using a similar approach as for the DSI[13, 20]. The sensitivity index for LSA is simply the STD of the log of the respective RC times the corresponding partial derivative, which we have obtained using Coupled 10
Finite Differences[9, 24]. Both, DSI and ME, provide almost identical results supporting the interpretation of the DSI as induced variability. But it should be noted that the main effect might assign zero sensitivity to important parameters, which is not the case for the DSI. Notably, the importance of the cus sites and the five most important steps have also been identified by LSA, albeit with significantly overestimated importance (for a linear dependence, the ME should agree with the defined LSA measure). Unlike GSA, where we have Si ≈ 0.4 for Obr/COcus, and Si > 0.1 for O ad. br/cus and CO ad. br, the LSA assigns zero importance to these three reactions. These can be found by LSA for increased CO partial pressures[9, 13]. However, these local sensitivities stay comparatively small and appear only in a narrow range of CO partial pressures. The significant global sensitivity of the Obr/COcus can not be explained from this LSA. We conclude that this reaction must have a significant local sensitivity over a wide range of parameter settings, albeit not in the vicinity of the nominal settings. Also, the small importance of the oxygen diffusion on br sites can not be identified by LSA even for changed CO partial pressure. Although we can not even predict whether catalyst is active or not, many of the qualitative findings from previous studies[5, 9] about important reaction steps therefore seem to be rather robust for the considered model and we can still draw conclusions on the atomistic driving forces. We explain this rather good agreement between LSA and GSA by the dominance of the relative binding stability of CO and oxygen on the surface and we can test this relative stability by varying the partial pressures. Also binding between br and cus sites differs significantly, such that the dominance of the cus sites survives the parameter variation. In general, we can not expect LSA to always be sufficient, especially not for more complex reaction networks. It must be noted, that we have considered the worst case scenario, with evenly distributed and completely uncorrelated errors. In reality, there is always correlation between parameters. For instance, there are over- or underbinding density functionals, such that errors in the binding energy point rather into the same direction. Other correlations, such as the Brønsted-Evans-Polanyi 11
principle (BEP), correlate adsorption energies with barriers. Creating a PDF, which represents these correlations, is a challenging task[21] and beyond the scope of this study. However, evenly distributed PDFs might serve as a starting point for a GSA based identification of those uncertainties and correlations, which are most promising for spending computational or experimental resources for estimating them. This could be done in a hierarchical setting, starting with a rather inaccurate semi-empirical parameter set and spending expensive first-principles based parameter and correlation estimation only for the most important factors. For the problem at hand, our analysis indicates that resources should be directed to the investigation of the CO ad/desorption on cus sites, if we want to improve the accuracy of the model for the considered reaction conditions. This could be done by higher level electronic structure methods, or by dedicated experiments. Once more accurate energies are available a new sensitivity analysis should be performed to be able to select those energies, which should be improved. This will most likely affect the oxygen adsorption on cus sites, the Obr/COcus or the Ocus/COcus. Conclusions and Outlook. We have presented a distribution based approach for global sensitivity analysis, which allows to quantify the error propagation in 1p-kMC models with reasonable costs. Its advantages are the interpretation as induced variability and the possibility to estimate it from a single given data set. Using the CO oxidation on RuO2 (110) as an example, we demonstrated which reliable conclusions can be drawn from an uncertain model. Our findings suggest, that one can not always trust the quantitative findings of 1p-kMC. On the other hand, trends and qualitative insights into the mechanism might still be accurate enough. For a well-founded discussion, a way for an error-aware construction and analysis of multiscale models is needed and, in our opinion, the presented approach is an important step into this direction. While we employed it to quantify error propagation, the approach is ready
12
to be used in the context of computational materials screening[2, 25]. For this, the parametric dependence is expressed in terms of semi-empirical theories, with computationally cheap, material specific parameters, some fixed, universal parameters and the computationally expensive deviation from the “truth”. GSA could serve as a mean to identify the semi-empirical theory, where these expensive deviations have the lowest impact, and the corresponding most important (cheap) indicators. As for a proper identification of the error propagation, a reasonable PDF, ideally incorporating the DFT error, for the parameter distribution is crucial for the approach to work. Candidates can then be selected based on probabilistic measures using an updated, material specific PDF. Future development of the approach will focus on improving the convergence properties and also the use of randomized QMC points, in order to be able to estimate the accuracy of the sampled DSI[16]. Also more work needs to be directed to make the results better interpretable. Currently, the sensitivity measures provide the global importance, but only limited insight into the ongoing chemistry. Supplementary Material. See supplementary material for details on the CO oxidation model and the distribution based sensitivity analysis as well as results for different parameter settings and local sensitivity analysis. Acknowledgment. This research was carried out in the framework of Matheon supported by Einstein Foundation Berlin. References [1] M.
Stamatakis,
catalytic
D.
reactions
G. via
Vlachos, kinetic
Unraveling
Monte
the
Carlo
complexity
simulation:
of
Cur-
rent status and frontiers, ACS Catal. 2 (12) (2012) 2648–2663. doi:http://dx.doi.org/10.1021/cs3005709. [2] J. E. Sutton, ing
relations
D. G. Vlachos, and
Effect of errors in linear scal-
Brønsted-Evans-Polanyi
13
relations
on
activ-
ity
and
selectivity
maps,
J.
Catal.
338
(2016)
273
–
283.
doi:http://dx.doi.org/10.1016/j.jcat.2016.03.013. [3] K.
Reuter,
M.
simulations
for
Scheffler,
First-principles
heterogeneous
CO oxidation at RuO2 (110),
catalysis:
kinetic
Monte
Application
to
Carlo the
Phys. Rev. B 73 (2006) 045433.
doi:http://dx.doi.org/10.1103/PhysRevB.73.045433. [4] B. Temel, H. Meskine, K. Reuter, M. Scheffler, H. Metiu, Does phenomenological kinetics provide an adequate description of heterogeneous catalytic reactions?, J. Chem. Phys. 126 (20) (2007) 204711. doi:http://dx.doi.org/10.1063/1.2741556. [5] H. Meskine, S. Matera, M. Scheffler, K. Reuter, H. Metiu, Examination
of
the
concept
of
degree
of
rate
control
by
first-
principles kinetic Monte Carlo simulations, Surf. Sci. 603 (2009) 1724. doi:http://dx.doi.org/10.1016/j.susc.2008.08.036. [6] G. J. Herschlag, S. Mitran, G. Lin, A consistent hierarchy of generalized kinetic equation approximations to the master equation applied to surface catalysis, J Chem. Phys. 142 (23) (2015) 234703. doi:http://dx.doi.org/10.1063/1.4922515. [7] D.-J. Liu, and
J. W. Evans, Transitions between strongly correlated
random
steady-states
for
catalytic
CO-oxidation
on
sur-
faces at high-pressure, J. Chem. Phys. 142 (13) (2015) 134703. doi:http://dx.doi.org/10.1063/1.4916380. [8] P. Gelß,
S. Matera,
C. Sch¨ utte,
without kinetic Monte Carlo: CO oxidation model,
Solving the master equation
Tensor train approximations for a
J. Comp. Phys. 314 (2016) 489 – 502.
doi:http://dx.doi.org/10.1016/j.jcp.2016.03.025. [9] M. J. Hoffmann, F. Engelmann, S. Matera, A practical approach to the sensitivity analysis for kinetic monte carlo simulation of 14
heterogeneous catalysis,
J. Chem. Phys. 146 (4) (2017) 044118.
doi:http://dx.doi.org/10.1063/1.4974261. [10] K.
S.
Exner,
F.
Heß, and
H.
bined
experiment
istry:
Stairway to heaven?,
Over,
theory
A.
approach
P.
Seitsonen,
Com-
surface
chem-
in
Surf. Sci. 640 (2015) 165 – 180.
doi:http://dx.doi.org/10.1016/j.susc.2015.01.006. [11] A. Farkas, F. Hess, H. Over, Experiment-based kinetic monte carlo simulations: Co oxidation over ruo2(110), J. Phys. Chem. C 116 (1) (2012) 581–591. doi:http://dx.doi.org/10.1021/jp204703p. [12] S. Pogodin, N. L´ opez, A more accurate kinetic monte carlo approach to a monodimensional surface reaction: The interaction of oxygen with the ruo2(110) surface, ACS Catal. 4 (7) (2014) 2328–2332. doi:http://dx.doi.org/10.1021/cs500414p. [13] See supplemental material at [URL should be inserted by elsevier] for details. [14] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, S. Tarantola, Global sensitivity analysis: the primer, John Wiley & Sons, 2008. [15] I. Sobol’, On the distribution of points in a cube and the approximate evaluation of integrals, USSR Comp. Math. Math. Phys. 7 (4) (1967) 86 – 112. doi:http://dx.doi.org/10.1016/0041-5553(67)90144-9. [16] P. L’Ecuyer, C. Lemieux, Recent advances in randomized quasi-monte carlo methods, in:
M. Dror, P. L’Ecuyer, F. Szidarovszky (Eds.),
Modeling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications, Springer US, Boston, MA, 2005, pp. 419–474. doi:http://dx.doi.org/10.1007/0-306-48102-2_20.
15
[17] M. J. Hoffmann, S. Matera, K. Reuter, kmos: A lattice kinetic Monte Carlo framework, Comp. Phys. Comm. 185 (7) (2014) 2138 – 2150. doi:http://dx.doi.org/10.1016/j.cpc.2014.04.003. [18] P.
J.
for
Bickel,
the
D.
bootstrap,
A.
Freedman,
Ann.
Some
Statist.
9
asymptotic
(6)
(1981)
theory
1196–1217.
doi:http://dx.doi.org/10.1214/aos/1176345637. [19] F.
Perez-Cruz,
Kullback-Leibler
tinuous
distributions,
posium
on
in:
divergence
2008
Information
estimation
IEEE
Theory,
of
International
2008,
pp.
conSym-
1666–1670.
doi:http://dx.doi.org/10.1109/ISIT.2008.4595271. [20] E. Plischke, E. Borgonovo, C. L. Smith, Global sensitivity measures from given data, Eur. J. Oper. Res. 226 (3) (2013) 536 – 550. doi:http://dx.doi.org/10.1016/j.ejor.2012.11.047. [21] J. E. Sutton, W. Guo, M. A. Katsoulakis, D. G. Vlachos, Effects of
correlated
based
parameters
chemical
kinetic
and
uncertainty
modelling,
in
Nat.
electronic-structure-
Chem.
(2016)
331
–
P.
Maˆıtre,
337doi:dx.doi.org/10.1038/nchem.2454. [22] A.
Alexanderian,
O.
M.
Knio,
tic
chemical
F.
Rizzi,
M.
Preconditioned kinetics,
J.
Sci.
Rathinam,
Bayesian Comp.
O.
regression 58
(3)
Le for
stochas-
(2014)
592–626.
doi:http://dx.doi.org/10.1007/s10915-013-9745-5. [23] C. T. Campbell, Future directions and industrial perspectives micro-and macro-kinetics: their relationship in heterogeneous catalysis, Top. Catal. 1 (1994) 353. doi:10.1007/BF01492288. [24] D. F. Anderson, An efficient finite difference method for parameter sensitivities of continuous time markov chains, SIAM J. Numer. Anal. 50 (5) (2012) 2237–2258. doi:http://dx.doi.org/10.1137/110849079.
16
[25] M. P. Andersson, T. Bligaard, A. Kustov, K. E. Larsen, J. Greeley,
T.
Johannessen,
C.
H.
Christensen,
J. K. Nørskov,
ward computational screening in heterogeneous catalysis:
To-
Pareto-
optimal methanation catalysts, J. Catal. 239 (2) (2006) 501–506. doi:http://dx.doi.org/10.1016/j.jcat.2006.02.016.
17
Highlights
First-principles kinetic Monte Carlo for the CO oxidation on Ru2(110) Impact of DFT errors on turn-over frequency (TOF) and coverages Global sensitivity analysis Identification of most important errors/elementary steps
Uncertainties due to DFT error
2
1
1
0.5
0
0 TOF
CObr
COcus
Obr
Ocus
ebr
ecus
coverage
log10 TOF (s-1cell-1)
Graphical Abstract
S0009-2614(17)30161-6 http://dx.doi.org/10.1016/j.cplett.2017.02.043 CPLETT 34556
To appear in:
Chemical Physics Letters
Please cite this article as: S. Döpking, S. Matera, Error propagation in first-principles kinetic Monte Carlo simulation, Chemical Physics Letters (2017), doi: http://dx.doi.org/10.1016/j.cplett.2017.02.043
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Error propagation in first-principles kinetic Monte Carlo simulation Sandra D¨ opking, Sebastian Matera∗ Institute for Mathematics, Freie Universit¨ at Berlin, Arnimallee 6, 14195 Berlin, Germany
Abstract First-principles kinetic Monte Carlo models allow for the modeling of catalytic surfaces with predictive quality. This comes at the price of non-negligible errors induced by the underlying approximate density functional calculation. On the example of CO oxidation on RuO2 (110), we demonstrate a novel, efficient approach to global sensitivity analysis, with which we address the error propagation in these multiscale models. We find, that we can still derive the most important atomistic factors for reactivity, albeit the errors in the simulation results are sizable. The presented approach might also be applied in the hierarchical model construction or computational catalyst screening. Keywords: stochastic simulation, sensitivity and uncertainty analysis, Quasi Monte Carlo, heterogeneous catalysis, first-principles kinetic Monte Carlo 2010 MSC: 65C05, 65C40, 82C99, 62-07 First-principles kinetic Monte Carlo (1p-kMC) have become popular for modeling the chemical kinetics in heterogeneous catalysis[1]. This popularity results from the predictive quality of the approach, enabling an understanding of the detailed surface chemistry without the possible masking by the fitting of parameters. In 1p-kMC, the mechanism and corresponding energy barriers are instead obtained from first-principles electronic structure methods and kinetic Monte ∗ Corresponding
author Email addresses: [email protected] (Sandra D¨ opking), [email protected] (Sebastian Matera)
Preprint submitted to Chemical Physics Letters
February 15, 2017
Carlo simulation then provides unbiased estimates of the solution of the resulting master equation. Albeit being predictive, practically applicable electronic structure methods, such as Density Functional Theory (DFT), provide only approximate results. Commonly reported uncertainties for DFT barriers are in the range of 0.2 − 0.3 eV. At low to medium temperatures, the corresponding rate constants carry then an uncertainty of some orders of magnitude. Hitherto, the propagation of these errors to the 1p-kMC simulations results has not been systematically addressed, only in the context of the much cheaper but also less accurate meanfield microkinetic modeling[2]. With the focus on the former, we employ Quasi Monte Carlo (QMC) based global uncertainty quantification (GUQ) to estimate the errors in the final 1p-kMC results. We further demonstrate a novel approach to global sensitivity analysis (GSA), which allows to quantify the error induced by each uncertain energy barrier using the same set of QMC points. This allows a global view on the main driving forces in the mechanism and generally provides a route to an error-aware analysis of such multiscale models. We exemplify this approach on the CO oxidation on the RuO2 (110) surface using the 1p-kMC model by Reuter[3], which is one of the most thoroughly investigated 1p-kMC models[4–9]. It is therefore a good test example, also because limitations of this model have been discussed in the literature[10–12]. For this model, we indeed find large possible errors, but also that many qualitative mechanistic conclusions seem to be not much affected by the DFT inaccuracies. Uncertainty Analysis. For an uncertainty analysis, the input parameters xi , the barriers or the rate constants in our case, must be considered as random variables, with a given joint probability density function (PDF). The output of a DFT code is deterministic, so this PDF must be interpreted in a information theoretical setting and reflects our knowledge about these parameters. We consider the CO oxidation on RuO2 (110) for CO and O2 partial pressures of 1 bar and a temperature of 600 K, in order to make contact with existing theoretical studies[4–6, 9], which investigate the behavior around these reaction
2
1
1
0.5
0
0
coverage
log10 TOF (s-1cell-1)
2
average nominal -1 TOF
CObr
COcus
Obr
Ocus
ebr
ecus
Figure 1: Turnover-frequencies (TOF) and coverages for CO occupied, O occupied and and empty bridge and cus sites for the CO oxidation on RuO2 . Grey bars represent the parameter averaged values and the error bars represent the standard deviation resulting from parameter uncertainty. For comparison, the values at the nominal parameter settings are displayed as green bars.
conditions. For these conditions, we also expect a rapid change of the reactivity when varying the energetic parameters within the DFT error; something which is particularly challenging for GUQ and GSA in high dimensions and therefore a good test for the methodology. We provide additional results for a CO partial pressure of one mbar in the supplement. The model considers two types of adsorption sites, bridge (br) and coordinatively unsaturated (cus) sites, which are arranged on a square lattice as alternating columns. The mechanism includes unimolecular adsorptions of CO on both types of sites, dissociative oxygen adsorption on nearest neighboring sites and the corresponding reverse elementary steps. Further, both adsorbates can make diffusive hops to neighboring empty sites and adsorbed CO and O on neighboring sites can react to gaseous CO2 . Altogether, this makes 22 uncertain parameters, for which rate constants (RC) and energy barriers are summarized in the supplementary material[13]. For simplicity, we assume that errors appear uncorrelated and the logarithm (to the base 10) of every RC is uniformly dis-
3
tributed between ±2 around its nominal value. This corresponds to uncertain DFT barriers with a uniformly distributed error in the bound [−0.25, 0.25] eV. The only exception from this rule are the adsorption steps, which have no barriers in the nominal setting. For these, we assume the nominal value to be the upper bound and allow variation only to smaller values. The lower bounds for these steps correspond then to a by a factor hundred smaller RCs. By the linear dependence of these RCs on the partial pressures, we therefore also address the dependence on these (at least to some extend). In order to access the impact of parameter uncertainty, we need to sample the parameter space and perform one kMC simulation for each sampling point. The results for these points are then used to calculate parameter-averaged key quantities, like the mean or the variance. The later then serves as a measure of the uncertainty of the simulation results[14]. We employ uniformly distributed Quasi Monte Carlo (QMC) points for this sampling of parameter space. This choice was motivated by the more even distribution of the points and a faster convergence rate in medium high dimensions compared to classical Monte Carlo sampling. In detail, we used 218 ≈260, 000 points from Sobol’s sequence for QMC[15, 16] and each kMC simulation was performed for a lattice of 20 × 20 units cells using the kmos kMC package[17]. Focusing on stationary operation, we employed 107 kMC steps for initial relaxation to reach the steady state and further 107 steps for sampling. This amounts to a total cost of ∼1000 CPU-hrs for the whole data set. The results from this parameter averaging are depicted in Figure 1. Averages are represented as gray bars and the error bars indicate the corresponding standard deviation (STD). For comparison, we have added the nominal values as green bars (i.e. those values which are obtained with the DFT barriers from Reuter[3]). The most important output from a 1p-kMC simulation is the turn-over frequency (TOF), i.e. the number of CO to CO2 reactions per second and surface unit cell. As for the RCs, we consider the logarithm of the TOF, i.e. what is shown is the average and the STD of log10 [TOF/(s−1 cell−1 )]. The nominal and averaged values for this quantity differ by 0.5 and the STD is 2.5, 4
1.5
CO des. cus CO ad. cus O2 des. cus O2 ad. cus Ocus/COcus Obr/COcus O ad. br/cus CO ad. br O diff. br CO des. br Ocus/CObr O2 des. br/cus O2 ad. br O2 des. br O diff. br/cus O diff. cus/br O diff. cus CO diff. br CO diff. br/cus CO diff. cus/br CO diff. cus Obr/CObr
1.25
Si
1
0.75
0.5
0.25
0 1000
100000
10000
sample size Figure 2: Convergence of the sensitivity measure Si with the number of Quasi-Monte-Carlo points. Shown are the sensitivities of the log10 [TOF/(s−1 cell−1 )] with respect to the logarithms of the rate constants of all elementary steps for the model for the CO oxidation on RuO2 (110).
which corresponds to a factor three respectively 300 on a linear scale. Additionally, Fig. 1 shows the results for the coverages on both site types, i.e. the relative number of empty (ecus/ebr), CO-covered (COcus/CObr) or O-covered (Ocus/Obr) cus and br sites. We find a similar picture for the CO and oxygen coverages on cus sites as for the TOF. For these the STD is not much smaller than one and the difference between nominal and average is of similar size. For the coverages on the br sites, the STD is a little smaller, but still sizable. Only for the coverages of empty sites, we find a good agreement between nominal and averaged values and also the STD is rather small. So, we can only safely predict that the surface will be not empty. Whether the surface is active or inactive, CO or oxygen covered or in a mixed phase is outside the accuracy of the model.
5
Sensitivity Analysis. Having large uncertainties, the question naturally arises, which microscopic errors are responsible for these. Derivative based local sensitivity analysis (LSA)[14] applies only to models which behave at least close to linear over the whole range of parameters. With RCs varying by orders of magnitude, we can not expect the model to behave linear and we need a global approach, which can properly deal with the existing non-linearities. We therefore employ a probability distribution based global sensitivity analysis (GSA), which allows to extract sensitivity information from the data for GUQ. For this, we consider the distribution based sensitivity index (DSI) Si = dxi f (xi )D(xi ) = dxi f (xi ) dy|F (y) − F (y|xi )|,
(1)
where f (xi ) is the marginal PDF of the input parameter xi . F (y) is the cumulative distribution function (CDF) of the output y and F (y|xi ) is the CDF, which arises when we fix xi . D(xi ) is the distance between F (y) and F (y|xi ) measured in the L1 -norm, which equals a Wasserstein metric in one dimension[18]. If y is independent of xi both CDFs would agree and Si = 0. If y depends sensitively on xi , we would see large differences, when we vary xi and Si would be large. To estimate Si from data, we decompose the data set along xi -axis into M subsets, such that we have the same number sampling points in each subset. We approximate Si by M 1 emp. (y)|, dy|F emp. (y) − Fm Si ≈ M m=1
(2)
where F emp. (y) is the linearly interpolated empirical CDF which results from emp. the whole data set and Fm (y) is the linearly interpolated empirical CDF[19]
resulting from the data in the mth subset. The above approach is similar to the approach proposed by Plischke et al.[20] with the difference of CDF instead of PDF based distance measures D(xi ) and the way to estimate it. The DSI Si can be interpreted as the uncertainty, which is caused by the uncertainty of the parameter xi . This interpretation can be motivated by the metric D(xi ), which measures the distance between two distribution in terms of an effective shift on the y-axis. Indeed, if both distributions have the same shape 6
but differ only in the mean, this distance is the absolute difference between both means. As F (y) is the average of F (y|xi ), the DSI Si measures the fluctuations of these shifts, i.e. the variability induced by the uncertainty of the parameter xi [13]. We want to investigate the logarithm to the base ten of the TOF, and the DSI must be interpreted in this setting. A DSI of one therefore corresponds to an induced uncertainty of one order of magnitude in the TOF and a DSI of two corresponds to two orders of magnitude, etc.. The DSI shares this interpretation with the so-called total sensitivity index[14]. However, in sampling approaches for the latter, the number of required data points scales linearly with the number of parameters[20]. Another GSA approach consist of averaging local derivative information[21], which has the drawback that derivatives are extremely cumbersome to estimate from 1p-kMC[9]. Polynomial Chaos based approaches might overcome both limitations[22], but we need to be able to represent the output as a low order polynomial of the input variables. The sensitivity indices for the TOF (precisely log10 [TOF]) with respect to all RCs for M = 16 are shown in Figure 2 for different sample sizes. Adsorptions are represented by solid lines and crosses, desorptions by dashed lines and open circles. Already at rather small sample sizes, the six most important elementary steps are properly identified with well enough converged DSIs for most purposes. The five most important steps are the ad- and desorptions of CO (blue) and O2 (red) on cus sites as well as the reaction of both adsorbates from adjacent cus sites (solid black). A little less important is the Obr + COcus reaction (Obr/COcus, solid gray). At larger sample sizes above 20, 000, we find that the O2 adsorption on two adjacent bridge and cus sites (O ad. br/cus, magenta), the CO adsorption on br sites (CO ad. br, turquoise) and the oxygen diffusion between two bridge sites (O diff. br, brown) have DSIs of 0.1 or higher. Also, that other reactions seem to have a DSI below 0.1, is already established at this relatively small sample size. Notably, those smallest Si seem to converge the slowest. However, with the variability interpretation of the DSI, these correspond to an induced uncertainty of at most 20% for the non-logarithmic 7
TOF. Considering the two order of magnitude total uncertainty, they seem to be sufficiently well converged. Not only the convergence with respect to the number of QMC points is important, but also the choice of M . In the supplementary information[13], we provide such convergence plots as Figure 2 for the choices M = 4 and M = 64. All three choices of M converge to the almost the same values for the DSIs when increasing the the sample size. We therefore concluded that M = 16 is a safe choice and introduces only a small bias. Also this allows us to judge to some extend about the accuracy of those DSIs, which converge slowly with increasing sample size. At smaller choices of M these converge significantly faster and they all converge to the same values as for M = 16. Thus our choice of M and of the number of QMC points provides a sufficient accuracy to estimate the DSIs. Another important aspect is the impact of the necessarily approximate kMC simulations from which we generated the data set. Here, we have two sources of errors: I) a systematic error due to an insufficient number of relaxation steps and II) a statistical error due to the estimation of the stationary TOF from a finite number of kMC steps. We have performed tests for 105 /105 , 106 /106 and 107 /107 relaxation/sampling steps, see the supplementary information[13]. We still see significant quantitative differences between 106 /106 and 107 /107 relaxation/sampling steps. We attribute this difference to insufficient relaxation to steady state, as a study with 107 relaxation but only 106 sampling steps show only very small quantitative deviations from the 107 /107 case. To be sure we recalculated the first 8000 QMC points with 108 /108 relaxation/sampling steps and the corresponding DSI estimates agree well with those for the afore mentioned cases with only 107 relaxation steps. We conclude that the approach seems to be rather robust against sampling errors. However, the systematic error due to insufficient relaxation has a significant impact. Discussion. The huge uncertainty of 1p-kMC results could have been expected. The TOF is a homogeneous function of grade one of the RCs, i.e. if we multiply all RCs with the same number the TOF changes by this number[9]. If we want
8
2
distribution based main effect local
sensitivity index
1.5
1
0.5
0
ff. br
us
s
cu
O
r
.b
ad
r/c
.b s
cu
O
s/
/C
di
br
O
O
O
C
ad O2
us
s
cu
.c
cu
O
C
s.
ad O2
us
s
cu
.c
ad
s.
de
de O2
O
C
O
C
Figure 3: Comparison of the distribution based sensitivity index (gray) with the main effect sensitivity index (green) and local sensitivity analysis (blue) for the CO oxidation on RuO2 (110). The two global sensitivity measures give almost identical results, while local sensitivity analysis assign zero sensitivity to the last three reactions.
9
the TOF to be accurate by a factor of two, the important energetic parameters need to be accurate to at least 40 meV. Also, our variation means that we change the relative binding energy of both adsorbates by up to 0.5 eV. We can not expect that the surface stays oxygen covered as at the nominal settings. Indeed, we expect that the surface gets CO covered when increasing the relative binding energy of CO; this happens at nominal parameter settings, when playing around with the partial pressures[3]. That the coverage of empty sites is relatively unaffected by the uncertainty, can be explained by the strong binding of both adsorbates. Regardless, how we change the binding energy, both like to bind to the surface and the coverage of empty sites stays small. The DSI indicates the impact of a rate constant’s variability and therefore connect to the notion of rate-determining steps, usually identified by a LSA[5, 9, 23]: When the local sensitivity is low for the whole parameter space, the corresponding DSI will be small. When the local sensitivity is large over a wide range of parameter settings the corresponding DSI will be large. In contrast to the classical approach, GSA therefore provides a more complete picture, identifying those steps which are most important over a large range of parameter settings. For the CO oxidation on RuO2 (110), we find a clear hierarchy in the importance of different elementary steps, with the strongest impact of reactions involving cus-sites. Among these the uncertainty of the COcus desorption has the highest impact with a respective effective factor of ∼20 of induced uncertainty for the TOF. Oxygen desorption and CO adsorption on cus are responsible for a factor ∼6, each, and Ocus/COcus and O ad. cus for ∼4 and ∼3, respectively. The Obr/COcus is the only reaction with an effective factor > 2, which involves br sites. For a more detailed discussion, figure 3 shows a comparison of the DSI (gray bars) with the variance based main effect (ME, unnormalized, green bars))[14] and local sensitivity analysis (blue) for the nine most important RCs[13]. The ME can be estimated using a similar approach as for the DSI[13, 20]. The sensitivity index for LSA is simply the STD of the log of the respective RC times the corresponding partial derivative, which we have obtained using Coupled 10
Finite Differences[9, 24]. Both, DSI and ME, provide almost identical results supporting the interpretation of the DSI as induced variability. But it should be noted that the main effect might assign zero sensitivity to important parameters, which is not the case for the DSI. Notably, the importance of the cus sites and the five most important steps have also been identified by LSA, albeit with significantly overestimated importance (for a linear dependence, the ME should agree with the defined LSA measure). Unlike GSA, where we have Si ≈ 0.4 for Obr/COcus, and Si > 0.1 for O ad. br/cus and CO ad. br, the LSA assigns zero importance to these three reactions. These can be found by LSA for increased CO partial pressures[9, 13]. However, these local sensitivities stay comparatively small and appear only in a narrow range of CO partial pressures. The significant global sensitivity of the Obr/COcus can not be explained from this LSA. We conclude that this reaction must have a significant local sensitivity over a wide range of parameter settings, albeit not in the vicinity of the nominal settings. Also, the small importance of the oxygen diffusion on br sites can not be identified by LSA even for changed CO partial pressure. Although we can not even predict whether catalyst is active or not, many of the qualitative findings from previous studies[5, 9] about important reaction steps therefore seem to be rather robust for the considered model and we can still draw conclusions on the atomistic driving forces. We explain this rather good agreement between LSA and GSA by the dominance of the relative binding stability of CO and oxygen on the surface and we can test this relative stability by varying the partial pressures. Also binding between br and cus sites differs significantly, such that the dominance of the cus sites survives the parameter variation. In general, we can not expect LSA to always be sufficient, especially not for more complex reaction networks. It must be noted, that we have considered the worst case scenario, with evenly distributed and completely uncorrelated errors. In reality, there is always correlation between parameters. For instance, there are over- or underbinding density functionals, such that errors in the binding energy point rather into the same direction. Other correlations, such as the Brønsted-Evans-Polanyi 11
principle (BEP), correlate adsorption energies with barriers. Creating a PDF, which represents these correlations, is a challenging task[21] and beyond the scope of this study. However, evenly distributed PDFs might serve as a starting point for a GSA based identification of those uncertainties and correlations, which are most promising for spending computational or experimental resources for estimating them. This could be done in a hierarchical setting, starting with a rather inaccurate semi-empirical parameter set and spending expensive first-principles based parameter and correlation estimation only for the most important factors. For the problem at hand, our analysis indicates that resources should be directed to the investigation of the CO ad/desorption on cus sites, if we want to improve the accuracy of the model for the considered reaction conditions. This could be done by higher level electronic structure methods, or by dedicated experiments. Once more accurate energies are available a new sensitivity analysis should be performed to be able to select those energies, which should be improved. This will most likely affect the oxygen adsorption on cus sites, the Obr/COcus or the Ocus/COcus. Conclusions and Outlook. We have presented a distribution based approach for global sensitivity analysis, which allows to quantify the error propagation in 1p-kMC models with reasonable costs. Its advantages are the interpretation as induced variability and the possibility to estimate it from a single given data set. Using the CO oxidation on RuO2 (110) as an example, we demonstrated which reliable conclusions can be drawn from an uncertain model. Our findings suggest, that one can not always trust the quantitative findings of 1p-kMC. On the other hand, trends and qualitative insights into the mechanism might still be accurate enough. For a well-founded discussion, a way for an error-aware construction and analysis of multiscale models is needed and, in our opinion, the presented approach is an important step into this direction. While we employed it to quantify error propagation, the approach is ready
12
to be used in the context of computational materials screening[2, 25]. For this, the parametric dependence is expressed in terms of semi-empirical theories, with computationally cheap, material specific parameters, some fixed, universal parameters and the computationally expensive deviation from the “truth”. GSA could serve as a mean to identify the semi-empirical theory, where these expensive deviations have the lowest impact, and the corresponding most important (cheap) indicators. As for a proper identification of the error propagation, a reasonable PDF, ideally incorporating the DFT error, for the parameter distribution is crucial for the approach to work. Candidates can then be selected based on probabilistic measures using an updated, material specific PDF. Future development of the approach will focus on improving the convergence properties and also the use of randomized QMC points, in order to be able to estimate the accuracy of the sampled DSI[16]. Also more work needs to be directed to make the results better interpretable. Currently, the sensitivity measures provide the global importance, but only limited insight into the ongoing chemistry. Supplementary Material. See supplementary material for details on the CO oxidation model and the distribution based sensitivity analysis as well as results for different parameter settings and local sensitivity analysis. Acknowledgment. This research was carried out in the framework of Matheon supported by Einstein Foundation Berlin. References [1] M.
Stamatakis,
catalytic
D.
reactions
G. via
Vlachos, kinetic
Unraveling
Monte
the
Carlo
complexity
simulation:
of
Cur-
rent status and frontiers, ACS Catal. 2 (12) (2012) 2648–2663. doi:http://dx.doi.org/10.1021/cs3005709. [2] J. E. Sutton, ing
relations
D. G. Vlachos, and
Effect of errors in linear scal-
Brønsted-Evans-Polanyi
13
relations
on
activ-
ity
and
selectivity
maps,
J.
Catal.
338
(2016)
273
–
283.
doi:http://dx.doi.org/10.1016/j.jcat.2016.03.013. [3] K.
Reuter,
M.
simulations
for
Scheffler,
First-principles
heterogeneous
CO oxidation at RuO2 (110),
catalysis:
kinetic
Monte
Application
to
Carlo the
Phys. Rev. B 73 (2006) 045433.
doi:http://dx.doi.org/10.1103/PhysRevB.73.045433. [4] B. Temel, H. Meskine, K. Reuter, M. Scheffler, H. Metiu, Does phenomenological kinetics provide an adequate description of heterogeneous catalytic reactions?, J. Chem. Phys. 126 (20) (2007) 204711. doi:http://dx.doi.org/10.1063/1.2741556. [5] H. Meskine, S. Matera, M. Scheffler, K. Reuter, H. Metiu, Examination
of
the
concept
of
degree
of
rate
control
by
first-
principles kinetic Monte Carlo simulations, Surf. Sci. 603 (2009) 1724. doi:http://dx.doi.org/10.1016/j.susc.2008.08.036. [6] G. J. Herschlag, S. Mitran, G. Lin, A consistent hierarchy of generalized kinetic equation approximations to the master equation applied to surface catalysis, J Chem. Phys. 142 (23) (2015) 234703. doi:http://dx.doi.org/10.1063/1.4922515. [7] D.-J. Liu, and
J. W. Evans, Transitions between strongly correlated
random
steady-states
for
catalytic
CO-oxidation
on
sur-
faces at high-pressure, J. Chem. Phys. 142 (13) (2015) 134703. doi:http://dx.doi.org/10.1063/1.4916380. [8] P. Gelß,
S. Matera,
C. Sch¨ utte,
without kinetic Monte Carlo: CO oxidation model,
Solving the master equation
Tensor train approximations for a
J. Comp. Phys. 314 (2016) 489 – 502.
doi:http://dx.doi.org/10.1016/j.jcp.2016.03.025. [9] M. J. Hoffmann, F. Engelmann, S. Matera, A practical approach to the sensitivity analysis for kinetic monte carlo simulation of 14
heterogeneous catalysis,
J. Chem. Phys. 146 (4) (2017) 044118.
doi:http://dx.doi.org/10.1063/1.4974261. [10] K.
S.
Exner,
F.
Heß, and
H.
bined
experiment
istry:
Stairway to heaven?,
Over,
theory
A.
approach
P.
Seitsonen,
Com-
surface
chem-
in
Surf. Sci. 640 (2015) 165 – 180.
doi:http://dx.doi.org/10.1016/j.susc.2015.01.006. [11] A. Farkas, F. Hess, H. Over, Experiment-based kinetic monte carlo simulations: Co oxidation over ruo2(110), J. Phys. Chem. C 116 (1) (2012) 581–591. doi:http://dx.doi.org/10.1021/jp204703p. [12] S. Pogodin, N. L´ opez, A more accurate kinetic monte carlo approach to a monodimensional surface reaction: The interaction of oxygen with the ruo2(110) surface, ACS Catal. 4 (7) (2014) 2328–2332. doi:http://dx.doi.org/10.1021/cs500414p. [13] See supplemental material at [URL should be inserted by elsevier] for details. [14] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, S. Tarantola, Global sensitivity analysis: the primer, John Wiley & Sons, 2008. [15] I. Sobol’, On the distribution of points in a cube and the approximate evaluation of integrals, USSR Comp. Math. Math. Phys. 7 (4) (1967) 86 – 112. doi:http://dx.doi.org/10.1016/0041-5553(67)90144-9. [16] P. L’Ecuyer, C. Lemieux, Recent advances in randomized quasi-monte carlo methods, in:
M. Dror, P. L’Ecuyer, F. Szidarovszky (Eds.),
Modeling Uncertainty: An Examination of Stochastic Theory, Methods, and Applications, Springer US, Boston, MA, 2005, pp. 419–474. doi:http://dx.doi.org/10.1007/0-306-48102-2_20.
15
[17] M. J. Hoffmann, S. Matera, K. Reuter, kmos: A lattice kinetic Monte Carlo framework, Comp. Phys. Comm. 185 (7) (2014) 2138 – 2150. doi:http://dx.doi.org/10.1016/j.cpc.2014.04.003. [18] P.
J.
for
Bickel,
the
D.
bootstrap,
A.
Freedman,
Ann.
Some
Statist.
9
asymptotic
(6)
(1981)
theory
1196–1217.
doi:http://dx.doi.org/10.1214/aos/1176345637. [19] F.
Perez-Cruz,
Kullback-Leibler
tinuous
distributions,
posium
on
in:
divergence
2008
Information
estimation
IEEE
Theory,
of
International
2008,
pp.
conSym-
1666–1670.
doi:http://dx.doi.org/10.1109/ISIT.2008.4595271. [20] E. Plischke, E. Borgonovo, C. L. Smith, Global sensitivity measures from given data, Eur. J. Oper. Res. 226 (3) (2013) 536 – 550. doi:http://dx.doi.org/10.1016/j.ejor.2012.11.047. [21] J. E. Sutton, W. Guo, M. A. Katsoulakis, D. G. Vlachos, Effects of
correlated
based
parameters
chemical
kinetic
and
uncertainty
modelling,
in
Nat.
electronic-structure-
Chem.
(2016)
331
–
P.
Maˆıtre,
337doi:dx.doi.org/10.1038/nchem.2454. [22] A.
Alexanderian,
O.
M.
Knio,
tic
chemical
F.
Rizzi,
M.
Preconditioned kinetics,
J.
Sci.
Rathinam,
Bayesian Comp.
O.
regression 58
(3)
Le for
stochas-
(2014)
592–626.
doi:http://dx.doi.org/10.1007/s10915-013-9745-5. [23] C. T. Campbell, Future directions and industrial perspectives micro-and macro-kinetics: their relationship in heterogeneous catalysis, Top. Catal. 1 (1994) 353. doi:10.1007/BF01492288. [24] D. F. Anderson, An efficient finite difference method for parameter sensitivities of continuous time markov chains, SIAM J. Numer. Anal. 50 (5) (2012) 2237–2258. doi:http://dx.doi.org/10.1137/110849079.
16
[25] M. P. Andersson, T. Bligaard, A. Kustov, K. E. Larsen, J. Greeley,
T.
Johannessen,
C.
H.
Christensen,
J. K. Nørskov,
ward computational screening in heterogeneous catalysis:
To-
Pareto-
optimal methanation catalysts, J. Catal. 239 (2) (2006) 501–506. doi:http://dx.doi.org/10.1016/j.jcat.2006.02.016.
17
Highlights
First-principles kinetic Monte Carlo for the CO oxidation on Ru2(110) Impact of DFT errors on turn-over frequency (TOF) and coverages Global sensitivity analysis Identification of most important errors/elementary steps
Uncertainties due to DFT error
2
1
1
0.5
0
0 TOF
CObr
COcus
Obr
Ocus
ebr
ecus
coverage
log10 TOF (s-1cell-1)
Graphical Abstract