Faraday Discuss., 1997, 106, 79»92 Computer simulation of zeolite structure and reactivity using embedded cluster methods Paul Sherwood,a* Alex H. de Vries,a Simon J. Collins,a Stephen P. Greatbanks,b Neil A. Burton,b Mark A. Vincentb and Ian H. Hillierb* a CCL RC Daresbury L aboratory, Daresbury, W arrington, UK WA4 4AD b Department of Chemistry, University of Manchester, Manchester, UK M13 9PL The use of bare cluster models to understand the nature of zeolite»substrate interactions may be improved to take account of the environment of the acid site.We consider two models for introducing the electrostatic Br‘nsted eÜects of the zeolite lattice. The –rst involves generating a specialised correction potential by –tting a non-periodic array of ca. 60 point charges to the diÜerence between the bare cluster and periodic potentials.The second starts by –tting a periodic array of atomic charges to the potential of the in–nite lattice and then builds up a classical cluster of ca. 2000 atoms into which the QM cluster is embedded. Such embedded cluster calculations, employing a T3 cluster, with electron correlation at the density functional theory level, are described, to model the interaction of water at a acid site.Br‘nsted Structures of the water»zeolite complex, and associated vibrational frequencies and 1H NMR shifts are calculated and compared with calculations of bare clusters of varying size and with experimental data. We then describe a mixed quantum mechanical»molecular mechanical (QM»MM) model derived by combining charges from the second model with a standard aluminosilicate force –eld.We report preliminary results on the eÜect of embedding on the energetics of a prototypical hydrocarbon cracking reaction ; the methyl-shift reaction of a propenium ion coordinated to the acid site. The importance of zeolites in many industrial catalytic processes has prompted continuing studies to understand the relationship between their structure and reactivity.A particular area of study has been the interaction of substrates with a acid site, Br‘nsted central to many reactions catalysed by zeolitic aluminosilicates, such as the cracking of hydrocarbons and the conversion of methanol to gasoline.1 Recently, there have been many examples of the functionalisation of zeolites by the exchange of protons by metal ions, leading to a range of catalytic applications such as the decomposition of nitrogen oxides over Cu-exchange zeolites.2 Interactions occurring at acid sites have been explored by a range of Br‘nsted experimental techniques, particularly NMR and IR spectroscopy, in an attempt to clarify the nature of the adsorbed substrate, with particular emphasis on the degree of proton transfer from the zeolite to the substrate.In parallel with these experimental studies, there have been continuing theoretical and computational developments in order to achieve the necessary realism to model the often quite subtle features of the zeolite»substrate interaction. The accurate modelling of these interactions is particularly challenging as the catalytic activity and selectivity arise from a number of diÜerent factors, including local 7980 Simulation of zeolites using embedded cluster methods structural eÜects on acidity, the level of aluminium substitution, the zeolite pore structure, chemical impurities and extra-framework ions.The most commonly used computational strategy is to employ clusters of varying sizes to represent the active site of the zeolite,3h5 and to use standard electronic structure methods, available in a number of commercial packages, to model the adsorbate»zeolite system.These models typically include two or three tetrahedral (T) atoms (Si, Al) and the bonds cleaved to generate the cluster are terminated by adding additional monovalent atoms, typically hydrogen.A drawback of these bare cluster models is that they omit long-range electrostatic eÜects and fail to include the geometric and electronic environment of the site speci–c Br‘nsted to a particular zeolite. These problems can, in principle, be addressed by the use of substantially larger clusters and studies using clusters up to T46 have been carried out.6,7 However, the routine use of ab initio quantum chemical methods is computationally prohibitive for systems of this size.A more realistic way of modelling the catalytic reaction, which avoids these problems, is to employ methods that can treat periodic systems, and to use these to model periodic supercells involving the active site and substrate. Such a strategy will often necessitate the use of quite large unit cells and hence may also be computationally prohibitive.However, such a strategy has been employed using periodic calculations based upon both a gaussian basis8 and a plane wave expansion within a density functional theory (DFT) formalism.9 The latter approach has been particularly useful in view of the ease of calculation of energy gradients allowing, in principle, both stationary structures to be obtained and molecular dynamics (MD) simulations to be performed.Embedded cluster models In other areas of condensed-phase modelling, particularly solvation studies and enzyme catalysis, hybrid methods are widely employed.10h12 Here, the reactive centre is described using a full QM treatment, and the more distant environment, which takes no direct part in the chemical reaction, is described at a much lower level of theory, often with a quite simple force –eld.The success of such hybrid methods naturally depends upon the use of an appropriate level of theory to describe the reactive centre, and the accurate description of both the electrostatic and other non-bonded interactions between the two regions, as well as the proper treatment of the junction between the two regions.In view of the polar nature of zeolites and other micro- and meso-porous materials, it is particularly important to model correctly the electrostatic environment of the active site resulting from the surrounding three-dimensional structure. We have suggested that a —exible and computationally economic way of achieving this is to use an embedded cluster model to study zeolite»substrate interactions.13 Here the cluster is embedded in an array of point charges which are chosen to model the electrostatic environment of the periodic lattice.However, since the formal atomic charges used in aluminosilicate force –elds are chosen to reproduce structural rather than electrostatic data, such charges may be inappropriate for use in embedded cluster models of zeolites.Indeed, oxygen charges vary from [0.26 for the CFF force –eld of Hill and Sauer,14 and [1.1 in the rigid ion model of van Beest et al.,15 to the formal [2.0 charge of the shell model of Catlow et al.16 One approach that we have used both to model the environment of the cluster and to correct for termination eÜects, is to derive a set of formal point charges by a leastsquares –t involving the electrostatic potential (ESP) of the in–nite 3D lattice, calculated at the Hartree»Fock (HF) level.13,17 Here the diÜerence between the ESP of the cluster and of the in–nite lattice in the region of the acid site is –tted to a point charge Br‘nsted potential to obtain point charges (potential-derived charges, PDC).This approach, which we henceforth label as model 1, has been validated by comparison with adsorp-P.Sherwood et al. 81 tion energies calculated at the full, periodic supercell level. Furthermore, when compared to full periodic results, indications are that our simpler description is competitive with more sophisticated embedding techniques.18 Although these charges which are derived from such a –tting technique, do re—ect the periodic nature of the lattice they need not be atom centred, and the individual values have no physical signi–cance.Indeed, even if they are atom centred, atoms of the same type will have diÜerent formal charges associated with them, depending upon the actual cluster chosen in the –tting procedure. Thus, they are unsuitable for inclusion in a general hybrid QM»MM method, where a complete force –eld is used to model the non-QM region.Furthermore, in the spirit of using general force –elds, it may be advantageous if such atomic charges were transferable between zeolites. We have, therefore, explored the use of an alternative strategy for obtaining PDC.19 Here formal atomic charges are derived for a periodic siliceous structure which reproduce the ESP for the same periodic structure calculated at the HF level.We have derived PDC for eight periodic siliceous zeolite systems by reference to ESPs from periodic HF calculations. Using a separate charge for each crystallographically unique atom, it was possible to –t the potential outside the van der Waals surface to within an RMS deviation of ca. 4 kJ mol~1. Fitted charges spanning [0.8 to [1.0 for oxygen, and 1.6 to 2.0 for silicon, were obtained when –tting to STO 3G calculations. Subject to the constraint that a common charge be assigned to all atoms of a given type, the RMS error in the point charge potential may rise as high as ca. 10 kJ mol~1 and a signi–cant reduction in the magnitude of the charges, relative to the unconstrained case, is observed.For –tting to an STO 3G ESP, a single assignment of an oxygen charge of [0.78 is suggested, whereas for a larger basis (6-21G*) a reduction in the compromise oxygen charge to ca. [0.6 is found. Once these periodic atomic charges, capable of reproducing the quantum periodic potential have been obtained, they may be used to construct a –nite set of point charges which, after treatment for termination eÜects, form the basis for embedded QM cluster models.We will describe a simple scheme of this type, referred to below as model 2. The advantage of both such embedded models is that the electrostatic eÜect of the PDC may be directly included in the one-electron Hamiltonian of conventional electronic structure codes, allowing for the location and characterisation of stationary structures, and for the calculation of molecular properties.However, there are some additional complications that need to be addressed. In view of the lack of van der Waals interactions with the non-QM region, geometric constraints must be imposed on the QM cluster to prevent the generation of unrealistic stationary structures.This has implications for the subsequent optimisations and calculation of vibrational frequencies, to be discussed later. Here, we describe embedded cluster calculations, based upon both approaches, to obtain PDCs. We discuss their use to model the interaction of water with a acid site, focusing on the prediction of substrate binding energies, vibrational Br‘nsted frequencies, and 1H chemical shifts. Hybrid QMñMM models Since the atomic charges from model 2 have been chosen to generate accurate electrostatic interactions, it is reasonable to assume that they may serve as a suitable atomic charge model for use in a classical force –eld for zeolite modelling, allowing the development of QM»MM models that can treat both the electrostatic interactions and the mechanical —exibility of the zeolite lattice.To date, no zeolite force –eld based on such a PDC charge model has been derived. We discuss mixed QM»MM results using the following three models for the QM»MM coupling, similar to those of Bakowies and Thiel.2082 Simulation of zeolites using embedded cluster methods (1) Mechanical embedding, where the QM»MM interactions are the classical bonded and non-bonded forces. (2) Electrostatic embedding, where the electrostatic potential due to the MM region is included in the electronic Hamiltonian of the QM region, allowing QM polarization.(3) MM polarization, where the electronic response of the MM region is modelled by coupled atomic polarizabilities in part of the MM region. We have explored these schemes based on combining –tted atomic charges with a standard aluminosilicate force –eld (CFF),21 and used them in a preliminary study reported here, of the so-called methyl»shift reaction in a propenium ion over an acid site in zeolite Y.The non-classical, cyclic propenium transition state structure in this reaction has been proposed by Sie22 as an alternative to earlier suggested pathways in hydrocarbon conversion reactions over zeolite sites.Computational details In the models used here, H-terminated T3 clusters based upon zeolite Y were chosen in view of the number of experimental studies involving the interaction of small substrate molecules with this zeolite. All PDC models were obtained using a periodic HF calculation, at the 3-21G level, of zeolite Y in which the aluminium sites were replaced by silicon, allowing the sodium counterions to be neglected.PDC model 1 A hydrogen-terminated T3 cluster, was excised from the zeolite Y structure, Si3O4H8 , the diÜerence in the ESP between the periodic structure and the cluster was –tted to obtain the PDC, as described in ref. 13. PDC model 2 The PDC used in model 2 were again obtained using an HF 3-21G periodic calculation on a purely siliceous zeolite Y structure.The oxygen and silicon charges which best reproduce the quantum mechanical ESP were found to be [0.92 and 1.84, respectively. However, these charges, which are appropriate for an in–nite periodic system must be modi–ed for use in our cluster model. This is accomplished by introducing, and treating, two boundary regions, the inner QM»MM boundary, and the outer boundary of the MM cluster.Initial modi–cations to the atomic charges are obtained by considering the charge on a given centre to arise from a sum of terms, one for each bond formed. The atomic charges are then assigned automatically, in a bond increment manner, based on connectivity, such that a silicon with only two bonds would have a charge of 1.84/2, i.e. 0.92. To correct for the loss of the long-range Madelung potential, we consider a wholly classical cluster, compare the ESP in the interior of the cluster with that from the in–nite array of periodic charges, and correct the former by addition of further charges outside the cluster. The values of these additional charges are derived by least-squares –tting to the diÜerence potential, by direct analogy with the PDC scheme of model 1.Fig. 1 illustrates the steps involved, comparing the periodic HF results [Fig. 1(a)] with the potential from the periodic point charges [Fig. 1(b)]. When the –nite cluster is terminated as described above the resulting potential in the interior of the cluster [Fig. 1(c)] is essentially identical to that of the periodic point-charge array ; the residual error is shown in Fig 1(d). To treat the boundary between the QM and MM regions, we –rst adjust the charges on the MM junction atoms using the bond charge-increment model described above.However, the problem of unrealistic electrostatic interactions between the MM junction and QM termination atoms remains. To overcome this we have chosen to redistributeP.Sherwood et al. 83 Fig. 1 Electrostatic potentials of the silicate : (a) periodic HF 3-21G calculation ; (b) –tted periodic point-charge potential ; (c) potential due to the corrected classical cluster (model 2) ; (d) diÜerence between the two potentials (b) and (c). Contours for (a), (b) and (c) at 12.5 kJ mol~1, (d) at 2.1 kJ mol~1. the charge of the junction atom amongst the adjacent MM sites, and to correct for the dipole of these charge shifts by adding a point dipole, implemented in this study by a pair of point charges, at the MM sites that receive the charge.The result of this procedure is illustrated in Fig. 2(a), which should be compared with the periodic HF result [Fig. 1(a)]. For comparison, the potential due to the bare cluster is illustrated in [Fig. 2(b)]. Energetics and spectroscopy of water adsorption Electronic structure calculations were carried out on the T3 cluster shown in Fig. 3 using these two sets of PDC. A 6-31G** and a larger basis [6- were 311]]G(2d,2p)] employed and electron correlation was included using a DFT method employing a hybrid functional (B3LYP).23 In this study, we utilise a modi–ed form of the redundant coordinate optimisation technique, implemented within the Gaussian9424 package to perform constrained optimisations of the cluster within the point-charge –eld.In particular, we freeze the positions of the cluster-termination hydrogen atoms to lie along the84 Simulation of zeolites using embedded cluster methods Fig. 2 (a) Model 2, electrostatic potential for the QM»MM cluster ; (b) potential for the bare cluster.Contours at 12.5 kJ mol~1. crystallographic SiwO bonds at 1.48 from the Si. Furthermore, the two Si atoms of ” the cluster are also kept –xed in their crystallographic positions. Thus, the unit, AlO4 acidic proton and substrate are allowed to relax fully within a rigid zeolite framework. These constraints included modi–cation to the cartesian to internal coordinate transformations by increasing the relative mass of the frozen atoms.Thus the vibrational frequencies re—ect an immobile framework and the constrained optimisations are, indeed, true minima in the reduced coordinate space. For PDC model 2, it was also necessary to freeze the two oxygen atoms at their crystallographic positions to maintain a reasonable SiO orientation and prevent collapse to We found that this AlO3]SiOH.further constraint results in a very similar bare cluster geometry. NMR shieldings were evaluated at the optimised geometries using the HF» GIAO25,26 approach with the TZP basis set of Ahlrichs,27 and PDC, where appropriate. Absolute shieldings can be compared with experimental shifts using the equation dTMS(complex)\dTMS(internal reference)]p(internal reference)[p(complex) Fig. 3 Acidic cluster Si2AlO4H9P. Sherwood et al. 85 Fig. 4 Methyl-shift reactions over a zeolite acid site As is usual, the substrate molecule is the internal reference from which the hydrogen shifts are measured. We use the experimental gas-phase shift of 0.73 ppm for water.28 Methyl-shift reaction We use a T3 cluster, as in the embedding calculation, but with hydrogen atoms terminating the aluminium atom (see Fig. 4) and without the geometric constraints. These calculations were carried out at a fairly low level, HF/3-21G, to explore the important features of the model. At the mechanical embedding level, we have used the aluminosilicate consistent force-–eld (CFF) due to Hill and Sauer21 for all zeolite interactions. The hydrocarbon fragment was uncharged, and all atoms within the QM fragment were given zero charge.Thus, there is no electrostatic interaction between the excised QM system and the MM system. The QM fragment includes terminating H atoms to replace the O atoms, bonded to the junction Si and Al atoms, constrained to lie along the QM»MM bond.The OwSiwO angle bend is assumed to be described correctly by the QM OwSiwH force, and the Si/AlwO(MM) stretches are described by the CFF terms, constraining the Si/AlwH distances to a –xed length. During optimisation of the QM cluster, the surrounding MM region was allowed to relax up to three bonds from the junction region. Non-bonded framework»hydrocarbon interaction parameters were the MM2 Lennard- Jones parameters taken from an MD study of absorption of propane in all-silica zeolites by Nicholas et al.29 In calculations at the electrostatic embedding level we utilised PDC (model 2) rather than the CFF charges, with cluster-termination eÜects on the potential at the boundary of the QM region corrected using the procedures described above for model 2.Computation of the eÜects of MM polarization were computed using the wavefunction and geometry from the electrostatic embedding calculations. The representation of the MM polarizability consisted of distributed polarizabilities at both Si and O atoms. (6.434 for Si, 6.074 for O).30 The number of MM atoms bearing a polarizability a0 3 a0 3 was limited to ca. 200 to save computer time; however, enough polarizabilities are taken into account to have converged the polarization energy (using 500 polarizabilities yielded less than 1 kJ mol~1 extra stabilization). The calculations were performed using a hybrid QM»MM code (ChemShell) coupling the GAMESS-UK QM code31 with the MM code, DL-POLY32, modi–ed to support the CFF force –eld. The MM polarisation calculations were carried out at –xed geometries utilising the direct reaction –eld version of HONDO8.1.33 We have investigated the eÜects on the reaction barrier of the methyl-shift reaction in the propenium ion of a single, isolated acid site in zeolite Y at the three progressive levels of QM»MM embedding described above.Hydrocarbon zeolite clusters of 35 radius from the single a0 Al-substituted T site were employed.86 Simulation of zeolites using embedded cluster methods Computational results Zeoliteñwater clusters Structures and binding energies.The calculated structures are summarised in Table 1. Both PDC models and diÜerent basis sets reveal the same trends. For the bare cluster there is a small increase in the OwH length due to the PDC.All calculations show a similar structure for the 1:1 complex with water. The water is bound as the neutral species, with a quite short hydrogen bond involving the acid site, and a longer Br‘nsted secondary hydrogen bond involving a framework oxygen atom. The eÜect of the PDC is to increase the zeolite OwH bond length, increasing the primary hydrogen-bond interaction and decreasing the secondary interaction. A similar eÜect has been found from bare cluster models, on increasing the size of the cluster.Thus, Krossner and Sauer34 studied both T3 clusters and a substantially larger model which is typical of the faujasite lattice. They found that, as the size of the model increases, the larger intermolecular distance [O(7)wH(11)] becomes longer and the shorter one [O(9)wH(3)] shortens. A comparison of the two PDC models at the 6-31G** level reveals an encouraging level of agreement between the predicted structures, with the intermolecular hydrogen bond lengths being well within 0.1 The eÜect of increasing the basis set size, investi- ”.gated for PDC model 1, is to increase the hydrogen-bond lengths, an eÜect observed previously in calculations on simpler hydrogen-bonded systems.The calculated binding energies of the water molecule to the acid site are Br‘nsted given in Table 2. These values re—ect the diÜerent intermolecular distances predicted by the various treatments. Thus, the eÜect of the PDC is to increase the binding energy by ca. 15 kJ mol~1 for all models considered. The binding energies are considerably smaller for the large basis set.It is of interest to note that our value for the large basis set, in the Table 1 Cluster bond lengths and 1:1 complexes with watera (”) cluster complex O(2)wH(3) O(2)wH(3) O(9)wH(3) H(11)wO(7) O(9)wH(10) O(9)wH(11) model 1 6-31G**/noneb 0.970 1.039 1.482 1.717 0.967 0.994 6-31G**/PDC 0.974 1.089 1.352 1.924 0.968 0.979 6-311]]G(2d,2p) 0.965 1.019 1.543 1.753 0.962 0.987 /none 6-311]]G(2d,2p) 0.970 1.056 1.424 2.164 0.964 0.968 /PDC model 2 6-31G**/none 0.969 1.042 1.479 1.713 0.967 0.995 6-31G**/PDC 0.973 1.098 1.350 1.835 0.968 0.987 a See Fig. 3 for atom labelling.b Bare clusters with no point-charge –eld. Table 2 Binding energies (kJ mol~1) of water model 1 6-31G**/none 92.6 6-31G**/PDC 105.1 6-311]]G(2d,2p)/none 64.8 6-311]]G(2d,2p)/PDC 79.4 model 2 6-31G**/none 102.0 6-31G**/PDC 114.5P.Sherwood et al. 87 presence of the PDC (79.4 kJ mol~1) is close the value found for the faujasite model of Krossner and Sauer34 (73.3 kJ mol~1). Vibrational frequencies. There have been a number of calculations of vibrational frequencies of bare water»zeolite clusters.35h37 The vibrational frequencies that we have calculated at the various levels are shown in Table 3.The changes in these values follow the structural changes already discussed. For the bare cluster, the OwH-stretching frequency is noticeably reduced by the presence of the PDC. The eÜect of the PDC on other hydrogen atom motions (d, c) is less readily explained. Considerable changes in the vibrational frequencies of both the zeolite and substrate water are found on binding.The OwH stretching frequency of the zeolite is reduced and the other modes (d, c) are at higher energy, in line with the increase in the cluster OwH bond length and restrictions on OwH deformation upon hydrogen bonding. The eÜect of the PDC is to increase all these vibrational shifts, particularly that of the OwH-stretching frequency, which correlates with the increased binding energy predicted by the PDC models.The energies of the OwH deformation modes are relatively insensitive, both to the basis set level and to the inclusion of the PDC, being at ca. 1100 cm~1 [c(OH)] and ca. 1500 cm~1 [d(OH)]. The predicted value of the OwH-stretching frequency in the complex diÜers by nearly 1000 cm~1 for the various models used.All models predict a substantial lowering of this frequency on the inclusion of the PDC, in line with the increased hydrogen-bond strength. It is probable that the values at the 6-31G** level are too small, owing to the overestimation of the intermolecular interactions. Thus, the most realistic value that we have obtained is probably 2192 cm~1 at the level.It has been 6-311]]G(2d,2p)/PDC observed that a number of functionals used in DFT calculations predict too great a red shift in the OwH-stretching frequency of a hydrogen-bond donor on hydrogen bond formation.38 However, the functional that we have used here, B3LYP, gives a value for this frequency close to the MP2 value for the case of the water dimer,38 so that our calculation may not give a value that is unduly low.There has been considerable discussion on the interpretation of the IR spectra of water»zeolite complexes. The adsorption of water onto H-ZSM-539 yields two new bands at 3700 and 3600 cm~1 and three bands at 2885, 2457 and 1630 cm~1. These bands were previously interpreted in terms of an ion-pair complex involving a hydroxonium cation.39 An alternative assignment, consistent with the general theoretical conclusion that the energetically most favourable complex involves neutral water,1,34h36 with no proton transfer from the zeolite, explains the observed bands in terms of the so»called A»B»C pattern found for strong H»bonded complexes.40,41 This characteristic triplet of bands is caused by resonant interactions between the OH-stretching vibration, strongly broadened by interaction with low-frequency O… … …O modes, and overtones of the bending modes of the strongly perturbed bridging OH group.The strong double maximum (A»B) observed in the experimental spectrum39 at 2885 and 2457 cm~1 may be attributed to the zeolitic OH vibration, downshifted by the H-bond interaction and strongly broadened by coupling O… … …O modes.The splitting of this band at 2675 cm~1 arises by interaction with the second overtone of the dOH bending mode, placing the fundamental at 1350 cm~1. This can be compared with our highest level result of 1513 cm~1. The considerable red shift of the OH-stretching mode that we predict, giving a frequency of 2192 cm~1, is perhaps a little low when compared with the experimental estimate in the region 2700»2300 cm~1.More controversy, however, concerns the assignment of the experimental peak at 1630 cm~1, viewed by some workers34,35,37 as a third part of the A»B»C triplet. The splitting at 2000 cm~1 may be due to interaction with the second overtone of the cOH bending mode or with the combination. However, an alternative cOH]dOH explanation36 considers such a resonance process unlikely and assigns this peak to the HwOwH deformation mode of the absorbed water. Our highest level results place the88 Simulation of zeolites using embedded cluster methods Table 3 Harmonic vibrational frequencies for cluster and 1:1 complex with water (cm~1) cluster complex l[O(2) H(3)] d[O(2) H(3)] c[O(2) H(3)] l[O(2) H(3)] d[O(2) H(3)] c[O(2) H(3)] l[O(9) H(10)] l[O(9) H(11)] model 1 6-31G**/none 3803 1057 321»334 2516 1516 1113 3845 3346 6-31G**/PDC 3736 1051 308»421 1861 1551 1167 3843 3629 6-311]]G(2d,2p)/none 3790 1045 276»378 2735 1470 1010»1019 3863 3391 6-311]]G(2d,2p)/PDC 3725 1088 313»418 2192 1513 1101»1109 3852 3738 model 2 6-31G**/none 3819 1050 245»349 2477 1490 1108 3851 3324 6-31G**/PDC 3766 1110 231»394 1823 1558 1215 3846 3494P.Sherwood et al. 89 mode at 1105 cm~1 and the HwOwH deformation at 1634 cm~1, consistent with cOH both explanations. The two high-energy bands, observed at 3700 and 3600 cm~1 can be assigned to the OH stretch of the free proton and the H-bonded proton, respectively, calculated to occur at 3852 and 3738 cm~1 at our highest level.The corresponding values calculated in the absence of the PDC, 3863 and 3391 cm~1 are quite close to the MP2 values of Krossner and Sauer34 (3740, 3328 cm~1) and show the improved agreement obtained by inclusion of the PDC. 1H NMR shifts. The experimental NMR data, for loadings of one water molecule per bridging hydroxy group, are reported to be 6.2 ppm (H-Y),42 4.3»5.8 ppm (H-Y),43 7.1 ppm (H-ZSM-5)44 and a broad peak at 7 ppm (H-rho).45 For model 1, all calculations predict an average shift in the range 7»9 ppm, with the value being larger when the PDC is included (see Table 4).Thus, on inclusion of the PDC, the shift of the zeolitic hydrogen increases, and that of the water proton involved in H bonding decreases, in line with the increase in the primary H-bonding interaction, and the decrease in the secondary H bonding upon inclusion of the PDC.A similar though considerably smaller eÜect was observed by Krossner and Sauer,34 when bare clusters of diÜerent sizes were studied. Thus, since comparison with experiment is not helpful in assessing the relative quality of the results from model 1, results for model 2 are not presented.Methyl-shift reaction Table 5 summarises the results for the methyl-shift reaction on zeolite Y at the three diÜerent levels of embedding. Activation barriers were obtained for the propoxide(I)» transition state (TS) (barrier I) and propoxide(II)»TS (barrier II) processes in zeolite Y, propoxide(I) and propoxide(II) being the minima in which a CwO bond is formed on either side of the Al (see Fig. 4). The O atoms involved are the O(1) and O(4) zeolite sites. For the gas-phase calculations, the barrier is 316 kJ mol~1, and no distinction can be made between the two minima. The eÜects of mechanical embedding on the geometries of the minima and the TS connecting them (see Fig. 5) are shown in Table 6. Looking closely at these gas-phase and embedded structures it can be seen that the con–ning of the acid site reduces the OwAlwO angle and pushes the propenium ion into the cage.The propenium ion is further encouraged into the cage by the Van der Waals forces of the surrounding cage (zeolite hydrocarbon non-bond interactions). The eÜects of the relaxation of the zeolite cage (MM energy) is seen to also reduce the barrier. The p-bond in the minima puts more constraints on these structures than on the loose TS structure, thus putting the minima at a disadvantage.If the QM contribution to the energy of the embedded clusters is compared with the energy of the corresponding gas-phase systems it can be seen that the energy cost of deforming the QM cluster to –t the zeolite backbone is greater for the minima than for the TS.The resulting stabilisation of the TS (by 55 kJ mol~1 relative to minimum I and 93 kJ mol~1 relative to 1H NMR shieldings Table 4 H2O d(ZOH) da(H2O) db(H2O) ave. model 1 6-31G**/none 14.0 8.1 1.6 7.9 6-31G**/PDC 17.9 6.0 3.3 9.1 6-311]]G(2d,2p)/none 12.4 7.5 1.4 7.1 6-311]]G(2d,2p)/PDC 16.0 4.4 3.2 7.9 a Hydrogen bonded. b Free hydrogen.90 Simulation of zeolites using embedded cluster methods Table 5 Energies, activation barriers (kJ mol~1) and dipole moments (D) of the QM cluster for methyl-shift reaction in propene on a zeolite acid site at diÜerent embedding levels model energy propoxide I TS propoxide II barrier I barrier II gas phase totala 0 316 0 316 316 mechanical totala 0 247 55 247 192 QMa 0 261 38 261 223 MMa,b 0 [6 9 [6 [15 Z»(C,H) non-bond [12 [20 [4 [8 [16 l (D) 6.3 10.9 6.1 electrostatic totala 0 253 68 253 185 QM»MM elecc [93 [103 [100 [10 [3 polarised QM»MM pold [30 [45 [33 [15 [12 a Energies for propoxide I, propoxide II and TS given relative to propoxide I.b Excluding the zeolite»hydrocarbon non-bonded interaction [Z»(C,H) non-bond]. c Electrostatic interaction between the DPC representation of the embedded QM cluster and the MM point charges.d Interaction of the DPC representation of the embedded QM cluster with the distributed polarizabilities on the MM atoms. minimum II) is the main contribution to the lowering of the barrier in the embedded case. To get an impression of the eÜects to be expected by electrostatic embedding, the dipole moments of the clusters at the mechanical embedding level are also presented in Table 5.Analysis of the diÜerent contributions to the barrier at the electrostatic embed- Fig. 5 Structural parameters for the methyl-shift reaction Table 6 Optimised geometries (in and degrees) for the methyl-shift ” reaction (see Fig. 5) using the mechanical embedding scheme. propoxide TS structural coordinate gas phase embedded gas phase embedded O(1)AlO(4) 110.6 92.4 106.4 99.5 SiO(4)Al 170.6 133.5 139.6 133.8 SiO(1)Al 114.2 124.7 139.2 130.1 hCXC 130.1 128.5 90.1 88.6 rXAl » » 3.37 3.44 rXC 1.95 1.90 1.79 1.80P.Sherwood et al. 91 ding level is not as straightforward as at the mechanical embedding level. Here, we give the interaction between the QM system and the surrounding point charges by representing the QM charge distribution by the so-called dipole preserving charges (DPCs).46 The computed interactions are all attractive (by 93, 100 and 103 kJ mol~1 for the two minima and the TS, respectively) but do not diÜer strongly for the diÜerent points along the reaction path.The additional eÜects of electrostatic embedding on the barrier are, therefore, not nearly as large as those resulting from the mechanical constraints, lowering barrier I by 10 kJ mol~1, and lowering barrier II by 3 kJ mol~1.This barrier change is relatively insensitive to changes in the MM charges, changing by no more than 5 kJ mol~1 when the bond contribution to the oxygen charge is changed from ^0.46 to ^0.30. The build-up of dipole moment on going from minimum to TS is also expected to aÜect the MM polarization contribution to the barrier.We can estimate this eÜect by computing the interaction energy for the DPC representation of the QM region and distributed atomic polarisabilities at the MM atom sites.47 The eÜect on the barrier of the MM polarization is found to be at least as big as that of the electrostatic interaction, even though the magnitude of the interaction is signi–cantly smaller (30»45 kJ mol~1).Discussion We have described two alternative strategies for obtaining PDC, which can be used in a embedded cluster model of zeolite reactivity to include the electrostatic environment and to take account of termination eÜects. When used in a model of water bound to a hydrogen-terminated T3 cluster, both give very similar eÜects, as judged by predicted structures, substrate binding energies and vibrational frequencies.The inclusion of the PDC results in substantial changes in all of these three properties compared with the bare cluster results. The trends we –nd on inclusion of the PDC mirror those found in bare cluster calculations on increasing the size of the cluster, although the changes are in general larger than those from the cluster calculations.It is clearly computationally prohibitive to use sufficiently large bare clusters that eÜectively reproduce the results of a full periodic calculation. For example, although quite large clusters have been used to predict structures and substrate binding energies, the determination of vibrational frequencies is more computationally demanding.A computational strategy that involves a hybrid QM»MM approach has obvious advantages and the consideration of the electrostatic potential for such a model, described in this paper, is a necessary step towards the construction of such a model. The water absorption and methyl-shift reaction show an interesting contrast in the in—uence of the charge –eld in the zeolite pore.Whereas the structural changes of the water adsorbed on the site induce large vibrational frequency shifts, the structures and barriers in the methyl shift appear rather insensitive to the electrostatic –eld. The latter do show, however, that electronic polarization of the surroundings must also be taken into account to complete a realistic modelling of zeolite chemistry. thank EPSRC and Shell Research and Technology Center in Amsterdam for We support of this research.We would like to thank A. M. Rigby and M. V. Frash for providing gas-phase structures for, and discussions on, the methyl-shift reaction. References 1 R. A. van Santen and G. J. Kramer, Chem. Rev., 1995, 95, 637. 2 M. Iwamoto, H. Yahiro, N. Mizuno, W. X. Zhang, Y. Mine, H. Furukawa and S.Kagawa, J. Phys. Chem., 1992, 96, 9360. 3 R. A. van Santen, Stud. Surf. Sci. Catal., 1994, 85, 273. 4 U. Fleischer, K. Kutzelnigg, A. Bleiber and J. Sauer, J. Am. Chem. Soc., 1993, 115, 7833.92 Simulation of zeolites using embedded cluster methods 5 L. A. Curtiss, H. Brand, J. B. Nicholas and L. E. Iton, Chem. Phys. L ett., 1991, 184, 215. 6 H. V. Brand, L. A. Curtiss and L.E. Iton, J. Phys. Chem., 1992, 96, 7725. 7 M. S. Stave and J. B. Nicholas, J. Phys. Chem., 1993, 97, 9630. 8 E. H. Teunissen, C. Roetti, C. Pisani, A. J. H. de Man, A. J. P. Jansen, R. Orlando, R. A. van Santen and R. Dovesi, Model. Simul. Mater. Sci. Eng., 1994, 2, 921. 9 R. Shah, J. D. Gale and M. C. Payne, J. Phys. Chem., 1996, 100, 11688. 10 J. Gao, in Reviews of Computational Chemistry, ed.K. B. Lipkowitz and D. B. Boyd, VCH, New York, 1995, vol. 7, p. 119. 11 M. J. Field, P. A. Bash and M. Karplus, J. Comput. Chem., 1990, 11, 700. 12 U. C. Singh and P. A. Kollman, J. Comput. Chem., 1986, 7, 718. 13 S. P. Greatbanks, P. Sherwood and I. H. Hillier, J. Phys. Chem., 1994, 98, 8134. 14 J. R. Hill and J. Sauer, J. Phys. Chem., 1994, 98, 1238. 15 B. W.H. van Beest, G. J. Kramer and R. A. van Santen, Phys. Rev. L ett., 1990, 64, 1955. 16 K. P. Schroé der, J. Sauer, M. Leslie, C. R. A. Catlow and J. M. Thomas, Chem. Phys. L ett., 1992, 188, 320. 17 S. P. Greatbanks, I. H. Hillier, N. A. Burton and P. Sherwood, J. Chem. Phys., 1996, 150, 3770. 18 S. P. Greatbanks, I. H. Hillier and P. Sherwood, J. Comput. Chem., 1997, 18, 562. 19 S.J. Collins, A. H. de Vries, P. Sherwood, S. P. Greatbanks and I. H. Hillier, J. Comput. Chem., submitted. 20 D. Bakowies and W. Thiel, J. Phys. Chem., 1996, 100, 10580. 21 J. R. Hill and J. Sauer, J. Phys. Chem., 1995, 99, 9536. 22 S. T. Sie, Ind. Eng. Chem. Res., 1992, 31, 1881. 23 A. D. Becke, J. Chem. Phys., 1993, 98, 1372; P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M.J. Frisch, J. Phys. Chem., 1994, 98, 11623. 24 M. J. Frisch, G. W. Trucks, H. B. Schlegel, P. M. W. Gill, B. G. Johnson, M. A. Robb, J. R. Cheeseman, T. A. Keith, G. A. Petersson, J. A. Montgomery, K. Raghavachari, M. A. Al-Laham, V. G. Zakrzewski, J. V. Ortiz, J. B. Foresman, J. Cioslowski, B. B. Stefanov, A. Nanayakkara, M. Challacombe, C. Y. Peng, P. Y. Ayala, W. Chen, M. W.Wong, J. L. Andres, E. S. Replogle, R. Gomperts, R. L. Martin, D. J. Fox, J. S. Binkley, D. J. Defrees, J. Baker, J. P. Stewart, M. Head-Gordon, C. Gonzalez and J. A. Pople, GAUSSIAN 94, Revision A.1, Gaussian Inc., Pittsburgh, PA, 1995. 25 R. Ditch–eld, Mol. Phys., 1974, 27, 789. 26 K. Wolinski, J. F. Hilton and P. Pulay, J. Am. Chem. Soc., 1990, 112, 8251. 27 A. Schaé fer, H.Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571. 28 W. G. Schneider, H. J. Bernstein and J. A. Pople, J. Chem. Phys., 1958, 28, 601. 29 J. B. Nicholas, F. R. Trouw, J. E. Mertz, L. E. Iton and A. J. Hop–nger, J. Phys. Chem., 1993, 97, 4149. 30 E. de Vos-Burghart, V. A. Verheij, H. van Bekkum and B. van de Graaf, Zeolites, 1992, 12, 183. 31 GAMESS-UK is a package of ab initio programs written by M. F. Guest, J. H. van Lenthe, J. Kendrick, K. SchoeÜel and P. Sherwood, with contributions from R. D. Amos, R. J. Buenker, M. Dupuis, N. C. Handy, I. H. Hillier, P. J. Knowles, V. Bonacic-Koutecky, W. von Niessen, A. P. Rendell, V. R. Saunders and A. J. Stone. The package was derived from the original GAMESS code due to M. Dupuis, D. Spangler and J. Wendoloski, NRCC Software Catalog, vol. 1, Program QG01 (GAMESS) 1980. 32 DL POLY is a package of molecular simulation routines written by W. Smith and T. R. Forester, copyright The Council for the Central Laboratory of the Research Councils, Daresbury Laboratory, nr. Warrington, 1996. 33 M. Dupuis, A. Farazdel, S. P. Karma and S. A. Maluendes, in MOT ECC-90, ed. E. Clementi, ESCOM, Leiden, 1990, p. 277. 34 M. Krossner and J. Sauer, J. Phys. Chem., 1996, 100, 6199. 35 A. G. Pelmenschikov and R. A. van Santen, J. Phys. Chem., 1993, 97, 10678. 36 S. A. Zygmunt, L. A. Curtiss, L. E. Iton and M. K. Erhardt, J. Phys. Chem., 1996, 100, 6663. 37 J. Sauer, P. Ugliengo, E. Garrone and V. R. Saunders, Chem. Rev., 1994, 94, 2095. 38 P. E. Sinclair and C. R. A. Catlow, J. Chem. Soc., Faraday T rans., 1997, 93, 333. 39 A. Jentys, G. Warecka, M. Derewinski and J. Lercher, J. Phys. Chem., 1989, 93, 4837. 40 M. F. Clydon and N. Sheppard, Chem. Commun., 1969, 1431. 41 A. G. Pelmenschikov, J. H. M. C. van Wolput, J. Janchen and R. A. van Santen, J. Phys. Chem., 1995, 99, 3612. 42 M. Hunger, D. Freude and H. Pfeifer, J. Chem. Soc., Faraday T rans., 1991, 87, 657. 43 P. Batamack, C. Doremieux-Morin and J. Fraissard, J. Chim. Phys., 1992, 89, 423. 44 P. Batamack, C. Doremeiux-Morin, J. Fraissard and D. Freude, J. Phys. Chem., 1991, 95, 3790. 45 Z. Luz and A. J. Vega, J. Phys. Chem., 1987, 91, 374. 46 B. T. Thole and P. Th. van Duijnen, T heoret. Chim. Acta, 1987, 63, 209. 47 B. T. Thole and P. Th. van Duijnen, T heoret. Chim. Acta, 1980, 55, 307. Paper 7/01790A; Received 14th March, 1997