QM/MM study of the reactivity of zeolite bound methoxy and carbene groups †

The conversion of methanol-to-hydrocarbons (MTH) is known to occur via an autocatalytic process in zeolites, where framework-bound methoxy species play a pivotal role, especially during catalyst induction. Recent NMR and FT-IR experimental studies suggest that methoxylated zeolites are able to produce hydrocarbons by a mechanism involving carbene migration and association. In order to understand these observations, we have performed QM/MM computational investigations on a range of reaction mechanisms for the reaction of zeolite bound methoxy and carbene groups, which are proposed to initiate hydrocarbon formation in the MTH process. Our simulations demonstrate that it is kinetically unfavourable for methyl species to form on the framework away from the zeolite acid site, and both kinetically and thermodynamically unfavourable for methyl groups to migrate through the framework and aggregate around an acid site. Formation of carbene moieties was considered as an alternative pathway to the formation of C–C bonds; however, the reaction energy for conversion of a methyl to a carbene is unfavourable. Metadynamics simulations help confirm further that methyl species at the framework acid sites would be more reactive towards formed C 2+ species, rather than inter-framework migration, and that the role of carbenes in the formation of the first C–C bond will be via a concerted type of mechanism rather than stepwise.


Introduction
5][6] The methanol-to-hydrocarbons (MTH) process is a promising technology for the conversion of methanol to fuel and light olefins, and is now used on industrial scale; 3,7 however, controlling the product selectivity and catalyst deactivation remains a challenge. 3,7,8In order to address these challenges, the catalytic mechanism must be understood.The MTH process is known to have an induction period that is characterized by the low reactivity of methanol until certain ''hydrocarbon pool'' (HP) species are formed, which themselves act as a co-catalyst to produce C-C species, [8][9][10] but even the initial steps during the induction period remain unclear.
Several mechanisms have been proposed for the formation of the first C-C bond, which involve stable compounds (dimethyl ether, methane, formaldehyde) or short life-time intermediates (trimethyl oxonium, carbene).Formaldehyde and methane are present in the reaction mixture and thus have been considered as reaction intermediates when forming the first C-C bond, particularly in the case of ethanol formation, but a very high reaction barrier (183 kJ mol À1 ) was calculated which led to the discounting of any mechanism involving these intermediates. 6][13] Another proposed mechanism involves a concerted process whereby abstraction of a hydrogen, from the methyl group of a methanol or DME molecule, by a basic zeolite oxygen, would allow a C-C bond formation with another methanol, DME, trimethyloxonium or zeolite bound methoxy group. 11omputational ONIOM techniques concluded that the energy barrier or the breaking of the covalent C-H bond by the available basic site is considerable (4200 kJ mol À1 ), which makes this mechanism seem energetically unlikely. 6,14Several other studies 15,16 report that pure methylated zeolite and SAPO frameworks (CH 3 -ZSM-5, CH 3 -Y, CH 3 -SAPO-34) have the potential to form a wide range of hydrocarbons (paraffins, olefins, aromatics) once heated over 573 K. Recent studies, 17 employing synchrotron based IR techniques, suggested strongly that carbene (:CH 2 ) moieties form from methyl groups, via the deprotonation of a zeolite bound methoxy, followed by either polymerisation to olefins or insertion in to a methanol or DME. 11Theoretical calculations 14,18 confirmed that the direct activation energy for carbene formation from methoxy is high (245; 14 326 kJ mol À1 18 ), whilst experimental H/D exchange studies showed that C-D bond breaking occurs in H-ZSM-5 17 but not in H-SAPO-34, 19 which does not completely rule out C-H bond cleavage but rather indicates a limited C-C interaction.Experimental studies involving CH 2 N 2 also showed that the presence of H-ZSM-5 enhances the conversion of CH 2 N 2 to ethene, further indicating that carbene could stabilise in H-ZSM-5. 20,21hese current observations provoke several important questions regarding the reaction mechanisms, especially with consideration of the reactivity of H-ZSM-5.A first question concerns the mobility of C 1 species within zeolites; any given acid can only convert one methanol, but the framework-bound moieties have to be in the vicinity of each other to react and form higher order hydrocarbons.In this study, we have used computational modelling to examine two potential routes for C 1 (methyl or carbene) mobility: (i) the methylation of an Si-O-Si basic site, followed by migration towards a methylated Al-O-Si site; and (ii) the direct migration of a methyl group away from the Al-O-Si active site to a second methylated active site.A second question concerns the conditions that may lead to carbene formation; to address this challenge, we have developed mechanistic models for stable carbene formation and migration on the zeolite framework, which is discussed following the summary of the computational methods employed.

QM/MM Model description
The hybrid quantum mechanical/molecular mechanical (QM/MM) approach, available in the Chemshell package, 22 is used to calculate energetics, having been used previously in modelling studies that investigated mechanistic processes for catalysis in zeolites. 23,24he approach ensures the correct long-range electrostatics for the active site of interest, whilst also removing the periodic boundary conditions that can hinder charge non-neutral models and/or high-level calculations of large unit cell systems, such as zeolites.Spherical embedded-cluster models of H-ZSM-5 were created from the experimental unit cell of siliceous MFI. 25 The embeddedcluster models were centred on a Si tetrahedral (T-)site of interest; MFI has 12 symmetry inequivalent T-sites, and we considered the Al substituent in the T12 position, so as to be at the intersection site of the sinusoidal and straight channels in MFI.The Al positions are difficult to characterise in experiment, and so their positioning in simulations introduces an aspect of selective bias, but some Al positions are commonly assumed more favourable than others, with the T12 site often reported to be substituted. 26here Al was substituted into the silicate framework, a chargecompensating hydrogen atom is added on an adjacent oxygen in the O26 position, in between two T12 sites, which is referenced as O Al throughout the manuscript; also considered was an O Al* site, where the H was added on the O8 position, between a T12 and a T3 site, with notations as per IZA database. 77In all cases, the Brønsted proton is directed towards the centre of the supercage as this configuration has the highest deprotonation energy, 27 i.e. most stable, and thus is assumed to be predominant.
The embedded-cluster models are partitioned with a central QM region, which is the chemically active part of our model, and an encapsulating outer MM region, where long-range structural and electrostatic effects are present to ensure correct bulk representation for the high-accuracy QM calculations.All QM atoms were unconstrained during geometry optimisation, whereas two concentric domains were employed in the MM region: an inner MM region where unconstrained atoms can move during a geometry optimisation; and an outer MM region that is frozen to ensure a bulk-like long-range structure, as shown in Fig. S1 of ESI.† In our calculations, the inner and outer MM regions extend from the central T-site to a radius of 10.58 Å (20 a 0 ) and 21.17 Å (40 a 0 ), respectively, and the QM region used herein includes atoms up to the fifth nearest neighbour (the third oxygen atom) from the central T-site.In its entirety, the total number of atoms in the cluster model of the active site is 2165, with 74 QM atoms and 197 inner MM atoms.During QM calculations, the terminal oxygen atoms at the edge of the QM region are saturated with hydrogen atoms: these artificial ''link'' atoms ensure correct electronic structure for the terminal atoms and do not inadvertently affect the electronic solution of the overall QM calculations, as a bonddipole correction is added at the boundary to ensure an accurate electrostatic embedding potential. 28hroughout, the QM energy has been calculated using hybrid-DFT with the dispersion corrected Becke 97-3 exchange-correlation (XC) functional, commonly referred to as B97-D, 29,30 as provided in the GAMESS-UK code. 31The atomic orbitals are represented using the Ahlrichs and Taylor TZVP Gaussian basis sets. 32The self-consistent field (SCF) convergence criteria was set to an energy change of less than 2.72 Â 10 À6 eV (1 Â 10 À7 Hartrees) between SCF iterations. 33,34he MM energy was calculated using DL_POLY, 35 employing the forcefield of Hill and Sauer, 36,37 with the coordination dependent charges in the original forcefield replaced with fixed point charges of 1.2 and À0.6 e for silicon and oxygen, respectively, as parameterised in the work of Sherwood et al. 28 As we have a neutrally charged system, we employed Restricted Hartree-Fock (RHF) conditions in our simulations, corresponding to all spins being paired and overall singlet spin multiplicity.
Geometry optimizations were performed using the Chem-Shell package, 22 in a Cartesian coordinate space, using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm.9][40][41] Vibrational frequencies were calculated with a task-farmed finite-difference approach, 31 allowing us to confirm that geometries correspond to local minima. 42,43The geometrically optimised models were then used to calculate the individual atomic charges using a Mulliken analysis. 44The transition state energies were determined by employing the Nudged Elastic Band (NEB) method, 45 in the task-farmed mode, with the reaction path represented by 5 images.For the NEB calculations, only the adsorbate and framework atoms up to the second nearest neighbour framework were free to move, as otherwise a prohibitively large computational cost is encountered; comparison of this approximation against displacement of all atoms in the QM region shows a small difference in the transition state energies (o10 kJ mol À1 ), which we consider as the error bar going forward.When comparing the energetics involved in the conversion of methanol to methoxy, two models were employed: one where methanol was adsorbed on the T12 site, with methyl pointing towards the active site, and a second with methyl pointing away from the active site.In the case of methanol adsorbed with methyl oriented away from the Al central T-site, the model was reconstructed to maintain the same configuration but with the Al placed at the T11 site, next to the T12 site, i.e. the methanol remained central to the embedded-cluster.This adjustment ensured the optimised methanol had exactly an equivalent QM region and electrostatics.A similar approach was taken for modelling the product state, which consists of methyl and water.
The energy of the deprotonated zeolite is corrected by the addition of the Jost correction, 46,47 which accounts for the truncation of the MM polarisation at the end of the first (flexible) MM region.Deprotonation leaves a charge in the QM region of the simulation model (À1), the electrostatic effects of which would be expected to be relatively longranging; however, our model prevents polarisation of atoms in the fixed MM region.The Jost correction provides a good estimate of the polarisation energy and is given by: where Q is the defect charge, R is the radius of the total cluster and e is the dielectric constant of the material, all in atomic units.The dielectric constant for MFI is taken in this work as 3.38, calculated using classical shell model methods. 48

Energy analysis
The adsorption energy (E ads ) of a sorbate to the zeolite framework is calculated as: where, E[ZeOH], E[Sorbate] and E[ZeOH + Sorbate] are the total energy of the zeolite sorbent, the gas-phase sorbate and the combined guest-host system, respectively, each in their optimised geometry.Due to our use of an atom-centred basis set, it is necessary to include a basis-set-superposition-error (BSSE) for the combined system, which is calculated thus: where the term in the first parentheses gives the BSSE (E BSSE ) for the framework when including the sorbate orbitals, and the second parentheses gives the E BSSE for the sorbate in the presence of the zeolite orbitals.Thus, in both parts the BSSE is calculated as the difference in energy of the system components (ZeOH and Sorbate) in an adsorbed geometry (denoted with ads ), with and without the basis functions (denoted as ''Basis'') for the second component of the complete system.e.g.E(ZeOH) is calculated with and without the basis functions of the sorbate present. 27E BSSE is included in all energies reported; generally, the error is r5 kJ mol À1 for a single adsorbed CH 3 OH.
In addition to the methylation reaction (E r ) and activation (E act ) energies reported here, three other thermodynamic energies are used.First the reaction energy for methyl or carbene to migrate from one bonding oxygen site (O A ) to another (O B ) is denoted as E mig , with E C1-Zeo , E C1 , and E Zeo-representing the absolute energies of the methyl or carbene bonded to the zeolite framework, the gas phase methyl (CH 3 + ) or carbene (CH 2 :) fragments, and the energy of deprotonated zeolite models, respectively.Thirdly, we calculate the carbene formation energy as: with E CH2-H and E CH3 representing the absolute energies of the carbene bonded to the protonated zeolite framework model, and methyl bonded to the zeolite framework, respectively.

Metadynamics simulations
1][52] Dynamic simulations provide a closer replica of the experimental conditions that lead to methoxy conversion to hydrocarbons in experiment.
Although the first hydrocarbons are detected after heating at 523 K, 16 the simulated temperature was set to 623 K to ensure no heat transfer effects influenced the reactivity.Ab initio molecular dynamics (MD) and metadynamical simulations were performed on the combined system, using 3-dimensional periodic boundary conditions, with the CP2K simulation package (version 6.1). 53The dynamics of the nuclei were governed by the Newtonian equations of motion, in which the potential from the Born-Oppenheimer electronic ground state is inserted.The self-consistent field (SCF) energy was evaluated with DFT using the revPBE functional 54 with Grimme D3 dispersion corrections 55 and the Gaussian Plane Waves method 56 that uses Gaussians as basis sets (DZVP-GTH 57 ) and planewaves (320 Ry cut-off) as an auxiliary basis.The SCF convergence criterion was set to an energy change of less than 1 Â 10 À5 Hartrees between SCF iterations.In both NPT ''production runs'' and the NVT metadynamics simulations, the integration time step was set to 0.5 fs.The initial geometry for the NVT metadynamics simulations was taken as the final snapshot of the NPT calibration models, themselves performed on empty zeolite models, with methoxy and carbene inserted into the zeolite pores.The cell parameters, presented in Table S1 of ESI, † were determined from a preliminary isothermal-isobaric (NPT) ensemble simulation of 50 ps on the zeolite cells containing methoxy or carbene, where the number of atoms, temperature (623 K) and external pressure (1 atm) are controlled by a chain of five Nose ´-Hoover thermostats 58,59 and pressure by the Martina-Tobias-Klein barostat. 602][63] Specifically, during the NVT MTD run, two geometric parameters selected to uniquely describe each reaction state, also named collective variables (CVs), were defined by coordination numbers (CN): 61 CNði; jÞ ¼ X in which r ij is the distance between bonded atoms i and j.
The parameters n and m were set to 6 and 12, respectively.
The reference distance, r 0 , was chosen to be similar to the transition state distance between atoms i and j (2.0 Å).The biased CV is a combination of two CVs, calculated as (CV2-CV1), which permits faster computation without compromising accuracy.For methyl migration, the width of the Gaussian hills is set to 0.02, and the height is initially 5.0 kJ mol À1 .Once the transition state has been identified and crossed twice, the height of the Gaussian hills is halved, which allows more accurate sampling of the activation barrier; this process is repeated until a final hill height of 0.65 kJ mol À1 is used, thus ensuring a refined representation of the energy surface.The metadynamics simulations were deemed complete when a change in the free energy barrier was equal or less than 5 kJ mol À1 between every 500 energy hills added, with the reported errors determined from the minimum and maximum barriers calculated at this point.
During initial MTD simulations, carbene commonly diffused away from the active site to the centre of the pore, or reacted with a neighbouring oxygen to form formaldehyde, and thus the sampling process was limited to Gaussian hills of 1.25 kJ mol À1 height and 0.02 width throughout, with the activation barrier being crossed once only.In addition, to ensure sampling of chemically relevant space with the MTD simulations, restrictions were applied to prevent the carbene from migrating to adjacent oxygen sites connected to Al.In particular, a series of singlesided energy ''walls'' were applied that extend from the barrier (B) towards larger values of the collective variable (CV3), represented by a quadratic potential k(CV-B) 2 , with k being a quadratic potential constant.The barrier is applied at CN(C CH 2 , O Al ) = 0.036 (k = 100 Ha), which corresponds to a bond distance of 2 Å.The FES of the reaction was reconstructed from a MTD calculation based on the sum of the spawned Gaussian hills; the free energies of reaction were then calculated as the difference between the lowest energy of the reactant state and the highest of the transition stage.Further details on the methodology and case studies are provided in previous studies. 61,64,65he metadynamics observables, specifically, free energy activation barriers for the forward (DF forward ) and reverse (DF reverse ) reactions are determined as: DF reverse = F product À F TS (9)  with F TS , F reactant and F product represent the highest energy of the transition state, and lowest energies of the reactant and products states, respectively, on the MTD mapped free energy surface.

Methyl formation and migration
Experimental IR and neutron spectroscopic evidence [66][67][68][69][70][71] suggest that methanol insertion in the zeolite catalyst and adsorption on the Brønsted acid site can lead to the disappearance of the zeolite O-H stretch vibrational band, and the appearance of the bands specific to the zeolite bound methyl group and water.This series of events is widely accepted as being properly described by the nucleophilic substitution framework methoxylation mechanism.Specifically, once multiple methanol molecules adsorb on the acid site, it is proposed 23,64,[72][73][74] that the Brønsted proton is spontaneously abstracted by the methanol oxygen, forming a methoxonium molecule (CH 3 OH 2

+
) or protonated cluster of methanol molecules and a negatively charged zeolite oxygen that acts as a Lewis basic site.Nucleophilic attack of the latter on the methyl group leads to the formation of water and a methylated framework oxygen on the oxygen neighbouring the Al site.
The mobility of methoxylated species is initially investigated considering the association of methanol, and the subsequent methoxylation reaction at an oxygen site coordinated to an Al site (O Al ), compared with the association of methanol at an oxygen atom with two silicon atom neighbours (O Si ), as shown in Fig. 1.
Initially, methanol was adsorbed on the Brønsted acid site with the methyl group oriented either towards the active site, centred on the Al T-site (Fig. 1, left), for the methyl to transfer on to O Al , and away from the active site, centred on the Si T-site (Fig. 1, right), for a methyl transfer on O Si .The optimised structures are similar for the two configurations, with the distances from the zeolite framework to the Brønsted proton and to the methanol molecule being very similar, as well as the Mulliken charge distribution (Table 1); furthermore, the methanol adsorption energies are similar (À120 and À126 kJ mol À1 , respectively, Table 1).Previous work 23 suggested that the methyl groups do not contribute substantially to the overall adsorption energy of methanol to the framework, with most of the energy attributable to the hydrogen bonding through the hydroxyl group, which the current results support.
The activation energy for methylating an Al-O-Si site was considered via the mechanism illustrated in Fig. 2, and is calculated to be 225 kJ mol À1 , with an overall reaction energy of 49 kJ mol À1 .The activation energy is three times higher to methylate a O Si site, i.e.Si-O-Si, with an energy barrier of 748 kJ mol À1 and an overall reaction energy of 118 kJ mol À1 (Table 1).These results make the direct methylation of an oxygen atom away from the Al site both thermodynamically and kinetically infeasible.
The methylation path (Fig. 2) requires a Brønsted proton to transfer completely from the zeolite active site to the methanol before the dissociation of the C MeOH -O MeOH bond and the formation of water.With the methanol adsorption energies and structures similar for both adsorption configurations considered, one expects the energy to remove the Brønsted proton from the active site to the methanol, in both cases, to be similar in magnitude; thus it may seem surprising that the energy barriers differ so significantly.The initial step of the methylation path, which is the deprotonation of the zeolite, appears not to contribute substantially to the difference in activation barrier between the two cases (which we refer to as methylation ''inside'' and ''outside'' of the active site when at the O Al and O Si sites, respectively); therefore, the significant difference in activation energy is due to energetic instability of the methyl fragment when ''outside'' the environment of the Al, i.e. at an oxygen atom surrounded by silicon, O Si , as opposed to Table 1 Summary of results for methanol adsorbed on active sites with methyl oriented towards aluminium (CH 3 OHÁO Al ) and silicon (CH 3 OHÁO Si ) centres, as shown in Fig. 1, alongside methoxy formation on oxygen closest to aluminium (CH 3 ÁH 2 OÁO Al ) and silicon (CH 3 ÁH 2 OÁO Si ).Adsorption energies (E ads ), reaction energies (E r ) and activation energies (E act ) are given in kJ mol À1 ; geometric distances (d) and Mulliken charge (q) data are presented in Ångstroms and electronic charge (e), respectively.Further structural and electronic data is provided in ''inside'' the environment (i.e.neighbouring) of the aluminium, i.e.O Al .We note, however, that the difference in activation barriers (523 kJ mol À1 ) is much greater than the stability of the methylated products (69 kJ mol À1 ), and therefore a more detailed investigation follows in the next section.

Methyl migration away from the active site
Here we consider the energetics of methyl migration both for isolated methyl groups and for multiple methyl groups at an active site.

Isolated methyl migration.
The energetics associated with the framework-adsorbed methoxy group migrating away from the active site, specifically from O Al to O Si , have been analysed.The bonding sites are presented in Fig. 3.
The methyl group bonds strongest to the zeolite framework through the C-O Al interaction (E bond = À635 kJ mol À1 ) with C-O Si being a slightly weaker interaction (E bond = À510 kJ mol À1 ).The reaction energy for methyl transfer (E mig ) from O Al to O Si , which is calculated as described in eqn (4), is 126 kJ mol À1 .Analysis of the atomic charges and bond lengths for the framework surrounding the methyl moiety are provided in Table 2 (and in more detail in Table S3 of the ESI †).The main structural differences that determine the endothermic reaction energy are over the first nearest neighbours of the methyl group, since the remaining extended structure (to the third nearest neighbour) was     3 and bonding energy (E bond ), given in kJ mol À1 , calculated using equations 4 and 5, respectively, from Section 2.2.Geometric and charge observations are presented in Ångstroms and electronic charge (e), respectively.The methyl bonding sites are highlighted in bold in the charge analysis to aid interpretation To understand how multiple methyl groups might interact, simulations were considered with two methyl groups on the same active site (Fig. 4).The bonding energies of the second methyl groups are B250 kJ mol À1 smaller per methyl added (i.e. less strongly bound) than when there is just a single methyl group, highlighting the relative instability of having two methyls in close vicinity.The addition of a methyl moiety at the active site led to a lower energy (81 kJ mol À1 ) of methyl transfer away from the active site, i.e. ''inside'' to ''outside'', compared to the calculations for an isolated methyl (126 kJ mol À1 ).This result is mainly attributable to the structural differences of the methyl active site of the single and double methyl models, given in Table 3; specifically, the C Al -O Al bond is slightly longer in the double methyl model, by 0.02 Å, which is due to charge-related electrostatic effects.The higher positive charge on aluminium in the double methyl case (1.07 e) compared to the single methyl case (0.94 e), shows that the addition of a second methyl draws more electronic density from aluminium to the neighbouring oxygen, further lengthening the C-O bonds, as highlighted in Table 3, with further details in the ESI † in Table S4.
The overall conclusion of our analysis is that thermodynamic factors will strongly inhibit the formation of isolated methyl groups at sites other than those surrounding the aluminium active site (i.e.O Al and O Al* ), and, as a consequence, migration of methyl groups is unlikely to be prominent in the mechanism for C-C bond formation.Our focus therefore changes to the potential role of carbene species.

Carbene migration
3.3.1.Carbene bonded on aluminium active site.Though carbene species have been proposed in experimental studies, 17 earlier theoretical studies 14,18 showed that the direct formation of carbene from methyl is highly energetically demanding (245; 14 326 18 kJ mol À1 ).In order to determine the conditions that could lead to carbene formation from methyl species, and to understand the possible role of carbenes in the production of hydrocarbons, several new models were analysed as illustrated in Fig. 5 and 6.
The calculated endothermic reaction energies of carbene formation, reported in Table 4, arise mainly from the difficulty of stabilising the carbene fragment on the zeolite framework in a configuration preventing the spontaneous conversion back to methyl (Fig. 5), which was also a challenge reported in other investigations. 18The carbene stabilises between the framework cation and oxygen, bonding to both atoms, in both models considered; the bond  6B).The carbene bonding energy is À174 kJ mol À1 for the carbene bonded closer to Al, illustrated in Fig. 6A, and À141 kJ mol À1 for carbene bonded towards Si (Fig. 6B); these bonding energies are substantially weaker than the methyl bonding discussed in Section 3.2 (from À510 to À635 kJ mol À1 ).
To determine the influence of the Brønsted proton on the strength of the carbene bond to the zeolite, the two models were optimised without the Brønsted proton present (i.e. for the carbene in isolation).The strength of the carbene-framework bond is calculated to be 20 kJ mol À1 weaker in the absence of the proton from the active site, showing how it aids stability.The carbene bonding energy is noted as exothermic throughout, which ensures the carbene remains stable on the zeolite framework; the energy for carbene migration (34 kJ mol À1 ) is also significantly less than for the methyl group in Section 3.1, due to  the weaker interaction of the carbene with the zeolite framework, leading us to conclude that carbenes are potentially mobile in the zeolite.
3.3.2.Carbene bonded on silicon active site.Although there is no experimental evidence for the position that carbene might occupy, for completeness, the formation of carbene from a methyl ''outside'' the active site, i.e. on a Si-O-Si site, was investigated (Fig. 7).
For these configurations, the carbene stabilises either ''outside'' the aluminium active site, or inserts itself in the Si-O bond (Fig. 7 and 8).Whilst the reaction energies are calculated to be endothermic, they are lower than for carbene coordinated ''inside'' the aluminium active site (Table 5).The carbene is most stable when bonded within the zeolite framework (Fig. 8A), with a reaction energy of 103 kJ mol À1 , as opposed to when carbene is positioned outside the lattice (Fig. 8B), with a formation energy of 239 kJ mol À1 (Table 5) similar to previous investigations (245 kJ mol À1 14 ).Since there is no charge neutralisation involved, the positioning of the carbene within the zeolite framework may be a consequence of the framework flexibility.Overall, whilst comparison of the carbene stabilised within the silicate structure (Fig. 8A) with that of carbene being placed closer to the active site (Fig. 8B) shows significant structural and energetic differences, the results are all considerably  Table 4 Summary of results for a carbene (CH 2 ) moiety bonded in the proximity of Al and Si, and also bonded in the proximity of Si and Si*, with reaction energy (E r ), migration energy from Si to Al (E mig ), and bonding energy (E bond ) given in kJ mol À1 , with geometric and charge data presented in Ångstroms and electronic charge (e), respectively.The specific bonding sites are highlighted in bold and illustrated in Fig. 6.A further detailed geometric and electronic set of observables is provided in Table S5    endothermic, which confirms that neither the Al-O-Si or the Si-O-Si coordination-sites enhance the stability of the carbene moiety.
Previous EPR analysis shows that triplet state species could be present in the ''hydrocarbon pool''. 11Thus, the conversion of methoxy to triplet state carbene was also analysed, which led to the determination of a considerably greater carbene reaction energy (526 kJ mol À1 ) than the singlet state carbene (283 kJ mol À1 ).As with the formation energies, the triplet carbene migration energy is almost double (79 kJ mol À1 ) that of the singlet state (34 kJ mol À1 ), which leads to the conclusion that the triplet carbene species is unlikely to play a role in the C-C formation mechanism.
3.3.3.Metadynamics simulations of free energy barriers for C 1 species migration.Since the migration energy for the singlet state carbene (34 kJ mol À1 ) was relatively small, metadynamics simulations were carried out, as detailed in Section 2.4, in order to determine the free energy barriers for translation inside the zeolite framework.For completeness, the methyl migration energy was investigated as well with the same procedure, even though the calculated migration energy from the potential energy surface was considerably large (126 kJ mol À1 ).The free energies of activation obtained for carbene migration (DF forward = 212 kJ mol À1 ; DF reverse = À125 kJ mol À1 ) from the unprotonated active site is considerably smaller than methyl migration (DF forward = 357 kJ mol À1 ; DF reverse = À281 kJ mol À1 ), although still very energetically demanding.By comparing our results to previous studies 56 using a similar technical setup to analyse methyl transfer from the active site to olefins, we find that methyl transferring from the active site to ethene (DF forward = 123 kJ mol À1 ), propene (DF forward = 112 kJ mol À1 ) or trans-2-butene (DF forward = 108 kJ mol À1 ), 61 is more favourable than framework migration, indicating that, in the presence of the hydrocarbon pool, methyl is unlikely to translate but rather to transfer to olefins.Previous studies proposed ethoxyde as an intermediate, 75,76 which could be formed here via a carbene insertion to methoxide, as shown in Fig. 9; however, this reaction route would involve carbene migration, which we have shown to be energetically unfavourable, and thus the formation of ethoxyde via such a carbene/methoxide reaction may be problematic.

Summary and conclusions
The migration and reaction of C 1 species around an active site within the H-ZSM-5 framework was considered by modelling of distinct reaction pathways with accurate QM/MM techniques.The formation of methoxy groups away from the active site is shown to be both kinetically and energetically demanding, leading us to conclude that these species are not present in significant concentrations under reaction conditions in the MTH process.Based on the methanol adsorption energy, we conclude that the rate determining step for methoxylation is the methyl transfer to the active site, rather than the Brønsted proton transfer to methanol.The results suggest that methyl migration would be unlikely within the zeolite; however, due to the small differences in the basicity of the oxygen atoms surrounding the active site needed to stabilise methyl, it could be speculated that a zeolite substituted with gallium or indium may enhance the chances of the methyl forming beyond the active site, due to the higher negative charge present on the adjacent sites.
The direct formation of carbene from methyl is also calculated with our QM/MM techniques, but is also highly energetically demanding; however, due to the strong bonding within the zeolite framework, in the absence of a Brønsted proton, the carbene moiety may have sufficient stability to take part in a concerted reaction with other species.Metadynamics simulations provide evidence for the reactivity of CH 3 under realistic MTH conditions, but our results suggest that the methyl moieties will react with formed C 2+ species rather than migrate along the zeolite framework, further suggesting that CH 3 is unlikely to migrate through the framework to form the initial C-C bond.
Table 5 Summary of results for a carbene (CH 2 ) moiety formed from the methyl conversion on O Si , bonded in the proximity of Si and Si* as illustrated in Fig. 8, with reaction energies given in kJ mol À1 .Geometric and charge data are presented in Ångstroms and electronic charge (e), respectively.For completeness, a further detailed geometric and electronic set of observables is provided in Table S6

Scheme 1
Scheme 1 Illustration of oxonium ylide mechanism via TMO to ethene.

Fig. 1
Fig. 1 Active site models, with methyl moiety of methanol oriented towards aluminium (left) and towards silicon (right); and geometric assignment, with Al, Si, O, C and H atoms represented as purple, yellow, red, green, and white.Methoxy bonding sites are denoted as O Al and O Si ; Si and Si* are also represented, as referenced in the structural analysis.Optimised bond lengths, marked with dashed lines, are given in Ångstroms.
similar in both cases.Analysis of the O Al and of O Si methyl bonded models shows that the C-O Si bond is slightly longer (1.52 Å) than the C-O Al distance of 1.48 Å.In addition, the charge on O Si (À0.45 e) is less than on O Al (À0.51 e), when each are bonded to the methyl, showing a reduced charge transfer between the methyl group and the Lewis basic site of the zeolite when the oxygen bonding site is surrounded by Si atoms; thus confirming the methyl is best stabilised adjacent to the Al T-site.To clarify further the factors influencing the stability of methyl at different framework positions, the energy of transferring a methyl from the O Al site to the most basic other site available, O Al* , was calculated.The methyl bonding energy calculated for O Al* (E bond = À621 kJ mol À1 ) is marginally weaker than on O Al (E bond = À635 kJ mol À1 ), with the C-O Al* bond of 1.48 Å similar in length to the C-O Al distance (1.47 Å).The reaction energy is much lower (14 kJ mol À1 ) compared with that

Fig. 2
Fig.2Methoxylation reaction path, with the colour scheme the same as Fig.1.The Brønsted proton is highlighted in blue for clarity.

Fig. 3
Fig. 3 Methoxy bonding mode, with colour scheme and labels the same as Fig. 1. O Al and O Al* are highlighted as two possible Al-neighbouring O sites ''inside'' the active site.

Fig. 4
Fig. 4 Models with an additional second methyl bonded on O Al (left) and O Si (right), alongside a methyl on the O Al* site.The colour scheme is as per Fig. 1.

Fig. 5
Fig. 5 Illustrations of reactions considered for methyl bonded on O Al when converted to hydrogen and a carbene bonded between O Al and Al (top) and O Al and Si (bottom).

Fig. 6
Fig. 6 Carbene migration from an Al-coordinated site (A) to an Sicoordinated site (B).The key is as per Fig. 1.

Fig. 7
Fig. 7 Carbene formation reactions from an O Si bonded methyl, leading to a framework integrated CH 2 and Brønsted proton (top) or a framework bonded CH 2 and silanol (bottom).

Fig. 8
Fig. 8 Models of carbene moiety bonded on the O Si site, in the proximity of Si (A) and Si* (B), with a key as per Fig. 1.

Fig. 9
Fig. 9 Suggested reaction mechanism involving C 1 species, specifically, carbene insertion to methoxy forming ethoxyde decomposing to ethene. 49 The (CV2-CV1) collective variable was biased by adding Gaussian hills every 25 fs.For the methoxy migration case, CV1 is defined by CN(C Me , O Z1 ), which describes the breaking of the C-O bond of the zeolite bound methyl; CV2 is then defined by CN(C Me , O Z2 ) to describe the subsequent formation of the C-O bond between the product methyl moiety and the bonding site external of the zeolite active site.As will be further discussed, carbene bonds with O and Al/Si, so CV1 is defined by CN(C CH 2 -O Z1 , Al), which describes the breaking of the C-O and C-Al bonds of the zeolite bound carbene; CV2 is then defined by CN(C CH 2 -O Z2 , Si) to describe carbene migration from the active site to the Si-O site.The definition of the collective variables is illustrated schematically in Fig.S2(ESI †).

Table 2
Summary of results for methyl (CH 3 ) group bonded on O Si , O Al , and O Al* sites, with migration energy (E mig ) for migration from O Al to O Si and from O Al to O Al* as illustrated in Fig.
Al to O Si (126 kJ mol À1 ), again illustrating that the methyl group is more stable on the Lewis basic sites that neighbour Al T-sites.In contrast, a more negative charge is present on the methyl carbon when bonded to O Al* (À0.30 e) than when bonded to O Al (À0.25 e), which may contribute to the reduced stability of the methoxy at the O Al* Lewis site.3.2.2.Neighbouring methyl migration.

Table 3
Summary of results when two methyl (CH 3 ) groups are bonded on O Al* and O Al (2CH 3 ÁO Al ) and on O Al* and O Si (2CH 3 ÁO Si ) sites as illustrated in Fig.3, with bonding energy (E bond ), and migration energy (E mig ) from O Si to O Al , given in kJ mol À1 .Geometric and charge observables are presented in Ångstroms and electronic charge (e), respectively.The methyl bonding sites are highlighted in bold This journal is © the 2021 Phys.Chem.Chem.Phys., 23, 17634-17644 | 17641 of ESI of ESI