Adsorption in zeolites using mechanically embedded ONIOM clusters

We have explored mechanically embedded three-layer QM/QM/MM ONIOM models for computational studies of binding in Al-substituted zeolites. In all the models considered, the high-level-theory layer consists of the adsorbate molecule and of the framework atoms within the first two coordination spheres of the Al atom and is treated at the M06-2X/6-311G(2df,p) level. For simplicity, flexibility and routine applicability, the outer, low-level-theory layer is treated with the UFF. We have modelled the intermediate-level layer quantum mechanically and investigated the performance of HF theory and of three DFT functionals, B3LYP, M06-2X and oB97x-D, for different layer sizes and various basis sets, with and without BSSE corrections. We have studied the binding of sixteen probe molecules in H-MFI and compared the computed adsorption enthalpies with published experimental data. We have demonstrated that HF and B3LYP are inadequate for the description of the interactions between the probe molecules and the framework surrounding the metal site of the zeolite on account of their inability to capture dispersion forces. Both M06-2X and oB97x-D on average converge within ca. 10% of the experimental values. We have further demonstrated transferability of the approach by computing the binding enthalpies of n-alkanes (C1–C8) in H-MFI, H-BEA and H-FAU, with very satisfactory agreement with experiment. The computed entropies of adsorption of n-alkanes in H-MFI are also found to be in good agreement with experimental data. Finally, we compare with published adsorption energies calculated by periodic-DFT for n-C3 to n-C6 alkanes, water and methanol in H-ZSM-5 and find very good agreement.


I. Introduction
Zeolites are among the most widely studied inorganic materials.7][8][9] Accurate description of sitespecific binding is critical to the calculation of catalytic pathways and the development of new catalysts but remains a challenge for theoretical models because the host-guest interactions are determined by long-range electrostatic interactions and the long-range electron correlation effects that give rise to dispersion forces.
Periodic density-functional theory calculations can adequately address the electrostatic problem, but LDA and GGA functionals fail to capture the rather weak dispersion forces, and advanced, meta-hybrid functionals (with demonstrated ability to capture dispersion, e.g., M06-2X or oB97X-D) are not easy to utilize because of algorithmic difficulties in computing the exact exchange.Although empirical dispersion corrections (DFT-D) 10 have been parameterized for various GGA functionals, and more accurate semilocal exchange functionals (vdW-DF2) are being developed designed to capture dispersion interactions in a nonempirical fashion, 11 periodic-DFT calculations remain computationally expensive for zeolites (cubic scaling with the number of atoms), especially if one wishes to map out catalytic pathways.
While it is widely recognized that the extended aluminosilicate framework around the active site of a zeolite-based catalyst is as important for the description of site-specific binding as the active site itself, it has also been argued that quantum mechanical description of atoms far from the active site and the substrate might not be critical.3][14] These layering strategies are in essence subtractive computational schemes, whereby the low-level-theory energy of an individual layer is subtracted from the energy of the real system and is substituted by a high-level-theory estimate.Typically, the reactive region is treated quantum mechanically (QM) and the surrounding environment is treated with a molecular mechanics force field (MM).However, it is not unusual to treat the environment quantum mechanically as well (QM/QM), either with low-level theories or with periodic-DFT, with concomitant increase in the computational cost.Notable in this respect is the QM/QM strategy of Tuma and Sauer (MP2/CBS:DFT + dispersion), whereby a small-size cluster is used to model the reactive region at the MP2 level, while the rest of the system is treated with periodic-DFT. 15,16This strategy has the advantage of handling long-ranged electrostatic interactions more accurately and of ameliorating the problem of GGA functionals not capturing dispersion interactions.One can surely envisage the development of variants with MP2 replaced by meta-GGA or hybrid-meta-GGA functionals specifically developed and parameterized to take into account dispersion and hydrogen bonding.
Hybrid QM/MM methodologies are inherently approximate as they primarily aim at routine applicability in large systems (e.g., enzymes, zeolites).The size of the QM domain, the theory level at which the QM domain is treated, the size of the overall system, the complexity of the MM force field or of the low-level quantum theory (e.g., HF or semi-empirical) in the environmental domain, and the interaction between the layers (e.g., mechanical versus electrostatic embedding) are all the result of a compromise between accuracy and computational efficiency.Choices that work universally remain a challenge.For example, what works for non-polar probe or substrate molecules might not work as well for polar ones, or when the adsorbate is basic enough to deprotonate proton-exchanged zeolites, leading to formation of an ion-pair structure.
8][19][20][21][22][23][24][25][26][27][28][29][30][31] ONIOM-(B3LYP/6-31++G(d,p):HF/3-21G) single point calculations on ONIOM(B3LYP/6-31++G(d,p):MNDO) optimized structures of a (12T:48T) cluster model of acidic CHA for the adsorption of H 2 O and NH 3 correctly predict proton donation to NH 3 and ion-pair formation. 25A comparison of computed binding enthalpies with experimental data was not made, but one expects that the deviation from experiment should be significant given that B3LYP does not capture dispersive and hydrogen-bonding interactions.The binding energy was, nonetheless, in satisfactory agreement with periodic-B3LYP calculations -albeit somewhat underestimated -which is remarkable considering the rather small size of the embedded cluster and the level of theory at which the outer-layer was modelled (HF/3-21G).Nevertheless, the reasonable agreement between periodic-B3LYP calculations and ONIOM(B3LYP/6-31++G(d,p): HF/3-21G//B3LYP/6-31++G(d,p):MNDO) suggests that electrostatic interactions are rather well captured, although, convergence with the outer-layer basis set size was not reported and thus we are unable to assess the role of error cancellation in the ONIOM estimates.Interestingly, when the same methodology was applied to a host of probe molecules within a wide range of proton affinities (i.e., gas-phase basicities), and upon comparison with experimental data, the binding was severely underestimated, especially that of ion-pairs. 24oekfa et al. have used a 5T:34T embedded model of H-MFI and computed binding energies for ethylene, benzene, ethylbenzene, and pyridine considering a combination of MP2 for the high-layer and various methods for the outer layer, including UFF (Universal Force Field), HF, B3LYP and M06-2X, with various combinations of triple-z-plus-diffuse and double-z quality basis sets, respectively. 28They have asserted that ONIOM(MP2/6-311+G(2df,2p):M06-2X/6-31G(d,p)//MP2/6-31G(d,p):M06-2X/6-31G(d,p)) gives good binding energies.This assertion was based on a comparison of the computed energies with experimental energies in H-FAU and not in H-MFI (the zeolite that was actually modelled) and thus it is difficult to ascertain the accuracy of the method.In fact, the rather acceptable error of ca.0][31] By investigating convergence with respect to the size of the framework QM domain, they have suggested that the QM region does not have to extend beyond the first coordination sphere and proposed the use of a 5T cluster which they model at the oB97X-D/6-311++G(3df,3pd)// oB97X-D/6-31G(d,p) theory level (namely, optimization at the oB97X-D/6-31G(d,p) level followed by a single-point energy calculation at the oB97X-D/6-311++G(3df,3pd) level).Thus the burden of capturing the interaction between the substrate molecule and the framework of the zeolite is on the force field that models the MM region, typically consisting of up to 44 tetrahedral atoms.They have computed binding enthalpies within 10% of the reported experimental values for species that are both physisorbed and chemisorbed in MFI, H-MFI, and H-BEA zeolites. 31This is a promising methodology which has demonstrated considerable improvements in accuracy, although it remains to be tested for a wider class of polar probe molecules in order to ascertain its accuracy and whether further re-parameterization of the CHARMM force field will be required.
Having in mind wide and routine applicability, in this article we explore a mechanically-embedded, three-layer QM/QM/MM ONIOM approach which retains the simplicity of the Universal Force Field (UFF) for the MM region.Philosophically, our decision to explore a QM/QM/MM approach stems from the understanding that, for polar and strongly binding substrate molecules, the binding energy greatly depends on QM-MM interactions, dispersive, and electrostatic.Electrostatic interactions, in particular, polarize the active site (always part of the QM region) and a reliable approach toward capturing such interactions more adequately would be to replace part of the MM region with a quantum region that would be treated at a medium theory level which, nevertheless, is capable of capturing some of the dispersion interactions as well.By considering a broad class of probe molecules, our intent is to understand This journal is © the Owner Societies 2016 the range of these interactions in the sense of how extended the intermediate layer needs to be and of the ''minimal'' theory level one needs to employ to describe it, while maintaining a good balance between reliability and computational cost.

A. Embedded cluster design
Zeolite clusters were cut out from crystals taken from the database of structures maintained by the IZA. 36Dangling silicon bonds from this structure were saturated with hydrogen atoms at a bond length of 1.47 Å along the corresponding Si-O bonds of the crystal.One silicon atom was replaced by an aluminum atom and an H + cation was introduced to the adjacent oxygen atom in the lowest-energy configuration.
8][39][40][41] We have investigated both two-and three-layer ONIOM models, whereby the total enthalpy of the system is given by where H, M, and L refer to the high, medium and low levels of theory and SL, IL, and RL refer to the small, intermediate, and real layers of the system. 37For this study, the small and intermediate ONIOM3 layers have been treated with quantum mechanical methods, while the real layer has been treated with the universal force field (UFF).The character of stationary points has been confirmed by vibrational frequency analysis.
The small layer was chosen to include the first two tetrahedral coordination spheres around the Al atom, shown for H-MFI in Fig. 1a and b.The small layer was mechanically embedded into a real layer which includes all the atoms up to the sixth tetrahedral coordination sphere of the Al atom, again shown for H-MFI in Fig. 1c-e.For the three-layer ONIOM models, an intermediate layer was investigated to include the third, fourth, and in the case of H-FAU, fifth tetrahedral coordination spheres around the Al atom.The real layer of the model was investigated to include the fourth, fifth, sixth, and in the case of H-FAU, seventh tetrahedral coordination spheres around the Al atom.Representations of H-MFI, H-BEA, and H-FAU can be seen in Fig. 2a-c, respectively.In all models, the small layer of the model was allowed to relax while the intermediate and real layers were frozen to maintain the integrity of the zeolite framework structure.We should note that, in preliminary studies, relaxation of the intermediate layer occasionally converged to structures with multiple imaginary frequencies.

B. Theory levels
In both the two-layer and three-layer models we treated the QM region with the hybrid-meta-GGA functional M06-2X, a global functional with 54% HF exchange, with demonstrated ability to capture dispersion interactions. 42We have considered two basis sets: 6-31G(d,p) and 6-311G(2df,p).We have not considered augmented basis sets, not only because diffuse functions increase the computational cost, but also because they can increase the basis set superposition error (BSSE).Extensive studies by Truhlar and co-workers have shown that better accuracy could be achieved if the extra cost were instead invested in a larger valence space (e.g., triple-zeta basis), or more polarization. 43,44n the two-layer model, the outer (low-level) layer has been simulated with the UFF.The UFF has enjoyed significant popularity in mechanically embedded QM/MM electronic structure calculations of large systems because of its simplicity.Additionally, because of the absence of partial charges on the UFF atoms, it affords us the flexibility to partition the system into QM and MM regions in a number of ways without worrying about the net charge of the MM region not being zero.]26,28,45 We have investigated the effect of the cluster size by keeping the QM region size constant and varying the MM region.So, in the case of H-MFI, we have considered the ONIOM2 models M06-2X/6-31G(d,p):UFF and M06-2X/6-311G(2df,p):UFF and the clusters 17T:120T and 17T:209T.
Calculated binding strengths have been benchmarked against available experimental values provided in Table 1, measured using a combination of TPD and microcalorimetry between 323 and 480 K. [46][47][48][49][50] Vibrational frequencies were used to calculate thermal corrections to the binding enthalpies.For the spurious soft vibrational modes of the framework and the wrong asymptotic behavior of the enthalpy at low frequencies, we have employed the quasi-rigid rotor harmonic oscillator (qRRHO) approximation proposed by Grimme 51 and Head-Gordon 31 (see the ESI † for more details).Binding enthalpies were calculated by eqn (3) and are reported only at atmospheric pressure and 25 1C, as changing the temperature within the 50-100 1C range results in enthalpy changes of less than 0.5 kcal mol À1 in all cases.
A full study of cluster and layer size, and of theory level effects on the calculated DH ads has been conducted using an H-MFI model because the greatest amount of experimental data was available for all sixteen probe molecules used here (top 16 molecules in Table 1).3][34][35] In order to investigate transferability across different frameworks, and because full sets of experimental DH ads for the sixteen probe molecules were not available, we also investigated the binding of a number of alkanes (see Table 1) in H-MFI, H-BEA and H-FAU and compared with experimental data, when available.204.0The calculations involving quantum layers are subject to the basis set superposition error (BSSE).There is no exact way to correct for BSSE in ONIOM calculations.The most obvious way for both two-and three-layer ONIOM models is to isolate the QM layer (capped with H atoms to saturate the dangling bonds) and perform a counterpoise correction on it and the substrate molecule; we follow this approach here as well. 21,23,24,26The adsorption schemes are summarized in eqn ( 4) and (5).
In those cases where the most stable adsorbed state is the one with the acidic proton on the adsorbate, the counterpoise corrections were performed on the protonated adsorbate and the zeolite conjugate base. 52

A. Quantum cluster and two-layer ONIOM
In Table 2, we present binding energies in H-MFI, modeled with a 17T quantum cluster and two two-layer ONIOM clusters, 17T:120T and 17T:209T.For the 17T cluster, the calculations are performed at the M06-2X/6-31G(d,p) and M06-2X/6-311G(2df,p) theory levels; the ONIOM calculations are performed at the M06-2X/6-31G(d,p): UFF and M06-2X/6-311G(2df,p):UFF levels.By the mean signed error (MSE) of 6.7 kcal mol À1 and the mean unsigned error (MUE) of 8.1 kcal mol À1 , we can see that the 17T quantum cluster under-binds almost systematically when we do not correct for the BSSE.Both the MSE and MUE for BSSE-corrected energies increase to 11 kcal mol À1 , indicating that, for the 17T quantum cluster, M06-2X/6-31G(d,p) under-binds systematically.This poor performance does not seem to be basis set related, as increasing the basis set size from 6-31G(d,p) to 6-311G(2df,p) increases the average BSSEcorrected error by 0.5 kcal mol À1 .The shortcomings of small or medium-size cluster calculations of adsorption energies are of course well known and here are presented mainly for comparison purposes.Small and medium-sized quantum clusters have been quite reliable for the calculation of activation energies of reactions in zeolites primarily because the long-range interactions essentially cancel out, leaving an energy controlled by local electronic interactions, provided that significant changes in geometry do not take place. 53However, long-range interactions and confinement phenomena remain an issue when we calculate adsorption energies, because similar error cancellation does not occur.
The two-layer models M06-2X/6-31G(d,p):UFF and M06-2X/ 6-311G(2df,p):UFF show a marked improvement (Table 2).The incorporation of the third, fourth, and fifth tetrahedral UFF coordination spheres around the Al atom improves the BSSE-corrected MUE to 6.1 and 6.0 kcal mol À1 for the 6-31G(d,p) and 6-311G(2df,p) small layers basis sets, respectively, but we still see systematic under-binding, despite the fact that the UFF is known to overestimate the van der Waals interactions.It is notable that, when we do not include BSSE corrections, the MUE is lower, 4.5 kcal mol À1 for both the 6-31G(d,p) and 6-311G(2df,p) small layers basis sets (Table 2).That notwithstanding, there is clear improvement over the 17T cluster which is most evident in the adsorbates with the larger proton affinities, which bind more strongly in the two-layer ONIOM model -even though the binding energies are still underestimated.Adsorbates with larger proton affinities are  able to accept the proton from the zeolite and the resulting ion pair seems to be stabilized through interactions with the pore.
Increasing the size of the real layer to also include the sixth tetrahedral coordination sphere reduces the BSSE-corrected MSE and MUE by ca.0.5 kcal mol À1 both in the case of M06-2X/ 6-31G(d,p):UFF and in the case of M06-2X/6-311G(2df,p): UFF (Table 2).Overall, we see that by extending the size of the real layer (i.e., by including more of the framework) we have additional stabilization of the adsorbed stated on account of attractive van der Waals interactions between the substrate molecule and the walls of the zeolite.However, both ONIOM2 models investigated here tend to under-bind on average, a behavior that becomes more pronounced and more systematic when we correct for the BSSE.

B. Three-layer ONIOM -method, basis set and layer size effects
In the following we present results for a number of three-layer QM/QM/MM ONIOM models using a combination of theories and basis sets with a mechanically embedded 17T:23T:97T cluster model of H-MFI.The 17T:23T:97T cluster is of the same size as the 17T:120T cluster used in the ONIOM2 calculations but with 23T atoms of the real layer reassigned to an intermediate quantum layer.In all the models considered, the small layer (17T) is treated at the M06-2X/6-311G(2df,p) level.For the intermediate layer (23T), we assess HF theory and the three functionals B3LYP, M06-2X and oB97x-D, and the basis sets 3-21G, 6-31G, 6-311G, 6-31G(d,p) and 6-31G(2df,p).The outer layer is modelled with the UFF.
At the ILT/3-21G, ILT/6-31G and ILT/6-311G (ILT = HF, B3LYP, M06-2X, oB97x-D) levels, we see a dramatic change compared to the ONIOM2 calculations as all four models now overbind (negative MSE) with MUE in the range 5-10.5 kcal mol À1 , without BSSE corrections (Tables 3 and 4).M06-2X and oB97x-D over-bind systematically and, in fact, we see an increase in the MUE and in the degree of over-binding as we move from the minimal basis set 3-21G (MUE ca.8.0 kcal mol À1 ) to 6-31G (MUE ca. 10 kcal mol À1 ).For the 6-311G basis set, however, the MUE drops to the 3-21G level -which is still significantly high and almost twice as large as in the ONIOM2 calculations.For the same basis sets, HF and B3LYP appear better balanced, with MUE in the range 5-6 kcal mol À1 .Correcting for the BSSE is more beneficial to M06-2X and oB97x-D than to HF and B3LYP.For example, the MUE at the M06-2X/6-311G and oB97x-D/6-311G intermediate levels drops by 2 kcal mol À1 (to 5.5 and 6.2 kcal mol À1 , respectively) when we correct for the BSSE.At the HF/6-311G and B3LYP/ 6-311G intermediate levels, the MUE drops by 0.4 kcal mol À1 (to 4.9 and 5.1 kcal mol À1 , respectively) when we correct for the BSSE.Interestingly, when one considers the BSSE-corrected MUE, all four models appear quite equivalent.
The rather large fluctuations in the MUE with the intermediate layer basis sets 3-21G, 6-31G and 6-311G have led us to investigate the effect of adding polarization (Tables 5 and 6).Specifically, we have performed single-point energy calculations with the intermediate layer basis sets 6-31G(d,p) and 6-31G(2df,p) on geometries optimized at the M06-2X/ 6-311G(2df,p):ILT/3-21G:UFF and M06-2X/6-311G(2df,p):ILT/ 6-31G:UFF (ILT = HF, B3LYP, M06-2X, oB97x-D).For M06-2X and oB97x-D, the MUE reaches a plateau value of ca. 4 kcal mol À1 , irrespective of the optimization geometry on which the singlepoint calculations were performed (Tables 5 and 6; see also Fig. 3a  and b).Further, correcting for the BSSE leaves the MUE practically unaffected.In contrast, the MUE for HF and B3LYP in the intermediate layer are largely unconverged with respect to the basis set (Fig. 3).Adding polarization seems less beneficial, as the MUE rises to ca. 9 kcal mol À1 accompanied by significant under-binding, with MSE of ca.6 kcal mol À1 without BSSE correction and ca. 9 kcal mol À1 with BSSE correction.Clearly, due their intrinsic shortcomings and despite the attendant computational savings, neither HF nor B3LYP are suitable for three-layer QM/QM/MM ONIOM calculations of binding energies; the dispersive interactions between the substrate molecule and the framework surrounding the active site play a significant role.Of all the probe molecules investigated for this study, pyridine and its chloro-fluoro-and methyl-derivatives and the amines posed a great challenge to the two-layer ONIOM methodology.The MSE calculated over the subset of the pyridines and amines is 8.1 and 5.0 kcal mol À1 , respectively.In contrast, the three-layer models perform significantly better.For example, at the M06-2X/6-311G(2df,p):M06-2X/6-31(d,p):UFF// M06-2X/6-311G(2df,p):M06-2X/6-31G:UFF level, the MSE dropped to 4.3 and À2.1 kcal mol À1 , showing the importance of the intermediate quantum mechanical layer for the binding of basic molecules.The decrease in the error relative to the ONIOM2 models -of the same size -and the almost similar behavior of the intermediate layer theories in the ONIOM3 models suggest the importance of electrostatic interactions and of charge polarization in particular.QM/MM methods with mechanically embedded QM regions do not consider the polarization of the quantum region (i.e. of the active site) by the framework and vice versa, force fields with fixed partial charges on the atoms do not react to changes in the electron density of the quantum region.QM/MM formulations with electrostatic embedding address the former problem because the Hamiltonian includes an extra term with the electrostatic potential, but they still suffer from the latter, which is of some significance when we calculate transition states, as there is no reason to assume that the polarization of the MM region will not be affected along the reaction pathway.However, by introducing a quantum mechanical intermediate layer, which replaces part of the MM layer, we partly ameliorate the lack of polarization; at the same time, we are also treating the dispersion interactions more accurately than with the UFF. 54e have also investigated the effects of the size of the intermediate and real layers on the binding energy (Table 7).We have considered three H-MFI clusters, 17T:23T:97T, 17T:61T:59T and 17T:23T:186, modelled at the M06-2X/ 6-311G(2df,p):M06-2X/6-31(d,p):UFF//M06-2X/6-311G(2df,p):M06-2X/3-21G:UFF level.In the 17T:61T:59T cluster, we have reassigned 38 atoms of the real layer of 17T:23T:97T to the intermediate layer, while keeping the total number of atoms fixed.In the 17T:23T:186T cluster, we have added atoms within the sixth coordination sphere of 17T:23T:97T.The results are presented in Table 7 and in Fig. 4. The MUE, with or without BSSE corrections, remains practically flat and independent of the size of the cluster, indicating that, overall, it is not affected in any significant way when we extend the MM layer beyond the sixth coordination sphere around the Al atoms of the active site, fluctuating around the value of ca.3.8 kcal mol À1 .We have seen similar behavior in the two-layer ONIOM models presented earlier.Extending the intermediate layer can lead to minor improvement.

C. Adsorption of alkanes in H-MFI, H-BEA and H-FAU
To test transferability, we have extended our calculations to include frameworks of different pore sizes, specifically H-BEA  and H-FAU.Experimentally determined adsorption strengths for the series of 16 adsorbates investigated above were not available for H-BEA and H-FAU.Numerous studies have, however, investigated the adsorption of linear alkanes in H-MFI, H-BEA and H-FAU.3][34][35] In Table 8, we present calculations with and without polarization in the intermediate layer using a 17T:23T:97T cluster for H-MFI, a 16T:18T:77T cluster of H-BEA and a 14T:16T:188T cluster for H-FAU.Two trends are observed, both in the experimental and computational results.First, as the alkane chain length increases, the adsorption strength increases.When comparing the adsorption strengths with increasing alkane chain length for a given zeolite system, the slope of the computational and experimental results are similar, suggesting that the model does capture this effect.Adding polarization to the intermediate layers in this model was very important, as can be seen in the case of H-MFI.Without polarization H-MFI has a MUE of 9.6 kcal mol À1 , which decreases to 5.8 and 5.3 kcal mol À1 for the 6-31G(d,p) and 6-311G(2df,p) basis sets.In the case of H-BEA, the MUE is equal to 3.1 kcal mol À1 , while in H-FAU the MUE is 1.7 kcal mol À1 .In the case of H-MFI, the error is somewhat larger than in the other two zeolites because of the smallest pore and the higher proximity of the adsorbate to the wall of the pore, which is treated with the UFF.As the pore size increases, the adsorbate does not interact as much with the UFF layer, and the error decreases.Nonetheless, the accuracy is quite satisfactory in all three cases, especially if one considers the uncertainty in the experimental enthalpies of adsorption.][34][35] It is well established that binding enthalpy increases (in absolute value) as the pore size decreases and becomes more confining.For these three frameworks, the pore size varies in the order H-FAU 4 H-BEA 4 H-MFI, 36 and thereby, the relative binding would be expected to follow the order H-FAU o H-BEA o H-MFI.

D. Relationship between proton affinity and enthalpy of adsorption
By making use of thermodynamic cycles, differences in the binding strengths of probe molecules in Brønsted acidic zeolites have been correlated with differences in the proton affinities (PA) of the probe molecules themselves.The thermodynamic cycle involves a state in which the acidic proton has been donated to the substrate molecule (see Fig. 5 inset).The result of this assumption is a linear correlation of unit slope between DH ads and the PA.For any given zeolite, this is not a strong correlation, because not all substrate molecules are strong enough bases to accept the acidic proton.][57] Adsorbates with gas phase PA greater than B200 kcal mol À1 demonstrate the ability to abstract the proton from the zeolite and form an ion pair structure.6][57] In this range of PA, the quantum intermediate layer has contributed to the electrostatic stabilization of the ion pair.At low PA of less than B180 kcal mol À1 , the adsorbate is unable to abstract the proton of the active site and binding is primarily affected through dispersion and hydrogen bonding interactions in the pore.6][57] In these cases, the use of a theory in the intermediate layer able to account for dispersion forces between the adsorbate and zeolite walls (M06-2X and oB97x-D), as well as the benefit of adding polarization to the basis set to allow directional interaction between the adsorbate and pore wall can be seen in the accuracy of these models.

IV. Discussion
Although detailed quantitative comparisons with experiment can at times prove quite challenging, as the experimental values are a statistical average over a distribution of acid sites, adsorbate-site geometries and adsorbate configurations whereas the calculated values refer to optimized binding geometries at a particular acid site and adsorbate configuration, comparisons between different computational strategies can, on the other hand, prove quite instructive.The handling of long-ranged electrostatic interactions by embedded cluster QM/MM strategies has been a source of concern to critics of this methodology.Below, we compare our results with published periodic-DFT calculations and find very good agreement, which should allay some of the reservations about hybrid QM/MM calculations.
Chiu et al. have also reported binding energies of water and methanol in H-ZSM-5 using periodic-DFT.For water, the PBE-D value of À22.2 kcal mol À1 fares very well upon comparison with the experimental value of À21.5 kcal mol À1 and is close to our value of À21.3 kcal mol À1 .On the other hand, the vdW-DF2 value of À17.2 kcal mol À1 somewhat underestimates water adsorption.For methanol, the PBE-D value of À24.4 is close to our value of À23.6 kcal mol À1 , whereas the vdW-DF2 functional somewhat underestimates binding, predicting an energy of À21.3 kcal mol À1 .
Even though more systematic benchmarking is required to fully access strengths and weaknesses in the performance of QM/QM/MM embedded cluster ONIOM models for predicting adsorption in zeolites, they do not seem to under-perform relative to periodic-DFT calculations.
In Table S1 of the ESI, † we show the calculated adsorption entropies of n-C1 to n-C8 alkanes in H-MFI, assuming a molecular surface area of 200 Â 600 (pm Â pm) 60 for 2-dimensional free translations of the guest molecule.In the temperature range of 25 to 200 1C, the computed values vary little with temperature.Upon comparison with experimental values reported by De Moor et al. for n-C3 to n-C6 alkanes, 60 the computed values systematically overestimate the entropic losses, but they are, nevertheless, in good agreement with experiment, with an error of no more than 10 J mol À1 K À1 .Despite the good agreement with experiment, we are of the opinion that the entropic error that is associated with the notoriously inaccurate low-frequencies in the harmonic approximation cannot be fully addressed by the qRRHO approximation.Sauer and co-workers, [61][62][63] in recent work that was inspired by earlier work by Njegic and Gordon, 64 have convincingly argued and demonstrated that correcting for the error in the adsorption entropy due to the soft modes requires not only consideration of the strong anharmonic character of these modes, but also that the diagonalization of the Hessian be performed in the curvilinear space of the intrinsic coordinates of the system.This approach is not entirely without problems, as it is not always easy to find appropriate intrinsic coordinates, nor is it entirely computationally inexpensive for large systems (factor of 10 computational overhead) 63 as the anaharmonicity calculation requires a number of single-point In the catalytic cycle, DPE Z is the deprotonation energy of the zeolite, PA A is the gas phase proton affinity of the adsorbate, and DH Z-AH is the interaction energy between the deprotonated zeolite and the protonated adsorbate.Dotted lines show second order fits of the experimental and computational data to guide the eye when compared to the blue line, showing the predicted binding strength using the gas phase proton affinity (PA) correlation with a slope of 1. 48,49,[55][56][57] This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 26094--26106 | 26105 energy calculations to sample the potential energy surface along a normal mode in curvilinear space.

V. Conclusions
We have explored mechanically embedded three-layer QM/QM/ MM ONIOM models for computational studies of binding in metal-substituted zeolites.In all the models considered, the high-level-theory layer consists of the adsorbate molecule and of the framework atoms within the first two coordination spheres of the metal atom and is treated at the M06-2X/6-311G(2df,p) level.For simplicity and flexibility in partitioning the system into QM and MM regions in a number of ways without worrying about the net charge of the MM region not being zero, the outer, low-level-theory layer is treated with the UFF.We have modelled the intermediate-level layer quantum mechanically and investigated the performance of HF theory and of three DFT functionals, B3LYP, M06-2X and oB97x-D, for different layer sizes and various basis sets, with and without BSSE corrections.We have studied the binding of sixteen probe molecules with a broad range of basicities in H-MFI and compared the computed adsorption enthalpies with published experimental data.We have demonstrated that HF and B3LYP are inadequate for the description of the interactions between the probe molecules and the framework surrounding the metal site of the zeolite, on account of being unable to capture dispersion forces.Both M06-2X and oB97x-D on average converge within ca.10% of the experimental values, at the M06-2X/ 6-311G(2df,p):M06-2X/6-31G(d,p):UFF//M06-2X/6-311G(2df,p): M06-2X/3-21:UFF and oB97x-D/6-311G(2df,p):oB97x-D/6-31G(d,p): UFF//oB97x-D/6-311G(2df,p):oB97x-D/3-21:UFF theory levels, respectively.In particular, the three-layer ONIOM models perform significantly better than the two-layer model M06-2X/ 6-311G(2df,p):UFF in stabilizing ion-pair structures of adsorbate molecules with high proton affinities.The mean errors over the subsets of the pyridines and of the amines drop significantly showing the importance of the intermediate quantum mechanical layer for the stabilization of ion-pairs.We have further computed the binding enthalpies of n-alkanes (C1-C8) in H-MFI, H-BEA and H-FAU at the M06-2X/6-311G(2df,p):M06-2X/6-31G(d,p):UFF//M06-2X/6-311G(2df,p):M06-2X/3-21:UFF theory level and found very good agreement with experiment, with mean unsigned errors of 5.3, 3.1 and 1.7 kcal mol À1 , respectively, demonstrating transferability of the model across zeolite frameworks and classes of molecules.The transferability is further demonstrated by the good agreement between computed entropies of adsorption and experimental values reported by De Moor et al. for the adsorption of n-C3 to n-C6 alkanes in H-ZSM-5.

Fig. 2
Fig. 2 ONIOM representations of (a) H-MFI, (b) H-BEA, and (c) H-FAU zeolites showing the layers of the standard models.The small layer is shown as a ball-and-stick representation, the intermediate layer is shown as a tubeframe representation, and the real layer is shown as a wireframe representation.
This journal is © the Owner Societies 2016

Fig. 4 a
Fig. 4 Behavior of mean unsigned error with intermediate and real layer sizes and intermediate layer theory.

Fig. 5
Fig. 5 Correlation between gas phase proton affinity and DH ads in H-MFI zeolites.Calculations performed using the 17T:61T:59T model of H-MFI at the M06-2X/6-311G(2df,p):M06-2X/6-31G(d,p):UFF//M06-2X/6-311G(2df,p): M06-2X/3-21G:UFF theory level with BSSE correction.In the catalytic cycle, DPE Z is the deprotonation energy of the zeolite, PA A is the gas phase proton affinity of the adsorbate, and DH Z-AH is the interaction energy between the deprotonated zeolite and the protonated adsorbate.Dotted lines show second order fits of the experimental and computational data to guide the eye when compared to the blue line, showing the predicted binding strength using the gas phase proton affinity (PA) correlation with a slope of 1.48,49,[55][56][57]

Table 5
Effect of intermediate layer basis set polarization on DH ads in H-MFI.The BSSE-uncorrected energies are shown in parentheses.Binding energies of probe molecules in a mechanically

Table 6
Effect of intermediate layer basis set polarization on DH ads in H-MFI.The BSSE-uncorrected energies are shown in parentheses.Binding energies of probe molecules in a mechanically