Activity Trends in the Propane Dehydrogenation Reaction Catalyzed by MIII Sites on an Amorphous SiO2 Model: A Theoretical Perspective

One class of particularly active catalysts for the Propane Dehydrogenation (PDH) reaction are well-defined M(III) sites on amorphous SiO2. In the present work, we focus on evaluating the catalytic trends of the PDH for four M(III) single-sites (Cr, Mo, Ga and In) on a realistic amorphous model of SiO2 using density functional theory-based calculations and the energetic span model. We considered a catalytic pathway spanned by three reaction steps taking place on selected MIII–O pair of the SiO2 model: σ-bond metathesis of propane on a MIII–O bond to form M-propyl and O–H group, a β-H transfer step forming M–H and propene, and the H–H coupling step producing H2 and regenerating the initial M–O bond. With the application of the energetic span model, we found that the calculated catalytic activity for Ga and Cr is comparable to the ones reported at the experimental level, enabling us to benchmark the model and the methodology used. Furthermore, results suggest that both In(III) and Mo(III) on SiO2 are potential active catalysts for PDH, provided they can be synthesized and are stable under PDH reaction conditions.


Introduction
Industrial processes involving the cleavage of strong C(sp 3 )-H bonds [1] are among the most challenging ones, being highly energy demanding processes [2][3][4]. Within alkane dehydrogenation, propane dehydrogenation (PDH) to propylene and hydrogen (Eq. 1) has gained major interest in the past years due to the discovery of shale gas [5] mainly in the US [6][7][8]. Moreover, the PDH reaction in combination with metathesis, i. e. propane to olefins (PTO) can provide a way to produce propene, ethene and butene in a singlestep [9].
Nevertheless, to overcome the thermodynamic limitations imposed by the highly endothermic PDH reaction (standard enthalpy ∆H°2 98 = +124 kJ mol −1 ), high temperatures are required to obtain significant conversions (typically above 550 °C). This leads to deactivation and poisoning of the catalysts via coke formation [10]. In industry, there are two historical heterogeneous catalysts used for the PDH reaction: the alumina-supported PtSn nanoparticles and the CrO x / Al 2 O 3 system, also known as the Houdry or Catofin catalyst [11], albeit both systems suffer from deactivation and need frequent regeneration. Therefore, the synthesis of new catalysts is crucial [12]. Among the first family of catalysts, i.e. based on Pt NPs, it was shown that silica-supported PtGabased catalysts were highly efficient for the PDH reaction [13] and their high performance and stability was recently rationalized via first principles calculations in combination with metadynamics [14].
Metal-based single-site catalysts on SiO 2 have also been found to display high activity towards the PDH reaction (1) [10]. For instance, Cr(III) centres on silica [15], analogous to the Houdry and Catofin catalysts in alumina, are also active catalysts in the reaction. Lewis acidic Ga(III) centres on SiO 2 also display high activity and remarkable selectivity and stability [16,17], Ga 2 O 3 is also catalytically active towards the PDH reaction [2,18] and recently Ga(III) sites have also been proposed as their active reaction centers [19].
In previous studies, we showed how the use of an amorphous SiO 2 model [24] is important to understand the reactivity of single-sites present on silica, being able to account for strain, an aspect which is not as straightforward to capture using cluster models [25]. This was exemplified by the olefin polymerization reaction catalyzed by Cr(III)/SiO 2 [26,27]. More recently, for the PDH reaction, we could also show by evaluating the reactivity of different Ga(III) sites that strain is a key parameter controlling the activity, showing that an intermediate strain of the reactive Ga-O pair leads to the highest catalytic activity in the PDH reaction [28].
The considered propane dehydrogenation pathway has the following reaction steps: (1) the C-H bond activation of propane on the M III site forming a M III -propyl intermediate and an OH group (Step 1), (2) the β-hydride transfer producing a hydride and propene (Step 2), (3) the H-H coupling or reverse σ-bond metathesis of H 2 (Step 3) following the de-coordination of propene to regenerate the initial M III site. This catalytic route was already considered for Cr [29] and Ga-single sites using cluster [30] and amorphous periodic models, respectively [28].
Given that Cr and Ga single-sites on SiO 2 are known to be active catalysts of the PDH reaction, together with the latter elements, we decided to evaluate two additional representative elements laying just below Cr and Ga in their respective groups of the periodic table, i.e. In and Mo. Therefore, in this work, we evaluated the reactivity of four different single metal(III) sites on silica (M = Cr, Mo, Ga and In) using the previously mentioned SiO 2 amorphous model [24]. For the present study, we selected the M-O pair that displayed the highest activity towards the PDH for the Ga(III)/SiO 2 system. Here, we extend our work on Ga to three additional metals and compare among them. Our study displays the catalytic differences among the evaluated metal sites, highlighting their differences in the key PDH steps, finally proposing which one should display the highest catalytic activity in PDH using the energetic span approach.

Construction of the M(III) single-sites on SiO 2
M III (Cr,Mo,Ga,In)/SiO 2 models are constr ucted using a previously developed amorphous silica model [24], which corresponds to a slab of dimensions 21.4 Å × 21.4 Å × 34.2 Å and contains 372 atoms. The silica model exposes five isolated silanol groups (SiOH) and has a surface silanol density (1.1 OH nm −2 ) close to the density found for silica partially dehydroxylated under vacuum at 700 °C (SiO 2-700 , 0.8 OH nm −2 ), which is for instance, the one used experimentally to prepare welldefined Ga III sites active in the propane dehydrogenation reaction and related systems [16]. The M III sites can be introduced into the silica model by substituting surface "SiOH" groups by a M 3+ i.e. turning (≡SiO) 3 SiOH sites into (≡SiO) 3 M III sites as previously carried out to build the corresponding Cr III and Ga III sites [26][27][28]. This model provided five types of Ga, Cr III /SiO 2 models, namely models I-V [26][27][28]. In our previous work on Ga, we probed the reactivity of all sites towards the propane dehydrogenation reaction [28]. We found that the most active site corresponds to the V-O3 site, in which V refers to the site label of the metal, while O3 refers to the specific oxygen involved in the PDH reaction. Therefore, in the current study, we selected V-O3 site to compare the reactivity of the single-site systems composed of the previously discussed metals: Cr, Mo, Ga and In, respectively. Figure 1 shows all the constructed tri-coordinated single-sites on the amorphous SiO 2 model together with the geometrical characteristics of each of the four metal sites (see Fig. 1

Evaluation of the Propane Dehydrogenation on the Selected M(III)/ SiO 2 Sites
We evaluated the propane dehydrogenation reaction catalyzed by four different single-site M III metals (Cr, Mo, Ga and In) considering the previously discussed catalytic cycle, as shown in Scheme 1.

Step 1: C-H Bond Activation of Propane
The reactivity of the four M(III)-O(3) site towards the C-H activation of propane was the first reaction step evaluated. Table 1  This site has a high degree of strain, and therefore is the most active site among all the previously evaluated sites for the case of the Ga single-sites on SiO 2 using the same amorphous model than the one used in the current study [28]. The energy barriers follow the trend: Ga < In < Cr ≪ Mo. The reaction energy barriers take values within 18.5 and 27.8 kcal mol −1 , while reaction energies range between − 2.0 and − 26.2 kcal mol −1 . In addition, for most metals, the energy barrier for the C-H activation decrease when the reaction energy is more favored. The only exception to this trend is between Ga and In, since the C-H activation energy barrier is slightly lower on Ga than on In (18.5 vs. 18.9 kcal mol −1 ) whereas the reaction is slightly less favored on Ga than on In (− 25.2 vs. − 26.2 kcal mol −1 ). In any case, the energies obtained for Ga and In are very alike, being the two most active metals for the direct propane activation via a σ-bond metathesis; forming a M-propyl group and a new O-H bond. The next most active metal for this reaction step is Cr with a slightly higher energy barrier (21.9 kcal mol −1 ) than In and Ga (18.9 and 18.5 kcal mol −1 , respectively). In contrast, Mo presents a significantly higher energy barrier than the other metals: 27.8 kcal mol −1 .
Similar to what we found when analyzing the C-H activation of propane on different Ga sites on silica, there is a Brønsted-Evans-Polanyi relationship followed to certain extent (R 2 = 0.816) between the electronic energy barrier and the reaction energy ( Figure S1 in electronic energies and Figure S2 and R 2 = 0.763 in Gibbs energies). There is actually a good correlation between the energy barrier and the reaction energies of the reverse step of the C-H activation, as reported in the ESI (R 2 = 0.956 in electronic energies, Figure S3

Step 2: β-H Transfer and Propene Decoordination
The next reaction step considered in the PDH reaction is the β-H transfer, in which one H atom in β-position of the propyl group is transferred to the metal atom forming an M-H bond, while propene remains coordinated to the M(III) center. Table 2 summarizes the energetics of this reaction step for all the evaluated metals. This step is endoenergetic and has rather high energy barriers for all metal atoms.
Moreover, the calculated energy barriers take significantly higher values than the previously evaluated C-H activation of propane reported in the former subsection (vide supra). The values of the energy barriers for the β-H transfer step follow the trend: Cr < Mo < Ga ≪ In. As expected, Cr and Mo, which are transition metals with partially filled d orbitals have lower energy barriers (35.4 and 37.2 kcal mol −1 , respectively) than Ga and In (41.7 and 53.7 kcal mol −1 , respectively), belonging to the group 13 of the periodic table.
The corresponding geometries for all the transition-state structures are shown in Fig. 3, which depicts their most characteristic geometrical features concerning the β-H transfer step. The geometry of the corresponding structure for Mo is different than for the other metals. In this case, the TS corresponds to the change of conformation of H from one side of the O-M-O plane to the other with propene barely interacting with the Mo (see Fig. 3c). In this case, the β-H transfer preceding this TS is slightly lower in energy than this H migration. Therefore, we took the latter process as the key transition state connecting M-propyl and the M-H + propene species.
In this case, the relation between ∆E ‡ (kcal mol −1 ) and ∆E does not follow a linear trend (R 2 = 0.584, Figure S5). On the other hand, in our previous study on five different Ga sites, this step followed a rather well transition-state scaling, Scheme 1 The proposed reaction mechanism for the dehydrogenation of propane on M III /SiO 2 sites. For some metals, an additional siloxane group coordinates to the M III centre which is shown in grey. In the present work, we report the reaction mechanism on four M(III) centers (M = Cr, Mo, Ga, and In)   Figure S6 and R 2 = 0.843 in Gibbs energies, Figure S7) with respect to the previously mentioned relation between ∆E ‡ and ∆E (R 2 = 0.57). The scaling trend is, however, not as high as the one we reported when comparing the five different Ga sites (R 2 = 0.997) [28]. Possibly, more data and complex functions, such as those provided by machine learning algorithms, are needed to predict the energy barriers of this reaction step with a high degree of confidence. The next step we considered in our mechanistic analysis was the decoordination of propene from the metal moiety. In all cases, the decoordination is endoenergetic taking values equal to 20.1, 27.2, 18.0 and 21.3 kcal mol −1 for Cr, Mo, Ga and In, respectively. The decoordination is preferred in the order of: Ga > In > Cr ≫ Mo. While Ga, In, and Cr present rather similar values ranging 18.0-21.3 kcal mol −1 for the decoordination of propene step, the decoordination of propene from Mo is much more energy demanding (27.2 kcal mol −1 ) since the bond between the propene and Mo is significantly stronger than for the rest. In any case, at the high temperatures needed for the PDH reaction, the propene decoordination is clearly exergonic (∆G < 0) for all metals as it is a highly entropically-driven step owing to the release of propene into the gas-phase.

Step 3: H-H Coupling and H 2 Formation
Finally, we calculated the last reaction step along the postulated mechanism for the PDH reaction, which is the coupling of M-H and O-H groups, i.e. the H-H coupling-or   Figure S8, and R 2 = 0.966 in Gibbs energies, Figure S9).
The transition-state structures corresponding to all H-H coupling steps are shown in Fig. 4, which depicts their most characteristic geometrical features.

Overall Catalytic Cycles for PDH
At this point, we can discuss the overall reactivity of the four evaluated metal sites based on the Gibbs energy profiles at 550 °C and 1 bar, which are shown in Fig. 4. The energy profile in electronic energies is given in the ESI ( Figure S10). Figure 5 provides the information to compare the reactivity among all four evaluated metal sites. As previously reported, the propane dehydrogenation reaction is endergonic by 7.4 kcal mol -1 , in good agreement with the experiment, since the PDH reaction has ca. 30% equilibrium conversion for propane at 550 °C and 1 bar [10].
Subsequently, we used the energetic span model [31] to compare the catalytic activity of the four evaluated systems. Within this approach, the turnover frequency (TOF) of a catalytic cycle is a function of the energetic span (∂E). This quantity (∂E) depends on the energy of the TOF-determining transition state (TDTS), and the TOF-determining intermediate (TDI). Simplifying the approach, the TDTS corresponds to the highest energy in the Gibbs energy profile, while the TDI is generally the most stable intermediate in the energy profile. When the TDTS of a catalytic cycle  appears after the TDI in the Gibbs energy profile, ∂E is the energy difference between these two steps. Conversely, when the TDI appears after the TDTS, the ΔG of the reaction (ΔG r ) is added to this difference. Overall, the energetic span model defines ∂E using the following equation.
Based on the energetic span E , the TOF is calculated at 550 °C for all systems using the expression: The former equation is in principle valid for exergonic reactions, which leads to a positive TOF when calculated. The TOF is understood as a catalytic flux via the span model, in analogy with Ohm's law in electric circuits [31]. Thus, within the flux analogy, it is positive and goes forward for exergonic reactions, while it is negative and goes backwards for endergonic reactions. Nevertheless, the TOF is defined differently at the experimental level and is always positive, since it is obtained based on the conversion of the reactants to products. The PDH reaction is endergonic under our evaluated conditions (the ΔG r term is positive at 550 °C), and therefore we would obtain a negative TOF, meaning the catalytic flux would go backwards. Thus, in order to compare the catalytic activity of the evaluated sites with the reported positive experimental TOF, we will use the above-mentioned equation to compare the catalytic activity of the four metal sites in a semiquantitative way, although the ΔG r term is positive in our case. We refer the reader to our previous work for further details and discussion about the specificities and caveats of the application of the span model to the PDH reaction [28]. For the Cr(III)/SiO 2 single-site, the TDTS corresponds to the β-H transfer step, being located at + 58.5 kcal mol -1 with respect to initial reactants. Since all intermediates are less stable than initial reactants, the TDI species within the energetic span model corresponds to the origin of energies of the Gibbs energy profile (catalyst + reactants). Therefore, the value for the energetic span is equal to + 58.5 kcal mol -1 , leading to a calculated TOF equal to 18.1 h −1 . Our calculated TOF agrees rather well with the initial turnover frequency (TOF) of 10.3 h −1 reported experimentally for Cr(III)/SiO 2 at 550 °C and 1.5 bar [15].
The TDTS of the Mo(III)/SiO 2 single-site is the TS of the β-H transfer step, which is located at 48.1 kcal mol −1 with respect to the initial reactants. Similar to Cr, the TDI in this case also corresponds to the origin of energies of the Gibbs energy profile. Therefore, taking the energetic span for Mo equal to 48.1 kcal mol -1 , leads to a rather high calculated TOF of 1.05·10 5 h −1 .
The TDTS of the catalytic cycle for the evaluated Ga site corresponds to the C-H activation step located at 31.7 kcal mol -1 in the Gibbs energy profile. In this case, Fig. 4 Transition-states corresponding to the initial H-H coupling steps for the four evaluated single-metal (III) sites on SiO 2 : a Cr, b Mo, c Ga and d In the TDI of the catalytic cycle appears after the TDTS. The TDI corresponds to the gallium hydride species, with a relative energy equal to − 19.8 kcal mol -1 with respect to initial reactants. Thus, the energetic span after summing up the ∆G r term (+ 7.4 kcal mol -1 ) is equal to 58.9 kcal mol -1 corresponding to a TOF of 14.2 h −1 . This result agrees well with the value reported experimentally for the Ga(III)/SiO 2 system, which took a value of 20.4 h −1 at 550 °C [32].
Finally, for In(III)/SiO 2 , the TDTS corresponds to the β-H transfer step, which is located at 53.5 kcal mol -1 above initial reactants. In this case, two species compete as the TDI, the In-propyl species formed after the activation of the C-H bond of propane and the In-H species after propene decoordination. Both species have similar energies, being equal to − 0.39 and − 1.27 kcal mol -1 , respectively. After the application of the energetic span model [31] the former species is the TDI, i.e. the In-propyl species, located at − 0.39 kcal mol -1 in the Gibbs energy profile. Making this assumption, the energetic span is equal to 53.9 kcal mol -1 and the calculated TOF takes a value of 3.02 × 10 2 h −1 .
Comparing the calculated TOF for all four metals on the evaluated sites, the obtained results suggest that the overall PDH activity follows the trend Ga ≈ Cr < In ≪ Mo. Therefore, the present study suggests In(III) and Mo(III) singlesites on SiO 2 could be potentially active catalysts for the PDH reaction, in particular Mo(III). One has to bear in mind, however, the following considerations: (i) the selectivity of the PDH reaction was not evaluated in the present work, but only the activity, (ii) the hypothetical In(III) and Mo(III) sites should be synthetized and remain stable under reactions, and (iii) surface heterogeneity was not evaluated in depth in the present study since we only evaluated one M-O site on the SiO 2 amorphous model.

Conclusions
In the present work, we evaluated at theoretical level the reactivity of four M(III) sites (Cr, Mo, Ga and In) on an amorphous model on silica towards the propane dehydrogenation reaction (PDH). We considered three activated steps for the reaction, taking place on a M-O reactive pair: namely the C-H activation of propane, the β-H transfer step from the M-propyl intermediate generating M-H and propene, and the H-H coupling step, which forms H 2 and regenerates the initial M-O bond.
A comparison in reactivity among the four evaluated metals based on electronic energies shows that the most active centres towards the σ-bond metathesis of propane are the Lewis acidic Ga(III) and In(III) centres, while the least active ones are the Cr(III) and Mo(III) sites. Conversely, for the β-H transfer step, Cr and Mo have much lower Gibbs energy barriers than Ga and In centres. The H-H coupling step, i.e. the inverse of the σ-bond metathesis step has also much lower Gibbs energy barrier for Cr and Mo than for Ga and In. When comparing the overall Gibbs energy profiles within the energetic span model for the four evaluated metals, the TDTS for Cr, Mo, and In is the β-H transfer step, while for Ga, it is the C-H activation step, which is located slightly higher than the β-H transfer in the Gibbs energy profile. Within the energetic span model, Cr and Ga are predicted to have similar catalytic activity with values in good agreement with the experimental results. Based on our results, In and Mo centers should have even higher activity, provided they are stable under reaction conditions and have high selectivity towards the desired products, i.e. propene and H 2 .

Computational Methods
DFT calculations based on the Gaussian and plane waves (GPW) formalism [33] were carried out using the Quickstep (QS) module [34] of the CP2K program package [35,36]. The functional chosen was PBE [37-39] with short range Gaussian double-ζ basis sets [40] optimized from molecular calculations. The energy cutoff of the auxiliary plane wave basis set was set to 500 Ry. Goedecker-Teter-Hutter pseudopotentials [41][42][43] were used. For the faster convergence of the wavefunction, the orbital transformation (OT) method was used [44,45]. The OT method is based on a minimization of the energy functional using a new set of variables to perform orbital transformations. This method is demonstrated to be very efficient in comparison to diagonalisation/DIIS based methods, especially for large systems and basis sets. A tetragonal simulation box of base area 21.4 Å × 21.4 Å and thickness 34.2 Å (ca. 24 Å of which correspond to a vacuum slab added in order to avoid interactions between images in the z direction) was used [24]. Ground state structures were obtained by energy minimization with the BFGS algorithm [46][47][48][49][50]. Initial transition state guesses were generally obtained from CI-NEB [51][52][53][54][55] band calculations. Transition state structure optimizations were performed using the dimer method [56,57] with the conjugate gradient optimizer and the two point based line search. It was further confirmed that TS structures actually corresponded to saddle points by frequency analysis via the finite difference method with a displacement equal to 0.01 bohrs. The Gibbs energies of the gas phase species was estimated by considering translational, rotational and vibrational degrees of freedom (DOF) at 550 °C and a pressure equal to 1 bar. All DOF of adsorbed species, including TSs, were treated as vibrational within the harmonic approximation under the same conditions to obtain the Gibbs energies. Thermal effects from the atoms of the catalysts were not considered.