In silico search for novel methane steam reforming catalysts

This paper demonstrates a method for screening transition metal and metal alloy catalysts based on their predicted rates and stabilities for a given catalytic reaction. This method involves combining reaction and activation energies (available to the public via a web-based application ‘CatApp’) with a microkinetic modeling technique to predict the rates and selectivities of a prospective material. This paper illustrates this screening technique using the steam reforming of methane to carbon monoxide and hydrogen as a test reaction. While catalysts are already commercially available for this process, the method demonstrated in this paper is very general and could be applied to a wide range of catalytic reactions. Following the steps outlined herein, such an analysis could potentially enable researchers to understand reaction mechanisms on a fundamental level and, on this basis, develop leads for new metal alloy catalysts.

of catalytic reactions. Following the steps outlined herein, such an analysis could potentially enable researchers to understand reaction mechanisms on a fundamental level and, on this basis, develop leads for new metal alloy catalysts.

Introduction
One of the ultimate goals in catalysis research is to help accelerate the search for novel active and selective catalysts for heterogeneous reactions [1]. Central to achieving this goal is a molecular level understanding of the interactions between the surface of a catalyst and the intermediates and transition states of a given reaction. Recent advances in electronic structure theory have made it possible to carry out reactions 'in silico' on large numbers of heterogeneous catalysts, and then screen these to find leads for new materials with high predicted rate, selectivity and stability. This method involves finding correlations between the energetics of the relevant reaction barriers and adsorption energies, finding appropriate descriptors for these correlations, and then utilizing these descriptors to simulate reaction kinetics in the reduced-dimensional descriptor space [2][3][4][5]. This approach has already assisted researchers in finding promising catalyst leads for heterogeneously catalyzed reactions [2,[6][7][8][9] and lends hope that a mapping of the catalyst genome could lead to general catalyst design strategies, where atomic-scale simulations would play a key role in an integrated approach to catalyst search [10].
In the present study we use the steam reforming of methane as a test reaction for demonstrating the utility of this strategy. CH 4 steam reforming is one of the most commonly used processes for the inexpensive production of H 2 and synthesis gas [11]. A number of transition metals have been reported to be active for this reaction. The highest CH 4 steam reforming rates have been observed for Rh and Ru-based catalysts, with Ni, Pt, Ir and Pd also showing high rates [12][13][14][15][16]. Commercial catalysts are typically formulated using Ni as the active phase, since it is the least expensive of these metals [17]. The Ni-catalyst is also very active for the steam and dry reforming of other reactants [18,19]. High reactor temperatures are required to achieve a high CH 4 conversion (from 750 to 1000 • C), because the reaction is strongly endothermic (equation (1)). Furthermore, high steam-to-carbon ratios (H 2 O/CH 4 > 3) are often needed to prevent catalyst deactivation through the deposition of graphitic carbon, which increases the capital and operating costs for H 2 generation: While the steam reforming of CH 4 using Ni-based catalysts is a relatively mature technology, important challenges still remain and are the subject of active research. Among these are the areas of improving catalyst rates while increasing the resistance to catalyst deactivation via sulfur poisoning, carbon formation and sintering [17,20]. One strategy in particular that has been investigated for increasing catalyst resistance to deactivation is through alloying the active metal (e.g. Ni), with another metal (e.g. Sn [21], Au [22,23], Ag [24], Ru [25,26]). Typically these experiments have resulted in an improvement in the long-term stability of the catalyst, albeit with a sacrifice in catalytic reaction rates.
In order to develop materials with high rates and stabilities under reaction conditions, it is essential to develop a better understanding of the steam reforming reaction mechanism. Previous studies, both experimental and theoretical, have proposed a range of reaction pathways for CH 4 steam reforming. Early theoretical work by Bengaard et al, for example, indicated that the rate-determining step for CH 4 steam reforming on Ni(111) and Ni(211) surfaces involves the breaking of the initial C-H bond of CH 4 [27]. They also reported that the CO formation proceeds through a direct C-O combination pathway. Other researchers suggest that CH 4 reforming on a Ni(111) surface, however, may involve a hydrogen insertion pathway involving CHO * and COH * , as opposed to the direct C-O formation pathway [28][29][30]. Experiments by Wei and Iglesia show that at high temperatures (823-1023 K) the rate-limiting step in both steam and dry reforming of CH 4 is the activation of the first C-H bond [31]. At somewhat lower temperatures (below 773 K), however, it has been suggested that the C-O coupling competes for the ratedetermining step [16]. Moreover, a study by Bhattacharjee et al is consistent with experimental results suggesting that on a Rh(211) surface at high temperatures the formation of CO proceeds through the direct C-O combination pathway. At lower temperatures, however, a hydrogen insertion mechanisms involving COH * may be favorable [32]. Shetty et al studied both the direct and H-assisted pathway for CO formation on Ru(1121), showing that the direct CO formation pathway has a lower overall barrier [33]. Chen et al studied the steam reforming pathways on Pt(111) and Pt(211), indicating that COH pathway is important for CO production [34]. Research by Inderwildi et al suggests that CHO * pathway is thermodynamically more favorable on Ru(0001), Pt(111), and Pd(111) surfaces [35].
As described elsewhere, microkinetic modeling calculations can incorporate multiple reaction routes, so it is unnecessary to know beforehand which pathways may become dominant under a given set of conditions [4,36]. By solving the model we can obtain reaction rates for each route and thus predict which elementary steps may become rate determining for different surfaces and under different reaction conditions. In a number of studies this method has been used to gain insight into heterogeneously catalyzed reactions [3,4,[37][38][39]. Herein we consider three potential CO formation pathways: direct formation via C * and O * coupling, as well as hydrogen insertion pathways involving CHO * or COH * . In principle many different parameters may affect the reaction rates and selectivities of a given catalyst. Nevertheless, it has been reported that correlations between various adsorption [40][41][42] and transition-state energies [27,[43][44][45] offer the possibility to reduce the number of these parameters to one or two [3,4,41]. In this paper, we coupled the microkinetic model with a descriptor-based approach in order to predict reaction rates in terms of just two parameters (i.e. the adsorption energies of atomic carbon and oxygen) [16]. Using this model, we screened a large number of transition metal alloys based on their predicted rates and stabilities for CH 4 steam reforming.

Methods
In order to develop catalyst leads for the steam reforming of CH 4 , a three-step approach was employed. Firstly, a microkinetic model of the elementary reactions considered was developed. Secondly, energies of reaction intermediates were scaled with descriptors and incorporated to the microkinetic model. Thirdly, using the results of this model, a large number of transition metal alloys were screened based on their predicted rates and stabilities under reaction conditions.

Microkinetic model method
This microkinetic model consisted of the following reaction steps for CH 4 and H 2 O dehydrogenation, as well as several possible pathways for CO formation via C * , CHO * and COH * intermediates: The parameterization of the microkinetic model was performed as described in literature [4,8,16], and can be found in section 2 of the supplementary material (available from stacks.iop.org/NJP/15/125021/mmedia).
In the microkinetic model, adsorbates may be adsorbed on one of four different sites: the 'step' ( * s) site corresponds to the upper part of a (211) step site. The 'four-fold' ( * f) site corresponds to the lower four-fold site of a (211) step. The 'terrace' ( * t) corresponds to a (211)  [16,48] shown to give a reasonable correlation with experiments in a number of previous studies. The energetics of the intermediates of SRM system are not strongly structure sensitive, thus we think this simplifying restriction gives a reasonable framework for the catalyst screening process. For more structure sensitive reactions, the impact of different facets and an estimation of their propensity should be considered for the screening to be reliable.
Interactions between adsorbates were also considered in the microkinetic model. Adsorption energies are described as a function of coverage using the model reported elsewhere [1,36,46,47]. Adsorption energies for the individual intermediates are calculated at 1 4 , 1 2 , 3 4 and 1 ML coverages for the lowest-energy adsorption site at low ( 1 4 ML) coverage. A brief summary of the interaction model, including the interaction parameters used in the microkinetic model, which can be found in section 3 of the supplementary material (available from stacks.iop.org/NJP/15/125021/mmedia).
Steady-state numerical solutions to the microkinetic model were calculated for one laboratory scale reactor conditions studied in the earlier work [16] and two model industrial reaction conditions, which were designed to simulate the inlet and outlet conditions of a typical steam reforming reactor designed to generate H 2 for ammonia synthesis [48]. The temperatures, pressures and gas compositions of these three conditions are described in table 1.

Scaling method
As described in previous studies, binding energies of reaction intermediates can be scaled with the binding energies of carbon (E C ) and oxygen (E O ) [40][41][42]. Similarly, the transition state energies for the various reaction steps can be scaled with the energies of the reaction products [4,[43][44][45]49]. Density functional theory (DFT)-calculated energies for the various reaction intermediates and transition states were accessed from CatApp, a web application that stores energetic information for elementary coupling reaction on transition metal surfaces [50], 7 and can be found in section 1 of the supplementary material (available from stacks.iop.org/NJP/15/125021/mmedia), along with the fitting parameters of all the scaling relations. Energies are given relative to the gas-phase energies of CH 4 , H 2 O and H 2 .
The DFT-calculated adsorption energies of surface intermediates and transition states for on the stepped (211) surfaces of Ag, Cu, Pd, Pt and Rh were accessed via CatApp [50,51]. Transition state energies for reactions (9) and (11) were calculated using the fixed bond length method in this work. Self-consistent, periodic DFT calculations for the surface intermediates and transition states were performed using the Dacapo plane wave code 8 , which represents the ionic cores using Vanderbilt ultrasoft pseudopotentials. Calculations were performed using the revised Perdew-Burke-Ernzerhof (RPBE) exchange-correlation functional [52], which uses a generalized gradient approximation, with a kinetic energy cutoff of 340 eV and a density cutoff of 500 eV. Self-consistent electron densities were determined by iterative diagonalization of the Kohn-Sham Hamiltonian, where the occupation of the Kohn-Sham states were smeared according to a Fermi-Dirac distribution with a smearing factor of k B T = 0.1 eV [53]. Energies were extrapolated to k B T = 0 eV. Stepped (211) surfaces were represented with nine-layered slabs (which correspond to three layers in the (111) direction) with (1 × 3) unit cells. All adsorbates and the top three layers were allowed to relax, while the remaining layers were kept fixed in their bulk-truncated positions. Brillouin zones were sampled using Monkhorst-Pack k-point meshes of 3 × 3 × 1. Successive slabs were separated with a vacuum of 15 Å.
Zero-point energies, entropies and internal energies of adsorbed intermediates were included due to the relatively high reaction temperatures. These values were calculated using the harmonic oscillator approximation from their vibrational frequencies on the Rh(211) surface (table S.5 of the supplementary material), and were assumed to be constant on all surfaces.

Alloy screening method
Based on this microkinetic model, a large number of transition metal alloys were screened based on their predicted rates, costs and stabilities for CH 4 steam reforming. Specifically, we thus imposed the following screening criteria for A 3 B and AB transition metal alloys, as depicted in figure 1: (1) The candidate bulk alloy should be stable or close to stable with respect to its bulk pure metal constituents. (2) The candidate bulk alloy should be stable or almost stable against oxidation at the given reaction conditions. (3) The candidate should not contain expensive metal constituents.
(4) The candidate should be reasonably active under the given reaction conditions.

2.3.1.
Alloy thermodynamic stability filter. The ability of two elemental metals to form a stable intermetallic compound is described by the alloy's formation energy, E f . As a simple model for the stability of the binary intermetallic compounds, the formation energies of the binary alloys of the form A 3 B and AB in simple fcc (face-centered cubic)-like crystal structures were calculated using DFT. The formation energy of these simple alloys is then used to predict the ability of the metals to form stable intermetallic compounds with each other. The formation energy of a given alloy was defined as where E f is the formation energy, E alloy is the DFT-calculated energy of the metal alloy, and N i and E i(bulk) are the molar ratios and the DFT-calculated energies, respectively, of the individual metal, i. According to this definition, the more negative the value, the more stable the alloy. For the purposes of the present analysis, only alloys with formation energy less than +0.2 eV per unit cell were considered to be potentially stable. A complete list of the investigated alloys, as well as their formation energies, is given in table S.8 of the supplementary material (available from stacks.iop.org/NJP/15/125021/mmedia).

Oxidation stability filter.
In order to mitigate catalyst deactivation due to coking and achieve an improved high CH 4 conversion, a high steam-to-carbon ratio is used under industrial conditions. With high H 2 O partial pressure, however, metals and metal alloys may deactivate via the formation of an inactive bulk metal oxide. Candidate alloys, therefore, must be evaluated in terms of their stability with respect to oxidation. In order to estimate the stability of a given candidate alloy (A x B y ) against oxidation, we calculated the Gibbs free energy of oxidation for each of the constituent metal atoms under reaction conditions. Then we compared those energies with the formation energy of the alloy as a whole using the equation where E f is the energy of formation, G A-oxide and G B-oxide are the energies of oxidation of metal A and B, respectively, and G alloy oxide is the energy of oxidation of the alloy. The energy of oxidation for metal M, which has a ground state oxide M m O n , was calculated using the equation where H 0 reaction and S 0 reaction are calculated using reference values from literature [54]. A detailed description of the method used here, as well as several example calculations, are given in section 4 of the supplementary material.

Cost and activity filters.
In order to be relevant on an industrial scale, it is often necessary to limit the use of expensive noble metals. Therefore, we used a simple interpolation method to estimate the cost of the various transition metal alloys, assuming that the cost is independent of the preparation method but depends solely on the cost of the constituted metals. In order to exclude alloys containing noble metals (e.g. Pt, Pd, Ir, Ru, Rh, Os and Au), we considered that alloys with predicted costs larger than $5000 kg −1 should be excluded. Reference prices of pure metals were taken from chemicool 9 , and can be found in table S.7 of the supplementary material (available from stacks.iop.org/NJP/15/125021/mmedia). Alloys were also screened based on their predicted rates for CH 4 steam reforming. In order to concentrate on those alloys that have the highest predicted rates, we screened them based on turnover frequencies (TOFs) larger then 10 −7 and 10 −2 s −1 for the inlet and outlet reactor conditions, respectively. We note that the activity threshold for inlet condition is very low. However, instead of focusing on the TOFs of the metals, our purpose here is to investigate whether the activity trend changes with different reaction conditions, thus the threshold is decided to identify the most active metals in that specific reaction condition. The predicted costs and rates of all thermodynamically stable alloys are given in table S.8 of the supplementary material.

Microkinetic model
The reaction energetics for CH 4 steam reforming on a stepped Rh(211) surface are given in figure 2. Rh is chosen because it is considered as one of the most active and coke-resistant metal in SRM reaction. Potential energy diagrams for less reactive metals (e.g. Pt, Cu) can be found in figure S4 in the supplementary material. Three reaction paths for the formation of CO * are depicted-the direct C-O formation pathway, as well as two hydrogen insertion pathways involving CHO * and COH * . The highest barriers for steam reforming on a Rh(211) surface (figure 2) are those for the formation of CO * . This is consistent with earlier work that suggested the formation of C-O as the rate-limiting step at low temperatures [16]. Similarly, it has been suggested that, while the direct formation of CO * from C * and O * is the main pathway, CO * formation from COH * and CHO * can also be significant at lower temperatures [55]. From the free energy pathway shown in figure 2(b), however, it is difficult to elucidate which CO formation pathway will be the dominant.
Combining scaling relations with microkinetic model is one strategy for developing a more general picture of catalyst reactivity in terms of theoretical activity volcanoes. For example, figure 3 shows the predicted CH 4 steam reforming rates as a function of E C and E O , as calculated by microkinetic modeling, with and without adsorbate-adsorbate interactions, for reaction conditions similar to a laboratory scale CH 4 steam reforming reactor. The highest rates are predicted for transition metals such as Rh, Ru, Ni and Ir, while lower rates are predicted for (a) Potential energy diagram for CH 4 steam reforming to CO (g) and H 2(g) on a Rh(211) surface; (b) free energy diagram for CH 4 steam reforming to CO (g) and H 2(g) on a Rh(211) surface. Reaction conditions are: T = 773 K, P = 1 bar, with a gas composition of 40% CH 4 , 40% H 2 O, 5% CO and 15% H 2 (corresponding to 29% approach to equilibrium). Three pathways are depicted in red, black, and blue for the formation of CO (g) via C * + O * , C * + OH * and CH * + O * , respectively. The main pathway is depicted in red. Values for the adsorption energies of reaction intermediates were taken from CatApp [50,56]. Pd and Pt. These results are consistent with previous experimental work under same reaction conditions, showing that catalysts involving Rh, Ru, Ni, Ir, Pd and Pt are active for steam reforming [13,16,19]. These results are also very similar to the previously volcano plots [16], where the dissociative adsorption of CH 4 and the direct formation of CO * from C * and O * are assumed to be competing as rate determining steps.
The overall shape of the activity volcano is better understood in terms of the coverages of C * , O * and CO * , as depicted in figure 4. Note that adsorbed carbon (C * ) bonds on the four-fold sites of the stepped surface, while CO * and O * are bonded more strongly on the on-top steps sites [47]. Based on a comparison of figures 3 and 4, it is evident that the most active transition metals have coverages of these strongly bonded intermediates somewhat in the area of 0.5 ML on the on-top and four-fold sites. As discussed elsewhere, interactions between adsorbates can influence their adsorption energies at coverages above a certain threshold (0.25-0.5 ML) [46,47,57], which potentially can affect the predicted rates. Therefore, we used an interaction model similar to that reported earlier to describe interactions between adsorbed reaction intermediates [46,47]. The fitted interaction parameters are calculated on Pd, Pt, Rh and Ru and scaled with the sum of C * and O * adsorption energies, as described in section 3 of the supplementary material (available from stacks.iop.org/NJP/15/125021/mmedia).
As shown in figure 3(c), the inclusion of adsorbate-adsorbate interactions only affects regions of the activity volcano where coverages of reaction intermediates exceed 0.5 ML. The predicted rates in the upper right section are hence unaffected since their coverages do not exceed this threshold. Similarly, the overall shape of the volcano, as well as the position of the maximum rates, are essentially unaffected. The volcano with adsorbate-adsorbate interactions differs from the one without interactions primarily in that it increases or decreases the predicted rates of the most strongly binding metals. The lower section of the volcano, as well as the peak vicinity, are predicted to have lower rates. In these areas, the dissociative adsorption of CH 4 has the highest activation barrier, thus is the rate-determining step. With the interactions included, high coverages of adsorbed oxygen and carbon destabilize the adsorption of CH * 3 , which increases the methane activation barrier.
By incorporating multiple reaction paths for CO formation, microkinetic modeling can also determine which path is likely to be the most important under the given reaction conditions. Based on the predicted reaction rates for each of the three CO formation paths (figure 5), the direct formation of CO from C-O coupling has the highest predicted rates at the top of the volcano, as compared to the hydrogen-assisted routes involving COH * and CHO * . This result is consistant with the theoretical work on metal (211) surfaces, which are conducted by Bengaard [27], Shetty [33] and Bhattacharjee [32], but contradictory to the results reported by Inderwildi [58], Wang [28,29] and Blaylock [30] on flat surfaces. There might be two reasons for the disparateness. Firstly, our scaling relations are obtained using energies calculated on stepped metal (211) surfaces. Generally the steps are more reactive than flat (111) surfaces for dissociation reactions, while whether the steps could promote bond formation depends on the individual reaction as well as the particular metal [59]. Therefore, the different structure sensitivity of elementary reactions maybe one of the reasons for different results of the primary pathways. Secondly, most of the work on the reaction mechanism only considered the potential energy of each pathway, neglecting temperature and coverage effects. From figure 2 we can see that the free energy barrier for activating gas-phase species are higher due to the significant loss of entropy at high temperatures. As depicated in figure 4, C * is one of the dominate species on the surface, the coverage of C * might be 2-4 orders of magnitue larger than CH * in the active area. Even if the activation barrier might be a litter lower for CH-O pathway, the rate for direct CO formation pathway could still be significantly higher. For less active metals, however, such as Pt and Pd, the hydrogen-assisted route involving COH * will likely also be important because coverage effects only play a minor role as compared to the activation barrier, in accordance with earlier work by Chen [34] and Inderwildi [35]. We also analysised the reaction pathways for CO formation under two industrial conditions, and the result that direct C-O coupling is the main route still holds.
Finally, the microkinetic model also provides a way to predict changes in rates as a function of changing gas compositions, similar to what might be experienced in a packed-bed reactor with changing temperatures and extents of conversion. In figure 6, the predicted rates (including interaction effects) for catalysts for the inlet and outlet conditions of such a reactor are given. While the overall shape of the volcano is conserved, it is evident that higher rates are predicted for the reactor outlet, primarily due to the increase in reactor temperature.

Thermodynamic and oxidation filters.
Having established a volcano shaped relationship for methane steam reforming as a function of two descriptors, we will now proceed to screen a large number of transition metal alloys with respect to their predicted rates, stability, and cost. Alloys with compositions of A 3 B and AB (where A and B are transition metals) were considered. The full list of alloys is given in section 4 of the supplementary material (available from stacks.iop.org/NJP/15/125021/mmedia). Using DFT, the formation energies of these alloys were calculated in connection with an earlier study of methanol formation [8]. We here consider those alloys with formation energy less than 0.2 eV per unit cell to be 'potentially' stable. The carbon and oxygen adsorption energies on the stepped (211) surfaces of these thermodynamically stable alloy materials were also calculated using DFT in the study by Studt et al [8] and are available through the CatApp [50]. As depicted in figure 7, we used the activity volcanoes from the previous section to predict the CH 4 steam reforming rates of these materials. We also assessed the stabilities of these alloys in terms of their predicted energies  figure 7(a), where it is evident that many of the alloys can be oxidized under inlet conditions. As expected, the majority of the alloys that are predicted to be resistant to deactivation via bulk oxidation are located at the upper right portion of the volcano-that is, they have relatively low binding energies of carbon and oxygen.

Price and activity filters.
One of the main drivers for new catalyst development is the desire to reduce the use of rare and expensive metals. As depicted in figure 7(b), we screened potential alloy formulations based on the costs of their component metals. Alloys with estimated costs less than $5000 kg −1 are denoted in the volcano plot using white stars, while more expensive alloys are denoted using black stars. This cutoff value of $5000 kg −1 was chosen to eliminate alloys containing metals that are relatively expensive (e.g. Pt, Pd, Rh, Ru, Os, Au, Ir). The remaining alloys, as depicted in figure 8, are the most interesting candidates for CH 4 steam reforming-since they are predicted to be stable under reaction conditions, formed from earth abundant materials, and have high predicted reaction rates.
It can be seen from figure 8 that most of the active, stable and inexpensive binary alloys contain Ni, Fe and Co in the A 3 B composition. It is also observed that the inlet conditions are more demanding in terms of stability than the outlet condition. This is reasonable because at the inlet the partial pressure of H 2 O is significantly higher, and the temperature lower, both of which tend to stabilize oxide formation. For all metals their oxides becomes less stable when the temperature increases. Two different (211)-like surface structure terminations are    [61]. Other materials, such as Ni-Ge and Co-Ge, have not been investigated previously, and are therefore interesting candidates for CH 4 steam reforming. A short list of the potentially active catalysts by screening method can be found in table 2.
Although some of the alloys are located close to the top of the volcano, they might experience problems with coking. As illustrated in figure 4, materials located on the left side of the activity volcano will likely have high coverages of carbon, and thus may be more prone to coking, whereas materials on the right side with lower carbon binding energies may be more resistant to coking and show better performance under steam reforming conditions. A rough guide for this, the equilibrium line for CO * decomposition to C * , is plotted on the volcano in figure 8. Alloys that are close to or below this line will have a high driving force for the formation of atomic carbon, which is a precursor for the formation of graphitic carbon.

Conclusion
We have constructed a descriptor-based volcano rate plot for methane steam reforming that combines a microkinetic model with scaling relations across the periodic table while taking adsorbate-adsorbate interactions into account. The model predicts the reaction rates of transition metals and metal alloys for the CH 4 steam reforming reaction based on carbon and oxygen adsorption as the two descriptors. This approach allowed for a detailed description of the reaction mechanism for different catalytic surfaces under the given reaction conditions. In addition we have used the model to screen a large number of transition metal alloys with respect to their steam reforming activity. The combination with filters that consider stability and cost enables the screening for novel leads for active, stable, and low-cost steam reforming catalysts.
The methods outlined in this paper are very general, making their application possible for a wide range of reactions. Starting from data stored in the public online database CatApp, one can begin to understand very quickly what characteristics an optimal catalyst might have. The described framework also allows one to consider what effects might be expected due to changes in catalyst composition and reaction conditions. This is especially true given the ability of the model to provide insights into the expected reaction mechanisms for a given catalytic surface. Therefore, we suggest the presented method herein could become useful as a starting point in the search for technical catalysts in the future.