Unveiling the catalytic potential of the Fe(IV)oxo species for the oxidation of hydrocarbons in the solid state

Abstract

We present the computational study of the ferryl-catalysed oxidation of methane into methanol in a solid-state system, the metal–organic framework MOF-74 with Fe(IV)O moieties in its cavities. We use spin-polarised ab initio molecular dynamics at the hybrid HSE06 level of theory to simulate this process as three consecutive steps: the hydrogen abstraction from methane by Fe(IV)O, the rebound of the resulting CH₃˙ radical to form a methanol molecule, and the detachment of the product from the reactive site. Our computational approach accounts for both enthalpic and entropic effects at room temperature. The calculations indicate that the overall oxidation process occurs with a free energy barrier of 95.6 kJ mol⁻¹, with the detachment of methanol as the rate-determining step. For the abstraction step, we estimate a free energy barrier of 51.1 kJ mol⁻¹ at 300 K and an enthalpy barrier of 130.3 kJ mol⁻¹, which indicates the presence of a substantial entropic contribution. van der Waals dispersion interactions play also a significant role in the overall reaction energetics. Our study suggests the potential applicability of metal–organic frameworks in the industrial production of fuels from saturated hydrocarbons and indicates that it is necessary to further investigate whether other factors, such as stability and easy regeneration, favour these materials.

---

1 Introduction

The iron-oxo Fe(IV)O (or ferryl) species plays a central catalytic role in many biological and industrial processes including aerobic respiration, catabolism, and the in vivo oxidation and separation of hydrocarbons.^1–12 This variety of roles is related to the high reactivity of Fe(IV)O, whose understanding has been the subject of substantial theoretical and computational work.^{9–11,13–16} These studies have shed light onto the complexity of the interaction between a coordinated Fe(IV)O moiety and its substrate during oxidation reactions, as well as in the processes that drive the reduction of O₂ or H₂O₂ and the formation of high-valent Fe species.^{14,15,17–19} It has been also shown that one of the most complex characteristic of Fe(IV)O is the frequent presence of several competing spin states, and the consequent appearance of two- and multi-state reactivity.²⁰ Typically, quintet and triplet states can be close in energy,²¹ but they exhibit different chemical properties and reaction pathways. In most situations, quintet systems show higher electrophilicity and oxidative activity, and they can therefore act as powerful electron acceptors. The electrophilic character of Fe(IV)O is crucially related to the energy of low-lying virtual orbitals. In particular, extensive computational studies carried out on a variety of gas-phase and solvated Fe(IV)O complexes have shown that an anti-bonding 3σ* orbital is anomalously stabilised in the quintet state, which leads to the experimentally observed dramatic increase in reactivity compared to triplet systems.^22–24

So far, the vast majority of the Fe(IV)O systems have been studied in the gas phase or in solution, while little is known about the catalytic activity of the Fe(IV)O species in the solid state. Although some theoretical studies have been published very recently on the oxidation of ethane^10,25–29 in metal–organic frameworks (MOFs) and methane in zeolites,³⁰ and of methane storage in MOF-74,³¹ the shortage of experimental work is surprising given that several crystalline systems can be synthesised relatively straightforward with structures and coordination properties mimicking those observed in highly reactive Fe(IV)O molecular systems. Noticing this enormous knowledge gap and that the family of MOF-74 has been identified as a potential candidate to catalyse the methane activation,³² we have recently used density-functional theory (DFT) calculations to assess the ability of MOFs to act as suitable guests for Fe(IV)O catalytic centres and their reaction substrates.^21,33 We have shown that MOF-74 can indeed be used to stabilise quintet Fe(IV)O groups, which can act as catalytic centres for substrate molecules adsorbed in the MOF cavities. In MOF-74, the coordination geometry of Fe(IV)O is close to the ideal one postulated from gas-phase calculations,²³ with four oxygen atoms in the equatorial positions, a relatively unhindered axial coordination site and the Fe–O axis point toward the interior of the cavity and easily accessible to the substrate. We have also shown, based on total energy calculations on a crystalline fully periodic MOF-74/Fe(IV)O model, that, as expected, these Fe(IV)O groups exhibit a reactivity in the hydroxylation of methane similar to those observed in some gas-phase and solvated quintet molecular Fe(IV)O complexes. However, this recent work has been limited to static calculations with no inclusion of temperature effects. Furthermore only the first step of the hydroxylation of methane, the abstraction of an H atom from the substrate, has been studied in detail.^29,34–37

The goal of this work is therefore to model the whole oxidation process of methane into methanol in a MOF-74 cavity in the presence of Fe(IV)O sites. This process is modelled as three consecutive steps that are illustrated in Fig. 1, as in the classical rebound mechanism proposed by Groves,^38,39 starting with the abstraction of an H atom from methane by a reactive Fe(IV)O group, followed by a rebound step in which the CH₃˙ radical binds to the newly formed Fe(IV)OH group to form a methanol molecule, and finally the detachment of CH₃OH from the catalytic centre. The present work is carried out using spin-polarised ab initio molecular dynamics (AIMD) at the hybrid DFT level of theory, since we have shown in previous work that range-separated hybrid functional such as sc-BLYP⁴⁰ provide accuracy in reaction barrier for this system comparable to global hybrid functionals such as B3LYP^41–43 and that the use of hybrid exchange in the density functional is required to remove spurious self-interaction effects responsible for anomalous reaction profiles and unphysical reaction energy barriers.^21,33 Obtaining accurate reaction (free) energy barriers is crucial for identifying the reaction's rate-determining step (abstraction, rebound, or detachment). In all three steps, we also evaluate entropic effects at room temperature, using AIMD based methods. To our knowledge, this work is the first attempt to model a hydrocarbon oxidation reaction using AIMD methods with accurate exchange–correlation functionals within an MOF cavity, which we model as a fully periodic crystalline system.

Fig. 1 Representation of the three steps involved in the conversion of methane into methanol in this work: abstraction, rebound, and detachment. After each step, our study generates the energies for the enthalpy ΔH and free energy ΔG.

After this introduction, this manuscript is organised as follows. In section 2, we provide details of the DFT methodologies and AIMD simulations and the methods used to compute free energy barriers from these simulations. The abstraction, rebound, and detachment steps of methane in Fe(IV)O/MOF-74 are fully analysed in section 3. Section 4 summarises our main findings and indicates potential directions for future work.

2 Simulation methods

The initial structure of Fe(IV)O/MOF-74 used in our calculations is derived from the data concerning acetylene/MOF-74 (ref. 7) deposited in the Cambridge Crystallographic Database,⁴⁴ which was modified using the Materials Studio suite package.⁴⁵ The space group is R3̄ (148), with a unit cell of dimensions of 25.9 × 25.9 × 6.9 ų and angles α = β = 90° and γ = 120°. Fig. 2 shows that this cell contains 60 atoms and corresponds to 6 primitive cells. The primitive cell of Fe(IV)O/MOF-74 contains 10 atoms: one Fe atom, three framework O atoms, four C atoms, one H atom and one O(oxo) atom bonded to Fe. We substituted the acetylene molecule in the original structure with an oxygen atom O(oxo), at a distance of 1.68 Å from the Fe atom, which is slightly larger than typical Fe(IV)–O(oxo) distances determined for gas-phase complexes (1.60–1.62 Å).²³ We then relaxed the atomic positions with the COMPASS2 (ref. 46) force field. After optimisation, the Fe(IV)–O(oxo) bond length decreases to 1.65 Å. The resulting structure is then employed as input for our DFT geometry optimisations and ab initio molecular dynamics.

Fig. 2 Representations of the periodic Fe(IV)O(oxo)/MOF-74 cavity with a methane molecule in the centre, with the top (a) and side (b) views of the unit cell used in the CRYSTAL17 calculations and its replicated cell (c). Fe atoms are blue, O atoms red (notice the O(oxo) atoms oriented toward the cavity centre), C atoms green and H atoms grey. These views are generated with the OVITO package.⁴⁷

Both types of first-principles computations are employed to determine reaction barriers such as enthalpies (ΔH), with geometry optimisations, and free energies (ΔG) with ab initio molecular dynamics. The purpose of calculating both ΔH and ΔG is to estimate the entropic contribution (ΔS) to each reaction step. We next describe the computational parameters for the DFT codes that we have used in the simulations in both cases.

2.1 CP2k setup for ab initio molecular dynamics

We use the CP2k/QUICKSTEP code^48,49 to run AIMD simulations, with standard GTH pseudopotentials^50,51 and double-ζ basis sets with polarisation functions (DZVP-MOLOPT-SR-GTH). van der Waals forces are computed using Grimme's DFT-D3 correction.⁵² We use the default tolerance for the charge density residual of 1 × 10⁻⁵ e⁻ Bohr⁻³ for the self-consistent solution of the Kohn–Sham equations with an energy cutoff of 500 Ry, which is justified in Fig. S1 in the ESI.† All CP2k calculations are carried out in the Γ point approximation on a super-cell containing 365 atoms, with dimensions a = b = 25.9 Å and c = 13.8 Å with α = β = 90.0° and γ = 120.0° (Fig. 3). This supercell contains 36 Fe atoms, which are configured in the quintet spin state all along this work, given that preliminary calculations show that this state has a lower total energy than the triplet's by 10.19 eV and singlet's by 11.8 eV when the methane molecule is far away from the Fe(IV)Ooxo reactive site. These HSE06 results are in agreement with our previous findings with GGA and other hybrid levels of theory.³³

Fig. 3 Top (a) and perspective (b) views of the periodic Fe(IV)O(oxo)/MOF-74 system with a methane molecule used in the CP2k calculations. Fe atoms are blue, O atoms in red (notice the O(oxo) atoms oriented toward the cavity centre), C atoms in green, and H atoms in white. These views are generated with OVITO package.⁴⁷

Free-energy reaction barriers are estimated using the potential of mean force (PMF) method^53,54 from constrained AIMD NVT simulations carried out at constant volume and at a temperature of 300 K. A timestep of 0.5 fs was used to integrate the equations of motion and a Nosé–Hoover thermostat^55,56 with a chain time constant of 1 ps was used to control the temperature. Due to the large computational cost required to run AIMD with HSE06 on spin polarised systems (ca. 350 seconds per AIMD step using 144 cores with CP2k v.6.1), the number of time steps is varied depending on the distance between the methane molecule and the O(oxo) atom. We used at least 1500 AIMD steps for all the cases in which the molecule is far from the O(oxo) atom (3.5 Å) and up to 4500 steps for short separation (ca. 0.9 Å), as longer simulation times are required to obtain meaningful averages at short distances. Depending on the reaction step to be considered, we adopt a different constraint length ξ between two atoms. In the case of the H-abstraction from methane (the first and, according to the standard hydrocarbon hydroxylation reaction mechanism, the rate determining step³⁸), ξ corresponds to the distance between the O(oxo) atom and one of the H atoms of the methane molecule. For the rebound step, ξ corresponds to the distance between O(oxo) and the C atom of the CH₃˙ radical generated at the end of the H-abstraction. We also consider the detachment of the CH₃OH molecule from the reaction site at the end of the second step. In this case, ξ corresponds to the distance between Fe and the O atom of the CH₃OH molecule (originally O(oxo)). For each of the three steps, we carry out a series of AIMD simulations, in which ξ is progressively decreased or increased, depending on the process, and for each value of the constrained distance, ξ_i, we compute the mean force of the constraint f(ξ_i) from an unbiased time averaged value of λ(ξ_i), where λ(ξ_i) is the average value of the Lagrangian multiplier λ(ξ_i) which is used to maintain the constraint to its fixed value:^24,53,54

Here, k_B = 1.38 × 10⁻²³ J K⁻¹, and T = 300 K, and f₀ is the value of f(ξ_i) for the initial simulation. The free energies associated with each process are given byFor the calculation of free energy barriers, ξ_N is chosen as the value of ξ closest to a free energy maximum.

2.2 CRYSTAL17 setup for geometry relaxations

Preliminary geometry optimisations and total energy calculations are carried out using the all-electron ab initio CRYSTAL17 code,⁵⁷ with default convergence settings. The results from these calculations are compared with our AIMD results with CP2k results to evaluate entropic effects. The following convergence criteria are used: 1.5 × 10⁻³ Ha Bohr⁻¹ for the largest component of the force and 10⁻⁷ Ha for the energy change between optimisation steps. This set of convergence thresholds has been used successfully in related systems, for instance graphenylene.⁵⁸ The DFT calculations are carried out with periodic boundary conditions and the range-separated hybrid HSE06 exchange–correlation functional. Pure density functionals are typically preferred in solid state applications, because of their reduced computational cost, especially in plane-wave based calculations and AIMD AI simulations.^59–62 In previous work,^21,33 we showed that the hybrid functionals B3LYP and sc-BLYP, which provide a reasonably accurate description of high-spin Fe(IV)oxo states in the gas phase,^63,64 are also more adequate than common generalised-gradient approximations (e.g. PBE and BLYP) in the solid state. This is largely to be attributed to the (partial) cancellation of self-interaction effects in hybrid functionals.³³ We do however notice that, for gas-phase or solvated systems, the OPBE functional may also provide an accurate alternative to B3LYP (ref. 24 and 64). A standard all-electron 6-31G**⁶⁵ basis set is used to represent the local atomic orbitals in terms of primitive Cartesian Gaussian functions. Polarization functions (p-functions for hydrogen and d-functions for carbon, oxygen and silicon) are used to ensure that the orbitals can distort from their original atomic symmetry, to better adapt to the molecular surroundings. Accurate truncation thresholds for the tolerances of the Coulomb and exchange bielectronic integral series are used in all calculations⁵⁷ to improve the convergence rate during the self-consistent solution of the Kohn–Sham equations. Brillouin zone integrations are carried out using a Monkhorst–Pack net of 2 × 2 × 2 k-points, and a ground-state energy convergence is enforced of 1 × 10⁻⁵ Hartree. van der Waals dispersion interactions are included semi-empirically using Grimme's DF2 parameterisation.⁵²

3 Results and discussion

3.1 Hydrogen abstraction by ferryl active site

Free-energy reaction barriers at room temperature are estimated using thermodynamic integration (section 2.1) with the constraint ξ_i corresponding to a set of O(oxo)–H(methane) distances. Fig. 4 shows the evolution of the corresponding mean force of the constraint as the distance between O(oxo) and H(methane) is progressively reduced. At large O(oxo)–H(methane) separations (3.6–2.7 Å), the mean force has a stable value of 5.7 kJ (mol⁻¹ Å⁻¹), which indicates a weak repulsion between these atoms. For shorter distances, the mean force increases and reaches a maximum of 52.0 kJ (mol⁻¹ Å⁻¹) at 1.8 Å. At shorter distances, the mean force decreases in consequence to the progressive formation of a covalent bond with O(oxo). The corresponding free energy profile, obtained from eqn (2) by numerical integration of the mean force is also shown in Fig. 4. The value of d_O(oxo)–H at which the mean force crosses the y axis (corresponding to a mechanical equilibrium between attractive and repulsive forces on the H atom) is 1.4 Å. We determine the free energy for the abstraction ΔG_abs = 38.3 kJ mol⁻¹, which corresponds to the maximum value of ΔG found at d_O(oxo)–H = 1.8 Å starting from 3.6 to 1.4 Å. Next, we estimate the enthalpy barrier ΔH_abs = 30.9 kJ mol⁻¹, which corresponds to the difference between the minimum and maximum of the total energy plotted in Fig. 5 and computed with CRYSTAL17. Then, we estimate a reaction entropy change ΔS_abs, defined as ΔG_abs andof −24.4 J (mol⁻¹ K⁻¹). A similar approach has previously been applied to estimate reaction entropies for FeO-catalysed H-abstraction reactions in water solution, which typically yields values of ΔS_abs of the order of 20 J (mol⁻¹ K⁻¹).²⁴ By comparison, we find that entropic effects in the H-abstraction are lower in our MOF environment than in water solution, a behavior that is consistent with the more ordered and rigid structure of this solid-state system. We will show in the next subsection that this finding is consistently observed as the reaction entropy increases when the rigidity of the system decreases upon the formation and detachment of the methanol molecule.

Fig. 4 Mean force f_O(oxo)–H of constraint, abstraction free energy ΔG_abs, and distance between the methane's carbon and hydrogen d_C–Hvs. the distance between the O(oxo) and this hydrogen obtained with AIMD using CP2k.

Fig. 5 The difference of the total (ΔE_DFT+disp), DFT (ΔE_EDFT), and long-range dispersion (Δ_disp) energies with respect to their values at 3.6 Å, and the distance between the central carbon atom with the reactive hydrogen (d_C–H) vs. the distance between the O(oxo) and methane's hydrogen obtained with CRYSTAL17 with the functional HSE06 during the abstraction step.

It is also interesting to observe that the entropic contribution in the case of MOF-74 is slightly negative, which indicates that the H-abstraction could be kinetically more favourable at low temperatures, although a more detailed analysis of the reaction free energy barrier as a function of temperature would be required to identify an optimal temperature range.

3.2 Rebound step: formation of the methanol molecule

Once the hydrogen atom has been transferred to O(oxo) (d_O(oxo)–H = 1.07 Å), we proceed to model with AIMD simulations the rebound step, in which the CH₃ radical approaches the O(oxo)–H group to form a methanol molecule. As in the previous section, AIMD simulations are carried out at constant volume and at a temperature of 300 K, with the constraint ξ set equal to the distance between O(oxo) and the C atom of the CH₃˙ radical. ξ varies over the distance range 3.8–1.3 Å, to yield the corresponding force of the constraint f_O–C, starting from the initial position of 3.2 Å. The calculated mean force of the constraint is shown in Fig. 6(a) and the corresponding free energy in Fig. 6(b).

Fig. 6 Rebound step: averages of the (a) force of constraint f_O–C; (b) rebound free energy ΔG_reb; (c) distance between Fe and O(oxo) atoms and Fe–O(oxo)–H and Fe–O–C angles.

The rebound step involves the formation of a covalent C–O bond between two radical species, CH₃˙ and Fe(–OH˙), which is energetically favourable, and the simultaneous change of the hybridisation of the C atom in CH₃˙ from sp² in the radical species to sp³ in the product molecule CH₃–OH. The latter process is associated with a change in the geometry of the –CH₃ group from planar to (distorted) tetrahedral. The mean force profile (Fig. 6(a)) exhibits a non-trivial dependence on the C–O(oxo) distance. Starting from a value f_O–C ≃ 0 at large C–O(oxo) separations, the mean force initially increases, which indicates a weak repulsive interaction between C and O(oxo). This repulsion is caused by the initial orientation of the Fe–O–H axis, which is pointing towards the C atom of CH₃˙ (cf.Fig. 6(c)). In order for the latter to approach the O(oxo) atom, a rearrangement of the Fe–O–H bond has to occur, in which the angle θ(Fe–O(oxo)–H) decreases from 143.4° (at 2.8 Å) to 117.8° at 2.5 Å and the direction of CH₃˙ attack (the θ(Fe–O–C) angle) varies from 148.4° to 158.3°. This rearrangement also induces a lengthening of the Fe–O(oxo) bond from 1.8 to 1.9 Å. This process has an associated free energy barrier of 6.2 kJ mol⁻¹ (Fig. 6(b)). On the other hand, Fig. 7 shows that CRYSTAL17 predicts an enthalpy energy barrier ΔH_reb of 16.3 kJ mol⁻¹ at d_O–C = 2.3 Å given that for shorter distances the formation of the O–C bond is thermodynamically favourable. i.e., no extra energy needs to be supplied to the CH₃˙ radical to move closer to the reactive oxygen. Hence, the activation entropy ΔS associated with rebound is 33.8 J (mol⁻¹ K⁻¹).

Fig. 7 The difference of the total (ΔE_DFT+disp), DFT (ΔE_EDFT), and long-range dispersion (Δ_disp) energies with respect to their values at 4.8 Å, and the distance between the incoming reactive Ooxo and the Fe(IV) (d_FeIV–Ooxo) vs. the distance between the O(oxo) and radical's carbon obtained with geometry relaxations made with CRYSTAL17 and the functional HSE06 during the rebound step.

Fig. 8 shows the distribution of the down component of the spin density at large O–C separation (a) and after the formation of the C–O bond (b). Initially, a spin down electron is localised on the C atom of the CH₃˙ group. Once the bond is formed, this contribution disappears, as the electron is paired with a spin up electron in the C–O bond.

Fig. 8 Isosurfaces representing 95% of the total spin down density at distances d_O–C of 3.0 Å (a) and 1.4 Å (b) during the rebound step.

Further electronic structure changes are shown in Fig. S4,† which illustrates the evolution of the Mulliken populations of the ↑ and ↓ spin orbitals, net charges, and spin moments during the approach of CH₃˙ to O(oxo). All these quantities experience a sudden change upon the readjustment of the H atom. During this readjustment, the ↑ and ↓ orbital populations and spin moments of O(oxo) appear to switch with those of the incoming C atom. However, after this rearrangement, a spin distribution similar to that observed at large O–C distances is recovered. We attribute these sudden changes in the spin populations of the C and O(oxo) atoms to spin contamination, which may affect our DFT/HSE06 calculations. Indeed we find that the deviation of 〈S²〉 from the ideal single-determinantal value amounts 0.1%. Although this value is small, we believe it to be significant. Since our simulation cell contains 24 Fe atoms, it is conceivable that the largest contribution to the spin-contamination responsible for this value arises from the reactive Fe–O group, whereas all the other FeO(oxo) sites remain virtually unaffected. At variance with Fe(II) containing MOF-74, in which magnetic coupling is observed between first- and second-neighbour metal centres,⁶⁶ our previous work indicates the existence of negligible magnetic interactions between Fe(IV)O groups in MOF-74,²¹ at least within the B3LYP approximation to DFT in the absence of spin–orbit coupling. Based on these considerations, we believe that the sudden changes in the spin and Mulliken populations of C and O as the distance between the two atoms decreases is related to implicit limitations of DFT within the HSE06 approximation. We cannot however fully support this hypothesis, as broken-symmetry DFT calculations, which could be used to pinpoint the origin of the spin contamination, are unfeasible within the context of AIMD with a mixed Gaussian-planewave basis set approach. Under the assumption that spin contamination increases the value of the mean force, and, therefore, of the corresponding free energy, we consider the calculated free energy barrier of 6.2 kJ mol⁻¹ (Fig. 6) to be an upper bound for the rebound step free energy.

After the first complex is formed, the angle θ(Fe–O(oxo)–H) does not recover its initial value of ca. 140° as it increases to 130.3° at 2.4 Å. This indicates that the H atom bonded to O(oxo) moves away from the path of approach of the C atom of CH₃˙. Furthermore, this angle decreases as the O–C shortens. The O–CH₃ bond can be considered formed at ca. 1.5 Å. During the formation of the bond, the Mulliken analysis reveals that Fe(IV) and O(oxo) transfer part of their spin moment to C. In particular, the spin moments of O(oxo), C, and H approach zero as the product is formed. This is consistent with the formation of a neutral and spin-unpolarised CH₃OH molecule.

3.3 Detachment of the methanol molecule from the catalytic site: AIMD vs. static approximation

Fig. 9 shows the mean force computed via AIMD at 300 K by constraining the distance between Fe and the O atom of methanol (d_Fe–O(oxo)) within the range 2.1–4.1 Å. The mean force initially decreases, owing to the restoring attractive interaction between Fe and O, to reach a minimum at 2.5 Å. For larger separations, this force evolves linearly to approach zero at 3.7 Å. This point can be considered as the distance at which the Fe–O(oxo) attraction becomes insufficient to prevent the free diffusion of the CH₃OH molecule away from the reaction site. By integrating the mean force as in previous sections, we calculate the free energy profile for the CH₃OH detachment, which exhibits a plateau at around 3.5 Å. From the value of the constraint at which the mean force vanishes, we then estimate a free energy of detachment ΔG_det of 51.1 kJ mol⁻¹.

Fig. 9 Mean force of constrain f_Fe–O(oxo) and free energy profile for the methanol detachment from Fe as functions of the distance between the Fe atom and the methanol O atom (O(oxo)) computed by thermodynamic integration at 300 K.

Static calculations carried out with CRYSTAL17 display a very similar energy profile for the detachment of the methanol molecule from the catalytic site (Fig. 10). From these calculations, we estimate the enthalpy of detachment, ΔH_det, from the work required to induce the separation of the methanol molecule from Fe (given by the energy difference between the plateau and the initial minimum), which amounts to 130.3 kJ mol⁻¹. For this final step of the oxidation reaction of methane to free methanol within an MOF cavity, our calculations therefore indicate a crucial role of dynamic (entropic) effects at 300 K. These are responsible for reducing the desorption energy from ΔH_det = 130.3 kJ mol⁻¹ to ΔG_det = 51.1 kJ mol⁻¹. From this value, we estimate an entropic contribution to the methanol detachment of 263.9 J (mol⁻¹ K⁻¹), which is the highest value for all steps of the oxidation of methane examined in this work. This large entropic contribution is likely to be explained by the high mobility of the methanol molecule, which is bound to the reduced Fe center by electrostatic forces (e.g. charge–dipole interactions), and can therefore access a large number of rotational, vibrational and translational configurations at finite temperature.

Fig. 10 The difference of the total (ΔE_DFT+disp), DFT (ΔE_EDFT), and long-range dispersion (Δ_disp) energies with respect to their values at 2.1 Å, and the angle between the methanol O–H and O–C bonds (θ_C–O–H) vs. the distance between the O(oxo) and Fe(IV) obtained with geometry relaxations made with CRYSTAL17 and the functional HSE06 during the detachment step.

Recent studies suggest that the production of methanol in confined environments, such as MOFs or zeolites, is highly dependent on the polar behaviour of the contact of this molecule with the support walls, which affects the selectivity of the catalyst in the conversion of methane to methanol in preference to other products such as CHOOH and CO₂.^67–69 We believe that the relatively high (free) energy barrier for the detachment of methanol can be related to the polar behavior of the Fe(III)-MOF74 environment, which can explain the significant effects of the long-range electrostatic forces on the CH₃OH detachment. In this case, further oxidation to CHOOH and CO₂ could occur at the reactive sites. Further work is however required to understand the effect of the MOF cavity polarity on the product yields.

We summarise in Fig. 11 the energetic data for the three steps of the oxidation of methane to methanol. Entropic effects appear to be significant in all three stages and, as noted above, especially relevant during the detachment of the methanol molecule from the Fe site. In addition, Fig. 11 shows that the rate determining step of the overall reaction is the desorption of methanol. These findings indicate that temperature and the structure of the MOF cavity (which can both contribute to enhance the entropic contribution to the detachment energy) can be crucial factors in driving the oxidation process. They also influence the ability of oxidants reaching the reduced Fe(II) site (e.g. H₂O₂ or O₂) to regenerate the reactive Fe(IV)O species. The latter process is of particular interest, considering recent DFT results suggesting that the N₂O decomposition to form ferryl species is the limiting step in several metal–organic frameworks, with estimated activation enthalpies as large as 140 kJ mol⁻¹.³⁷

Fig. 11 Calculated enthalpies (ΔH), free energies (ΔG) and entropies (ΔS) for the H-abstraction, rebound, and detachment steps of the oxidation of methane to methanol.

In summary, according to the standard model of oxidation reactions catalysed by Fe(IV)oxo species, the oxidation of a methane molecule to methanol,

involves (a) a proton-coupled electron transfer reaction to abstract an H atom from methane to generate a methyl radical,

CH₄ + Fe(IV)O → CH₃˙ + Fe(III)O–H, (5)

which, according to our simulations, requires a free energy ΔG_abs = 38.3 kJ mol⁻¹ at 298 K; (b) a CH₃˙ rebound step

CH₃˙ + FeO(III)–H → CH₃OH–Fe(III) (6)

with a calculated free energy ΔG_reb = 6.2 kJ mol⁻¹ (as an upper bound); (c) the detachment of the CH₃OH molecule from the Fe catalytic site

CH₃OH–Fe(III) → CH₃OH + Fe(III) (7)

with ΔG_det = 51.1 kJ mol⁻¹. Therefore, we estimate a total free energy for the process (4) to occur at 298 K ΔG_tot = ΔG_abs + ΔG_reb + ΔG_det = 95.6 kJ mol⁻¹.

Finally, we compare our values with available data on previous works that have studied the methane oxidation to methanol by molecular Fe(IV)oxo complexes. For instance, reaction barriers for the hydrogen abstraction by [FeO(H₂O)₅]₂₊ was calculated with DFT 29 kJ mol⁻¹ in vacuum and 92 kJ mol⁻¹ when solvated in water.²² In addition, DFT work with gas-phase complexes of composition [FeO(H₂O)_n(L)_5−n]²⁺ (n = 4, 1, 0) by substitution of ligand water molecules with L = NH₃, CH₃CN, H₂S and BF₃ predicted enthalpy barriers of 54, 111 and 103 kJ mol⁻¹ for axial (n = 4), equatorial (n = 1) and fully-ammoniated (n = 0) orientations, respectively.²³ If we compare these values with predicted enthalpy barrier of 30.9 for the abstraction, we find that solid-state systems such as MOFs can be promising catalysts for the conversion of methanol into methanol with the advantage of not depending on other factors such as solvent and pH.

4 Conclusions

We have used hybrid ab initio molecular dynamics and DFT calculations to model the conversion of methane into methanol catalysed by Fe(IV)O moieties in MOF-74. The calculations show that this reaction occurs with modalities similar to analogue reactions in gas-phase and in water solution. The H-abstraction step occurs through a proton-coupled electron transfer from the substrate to the O(oxo) atom of the Fe(IV)O group, which we model using two methods: a sequence of geometry optimisations, to compute the enthalpy of reaction, and constrained molecular dynamics simulations, to compute the corresponding free energies. We estimate an enthalpy barrier of 30.9 kJ mol⁻¹ and a free energy barrier of 38.2 kJ mol⁻¹ at 300 K. These values indicate that entropic effects have a substantial contribution in the reaction, which remains however below the estimates for the same reaction carried out in water solution. Modelling the rebound step leading to the formation of a CH₃OH molecule with similar methods is complicated by spin-contamination effects. We could nonetheless estimate a free energy barrier not larger than ca. 6 kJ mol⁻¹ for this reaction and an enthalpy of 16.3 kJ mol⁻¹. Finally, the detachment of the product CH₃OH from the reactive site requires the highest free energy of 51.1 kJ mol⁻¹ and an enthalpy of 130.3 kJ mol⁻¹. Overall, free energy cost involved in the formation of a free CH₃OH inside a MOF-74 cavity from a CH₄ molecule accounts for 95.6 kJ mol⁻¹, a value that is obtained by summing all three free energies ΔG_abs, ΔG_reb, and ΔG_det. The corresponding activation entropies for all three steps are −24.3, 33.8, and 263.9 J (mol⁻¹ K⁻¹), revealing the influence of dynamic structural disorder in the desorption of the methanol molecule from the active site. Given the novelty and computational cost of our approach, we cannot at this point verify whether these changes in entropy, and hence in reaction free energies, are also observed in other systems. However, we argue that they can likely be a general phenomenon, observable also in the oxidation of other light alkanes, potentially with significant changes in the values of the entropy caused by the larger volume of the molecular product of the conversion.

Although the overall free energy of reaction remains higher that typical estimates for reactions assisted by solvent molecules in water solution, it indicates that FeO(oxo) moieties synthesised within MOF-74 can have a substantial potential as oxidation catalysts. Our work shows that the initial proton coupled electron transfer to Fe(IV)O and the final detachment of the product molecule from the reactive site are the stages that need to be optimised in order to further enhance the catalytic properties of FeO(oxo) moieties in MOF-74.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was supported in part by the University of Pittsburgh Center for Research Computing through the resources provided and the computing resources provided by STFC Scientific Computing Department's SCARF cluster.

References

1. Y. Dong, L. J. Que, K. Kauffmann and E. Muenck, J. Am. Chem. Soc., 1995, 117, 11377–11378.
2. B. Meunier, Biomimetic Oxidations Catalyzed by Transition Metal Complexes, Imperial College Press, World Scientific Publishing Co., 2000.
3. J. C. Price, E. W. Barr, B. Tirupati, J. M. Bollinger and C. Krebs, Biochemistry, 2003, 42, 7497–7508.
4. M. Costas, M. P. Mehn, M. P. Jensen and L. Que, Chem. Rev., 2004, 104, 939–986.
5. B. Meunier, S. P. De Visser and S. Shaik, Chem. Rev., 2004, 104, 3947–3980.
6. L. Que and W. B. Tolman, Nature, 2008, 455, 333–340.
7. E. D. Bloch, W. L. Queen, R. Krishna, J. M. Zadrozny, C. M. Brown and J. R. Long, Science, 2012, 335, 1606–1610.
8. P. Verma, X. Xu and D. G. Truhlar, J. Phys. Chem. C, 2013, 117, 12648–12660.
9. A. N. Biswas, M. Puri, K. K. Meier, W. N. Oloo, G. T. Rohde, E. L. Bominaar, E. Muenck and L. Que Jr, J. Am. Chem. Soc., 2015, 137, 2428–2431.
10. P. Verma, K. D. Vogiatzis, N. Planas, J. Borycz, D. J. Xiao, J. R. Long, L. Gagliardi and D. G. Truhlar, J. Am. Chem. Soc., 2015, 137, 5770–5781.
11. L. Bernasconi, A. Kazaryan, P. Belanzoni and E. J. Baerends, ACS Catal., 2017, 7, 4018–4025.
12. G. Olivo, O. Cussó, M. Borrell and M. Costas, JBIC, J. Biol. Inorg. Chem., 2017, 22, 425–452.
13. F. Neese and E. I. Solomon, J. Am. Chem. Soc., 1998, 120, 12829–12848.
14. B. Ensing, F. Buda, P. Blöchl and E. J. Baerends, Angew. Chem., Int. Ed., 2001, 40, 2893–2895.
15. B. Ensing, F. Buda, M. C. Gribnau and E. J. Baerends, J. Am. Chem. Soc., 2004, 126, 4355–4365.
16. F. Neese, J. Inorg. Biochem., 2006, 100, 716–726.
17. F. Buda, B. Ensing, M. C. Gribnau and E. J. Baerends, Chem. – Eur. J., 2001, 7, 2775–2783.
18. F. Buda, B. Ensing, M. C. Gribnau and E. J. Baerends, Chem. – Eur. J., 2003, 9, 3436–3444.
19. B. Ensing, F. Buda, P. E. Blöchl and E. J. Baerends, Phys. Chem. Chem. Phys., 2002, 4, 3619–3627.
20. S. Shaik, H. Hirao and D. Kumar, Acc. Chem. Res., 2007, 40, 532–542.
21. F. Saiz and L. Bernasconi, Phys. Chem. Chem. Phys., 2019, 21, 4965–4974.
22. M. J. Louwerse and E. J. Baerends, Phys. Chem. Chem. Phys., 2007, 9, 156–166.
23. L. Bernasconi, M. J. Louwerse and E. J. Baerends, Eur. J. Inorg. Chem., 2007, 2007, 3023–3033.
24. L. Bernasconi and E. J. Baerends, J. Am. Chem. Soc., 2013, 135, 8857–8867.
25. D. J. Xiao, E. D. Bloch, J. A. Mason, W. L. Queen, M. R. Hudson, N. Planas, J. Borycz, A. L. Dzubak, P. Verma and K. Lee, et al. , Nat. Chem., 2014, 6, 590.
26. H. Hirao, W. K. H. Ng, A. M. P. Moeljadi and S. Bureekaew, ACS Catal., 2015, 5, 3287–3291.
27. P. Liao, R. B. Getman and R. Q. Snurr, ACS Appl. Mater. Interfaces, 2017, 9, 33484–33492.
28. B. L. Suh and J. Kim, J. Phys. Chem. C, 2018, 122, 23078–23083.
29. M. Barona and R. Q. Snurr, ACS Appl. Mater. Interfaces, 2020, 12, 28217–28231.
30. F. Goltl, C. Michel, P. C. Andrikopoulos, A. M. Love, J. Hafner, I. Hermans and P. Sautet, ACS Catal., 2016, 6, 8404–8409.
31. S. Hyeon, Y.-C. Kim and J. Kim, Phys. Chem. Chem. Phys., 2017, 19, 21132–21139.
32. A. S. Rosen, J. M. Notestein and R. Q. Snurr, ACS Catal., 2019, 9, 3576–3587.
33. F. Saiz and L. Bernasconi, Phys. Chem. Chem. Phys., 2020, 22, 12821–12830.
34. F. Göltl and J. Hafner, J. Chem. Phys., 2012, 136, 064501.
35. F. Göltl and J. Hafner, J. Chem. Phys., 2012, 136, 064502.
36. F. Göltl and J. Hafner, J. Chem. Phys., 2012, 136, 064503.
37. J. G. Vitillo, A. Bhan, C. J. Cramer, C. C. Lu and L. Gagliardi, ACS Catal., 2019, 9, 2870–2879.
38. J. T. Groves and M. Van der Puy, J. Am. Chem. Soc., 1974, 96, 5274–5275.
39. J. T. Groves and G. A. McClusky, J. Am. Chem. Soc., 1976, 98, 859–861.
40. H. Iikura, T. Tsuneda, T. Yanai and K. Hirao, J. Chem. Phys., 2001, 115, 3540–3544.
41. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
42. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
43. S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys., 1980, 58, 1200–1211.
44. C. R. Groom, I. J. Bruno, M. P. Lightfoot and S. C. Ward, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 171–179.
45. D. S. BIOVIA, Materials Studio Suite, v.5.5.2, 2017.
46. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364.
47. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
48. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25.
49. J. VandeVondele and M. Sprik, Phys. Chem. Chem. Phys., 2005, 7, 1363–1367.
50. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
51. M. Krack, Theor. Chem. Acc., 2005, 114, 145–152.
52. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
53. E. Carter, G. Ciccotti, J. T. Hynes and R. Kapral, Chem. Phys. Lett., 1989, 156, 472–477.
54. M. Sprik and G. Ciccotti, J. Chem. Phys., 1998, 109, 7737–7744.
55. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
56. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695.
57. R. Dovesi, R. Orlando, A. Erba, C. M. Zicovich-Wilson, B. Civalleri, S. Casassa, L. Maschio, M. Ferrabone, M. De La Pierre, P. D'Arco, Y. Noël, M. Causà, M. Rérat and B. Kirtman, Int. J. Quantum Chem., 2014, 114, 1287–1317.
58. Y.-X. Yu, J. Mater. Chem. A, 2013, 1, 13559–13566.
59. D. Marx and J. Hutter, Ab initio molecular dynamics: basic theory and advanced methods, Cambridge University Press, 2009.
60. F. Saiz, N. Quirke, L. Bernasconi and D. Cubero, Chem. Phys. Lett., 2016, 664, 143–148.
61. F. Saiz and N. Quirke, Phys. Chem. Chem. Phys., 2018, 20, 27528–27538.
62. F. Saiz, D. Cubero and N. Quirke, Phys. Chem. Chem. Phys., 2018, 20, 25186–25194.
63. S. Choomwattana, T. Maihom, P. Khongpracha, M. Probst and J. Limtrakul, J. Phys. Chem. C, 2008, 112, 10855–10861.
64. M. Swart, A. R. Groenhof, A. W. Ehlers and K. Lammertsma, J. Phys. Chem. A, 2004, 108, 5479–5483.
65. G. Petersson, A. Bennett, T. G. Tensfeldt, M. A. Al-Laham, W. A. Shirley and J. Mantzaris, J. Chem. Phys., 1988, 89, 2193–2218.
66. R. Maurice, P. Verma, J. M. Zadrozny, S. Luo, J. Borycz, J. R. Long, D. G. Truhlar and L. Gagliardi, Inorg. Chem., 2013, 52, 9379–9389.
67. Z. Jin, L. Wang, E. Zuidema, K. Mondal, M. Zhang, J. Zhang, C. Wang, X. Meng, H. Yang and C. Mesters, et al. , Science, 2020, 367, 193–197.
68. J. Kang and E. D. Park, Catalysts, 2020, 10, 299.
69. X. Yu, B. Wu, M. Huang, Z. Lu, J. Li, L. Zhong and Y. Sun, Energy Fuels, 2021, 35, 4418–4427.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1cy00551k

---