Modelling metal centres, acid sites and reaction mechanisms in microporous catalysts

We discuss the role of QM/MM (embedded cluster) computational techniques in catalytic science, in particular their application to microporous catalysis. We describe the methodologies employed and illustrate their utility by brie ﬂ y summarising work on metal centres in zeolites. We then report a detailed investigation into the behaviour of methanol at acidic sites in zeolites H-ZSM-5 and H-Y in the context of the methanol-to-hydrocarbons/ole ﬁ ns process. Studying key initial steps of the reaction (the adsorption and subsequent methoxylation), we probe the e ﬀ ect of framework topology and Brønsted acid site location on the energetics of these initial processes. We ﬁ nd that although methoxylation is endothermic with respect to the adsorbed system (by 17 – 56 kJ mol (cid:1) 1 depending on the location), there are intriguing correlations between the adsorption/reaction energies and the geometries of the adsorbed species, of particular signi ﬁ cance being the coordination of methyl hydrogens. These observations emphasise the importance of adsorbate coordination with the framework in zeolite catalysed conversions, and how this may vary with framework topology and site location, particularly suited to investigation by QM/MM techniques.

Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s. S e e h t t p://o r c a . cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s. Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

I. Introduction
Computational methods are now very widely used in catalytic science and are increasingly powerful in obtaining an understanding of catalysis at the molecular level, where they yield models for both the structure and mechanism that assist and complement the interpretation of experimental data. The range and scope of the eld is illustrated by the recent review of their applications in zeolite catalytic science given in ref. 1.
The majority of contemporary computational studies in catalysis use electronic structure methods, especially Density Functional Theory (DFT) which is very widely applied, employing methods based on periodic boundary conditions (PBCs). Such methods have enjoyed considerable success and have indeed become increasingly predictive. An alternative approach is to use "QM/MM" or embedded cluster methods, which have been extensively used in computational chemistry and biomolecular sciences and in modelling localised states in solids. Such methods describe a nite region (which in the case of catalysis includes the active site, the reacting species and surrounding atoms) by a quantum mechanical method, while describing more distant regions of the solid by a more approximate method, oen based on interatomic potentials. They have a number of advantages which will be discussed in greater detail in the following section; but it can be argued that they are inherently more appropriate for modelling a localised state, such as an active site on a surface or in an enzyme, as they focus the computational effort on the region containing the active site rather than the whole molecule or solid. This paper will present recent applications of QM/MM techniques to the important and topical problem of catalytic methanol conversion in ZSM-5 and where we obtain new information on mechanistic aspects inspired by recent observations from neutron scattering studies.
In the next section we describe in greater detail the basis of the methodologies employed and illustrate briey their application by reference to earlier studies of metal centres in zeolites. The subsequent section presents our study of methanol/ ZSM-5 which demonstrates well the power of the techniques and yields new insight into the interaction of methanol with acidic sites in zeolites.

Methodology
In our work, we have employed a hybrid QM/MM model of zeolites developed and implemented by Paul Sherwood and colleagues 2-6 in the computational quantumchemical environment soware ChemShell 7-10 as illustrated in Fig. 1.
The model draws on the strengths of simple molecular cluster approaches that terminate a cluster carved out of the material with hydrogen, but also accounts for steric constraints and the elastic response of the framework. The cluster is treated with an appropriate level QM method, which is chosen as a compromise between the required accuracy and available compute resources. The ChemShell zeolite model thus has much in common with alternative subtractive mechanical embedding approaches by Sauer [11][12][13] (implemented in QM-POT) and Morokuma 14 (ONIOM). The ChemShell methodology, however, goes beyond such schemes and provides an accurate description of the electrostatic interactions between the QM region and the framework, including polarisation and embedding effects. More recently, similar approaches (although more focused on a careful termination of the QM region with specially calibrated link atoms) have been employed, for example, by Nasluzov et al. 15,16 (cov-EPE) and Sushko et al. 17,18 (GUESS). The native ChemShell zeolite model briey described in this section is complemented in ChemShell with support for both the QM-POT, cov-EPE and GUESS methods. The exible, modular approach adopted in ChemShell extends these methodologies to other classes of material, from metal oxides to enzymes. As argued above, these schemes serve as a computationally cheaper and at the same time chemically more accurate alternative to popular DFT PBC approaches.
To describe the long-range polarisation effects in these materials, the Chem-Shell model employs the Hill-Sauer molecular mechanical force-eld 19,20 with rigid ions described by modied constant point charges q Si ¼ 1.2e and q O ¼À0.6e, and atoms assumed to bind to each other via polar covalent bonds. These charges have been tted to reproduce the electrostatic eld in the pores of the zeolitic framework materials described using PBCs. When we dene a QM cluster embedded in a siliceous system, we have to cut through Si-O bonds that connect atoms, one in the QM and the other in the MM regions. As under-coordinated QM atoms on the QM region boundary are terminated with hydrogens, an articial strain on the system is introduced due to spurious interactions between the hydrogen and the nearest neighbours in the MM region. To remedy this artefact, a charge on the corresponding MM atoms is shied outwards, further onto the next nearest neighbours, whereas the dipole created by such a shi is compensated locally by a pair of point charges placed along the corresponding bonds. Further information on the method can be found in the original publications, and relevant model construction details are given below.
Next, we provide a short outline of our earlier work on metal centres in zeolites, which illustrate the power of our approach and its potential in the elucidation of the structure and properties of catalytically active sites in zeolites.

Applications to metal centre in zeolites
Incorporation of metal sites in zeolitic frameworks has been the subject of numerous investigations. In particular, Ti doped silicalite-1 (TS-1), 21 a siliceous counterpart of ZSM-5, has been of great interest, since it is used with aqueous H 2 O 2 as the primary oxidant in a large number of important reactions such as the epoxidation of linear olens, oxidation of linear alkenes to alcohols and ketones, oxidation of alcohols, hydroxylation of aromatics, oxidation of amines, and oxidation of sulphur compounds and ethers. 22 In our work, [23][24][25] we have proposed a number of structural models for a Ti species that would result from the framework hydrolysis and site inversion that can take place during the synthesis and post-synthetic treatment and/or, depending on the catalytic process, occur as a part of the catalytic transformation. The models of interest are shown in Fig. 2.
The geometrical parameters of the tripodal sites obtained with ChemShell calculations proved to be in excellent agreement with experimental EXAFS data, 26 and therefore could be usefully employed to rationalise the difference in the observed spectroscopic signature of metal centres in differently treated samples, as summarised in Table 1.
Our analysis of site transformation in the presence of water, illustrated in Fig. 3, 27 allowed us to draw the crucial conclusion that the catalytically active sites in TS-1 are hydrated tripodal rather than tetrapodal, which is strikingly similar to Ti doped mesoporous materials that had been previously explored using simple cluster models. 28 Further, we consider the site reactivity towards H 2 O 2 in anhydrous and hydrous conditions according to the reaction pathways drawn in Fig. 4 and 5. 29 Our calculations conrmed the close competition between the h 1 and h 2 complexes as candidates for the oxygen donating species in TS-1. Curiously, a much higher energy (by ca. 400-450 kJ mol À1 ) superoxide radical species has also been reported (from EPR studies) to emerge as a result of interaction of aqueous H 2 O 2 and non-aqueous urea H 2 O 2 . 30 Our calculations unambiguously identied the superoxide supported on tripodal Ti as the magnetically active      centre with the spectroscopic signature in close agreement with the experiment as shown in Table 2.
From this short outline, we conclude that the hybrid QM/MM methodology implemented in ChemShell provides us with a reliable tool to investigate both the structure and physical and chemical properties of active sites in zeolites. In the next section we employ ChemShell in a study of methoxylation of two zeolites with very different behaviour.

III. Methanol-active site interactions in zeolites
The zeolite catalysed conversion of methanol to hydrocarbons (MTH) has been researched signicantly since its initial discovery by Chang and Silvestri, 31 with the rst commercialization of the MTH process in New Zealand in 1985. 32 Demand for light olens has also lead to extensive research into the use of zeolitic catalysts for converting methanol to olens (MTO). 33 An early step of interest in the reaction mechanism is the formation of framework methoxy species aer initial physisorption of methanol through Hbonding to the zeolite Brønsted acid site. [34][35][36] However, studying such species experimentally is difficult, as secondary reactions dominate rapidly at higher temperatures. 37 Recent studies using neutron scattering 38 have observed the total conversion of methanol to framework methoxy at room temperature in commercial samples of H-ZSM-5 (in contrast to previous experiments discussed in ref. 38, which suggest that elevated temperatures are necessary), which is, however, not observed in commercial samples of H-Y with the same Si/Al ratio. This difference was attributed to the necessary dealumination procedure for high silica H-Y synthesis (creating defects which decrease the activity). However, detailed theoretical studies into the effect of framework topology on methanol behaviour in zeolite catalysts are necessary for improved understanding and further catalyst design. A short discussion of theoretical studies of methanol behaviour in zeolite catalysts will now follow.
Though the MTG/MTO process has received a signicant amount of attention from a theoretical perspective, 39,40 the aforementioned earlier steps are not understood in detail. Ab initio studies of methanol adsorption began with small aluminosilicalite clusters. 41,42 Numerous studies using three tetrahedral site clusters have concluded that the hydrogen bonding of the methanol oxygen with the Brønsted acid hydroxyl is exothermic, generally calculating adsorption energies of between 50 and 100 kJ mol À1 . 35,36,[43][44][45] It has also been found that there is an energetic difference between 'side-on' and 'end-on' adsorption, with the end-on geometry (where the methyl group protrudes almost perpendicularly from the aluminosilicate cluster) calculated to be more favourable, and methoxylation taking place via the rapid equilibration between the end-on and side-on geometry. 43,46 All studies into the methoxylation process using clusters conclude a barrier to methoxylation with values of $180-225 kJ mol À1 . 35,36,43,45,46 Methoxylation and adsorption has also been studied in periodic zeolite structures, with static calculations showing endothermic methoxylation in chabazite (with a barrier of 56 kJ mol À1 ), 47 and DFT-MD giving an adsorption energy of À94 kJ mol À1 . 48 Andzelm and Govind 49 studied the methoxylation process in a periodic ferrierite structure, nding a methoxylation barrier of 226 kJ mol À1 . However, it was found that when a second methanol molecule was added to the system, the barrier to methoxylation was then reduced by 42 kJ mol À1 due to the formation of hydrogen bonding networks spanning the width of the channel, which signicantly stabilise the methoxylation transition states. We note that the use of a discrete cluster model is less accurate than the use of a periodic zeolite structure, due to the lack of any long-range electrostatic potential inuencing the system energetics. However, as noted earlier, in addition to being computationally expensive in calculations involving large unit cells, the use of periodic boundary conditions can bring problems when calculating reaction barriers and adsorption energetics, such as incomplete relaxation of the system or the inuence of the interaction between sorbate mirror images. This problem can be especially signicant when comparing the effect of the framework structure on sorbate behaviour. A large zeolite unit cell would minimise the effe c to ft h e s es o r b a t e -sorbate interactions on the methoxylation process compared to, say, chabazite where they may become signicant. These complications can be bypassed using the quantum mechanical/molecular mechanical (QM/MM) embedded cluster technique 9,50 outlined in the previous section, where large framework structures may be modelled relatively cheaply, providing an accurate electrostatic potential that will affect the behaviour and reactivity at an active site.
We present a comparison of the deprotonation, methanol adsorption, and methoxylation energies in H-ZSM-5 and H-Y zeolite frameworks, in an attempt to investigate the effect of framework topology and active site location on reactive sorbate behaviour. The deprotonation energy is rst investigated as it is a major factor controlling the Brønsted acidity of the differing frameworks; our calculations include a comparison of three different sites in the ZSM-5 structure. The adsorption energy of methanol at these sites is then compared to quantify the differences in interaction strength between methanol and the differing acidic hydroxyl sites, aer which the methoxylation energy is then calculated to quantify the effect of these differing environments on the reaction.

Methodology
We rst discuss the construction of the embedded cluster models for the H-Y and H-ZSM-5. We begin with the purely siliceous spherical structure generated from the periodic experimental unit cell structure of silicalite 51 and siliceous faujasite. 52 The use of the experimental unit cell is assumed to have minimal effects on the geometric and electronic structure. Each spherical cluster centres on a Si T-site. Only one cluster is generated for the H-Y structure as the tetrahedral sites are equivalent, which, however, is not the case for H-ZSM-5. Therefore, three embedded clusters are generated based on three sites in the ZSM-5 structure: 53 a T-site at the straight channel (M7-T1), the sinusoidal channel (Z6-T4) and the channel intersections (I2-T12). The clusters, including their locations in the MFI structure are shown in Fig. 6. 54 T h eQ Mr e g i o ni se x p a n d e df r o mt h ec e n t r a la t o mt oi n c l u d et h eh neighbour, i.e. the third oxygen atom away from the T-site as shown in Fig. 6. The cluster is then terminated with hydrogen atoms. We also include two  concentric regions of molecular mechanics (MM) atoms spherically around this QM cluster (Fig. 7). The MM atoms in the rst concentric region are allowed to relax during the energy calculation; this region reaches a radius of 10.   The calculation of the energy and gradients of the QM region used the Gamess-UK package, 62 and the calculation of the energy and gradients of the MM region employed the DL_POLY package, 63 with the QM and MM calculations coupled in the ChemShell environment. 53 Electronic minimisation was deemed complete for the QM system when the energy variation was below a threshold of 2.72 Â 10 À6 eV was achieved (1 Â 10 À7 Ha). The geometries were optimised using the Limited-Memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) 64-67 algorithm with a convergence threshold of 0.015 eV A À1 .

Results and discussion
2.1. Calculation of deprotonation energies. In our rst comparison of H-Y and H-ZSM-5, we investigate the deprotonation energy differences between the structures. The deprotonation energy inuences acidity 68-70 as a lower deprotonation energy will favour protonation of any sorbate. We begin by determining the oxygen site around the central aluminium T-site, that has the lowest deprotonation energy. This site will then be used for further studies into the methanol adsorption, and the later methoxylation studies. We dene the deprotonation energy as: Here, E zeoH is the energy of the protonated zeolite and E zeo À is the energy of the deprotonated zeolite. Polarisation effects outside of our active region, where all atoms relax during geometry optimisation, are neglected in the following discussion as they represent a relatively constant shi in the deprotonation energies for all Brønsted sites considered (< 50 kJ mol À1 using Born and Jost's formulae 9,71 ). Fig. 8-11 show geometries of the four protonated systems around the central Al T-site in each zeolite system, labelled O1-O4. The deprotonation energies of each oxygen site were calculated using the PW91 exchange correlation functional and are listed in Table 3.
Using the PW91 functional, the most acidic oxygen sites are O1 for H-Y, O3 for H-ZSM-5 (I2), O3 for H-ZSM-5 (M7) and O2 for H-ZSM-5 (Z6). The difference between the lowest and highest deprotonation energies, which are respectively H-Y and H-ZSM-5 (M6), is 90.2 kJ mol À1 . H-Y has the lowest deprotonation energy of the series, followed by ZSM-5 (Z6) located in the sinusoidal channel, ZSM-5 (I2) in the intersections, and the T-site with the highest deprotonation energy is in the straight channel (M7).
We now study the deprotonation energies of the most acidic sites in each structure using the B3LYP and B97-2 hybrid exchange-correlation functionals, as listed in Table 4. The deprotonation energies calculated using the B3LYP functional were higher than when using the PW91 functional, and higher still using the B97-2 functional, which can be attributed to the superior representation of the kJ mol À1 calculated for ZSM-5 (I2) is 122 kJ mol À1 lower than obtained by ab initio molecular dynamics for the MFI structure, 70 using the same PW91 functional.
Possible reasons for this difference could be the aforementioned sorbate-sorbate interactions, the more complete representation of polarisation afforded using the QM/MM method, or the charge compensation parameters needed for periodic DFT calculations. We also note that ref. 70   Though the acidity of a zeolite is a typical indicator of the activity in a catalytic process, the magnitude of the sorbate interaction with the active site may also be a deciding factor in a methoxylation process. We note that although a favourable adsorption energy can facilitate the process, it may also hinder the methoxylation through increasing the barrier to the reaction. The adsorption energy of methanol is calculated as:

used a different intersectional T-site
where E MeOH is the calculated energy of a methanol molecule in vacuo, E ZeoH is the energy of the protonated embedded zeolite cluster, E MeOH+ZeoH is the energy of the methanol adsorbed onto the protonated embedded zeolite cluster, and E BSSE is the basis set superposition error calculated using the counterpoise correction method. 74 The geometry of the adsorbed methanol for each cluster, aer optimisation with the B3LYP functional, is shown in Fig. 12.
We have studied the 'side-on' adsorption in all zeolite systems. Though we note in previous studies that 'end-on' adsorption is preferred, 43,46 the side-on adsorption is more conducive to methoxylation, as is the emphasis of this study, and the barriers to this switching between the geometries are very small and may be overcome under experimental conditions. From Fig. 12 it is observed that the H-bond length between the hydroxyl site and the methanol hydroxyl is 1.476 A for the H-Y cluster and 1.452-1.478 A for the H-ZSM-5 clusters, with the shortest Hbond being in the M7-T1 cluster despite having the highest deprotonation energy. However, we note that the range of bond distances of 0.026 A is quite small and so a correlation with deprotonation energies is not necessarily meaningful.
We should note that the BSSE values were calculated with the B3LYP functional only. As an approximation, we have then taken the percentage error (5.5%) that these values contribute to each system and applied it to the adsorption energy calculations using the PW91 and B97-2 functionals. The calculated adsorption energies are listed in Table 5 aer subtraction of this percentage error.
For all functionals, the order of methanol adsorption energy is H-ZSM-5 (M7) < H-Y < H-ZSM-5 (I2) < H-ZSM-5 (Z6). We note that the difference in adsorption between the most favourable site, H-ZSM-5 (M7), and the least favourable site, H-ZSM-5 (Z6), is between 60 and 70 kJ mol À1 , depending on the functional used. The methanol in the H-ZSM-5 (M7) system is H-bonded to two framework oxygens by two C-H hydrogens (depicted in Fig. 12c), suggesting that this geometry is more favourable than a short, single C-H/O bond, as is demonstrated by the less favourable adsorption energy in the H-ZSM-5 (Z6) system. It may be signicant that the most favourable adsorption process also occurs at the site with the highest deprotonation energy, but we also see that the site with the lowest deprotonation energy (H-Y) has the second lowest adsorption energy. The calculated adsorption energy for methanol in ZSM-5 (M7) appears anomalous when considering the overall correlation between deprotonation and methanol adsorption energies for FAU and MFI, however it is perhaps more reasonable to consider that there should not be a correlation at all, and in fact the adsorption  energy is dominated by the H-bonding capability of the methyl group of the methanol. We note that previous calculations of Svelle et al. reported a weaker adsorption for methanol at the I2-T12 site than we present here: E ads was calculated as À86 and À114 kJmol À1 when using a periodic model with the PBE and PBE-D exchange-correlation functionals, respectively. 75 Preliminary comparative calculations indicate that the electrostatic embedding environment makes a stronger contribution to our adsorption energy, but the signicant contrast between our model and that of Svelle et al. necessitates further investigation before any certain conclusions can be drawn. 2.3. Calculation of methoxylation energies. We calculated the energy of methoxylation for each cluster. The energy is dened as: where E ZeoMe(aq) is the calculated energy of the methoxylated species with the neighbouring product water molecule. The energy change associated with this process should act as a marker for reactivity differences between the two frameworks and the differing sites within the MFI framework. The methoxylation geometries were based on those obtained by Andzelm and Govind 49 from the S N 2 mechanism followed from the original side-on adsorption geometry, where the resultant water molecule H-bonds to the formerly protonated oxygen, and the methyl group is attached to the adjacent oxygen across from the aluminium Tsite. The geometries were optimised and energies calculated using solely the B3LYP functional. The values of E methox are listed in Table 6, along with the difference between E ads and E methox . The optimised methoxylation geometries are shown in Fig. 13. Firstly, we note that all the methoxylation energies are exothermic overall but are endothermic with regard to the adsorbed methanol molecule, as observed in previous cluster simulations 35 and in the periodic chabazite structure. 42 We observe that the order of the methoxylation energy is H-Y < H-ZSM-5 (I2) < H-ZSM-5 (M7) < H-ZSM5-I2 (Z6), in agreement with deprotonation energies but in contrast to the order of methanol adsorption energies shown in the previous section. This suggests that the tendency of methoxylation is not related directly to the adsorption energy, at least when only one methanol molecule is present in the immediate environment.
Finally, we turn our attention to the geometries of the methoxylated system. The most exothermic methoxylation takes place in H-Y, which is also the more coordinated system: the water molecule is coordinated to three framework oxygen sites. Threefold coordination is also the case for the ZSM-5 (I2) system, which has the second most exothermic methoxylation energy, but for H-Y the water molecule is also coordinated to two of the methyl hydrogens. The H-bonds associated with the H-Y system are not the shortest of all the congurations: the shortest is 2.24 A, compared to 2.04 A in the ZSM-5 (Z6) system. However, the longer H-bonds in H-Y are compensated by their increased quantity, thus highlighting the importance of highly coordinated structures in stabilising the system when considering the interaction of sorbates such as methanol with the active sites of acidic zeolite catalysts. Using the QM/MM methodology we have, therefore, obtained new insights into the interaction, and reaction, of methanol with zeolite acidic sites and discussed how this may change with framework structure and aluminium location. Future work will use the embedded cluster models to investigate the barriers of methoxylation in zeolites using transition state searches. The effect of including multiple methanol molecules will also be investigated as the formation of Hbonded networks reduces the total methoxylation barrier for the FER framework, 49 thus further emphasising the importance of highly coordinated systems in controlling the barrier of such reactive processes.

IV. Summary and conclusions
QM/MM methods clearly have a major role to play in the study of catalytic systems. Earlier work demonstrated their efficacy in modelling active sites in zeolites. Here, we investigated the energetics of deprotonation, methanol adsorption and methoxylation in two acidic zeolite catalysts. Embedded cluster systems have been used to model zeolites H-Y and H-ZSM-5, with acidic sites centred in the straight channel, sinusoidal channel and intersection region of the latter system. The hierarchy of deprotonation energies was calculated as H-Y < H-ZSM-5 (Z6) < H-ZSM-5 (I2) < H-ZSM-5 (M7). The side-on adsorption energy of methanol in each structure was also calculated, giving a trend of H-ZSM-5 (M7) < H-Y < H-ZSM-5 (I2) < H-ZSM-5 (Z6). The most favourable adsorption site, which was H-ZSM-5 (M7), exhibited an adsorbed geometry where two methyl hydrogens were able to hydrogen bond to framework oxygen atoms, suggesting that the Hbonding behaviour of the methyl group may be dominant in controlling the adsorption energy. The methoxylation energy was calculated, giving a trend of H-Y < H-ZSM-5 (I2) < H-ZSM 5 (M7) < H-ZSM5-I2 (Z6). The ordering of our results suggests the total exothermicity of methoxylation is not dependent on the adsorption energy of methanol. However, the effect of increased framework coordination, including H-bonding of methyl groups, was emphasised as the most coordinated system gave the lowest energy conguration, similar to our observations for methanol adsorption. The present paper shows the role of the QM/MM modelling techniques in zeolite catalysis, as well as the limitations in using the deprotonation energy as a measure of catalytic reactivity for a material. We believe that QM/MM modelling in catalytic science can, however, be made much broader.