Trade-Off between Accuracy and Universality in ... - ACS Publications

To screen heterogeneous catalysts in silico, the linear energy relationships derived from the Brønsted–Evans–Polanyi principle are extremely usef...
2 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCC

Trade-Off between Accuracy and Universality in Linear Energy Relations for Alcohol Dehydrogenation on Transition Metals Jérémie Zaffran, Carine Michel, Françoise Delbecq, and Philippe Sautet* University of Lyon, CNRS, Laboratoire de Chimie, UMR5182, Ecole Normale Supérieure de Lyon, Lyon, CEDEX 07, France S Supporting Information *

ABSTRACT: To screen heterogeneous catalysts in silico, the linear energy relationships derived from the Brønsted−Evans−Polanyi principle are extremely useful. They connect the reaction energy of a given elementary step to its activation energy, hence providing data that can be fed to kinetics models at a minimal cost. However, to ensure reasonable predictions, it is essential to control the statistical error intrinsic to this approach. We derived several types of linear energy relations for a series of CH and OH bond scissions in simple alcohol molecules on compact facets of seven transition metals (Co, Ni, Ru, Rh, Pd, Ir, and Pt) aiming at a single but accurate relation. The quality of the relation depends on its nature and/or on the manner the data are split: a single linear relation can be constructed for all metals together on the basis of the original Brønsted−Evans−Polanyi formulation with a mean absolute error smaller than 0.1 eV, whereas the more recent transition state scaling approach requires considering each metal individually to reach an equivalent accuracy. In addition, a close statistical analysis demonstrates that errors stemming from such predictive models are not uniform along the set of metals and of chemical reactions that is considered opening the road to a better control of error propagation.



INTRODUCTION Nowadays, the in silico design of catalysts can be achieved in principle using a multiscale bottom-up approach.1−3 First, one needs to build an extensive set of ab initio data, generally based on the density functional theory (DFT), that provides the reaction energy and the activation energy for each elementary step under consideration. Then, those energies can be input in a kinetic model (mean field approach or kinetic Monte Carlo) providing yields and selectivity as a function of the external conditions (pressure, temperature, and nature of the reactor).4,5 Following this strategy, the generation of the DFT data set is currently a real bottleneck, especially for a large reaction network. For instance, in the fast-growing field of biomass conversion, the selective hydrogenolysis of monosaccharide compounds or furanic compounds is a timely challenge with a potential reaction network containing thousands of elementary steps.6 Linear energy relationships are in this respect essential to the analysis of complex reaction networks, reducing considerably the load of DFT computations.7−14 For instance, Daoutidis and co-workers designed an automated identification of feasible mechanisms on the basis of the relations previously proposed by other groups10,11 and applied this tool to the glycerol conversion on transition metals, generating a network of 3300 reactions and 500 species.15 We can distinguish two main linear energy relationships. The scaling relation correlates the adsorption energy of a complex molecule with the adsorption energy of simple fragment such as single atom or diatomic fragment.16,17 The second class of energy relationship is derived © XXXX American Chemical Society

from the classic correlation between kinetics and thermodynamics evidenced by Brønsted experimentally.18 It has been rationalized by Bell19 and then Evans and Polanyi20,21 in the BEP relationship that connects activation energies (E⧧) and reaction energies (ΔE) within a given family of elementary steps: E ⧧ = αΔE + β

The slope α characterizes the earliness of the transition state (the smaller the value, the earlier the transition state), and the intercept β is the reference activation energy for an athermic reaction. The use of this relationship was limited to organic chemistry for decades.22 In heterogeneous catalysis, the BEP relationship was extended into a transition state scaling (TSS) relationship where the transition state (TS) energy ETS is correlated with the initial state (IS) energy EIS or the final state (FS) energy EFS, the energies being defined with respect to a given reference, generally the energy of the free molecule in gas phase (Figure 1).23,24 This TSS relationship has been shown to be an approximation of the original BEP formulation by Sutton and Vlachos.25 Both the original BEP and the TSS relations are extensively used in heterogeneous catalysis, especially for metalcatalyzed reactions.26 They have been more recently used for other classes of materials such as metal carbides, metal oxides, Received: February 19, 2015 Revised: April 30, 2015

A

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

and each set has to be sufficiently large to be meaningful. So, it is important to assess the impact of dividing a large set into several subsets going from a generic relation to more specific ones. For instance, we have recently shown that on a set of CH and OH bond breakings in simple alcohols on Rh(111), splitting the complete set into two subsets according to the nature of the bond leads to a gain in quality (MAE and maximum error, MAX, divided by two at least).33 In the same study, we also demonstrated that the reactivity of polyalcohols such as glycerol could be predicted via the linear relations obtained on the simple alcohols with a limited mean signed error (MSE), lower than ± 0.13 eV. In that study, we have extensively compared the various possible BEP and TSS relationships and shown that in the case of Rh(111) their performances are similar. This is consistent with other studies in the literature.25,31 In the work presented here, we extend the approach to the close-packed surfaces of a set of seven transition metals (Co, Ni, Ru, Rh, Pd, Ir, and Pt). Many works can be found in the literature concerning the dehydrogenation reaction on metallic catalysts for various simple alcohols,12,25,34−36 but the presented linear relations are in many cases established for one unique metal catalyst. Besides, even when several metallic surfaces are considered,31,37 those relations are not specific to alcohols and bring together various molecules of different nature (water, hydrocarbons, diatomics, etc.) losing quality of prediction. In the current paper, our aim is to attempt the design of linear energy relations valid for a wide range of metals for CH and OH scission in monoalcohols, aiming at a fast prediction of the reactivity of families of alcohols and polyalcohols extracted from biomass on metal catalysts. A unique global relation could be very attractive compared to several metal-dependent relations and would be easier to apply. However, the committed error has to be assessed, applying a careful statistical analysis of the linear relation quality. The effect of dividing our set into several subsets according to the metal or to the nature of the reaction will be addressed for both BEP and TSS relationships. We will show that the original BEP is more robust than the TSSs when several metals are considered in the same set, and we will discuss the underlying factors controlling accuracy with pure geometric considerations.

Figure 1. Typical energy profile for a dehydrogenation elementary step (X−H splitting) on a metallic slab represented as a gray rectangle. Each dehydrogenation elementary step as well as the energy of initial state (EIS), final state (EFS), and transition state (ETS) are defined relative to the fragments in gas phase and the slab. The activation energy (E⧧) and the reaction energy (ΔE) are defined relative to the adsorbed initial state.

or perovskites.27,28 In combination with scaling relations, they give access to large reaction networks at the price of a small one,7−11,15 or they accelerate the prediction to other metals.12−14 Indeed, from the adsorption energies of a subset of intermediates, a scaling relationship can be derived and used to extrapolate the adsorption energy of the full set under consideration. Similarly, from a subset of elementary steps, a BEP (resp. TSS) relationship ensures the subsequent prediction of the activation energy E⧧ (resp. transition state energy ETS) based only on the reaction energy ΔE (resp. on the FS or IS energy). The predicted data come with a certain statistical error that can be quantified by the mean average error (MAE).25,29 The control of this error is essential for the performance of the kinetic model that relies on these data. Because of the exponential nature of the Eyring equation, an error as small as 0.05 eV on the activation barrier leads to an error of at least 1 order of magnitude on the corresponding kinetic constant. With a good knowledge of the committed error, the model can be eventually refined with a more accurate evaluation of energies of the most sensitive steps.9,30 In this work, we focus on the BEP and TSS relationships. The quality of those relationships is linked to the similarity of the elementary steps that are included in the set under consideration. Searching for a universal relationship will be necessarily at the price of an increased mean average error. Bligaard and co-workers have collected the DFT data for 249 hydrogenation/dehydrogenation reactions over close-packed and open-facets of various transition metals and including a variety of CH, OH, and NH bond scissions.31 The universal TSS relation leads to an MAE of 0.28 eV. This error can be reduced when using a limited set; with the subset including only the successive water splitting into OH and then O, the MAE decreases to 0.20 eV. The impact of the homogeneity of the set of reactions is even more clear when the original BEP is considered; the MAE decreases also from 0.27 eV for the global set to 0.17 eV for the limited set of water splitting. It can be further decreased when the set is restricted to a given reaction on a given facet; the MAE for water splitting into OH on open facets is lower than 0.1 eV. Similar observations can be done for CC, CO, CN, NO, NN, and OO bond scissions catalyzed by metal facets; the MAE grows with the inhomogeneity of the set of reactions.32 However, developing specific relations means computing several sets to tackle the targeted reaction network,



COMPUTATIONAL DETAILS

DFT Calculations. We performed total energy calculations in the framework of periodic DFT using the VASP 4.36 code.38,39 The catalyst surfaces are represented by a four layer slab that exposes the most stable facet, meaning the closepacked surface (for Ni, Rh, Pd, Ir, and Pt: (111); for Co and Ru: (0001)). Each species was modeled adsorbed on the upper face of the slab, with a vacuum in the z direction equivalent to five layers. During the ionic relaxation, only the two uppermost layers were released, the lower ones being frozen in a bulk-like geometry. The coverage ratio was fixed at 1/9 ML for all the reactions that we considered using a p(3 × 3) cell. Besides, gas molecules relaxations were carried out in a 20 × 20 × 20 Å3 box. Concerning the computational parameters, we used the PW91 exchange-correlation functional40 with the PAW formalism41 to describe electron−ion interaction. All the calculations were performed with an energy cutoff of 400 eV for the planewave basis set, and a (3 × 3 × 1) k-mesh in the Monkhorst− Pack scheme.42 The convergence criterion was fixed to 10−6 eV for electronic energy optimization and to 10−2 eV/Å for ionic B

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

“range of errors”. Any point occurring out of this range of data is thus considered as an outlier, meaning nonstatistically representative. Apart from this relevant visual tool, we can define other quantitative parameters to study the performance of a statistical model for a sample of data with N points: The mean average error (MAE), the maximum error (MAX), and the mean signed error (MSE).

relaxation. Besides, for Ni and Co catalysts we performed open shell calculations in relation with their magnetic properties. Relating to the minimum energy path (MEP) optimization, for each elementary dissociation reaction we selected the most stable IS and FS, the energy of the FS being considered with the dissociated H on a separate slab. The TS was searched using the nudged elastic band (NEB) method,43,44 in combination with the dimer45 and quasi-Newton algorithms and confirmed by a vibrational analysis. Concerning the linear energy relationships that we present in this paper, we tested various exchange-correlation functionals, also including dispersion effects, and the correlation parameters and the error distributions were not significantly impacted. Similar conclusions were observed increasing the k-mesh density, the energy cutoff, the slab thickness, and the supercell area. For instance, the slopes (and similarly for the intercepts) of the BEP relationship for OH dissociations are only different by 0.02 units between 3 × 3 × 1 and 7 × 7 × 1 k-points calculation, for both Ni and Pd. Besides, between the two kmesh densities, MAE only changes by approximately 0.01 eV for those metals. Statistics. Because our objective is to establish linear models to efficiently predict kinetic properties of a chemical process, statistics are a central concept all along this work.46,47 All the statistical analysis were performed using the R software.48 In our case, the main issue is that linear energy relationships often deal with restricted sample size, typically less than 10 points for each catalyst for a given reaction (OH or CH breaking) in the current paper. In such conditions, usual criteria to assess the quality of linear regressions, especially the coefficient of determination R2, are not relevant. Indeed, this parameter is not robust enough and may be highly impacted by any minor change in the data, thus leading to misleading conclusions. As a result, we need to present other statistical tools that are appropriate to this situation. Let us consider the BEP framework, with a sample of N reactions. In this definition, activation energies {E⧧1 , E⧧2 ,..., E⧧n } are correlated with reaction energies {ΔE1, ΔE2,..., ΔEn}. Once a linear regression is plotted in this set of points, we defined the statistical errors εi as εi = Ei⧧ − Ei⧧

⎧ ∑ |ε | ⎪ MAE = i i N ⎪ ⎪ ⎨ MAX = max|εi| i ⎪ ⎪ ∑ εi ⎪ MSE = i ⎩ N

Although MAE and MAX are useful when a linear regression is extracted from one given set of points, MSE is especially relevant to assess the quality of a statistical model when it is used to perform predictions on a different sample. Hence, to discuss the quality of a BEP relationship established for instance in the case of CH dissociation on one given metal, MAE and MAX will be sufficient. However, if the linear energy relation is established for a set of various different metals, then MSE is necessary to assess the quality of the prediction on each metal taken individually.



RESULTS We considered a set of various monoalcohols undergoing either an OH or a CH bond dissociation on the close-packed surface of seven transition metals: Ru(0001), Co(0001), Rh(111), Ir(111), Ni(111), Pd(111), and Pt(111). Among all the possible CH bonds in an alcohol, we focus on the CH in α position to the OH group (denoted “CHα” in the following) because it has a specific behavior.33 Indeed, although the CH bonds in β or γ positions to the OH group have a reactivity that is similar to hydrocarbon CH bonds, CHα is very different. To be chemically representative, we selected reactions leading to products of various chemical natures, namely, radical species, carbonyl derivatives, and enols. Note that when chemisorbed on the metal surface species that are radicals in gas phase turn to closed-shell from the electronic coupling with the metal surface. The list of chosen reactions is shown in Figure 3. Altogether, our statistical set includes 56 C−H and 49 O−H elementary dissociation steps. As systematically described in an earlier study,33 there are various choices of BEP and TSS relations, depending on the direction chosen for the reaction (dissociation or formation of the bond) or for the TSS type depending on the energies that are correlated and on the energy reference that is chosen. In the following, we will focus on three relations, all defined for a set of elementary dissociation steps. The first one follows the original BEP definition and connects activation energies E⧧ and reaction energies ΔE for the C−H or O−H bond dissociation on the surface; we called it “BEP.diss”. The two others belong to the TSS framework. They correlate the TS energy ETS with either the IS energy EIS or the FS energy EFS, the reaction being considered in the dissociation direction, both with the dissociated FS in gas phase as energy reference (Figure 1). We called them “TSS.diss.IS/FS” and “TSS.diss.FS/FS”, respectively.

̂

E⧧i is the DFT calculated activation energy of the ith point, and ̂ Ei⧧ is the BEP-predicted activation energy. A convenient way to display the analysis of errors is by the means of a box-plot diagram. A typical one is given in Figure 2. In such a scheme, representative errors appear directly between the whiskers of the box. This interval is called “error span” or

Figure 2. Box-plot analysis of the errors for a linear relation. C

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

eV) for such a wide set of metals.31,49,50 The nice behavior of those BEP linear relationships is probably related to the high homogeneity of the set of elementary steps under consideration here, all related to simple aliphatic alcohols OH or CHα bond scission on dense metallic facets. The slope of the linear relation for the OH bond dissociation is small (0.11), whereas it is much larger for that of the CH bond (0.60). This point will be commented further later. For OH scission, all the points seem to fit well to the linear model whatever the metal, but a closer look at Pd and Ir data shows that the OH scission barriers on Pd are generally underestimated by the BEP line, whereas barriers in Ir are overestimated. In addition, some points seem singular for CH dissociation, mainly related to Pt or Pd. Those observations are better quantified with a representation of the distribution of the deviation between the data and the resulting linear relation. To start with, we keep the global linear relation of Figure 4, and analyze the error distribution metal-by-metal, as represented in the left panels of Figure 5. Let us focus first on the OH bond. The error distribution is markedly dependent on the metal, Co exhibiting the highest error span (from ∼−0.20 eV to ∼+0.20 eV), followed by Ni, Pd, and Rh (from ∼−0.15 eV to ∼+0.15 eV), and then other metals (Ru, Ir, and Pt) that present tight error distributions. Besides, although the MSE is low for almost every metal, Pd and Ir have opposite MSEs with a significant magnitude (+0.06 and −0.08 eV, respectively), in line with our initial observations. Thus, using the global linear relationship may induce some difficulties when comparing the reactivity of two different catalysts, especially Ir and Pd. In a second step, it can be interesting to see if the quality of the prediction using BEP relationships can be improved further when this prediction is performed via metal-dependent models instead of a global linear relation. We hence established these individual linear relations, using for every metal only the elementary steps on that metal. The corresponding errors are represented on the right panel of Figure 5. For almost every metal, error spans are similar to the global model predictions. The range of errors is slightly decreased only in the case of Pd. Obviously, the systematic deviation vanishes, especially for Ir and Pd. Hence, it is easier to compare the reactivity of two different catalysts. Nevertheless, the global model for OH scission is satisfying for each metal taken individually, and the committed error can be assessed precisely metal-per-metal. If we now consider CH bond dissociations, then we still observe opposite systematic deviations for Pd and Ir (+0.10 and

Figure 3. List of elementary reactions used to establish BEP-type linear relations for CH and OH scissions. Triangles, crosses, and squares refer to radical, enol, or carbonyl products, respectively, in relation with the plots of the linear energy relations provided later.

We will specifically focus on the possibility to build a single universal and accurate relation or in contrast on the need to split the data as a function of either the metal involved or along the nature of the product obtained after the dehydrogenation reaction. Ideally, a robust relation would be universal, and the subdivision of the data set along the nature of the metal and of the type of product would not be necessary. BEP Linear Relationships. We start with the classical BEP relation in its most common form, the BEP.diss definition. A clear and good statistical quality linear relation is observed for both OH and CH bond dissociation processes when taking all the elementary reactions and metals together (Figure 4). The distributions of errors are shown as box plots in the insets. Concerning the OH scission, the quality of the BEP relationship is good with an MAE of 0.08 eV and errors ranging approximately between −0.15 and +0.20 eV. Concerning the CH dissociation, errors range between −0.20 and +0.20 eV, which is slightly larger than for OH bond, but the quality of the relation is still good (MAE = 0.09 eV). Besides, note the tightness of the confidence interval on the slope (±0.03) for both bond breakings, attesting to the robustness of the predictive model. The MAEs that were observed both for OH and CH breaking are fully satisfying and smaller than some values available in the literature (0.15−0.20

Figure 4. BEP.diss relations for OH and CH scissions in monoalcohols on various metallic catalysts. We present here only the global predictive model common to all metals and its corresponding error distributions in the bottom right corner. The color of the symbol denotes the metal, and the type of symbol is related to the nature of the product formed in the considered elementary step (square for carbonyls, crosses for enols, and triangles for radicals). D

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 5. Error distributions for the prediction of OH (top) and CHα (bottom) dissociations with the BEP.diss relation on the chosen transition metals. On the left panels, we presented the error distributions related to each metal when the prediction is performed using one global predictive model common to all catalysts. On the right panels, we presented the error distributions when the prediction is performed using metal-dependent models. The crosses denote the mean signed error (MSE) obtained for each metal, and the numbers below each box refer to the mean absolute errors (MAE) related to each metal as well.

−0.06 eV) in the global model predictions, as for OH dissociation. The error distribution is again strongly dependent on the metal but with a different order. Co still presents the widest error distributions, followed by Pt, Pd, and then Ir and Rh, which have the tightest error span. Splitting the global model into single-metal ones results in a lowering of error magnitudes for almost every metal, especially for Co, Pd, and Ir. However, let us highlight that even if the prediction quality is slightly improved by using metal-dependent relations, the global model still yields highly satisfying results, with an MAE under 0.1 eV that is comparable with the intrinsic error of the DFT energy (0.15 eV). As we did for the various metallic catalysts, it is also possible to split the set of elementary steps according to the nature of the products that are formed. It was already shown that splitting the data according to the reactions might lead to a large lowering of statistical errors (from up to 0.30−0.40 eV to lower than 0.10 eV in some cases).25,31,35 It stems from Figure 3 that two kinds of products may appear from alcohol dehydrogenation, either radical or (closed-shell) unsaturated species. Figure 6 shows that the distribution of errors is not equivalent for the two types of products. This is especially the case for CH dissociation where radical products give a smaller error (MAE = 0.06 eV) compared to that of closed-shell unsaturated species (MAE = 0.13 eV). The situation is similar for OH with however a smaller effect, but a non-negligible systematic deviation (∼−0.10 eV) appears in the case of unsaturated species treated with the global model. In a second step, we can design two product-dependent models (one to predict radical formation and one to predict unsaturated species formation) and see if the quality of the prediction is improved comparatively to the

global model. As shown in Figure 6, the error distributions are not significantly broadened when switching from the reactiondependent models to the global model. Hence, constructing different linear relations for the two types of products does not bring any advantage, and a single linear relation can be kept. It should only be acknowledged that the two types of products are associated with different errors and that the predictions for reactions yielding radical products would be more accurate than those giving closed-shell unsaturated molecules. TSS Linear Relationships. In this second part, we consider two TSS relations, which are also widely used in literature: the TSS.diss.FS/FS and the TSS.diss.IS/FS. They can constitute an interesting alternative to the BEP approach in some cases51 and were shown to yield predictions of a high statistical quality for OH and CH bond scissions on a given metal and facet.11,25,34−36 The first observation arising from Figure 7 is that in our case it is difficult to find a global linear energy relation common to every metal leading to predictions of a satisfying quality. Only the TSS-diss.IS/FS relation in the case of OH scission exhibits a tight error distribution, whereas the TSS.diss.FS/FS relation for the same OH bond scissions leads to large errors. Both formulations for the CH scissions lead to large MAX errors. Concerning the OH bond scissions, the TSS.diss.FS/FS gives a rather poor correlation with an MAE of 0.30 eV and a large MAX of 0.68 eV (Figure 7). A closer look to the repartition of the data according to the nature of the metal in Figure 7, top left panel, shows that metal-dependent models should give better results. Indeed, the various slopes of the linear relations are rather similar between the metals, but it is the different offsets that prevent the establishment of an accurate global E

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

TSS-diss.FS/FS relation for OH dissociation (see the Supporting Information for the numerical values). This can be seen also in Figure 8 where the systematic deviation for each metal is represented. The MSE is small for Rh (0.08 eV) but large and positive for oxophilic metals such as Co (+0.44 eV) and large and negative for non-oxophilic metals such as Pt (−0.54 eV). Hence, the MAE is at least divided by 2, i.e., lower than 0.15 eV, for all metals when the set of elementary reactions is split according to the metal for the TSS-diss.FS/FS relationship. Conversely, the TSS-diss.IS/FS relation provides a good correlation for OH bond splitting on the complete set, and the separation according to the nature of the metal does not improve its quality (Figure 8, bottom panels). Regarding the CH dissociations, average errors (0.15 and 0.11 eV for TSS.diss.IS/FS and TSS.diss.FS/FS, respectively) are only slightly higher than the one observed with the BEP definition (0.09 eV). However, MAEs are especially high in the global TSS relations (∼0.50 eV for both TSSs vs 0.30 eV for the BEP). In addition, we analyzed also the committed error per metal in Figure 9, left panels. The MSE is rather large and varies a lot from one metal to another for both TSS relations (from −0.17 eV to +0.22 eV). Here again, the amplitude of the MSE seems to be related to the oxophilicity of the metal. For instance, it is positive for the more oxophilic metals (Co and Ru) and negative for the less oxophilic metals (Pt and Pd) in the case of TS.diss.IS/FS. It is hence difficult to extract and use one global TSS relation common to all metals in the case of CH bond dissociation. However, if we consider the metals separately, accurate TSS relations can be constructed for each of them (Figure 9, right panels). For the TSS.diss.FS/FS, the MAE is nicely controlled under 0.06 eV for Pd and Rh and higher but still acceptable for the other metals (around 0.12 eV). Concerning the TSS.diss.IS/FS, errors as low as 0.04 eV are found for Ru, Rh, and Ir, whereas the relation is somewhat less accurate on Pd (MAE = 0.12 eV).

Figure 6. Error distributions for the prediction of OH (top) and CHα (bottom) dissociations with the BEP.diss relation. On the left side, we presented the error distributions related to each kind of species when the prediction is performed using one global model common to both radicals and unsaturated species. On the right side, we presented the error distributions when the prediction is performed using productdependent models (radical vs unsaturated). The numbers above each boxplot refer to the mean absolute error (MAE) in eV.

Figure 7. TSS relations for OH and CH scissions in monoalcohols on various metallic catalysts. We present here only the global predictive model common to all metals and its corresponding error distributions in the bottom right corner. The color of the symbol denotes the metal, whereas the type of symbol is related to the nature of the product formed in the considered elementary step (square for carbonyls, crosses for enols, and triangles for radicals). F

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 8. Error distribution for OH scission in the TSS framework. On the left panels, we present the error distributions related to each metal when the prediction is performed using one global predictive model common to every catalyst. On the right panels, we present the error distributions when the prediction is performed using metal-dependent models. The crosses denote the mean signed error (MSE) obtained for each metal, and the numbers below each box refer to the mean absolute errors (MAE) related to each metal as well.

Figure 9. Error distribution for CH scission in the TSS framework. On the left panels, we present the error distributions related to each metal when the prediction is performed using one global predictive model common to every catalyst. On the right panels, we present the error distributions when the prediction is performed using metal-dependent models. The crosses denote the mean signed error (MSE) obtained for each metal, and the numbers below each box refer to the mean absolute errors (MAE) related to each metal as well.

or unsaturated species (carbonyls and enols). At first glance, the data are indeed clustered according to the products, especially

An alternative is to split the set of elementary steps in function of the nature of the products that are formed: radicals G

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 10. Error distributions for the prediction of OH (top) and CHα (bottom) dissociations for the two TSS relationships. On the left side of each box, we present the error distributions related to each kind of species when the prediction is performed using one global predictive model common to both radicals and unsaturated species. On the right side of each box, we present the error distributions when the prediction is performed using product-dependent models (radical vs unsaturated). The crosses denote the mean signed error (MSE) obtained for product, and the numbers below each box refer to the mean absolute errors (MAE) related to each product as well.

Pt (leading to a deviation of 0.96 eV on the TS energy). This tendency to reach larger errors for TSS than for BEP relations has already been observed by Bligaard and co-workers,31 considering a wide set of dense and stepped facets of 13 transition metals. For the successive OH bond scissions of water, they found an MAE of 0.15 eV versus 0.20 eV for BEP and TSS.diss.FS, respectively. Similarly, for the successive CH scission in methane, ethane and propane, they obtained an MAE of 0.21 versus 0.25 eV for BEP and TSS.diss.FS, respectively. The quality of the TSS relation strongly depends on the nature of the bond (OH vs CH) and on the stable state (IS vs FS) that is chosen. For instance, the TSS-diss.IS/FS yields low deviations (MAE = 0.07 eV), whereas the TSS.diss.FS/FS leads to strong deviations (MAE = 0.30 eV) on the same global set of OH dissociations (Figure 7). This can be explained by examining the adsorption energy of some typical species. Indeed, if we consider the TSS-diss.FS/FS definition, the TS energy (ETS) is correlated with the FS energy (EFS). In the case of OH breaking, various alkoxy radicals and carbonyl derivatives can be formed as a FS. However, although adsorption energies of carbonyls are close to each other between all the metals, they are very different in the case of alkoxy radicals, ranging approximately between −1.67 eV (Pt) and −2.93 eV (Co), with the strength of the interaction increasing with the metal oxophilicity as reported in Figure 11. In contrast, the TS adsorption energies are not significantly impacted when changing the nature of surface. In Figure 7, one can notice that the energy of the TS is insensitive on the metal when the OH scission leads to a radical alkoxy. When the dissociation leads to a CO double bond, the energy of the TS is slightly

the data for the OH bond scission for the TSS.diss.FS/FS (Figure 7, squares for carbonyls, crosses for enols, and triangles for radicals). However, the analysis of the errors depending on the obtained product does not show large signed errors (Figure 10); for instance, an MSE of −0.04 eV for radicals and +0.09 eV for unsaturated species (carbonyls and enols together) is found using TSS.diss.FS/FS for OH bond scission. Accordingly, the usage of specific models does not bring any improvement of the committed error.



DISCUSSION The BEP relations derived from the global sets of OH and CH bond dissociations are almost as good as those established for each specific metal with a low MAE of 0.09 eV and a reasonable MAX of 0.3 eV. On the contrary, the TSS relations generally have to be obtained from specific subsets, one per metallic surface. Under this condition, the MAE and the MAX error can be limited to 0.1 and 0.25 eV, respectively, for most of the cases. However, the TSS-diss.IS/FS yields reasonable deviations even when built on the global set when OH scissions are only considered. In that case, the MAE of 0.07 eV and the MAX error of 0.23 eV can compete with the results obtained via the original BEP relation. To sum up, the BEP in its original formulation is more robust and universal than the more recent and popular TSS formulation upon the extension to other transition metal surfaces. In addition, the use of the TSS formulation developed for a given metal to predict the reactivity of other metals appears as rather risky and may lead to systematic error larger than 0.3 eV; for instance, the energy of the TS for the OH bond scission on Co would be strongly underestimated on the basis of a TSS-diss.FS/FS established on H

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 12. Isoenergetic curves in function of two geometrical parameters s1 and s2 in the region of a TS. s1 could be for instance the C−metal distance and s2 the C−H distance for a CH bond scission. The reaction coordinate is represented with a dashed curved line. The straight arrows correspond to DIS−TS and DTS−FS, the Euclidian distances computed between, on one hand, IS and TS and, on the other hand, TS and FS.

Figure 11. Adsorption energies (in eV) as a function of the metal. Six typical species involved in the set of elementary steps represented in Figure 1 are shown: isopropanol (iPrOH, □), acetone (×), hydrogen atom (○), enol (CH2CH−OH, *), the hydroxyl-isopropyl (MeCOHMe, △), and the isopropoxy (+) radicals.

modulated by the metal, being more stabilized on the less oxophilic metal such as Pt than on Co. In fact, the TS does not behave as the FS but rather as the IS. Indeed, when the FS is an alkoxy, the IS is an alcohol, and its energy follows the trend of the corresponding TS because it does not vary much with the nature of the metal (Figure 11). Similarly, when the FS is an aldehyde or a ketone, the IS for the OH scission is a hydroxyalkyl, which is more stabilized on Pt or Pd than on Ru and Co as is the corresponding TS (Figures 7 and 11). In a nutshell, for OH bond scission in alcohols at metallic compact surfaces, the TS is closer to the IS than the FS and that is why the global TSS.diss.FS/FS definition gives high errors, whereas TSS.diss.IS/FS definition leads to a relation of a high statistical quality. This is also why splitting the data into metal-specific subsets yields much more reasonable TSS.diss.FS/FS linear relationships for OH bond scission. The earliness of the OH scission TS is in agreement with the slope α close to zero (0.11) in the BEP correlation.26 Concerning the CH bond scission, the quality of the TSS relationship does not depend on the choice between the IS and the FS (see Figure 7). According to the slope of α = 0.60 of the corresponding BEP linear relationship, the TS is intermediate to late.26 This finding that TS related to OH dissociations are earlier than those related to CH breaking is confirmed by various papers33,35 because the transfer coefficient is generally lower for OH than for CH. In the following, we will present a way to correlate those observations on the basis of energetics with some geometrical considerations. There have been few attempts to determine the early or late character of a TS. Parameters to qualify the proximity between TS and reactants were built on the basis of the first-order density matrices or on atomic positions.52,53 We propose here a simple parameter ρ(%), exclusively based on geometrical considerations, to quantify the earliness (or lateness) of the TS. ρ is defined as

them, we determined the FS with the dissociated hydrogen in the same cell to evaluate ρ. According to Figure 13, it is clear that the nature of the TS shifts globally from early to late when the reaction energy moves from high exothermicity to high endothermicity, as expected. A good-quality linear relation between ρ and ΔE is shown for OH bond-cleavage steps, whereas for CH, the relation is much more qualitative. TSs related to the OH bond are globally early or intermediate, in agreement with the BEP transfer coefficient close to zero (α = 0.11), whereas those for the CH bond dissociation are distributed from early to late, explaining the intermediate value of the BEP transfer coefficient for this set of elementary steps (α = 0.60). The better quality of the linear relationships for the OH scission compared with the CH scission probably lies in the narrower range of ρ for the OH scission.



CONCLUSIONS In this work, we showed that it is possible to design a global BEP relation, aiming at predicting activation energies for monoalcohol dehydrogenation on a wide set of metallic catalysts. Such a relation leads to results of a high statistical quality with an MAE lower than 0.10 eV both for CH and OH breaking. We showed that the quality of the prediction depends on not only the nature of the metal but also the chemical reaction that is performed. The largest errors are obtained in the case of Co, Pd, and Pt metal surfaces and for unsaturated species formation. Although BEP relationships are easily extendable to metals of different nature, it is not obvious to design a global TSS relation valid for such a wide range of catalysts. Accurate TSS relations can only be obtained specifically for each metal separately. We evidenced that the origin of this phenomenon is related both to the discrepancies of adsorption energies of some key chemical species between the various surfaces and to the different nature of the TS for the OH (early) and the CH (intermediate) bond dissociation, which makes difficult the choice of the appropriate thermodynamic state for the TSS relation. As a result, the BEP-type linear relations are shown to offer the best compromise between accuracy and universality for the fast exploration of complex reaction systems. Relations in this work were established for the dense surfaces ((111) facet for fcc metals and (0001) facet for hcp ones), and it can be expected that different slope and offset parameters for BEP-type relationships

ρ(%) = 100 × D(IS − TS)/(D(IS − TS) + D(FS − TS))

where D(IS−TS) is the Euclidian distance (in Å) between the IS and the TS (and, respectively, FS for D(FS−TS); Figure 12). Hence, ρ is between 0 and 100%. We propose to classify the TS as early if ρ is between 0 and 40%, as late if ρ is between 60 and 100%, and as intermediate if ρ is between 40 and 60% (see Figure 13). In our data set of reactions and metals, we selected some reactions for each metal for both CH and OH scissions, leading to products of different chemical nature. For each of I

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 13. The ρ parameter in function of the reaction energy ΔE for selected CH and OH bond scissions. In general, early TS (ρ < 40%) are found for highly exothermic reactions and late TS (ρ > 60%) for highly endothermic reactions. (4) Kang, H. C.; Weinberg, W. H. Modeling the kinetics of heterogeneous catalysis. Chem. Rev. 1995, 95, 667. (5) Stamatakis, M.; Vlachos, D. G. Unraveling the complexity of catalytic reactions via kinetic Monte Carlo simulation: current status and frontiers. ACS Catal. 2012, 2, 2648. (6) Besson, M.; Gallezot, P.; Pinel, C. Conversion of biomass into chemicals over metal catalysts. Chem. Rev. 2014, 114, 1827. (7) Chen, Y.; Salciccioli, M.; Vlachos, D. G. An efficient reaction pathway search method applied to the decomposition of glycerol on platinum. J. Phys. Chem. C 2011, 115, 18707. (8) Guo, N.; Caratzoulas, S.; Doren, D. J.; Sandler, S. I.; Vlachos, D. G. A perspective on the modeling of biomass processing. Energy Environ. Sci. 2012, 5, 6703. (9) Salciccioli, M.; Stamatakis, M.; Caratzoulas, S.; Vlachos, D. G. A review of multiscale modeling of metal-catalyzed reactions: mechanism development for complexity and emergent behavior. Chem. Eng. Sci. 2011, 66, 4319. (10) Salciccioli, M.; Vlachos, D. G. Kinetic modeling of pt catalyzed and computation-driven catalyst discovery for ethylene glycol decomposition. ACS Catal. 2011, 1, 1246. (11) Liu, B.; Greeley, J. Decomposition Pathways of Glycerol via C− H, O−H, and C−C Bond Scission on Pt(111): A Density Functional Theory Study. J. Phys. Chem. C 2011, 115, 19702. (12) Liu, B.; Greeley, J. A density functional theory analysis of trends in glycerol decomposition on close-packed transition metal surfaces. Phys. Chem. Chem. Phys. 2013, 15, 6475. (13) Yoo, J. S.; Abild-Pedersen, F.; Nørskov, J. K.; Studt, F. Theoretical Analysis of Transition-Metal Catalysts for Formic Acid Decomposition. ACS Catal. 2014, 4, 1226. (14) Ferrin, P.; Simonetti, D.; Kandoi, S.; Kunkes, E.; Dumesic, J. A.; Nørskov, J. K.; Mavrikaki, M. Modeling ethanol decomposition on transition metals: a combined application of scaling and Brønsted− Evans−Polanyi relations. J. Am. Chem. Soc. 2009, 131, 5809. (15) Rangarajan, S.; Brydon, R. R. O.; Bhan, A.; Daoutidis, P. Automated identification of energetically feasible mechanisms of complex reaction networks in heterogeneous catalysis: application to glycerol conversion on transition metals. Green Chem. 2014, 16, 813. (16) Abild-Pedersen, F.; Greeley, J.; Studt, F.; Rossmeisl, J.; Munter, T. R.; Moses, P. G.; Skúlason, E.; Bligaard, T.; Nørskov, J. K. Scaling properties of adsorption energies for hydrogen-containing molecules on transition-metal surfaces. Phys. Rev. Lett. 2007, 99, 016105. (17) Montemore, M. M.; Medlin, J. W. Scaling relations between adsorption energies for computational screening and design of catalysts. Catal. Sci. Technol. 2014, 4, 3748. (18) Brønsted, J. N. Acid and Basic Catalysis. Chem. Rev. 1928, 5, 231. (19) Bell, R. P. The theory of reactions involving proton transfers. Proc. R. Soc. London, Ser. A 1936, 154, 414. (20) Evans, M. G.; Polanyi, M. On the introduction of thermodynamic variables into reaction kinetics. Trans. Faraday Soc. 1937, 33, 448. (21) Evans, M. G.; Polanyi, M. Inertia and driving force of chemical reactions. Trans. Faraday Soc. 1938, 34, 11.

may be found for other facets but still valid within a large set of metals, as already shown for some small molecules.16,23,36,50 Even if such relationships were designed for monofunctional molecules, we suggest applying them to more complex systems. We already demonstrated the validity of this concept for glycerol dehydrogenation on Rh(111).33



ASSOCIATED CONTENT

S Supporting Information *

Correlation parameters of the three kinds of relationships discussed here, for the metal dependent and reaction dependent predictive model, and the complete list of structures and energies of the various species used in this article to established BEP-type relationships. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b01703.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

J.Z.: Department of Materials Science and Engineering, Technion − Israel Institute of Technology, Haifa, Israel. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the ANR through the GALAC Project (ANR-10-CD2I-011) and was performed using HPC resources from GENCI-CINES/IDRIS (grant 2013080609) and PSMN at ENS Lyon.

■ ■

ABBREVIATIONS DFT, density functional theory; MAE, mean average error; MAX, maximum error REFERENCES

(1) Nørskov, J. K.; Bligaard, T.; Rossmeisl, J.; Christensen, C. H. Towards the computational design of solid catalysts. Nat. Chem. 2009, 1, 37. (2) Raimondeau, S.; Vlachos, D. G. Recent developments on multiscale, hierarchical modeling of chemical reactors. Chem. Eng. J. 2002, 90, 3. (3) Salciccioli, M.; Stamatakis, M.; Caratzoulas, S.; Vlachos, D. G. A review of multiscale modeling of metal-catalyzed reactions: Mechanism development for complexity and emergent behavior. Chem. Eng. Sci. 2011, 66, 4319. J

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (22) Shorter, J. Correlation Analysis in Organic Chemistry: An Introduction to Linear Free Energy Relationships; Clarendon Press: Oxford, U.K., 1973. (23) Nørskov, J.; Bligaard, T.; Logadottir, A.; Bahn, S.; Hansen, L.; Bollinger, M.; Bengaard, H.; Hammer, B.; Sljivancanin, Z.; Mavrikakis, M.; et al. Universality in heterogeneous catalysis. J. Catal. 2002, 209, 275. (24) Michaelides, A.; Liu, Z. P.; Zhang, C. J.; Alavi, A.; King, D. A.; Hu, P. Identification of General Linear Relationships between Activation Energies and Enthalpy Changes for Dissociation Reactions at Surfaces. J. Am. Chem. Soc. 2003, 125, 3704. (25) Sutton, J. E.; Vlachos, D. G. A Theoretical and Computational Analysis of Linear Free Energy Relations for the Estimation of Activation Energies. ACS Catal. 2012, 2, 1624. (26) Van Santen, R. A.; Neurock, M.; Shetty, S. G. Reactivity theory of transition-metal surfaces: a Brønsted−Evans−Polanyi linear activation energy−free-energy analysis. Chem. Rev. 2010, 110, 2005. (27) Li, H. Y.; Guo, Y. L.; Guo, Y.; Lu, G. Z.; Hu, P. C-H bond activation over metal oxides: a new insight into the dissociation kinetics from density functional theory. J. Chem. Phys. 2008, 128, 051101. (28) Vines, F.; Vojvodic, A.; Abild-Pedersen, F.; Illas, F. BrønstedEvans-Polanyi Relationship for Transition Metal Carbide and Transition Metal Oxide Surfaces. J. Phys. Chem. C 2013, 117, 4168. (29) Sutton, J. E.; Vlachos, D. G. Error estimates in semi-empirical estimation methods of surface reactions. J. Catal. 2013, 297, 202. (30) Stegelmann, C.; Andreasen, A.; Campbell, C. T. The degree of rate control: how much the energies of intermediates and transition states control rates. J. Am. Chem. Soc. 2009, 131, 8077. (31) Wang, S.; Petzold, V.; Tripkovic, V.; Kleis, J.; Howalt, J. G.; Skulason, E.; Fernandez, E. M.; Hvolbaek, B.; Jones, G.; Toftelund, A.; et al. Universal transition state scaling relations for (de)hydrogenation over transition metals. Phys. Chem. Chem. Phys. 2011, 13, 20760. (32) Wang, S.; Temel, B.; Shen, J.; Jones, G.; Grabow, L. C.; Studt, F.; Bligaard, T.; Abild-Pedersen, F.; Christensen, C. H.; Nørskov, J. K. Universal Brønsted-Evans-Polanyi relations for C−C, C−O, C−N, N− O, N−N, and O−O dissociation reactions. Catal. Lett. 2011, 141, 370. (33) Zaffran, J.; Michel, C.; Auneau, F.; Delbecq, F.; Sautet, P. Linear energy relations as predictive tools for polyalcohol catalytic reactivity. ACS Catal. 2014, 4, 464. (34) Greeley, J.; Mavrikakis, M. Competitive paths for methanol decomposition on Pt(111). J. Am. Chem. Soc. 2004, 126, 3910. (35) Wang, S.; Vorotnikov, V.; Sutton, J. E.; Vlachos, D. G. BrønstedEvans-Polanyi and transition state scaling relations of furan derivatives on Pd(111) and their relation to those of small molecule. ACS Catal. 2014, 4, 604. (36) Wang, H. F.; Liu, Z. P. Comprehensive mechanism and structure-sensitivity of ethanol oxidation on platinum: new transitionstate searching method for resolving the complex reaction network. J. Am. Chem. Soc. 2008, 130, 10996. (37) Michaelides, A.; Liu, Z. P.; Zhang, C. J.; Alavi, A.; King, D. A.; Hu, P. Identification of general linear relationships between activation energies and enthalpy changes for dissociation reactions at surfaces. J. Am. Chem. Soc. 2003, 125, 3704. (38) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B 1993, 47, 558. (39) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169. (40) Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B 1992, 45, 13244. (41) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 1999, 59, 1758. (42) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188. (43) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 2000, 113, 9901.

(44) Henkelman, G.; Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 2000, 113, 9978. (45) Heyden, A.; Bell, A. T.; Keil, F. J. Efficient methods for finding transition states in chemical reactions: Comparison of improved dimer method and partitioned rational function optimization method. J. Chem. Phys. 2005, 123, 224101. (46) Spatz, C. Basic Statistics: Tales of Distributions, 10th ed.; Wadsworth Cengage Learning: Belmont, CA, 2010. (47) Utts, J. M.; Heckard, R. F. Statistical Ideas and Methods; Thomson Brooks/Cole: Belmont, CA, 2006. (48) Crawley, M. J. The R book, 2nd ed.; John Willey & Sons: West Sussex, 2013. (49) Michel, C.; Goeltl, F.; Sautet, P. Early stages of water/hydroxyl phase generation at transition metal surfaces - synergetic adsorption and O-H bond dissociation assistance. Phys. Chem. Chem. Phys. 2012, 14, 15286. (50) Fajin, J. L.; Cordeiro, M. N. D.; Illas, F.; Gomes, J. R. Descriptors controlling the catalytic activity of metallic surfaces toward water splitting. J. Catal. 2010, 276, 92. (51) Loffreda, D.; Delbecq, F.; Vigné, F.; Sautet, P. Fast prediction of selectivity in heterogeneous catalysis from extended Brønsted−Evans− Polanyi relations: a theoretical insight. Angew. Chem., Int. Ed. 2009, 48, 8978. (52) Ciolowsky, J. Quantifying the Hammond postulate: intramolecular proton transfer in substituted hydrogen catecholate anions. J. Am. Chem. Soc. 1991, 113, 6756. (53) Manz, T. A.; Sholl, D. S. A dimensionless reaction coordinate for quantifying the lateness of transition states. J. Comput. Chem. 2009, 31, 1528.

K

DOI: 10.1021/acs.jpcc.5b01703 J. Phys. Chem. C XXXX, XXX, XXX−XXX