Faraday Discuss., 1997, 106, 41»62 Structure and reactivity of silica and zeolite catalysts by a combined quantum mechanicsñshell-model potential approach based on DFT Marek Sierka and Joachim Sauer*§ Humboldt zu Berlin, Institut Chemie, Arbeitsgruppe Quantenchemie, Universitaé t fué r 10/11, D-10117 Berlin, Germany Jaé gerstr. An ion-pair shell-model potential with functional parameters derived from the results of quantum mechanical density functional theory (DFT) calculations on small molecular models is presented.It is used to predict the structure and properties of diÜerent silica and zeolite catalysts. Characteristic diÜerences between the Hartree»Fock and DFT structures of quartz, silica sodalite and silicalite are revealed. A combined quantum mechanics» ion-pair shell-model potential scheme is presented and applied to embedded cluster calculations on catalytically active sites in periodic framework structures. 1 Introduction Understanding the catalytic activity of solids at the molecular level requires assistance from accurate quantum chemical calculations. Advances in computer technology, e.g. the advent of parallel machines and recent developments in quantum chemical codes, permit more complex calculations on more realistic model systems, often involving more than 1000 basis functions and more than 100 atoms, depending on symmetry.Zeolites appear as particularly interesting systems for such studies because of their well de–ned framework structures, the wealth of experimental studies, and their technological importance.Quantum chemical calculations for solids, which, microscopically, are threedimensional in–nite systems, are computationally very demanding, even when applying periodic boundary conditions. Active sites or defects introduce further complications owing to broken space and translational symmetry, which require larger pseudo-unit cells. One solution to this problem is the –nite cluster approach. The treatment is limited to a part of the system and the in—uence of the crystalline environment is neglected or only included approximately.Thus, cluster models are particularly suited to describe local phenomena, e.g. catalytically active sites. There exist, however, classes of problems which require theoretical predictions of the structure of the whole system.The interatomic potential functions (also known as force –elds) serve these needs. They provide an analytical approximation to the potential energy surface of the system. Traditionally, empirical interatomic potentials are used. However, for many but the simplest systems, there is incomplete information for parameter determination. Therefore, the following procedure has been suggested :1h3 quantum chemical calculations are performed for –nite molecular models of the periodic zeolite structures. The results obtained are used as a data base to –t the parameters of the interatomic potential functions.These parameters are assumed to be transferable and § E-mail: js=qc-ag-berlin.mpg.de 4142 Structure and reactivity of silica and zeolite catalysts are used for calculating the energies of the periodic structures and the forces on their atoms.This approach has several advantages. The set of data that may be included in the –t is, in principle, unlimited and there is a one-to-one correspondence between the ìcalculatedœ (by the potential function) and ìobservedœ (by the quantum chemical calculations) data. As ìobservablesœ serve the energy and its –rst and second derivatives with respect to the displacement of all nuclei, ìobservablesœ can be generated, not only for the energy minimum structure of the models but also for an, in principle, unlimited number of distorted structures.A similar procedure has been applied by other authors.4 The disadvantages of interatomic potential functions are that no wavefunction is obtained and that reactions are difficult to model.Therefore, the idea emerges to combine the virtues of the quantum mechanical cluster calculation (QM) with a lattice energy minimization of the periodic solid using interatomic potential functions (Pot). Such a combined scheme, QM»Pot, is particularly useful for studying the properties and reactivity of an active site in diÜerent environments.Our implementation uses the ion-pair shell-model potential,5 which takes into account the polarization of the environment by the active site, to describe the periodic crystal. It accounts for both the long-range crystal potential and the ìmechanicalœ constraints of the position of the atoms of the cluster by the periodic lattice. In principle, any of the quantum chemical methods available can be combined with interatomic potential functions.However, for an accurate description of the active acidic sites in zeolites and of reactions involving such sites, inclusion of electron correlation proved to be important. For example, inclusion of electron correlation was found to be essential in order to correctly describe the interaction of methanol with zeolite Br‘nsted sites.6 Among quantum chemical methods incorporating correlation eÜects, for instance MP2, MP4 or CI, density functional theory (DFT) is particularly attractive since, for large systems, the computational eÜort is much less than for even the simple MP2 approximation.In practice, the computational cost of DFT calculations grows as nF2.6 with the number of basis functions, compared to the nF4.0 dependence of efficient MP2 codes.7 In this paper we –rst derive an ion-pair shell-model potential for silica and zeolites based on DFT calculations using the B3-LYP exchange-correlation functional.8 We then use this DFT parametrized potential to predict the structures and vibrational properties of dense and microporous silica polymorphs.A systematic overestimation of SiwO bond lengths is observed, which is shown to be an inherent feature of the DFT B3-LYP method for the basis set used.This conclusion is based on DFT calculations for octahydridosilasesquioxane, and disiloxane, We then reveal charac- H8Si8O12 H6Si2O. teristic diÜerences between the HF and the DFT potential energy surfaces for quartz, sodalite and silicalite.The former tends to predict the higher symmetry structures incorrectly as the most stable ones. Finally, the performance of the potential for protonated zeolites is discussed and the embedded cluster approach (QM»Pot scheme) is applied to sites in faujasite. Br‘nsted 2 Methods and computational details 2.1 Quantum chemical calculations The cluster calculations use the DFT together with the B3-LYP8 exchange-correlation functional.For the embedded cluster calculations we use the TURBODFT9 program. Unless speci–ed explicitly, we apply the fully optimized basis sets from Ahlrichsœs group,10 double zeta (DZ) for Si, Al and H and the triple zeta (TZ) for O. Polarization functions with the following exponents are added: 0.4 (Si), 0.35 (Al), 1.2 (O) and 0.8 (H).This combination is designated as T(O)DZP basis set. The quantum chemical program package GAUSSIAN9411 is employed for the free cluster calculations. For numericalM. Sierka and J. Sauer 43 integration in GAUSSIAN94 we use the option ì–negridœ which produces a grid formally consisting of 75 radial and 302 spherical points per atom. The DFT calculations in TURBODFT are made with the following grid (see ref. 9 for details) : for elements of the –rst row 5890 points (35 radial and 302 spherical), for the second row 9212 points (40 radial and 434 spherical), and for the third row 10 542 points (45 radial and 434 spherical). The calculations of absolute shielding constants p for the optimum structures of the embedded clusters obtained using the combined QM»Pot scheme are performed within the coupled Hartree»Fock approach using the gauge including atomic orbitals (GIAO)12 implemented in the SHEILA13 module of the TURBOMOLE14 program.The TZP basis sets from Ahlrichsœs group10 are used on all the atoms. The relative 1H chemical shifts are calculated using methanol as internal reference (secondary dTMS standard) :15 dTMS(cluster)\dTMS(CH3OH)]p(CH3OH)[p(cluster) The calculated absolute shielding constants were 32.71 and 31.79 ppm at the p(CH3OH) HF and DFT B3-LYP optimized structure, respectively, while the experimental gas phase value is 0.02 ppm.dTMS(CH3OH) 2.2 Database generation Since four-, –ve- and six-membered aluminosilicate rings are characteristic secondary building units of zeolites, they are adopted as models and saturated with hydrogen atoms (Fig. 1). For a proper description of the aluminium-rich zeolites, we consider all possible Al for Si substitutions permitted by Loé wensteinœs rule, which forbids that Al atoms occupy neighbouring tetrahedral framework sites. The equilibrium structures of Fig. 1 Molecular models used for generating the database using the DFT B3-LYP method and the symmetry constraints applied in structure optimizations.For simplicity the oxygen atoms and terminal OwH groups are omitted.44 Structure and reactivity of silica and zeolite catalysts all the models are determined by imposing symmetry restrictions. This speeds up the calculations and avoids the formation of non-realistic intramolecular hydrogen bonds involving terminal H atoms.The point group chosen is a compromise between a small number of imaginary vibrational frequencies and a high order of the symmetry group. In two cases, cyclotetrasilicic and cyclopentasilicic acid the optimized structures (S4) (D5), prove true minima by frequency calculations. For all the equilibrium structures, the force constant matrix is calculated and used as data in the potential-–tting procedure.The only exception is Al-cyclohexasilicic acid as its low symmetry and large size (36 (Cs) atoms) do not allow computation of force constants analytically. In addition, a set of ca. 60 distorted structures is generated by distorting the internal coordinates, starting from the equilibrium structures. The magnitude of typical distortions is ^10 pm for bond distances and ^4° for bond angles.For the resulting structures the DFT B3-LYP gradients are computed. The entire database used for the –tting procedure consists of a set of 10 equilibrium structures (coordinates and zero gradients) together with their harmonic force constants and a set of 60 distorted structures (coordinates) along with their gradients. 2.3 Parameter –tting We apply the same procedure as used previously by Schroé der and Sauer.1 DiÜerent ion types are assumed for the hydrogen and oxygen ions of bridging and terminal hydroxy groups and O, respectively) and diÜerent short-range parameters are –tted (Hb, Ob Ht , for the interactions between them.To obtain a well balanced parametrization, diÜerent weights are given to diÜerent types of data.The DFT force constants are weighted with a factor 5]10~4 since their number is ca. 2000 times larger than the number of corresponding gradient components. This results in parameters which combine good predictions of structures with good predictions of vibrational spectra. To limit the in—uence of the terminal hydroxy groups, which are not present in the periodic structures, all data connected with them are weighted with a factor 0.1.All potential parameters are –tted simultaneously by the least-squares method. The –nal values are listed in Table 1. 2.4 Lattice-energy minimization Constant pressure minimizations of the lattice energy are performed with respect to positions of the ions, shells and the unit cell parameters. A cut-oÜ radius of 10 is ” chosen for the summation of the short-range interactions.For some of the equilibrium structures we determine properties such as static constants and relative permittivities, as well as the harmonic phonon frequencies. All calculations are performed using the GULP16 and METAPOCS17 codes. 2.5 Combined quantum chemicalñinteratomic potential approach The embedding scheme used18 partitions the entire system (S) into two parts : the inner part (I) usually containing the site in question e.g.the site, and the outer part Br‘nsted (O). The energy of the combined systems can be decomposed as : E(S)\E(I)]E(O)]E(I»O) (1) where E(I»O) is the interaction between the inner and the outer region. The energy of the outer part, E(O), and the interaction between the inner part and the outer part, E(I»O), are approximated by the potential : E(S)BEQM(I)]EPot(O)]EPot(I»O) \EQM(I)]EPot(S)[EPot(I) (2)M.Sierka and J. Sauer 45 Table 1 Parameters of the DFT B3-LYPderived shell-model potential charges/e core shell Si 4.0a » Al 3.0a » O 1.228 58 [3.228 58 Ob 0.817 53 [2.817 53 Hb 1.0a » Ht 1.0a » short-range repulsion A/eV o/” SiwO 1612.459 20 0.299 55 SiwOb 997.880 97 0.332 12 AlwO 1395.774 63 0.304 49 AlwOb 1644.881 77 0.291 39 ObwHb 368.648 03 0.225 11 OwHb 7614.580 03 0.199 13 OwHt 772.068 14 0.185 24 core»shell interaction k/eV ”~2 O 122.478 53 Ob 70.151 23 three-body interaction kb/eV rad~2 H0/rad OwSiwO 0.144 703 109.47a ObwSiwO 0.384 711 109.47a OwAlwO 0.893 930 109.47a ObwAlwO 0.686 678 109.47a a Not adjusted in the –tting procedure.Hence, the energy of the total system is obtained approximately from the QM energy of the internal part and the diÜerence in the energies of the total system and the internal part, calculated by the interatomic potential function. For partially covalent solids, such as zeolites, the partitioning creates dangling bonds which need to be saturated by hydrogen atoms, called link atoms.The inner part (I) together with link atoms (L) form the cluster (C). Now we have: E(S)\EQM(C)]EPot(S)[EPot(C) (3) only if : D\EQM(L)]EQM(I»L)[EPot(L)[EPot(I»L)\0 (4) where and are the QM and potential function interaction energies EQM(I»L) EPot(I»L) between the inner part and link atoms. Hence, the energy of the total system can only be obtained from the subtraction scheme, eqn.(3), if D is small and negligible. This is46 Structure and reactivity of silica and zeolite catalysts expected to be the case if the embedding potential is parametrized using the same quantum mechanical method as used for the cluster. The diÜerentiation of the above expression leads to the following formulae for the energy gradients : Fa(S)\Fa, QM(C)]Fa, Pot(S)[Fa, Pot(C) ; a ½ I (5) Fb(S)\Fb, Pot(S) ; b ½ O (6) For the atoms present in the inner region [a, eqn.(5)] all three components contribute, while forces acting on the atoms of outer region [b, eqn. (6)] are calculated from the interatomic potential function only. The link atoms are not independent variables in the structure optimization but are kept at a –xed distance from atom X on the bond XwY they terminate : E\E(rX , rH)\E[rX , rH(rX)] (7) where the atoms X and Y belong to the inner and outer part, respectively. Hence, the forces on the atoms X and Y have to be modi–ed:18 dE drX, m \ LE LrX, m ] LE LrH, m LrH, m LrX, m (8) The present implementation18 of the embedding scheme describes the outer part by a shell-model ion-pair potential.It couples the GULP program16 for the lattice energy minimization with TURBOMOLE14 and TURBODFT9 programs for the quantum mechanical calculations.When the HF method is applied to the cluster the HF-derived potential is used for the outer part, for the DFT calculations we apply our DFT-based shell-model potential. The hydrogen link atoms of the cluster models are kept at a –xed distance from the oxygen atoms they terminate.These distances are the optimum distances found in free cluster calculations. The (Si)OwH distances are 94.50 pm for the HF and 96.66 pm for the DFT calculations, and the (Al)OwH distances are 94.00 pm for HF and 96.28 pm for DFT calculations. All combined QM»Pot optimizations are performed without any symmetry constraints, i.e. the space group is assumed. P1 However, the unit cell parameters are –xed to values found by a lattice-energy minimization using the embedding shell-model potential alone.In the combined QM»Pot scheme the atoms are relaxed until the largest cartesian gradient component and the RMS gradient are less than 0.0001 Eh a0~1. The OwH stretching frequency is calculated using 5 points (step size 0.3 pm) on the combined QM»Pot potential energy surface along the bridged OwH bond.19 2.6 Deprotonation energy The calculation of deprotonation energies, i.e.the energy of the reaction : ZwOH]ZwO~]H` requires application of the combined QM»Pot scheme to the protonated and deprotonated acidic sites. If the QM cluster is large enough then all short-range terms resulting from the interatomic potential cancel approximately and the total reaction energy can be written as :20 *EQMhPot\*EQM@@QMhPot]*ELR@@QMhPot (9) The –rst term describes the ìmechanicalœ part arising from the fact that the structures of the embedded clusters are diÜerent from free space clusters.The second term describes the correction due to the crystal potential.M. Sierka and J. Sauer 47 Two problems arise when applying periodic boundary conditions to a deprotonation reaction.19 First, the calculation of the energy of the negatively charged faujasite framework is not possible using the traditional Ewald summation technique because of the in–nite Coulomb repulsion of charged unit cells.This can be avoided by applying a neutralizing homogeneous background charge distribution.21,22 Second, by taking the diÜerence between QM»Pot energies of the neutral and deprotonated zeolites, one obtains the energy for removing the proton from every unit cell.We are interested in the proton removal from a single acidic site. This can be obtained by removing the interaction between the charged defects in the anionic unit cells. We follow the method proposed by Leslie and Gillan.21 Details about the implementation are given elsewhere. 19 3 Results and Discussion 3.1 Structure and properties of silicates To test the quality of the DFT-derived potential we determined the equilibrium structures of a selected set of silica modi–cations. Table 2 shows the unit cell parameters obtained and compares them with predictions by two parameter sets derived previously from empirical5 and ab initio HF data.1 The cell parameters predicted by the present DFT-derived potential deviate by 1.4% on average from the observed data.This is a slight improvement compared with the HF-derived shell-model potential (1.9%). However, for the former, the deviations from the observed values are less consistent. In the case of a-quartz, the –rst cell parameter is smaller than the experimental one, while Table 2 Cell parameters of microporous and dense modi–cationsa,b SiO2 empirical HF shell DFT shell modi–cation shell modelc modeld model obsd.e faujasite A\B\C 2423 2463(1.5) 2466(1.6) 2426 sodalite A\B\C 882 895(1.4) 885(0.2) 883 theta-1 A 1382 1414(2.0) 1420(2.5) 1386 B 1739 1777(2.0) 1768(1.5) 1742 C 500 516(2.4) 502(0.4) 504 mordenite A 1802 1830(1.1) 1812(0.1) 1810 B 2004 2049(0.5) 2069(1.5) 2038 C 743 758(1.2) 762(1.7) 749 ZSM-5 A 1998 2043(1.6) 2020(0.4) 2011 B 1974 2021(1.7) 1996(0.4) 1988 C 1332 1363(1.9) 1348(0.8) 1337 90.8 90 90.7 90.7 zeolite RHO A\B\C 1477 1501(1.1) 1505(1.3) 1485 a-quartz A\B 484 499(1.4) 489(0.6) 492 C 535 551(1.8) 545(0.7) 541 b-quartz A\B 500 510(2.0) 511(2.2) 500 C 550 562(2.9) 560(2.6) 546 a-cristobalite A\B 497 513(3.0) 498(0.0) 498 C 701 727(4.6) 717(3.2) 695 b-cristobalite A\B\C 731 739(3.1) 746(4.0) 717 mean dev.f 0.7 1.9 1.4 a Units: pm.b Percentage deviation from observed data is given in parentheses. c Ref. 5. d Ref. 1. e Faujasite : ref. 23, sodalite : ref. 24, theta-1 : ref. 25, mordenite: ref. 26, ZSM-5: ref. 27, zeolite rho: ref. 28, a-quartz : ref. 29, b-quartz : ref. 30, a- and bcristobalite : ref. 31. f Mean percentage deviation from observed data.48 Structure and reactivity of silica and zeolite catalysts the HF-derived potential gives cell parameters slightly larger for all the structures. The predictions of the empirical shell-model potential are in better agreement with experimental data (only 0.7% average deviation).A further test of the potential is the comparison of calculated and measured elastic constants and relative permittivities of a-quartz (Table 3). The calculated elastic constants are of the same quality as those obtained with the empirical potential and the potential –tted to the HF data set. However, relative permittivities seem to be reproduced better with the present potential than with the HF parametrized one.As an example of the structure prediction using the DFT parametrized potential Table 4 compares observed and calculated bond lengths and angles for silica faujasite. The averaged calculated SiwO bond length is ca. 2.7 pm longer than the observed one. Table 3 Elastic constants and relative permittivities of a-quartz at 0 Ka empirical shell HF shell DFT shell modelb modelc model obsd.d C11 9.47 8.45 7.59 8.69 C33 11.61 9.63 10.56 10.60 C44 5.01 4.11 3.81 5.83 C66 3.82 3.65 2.61 3.99 C14 [1.45 [1.37 [0.79 [1.81 C13 1.97 2.39 2.89 1.19 C12 1.84 1.15 2.36 0.70 e11 4.74 4.07 4.47 4.52 e33 5.01 4.42 4.65 4.64 e= 2.12 1.76 1.86 2.40 a Units: elastic constants in 1010 N m~2.b Ref. 5. c Ref. 1. d Ref. 29. Table 4 Bond lengths and SiwOwSi bending angles in pure silica faujasite optimized with present DFT B3-LYP-derived potentiala,b empirical HF shell DFT shell shell modelc modeld model obsd.e r(SiwO1) 161.4(0.7) 162.2(0.9) 163.7(1.9) 160.7 r(SiwO2) 159.9(0.2) 160.9(0.8) 162.8(1.9) 159.7 r(SiwO3) 160.9(0.5) 161.7(0.8) 163.7(2.1) 160.4 r(SiwO4) 160.8(0.6) 161.6(0.1) 163.1(1.1) 161.4 average 160.8 161.6 163.3 160.6 mean dev.f 0.50 0.65 1.8 n(SiwO(1)wSi) 138 141(2.2) 139.1(0.6) 138 n(SiwO(2)wSi) 149 154(3.4) 146.9(0.9) 149 n(SiwO(3)wSi) 146 152(4.1) 144.2(0.8) 146 n(SiwO(4)wSi) 141 145(2.8) 143.2(1.0) 141 mean dev.f 0.0 3.1 0.8 a Units: distances in pm, angles in degrees.b Percentage deviation from observed data is given in parentheses. c Ref. 5. d Ref. 1. e Ref. 23. f Mean percentage deviation from observed data.M. Sierka and J. Sauer 49 The calculated bond angles agree very well with the experimental values. The average deviation, 0.8%, is lower than for the HF-derived potential, which overestimates the angles by more than 3%. Thus, the DFT-parametrized shell-model potential provides an excellent description of the SiwOwSi bond angles, but systematically overestimates SiwO bond distances.In contrast, the HF-parametrized shell-model potential gives good bond distances, but systematically too large bond angles. Do these diÜerences re—ect the eÜect of electron correlation on the potential-energy surface (PES) or are they a result of the representation of the PES by a parametrized shell-model potential ? To answer this question we performed test calculations on the disiloxane, and H6Si2O, octahydridosilasesquioxane molecules, because they contain the character- H8Si8O12 , istic SiwOwSi bond and both experimental and computational data are available for comparison. 3.2 DFT description of the SiwO bond 3.2.1 Octahydridosilasesquioxane. The highly symmetric molecule of octahydridosilasesquioxane (Fig. 2) is a model of the double four ring (D4R) H8Si8O12 found, e.g. in zeolite A. Its structure has been determined experimentally32 in the molecular crystal and the vibrational spectrum is also known.33 The equilibrium structure belongs to the point group32 but shows only small distortions from the symmetry Th Oh due to the crystalline environment. This allows one to perform calculations for the higher symmetry structure and thus to simplify the computational task.Structural (Oh) data predicted at the MP2 and HF level are available34 that used the same basis set as the present DFT study. Appropriate structural parameters are listed in Table 5. The SiwO distance predicted by DFT is ca. 3 pm too long, while the SiwOwSi bond angle deviates by only 1.5° from the observed one.The HF calculations yield an SiwO distance which deviates by less than 1 pm from the observed one, but the predicted angles are too large by ca. 3°. This is the same pattern as observed for the structures of the silica polymorphs predicted with the DFT- and HF-parametrized potentials. This points to an intrinsic eÜect of electron correlation. This is further supported by the MP2 results, which are virtually identical with the B3-LYP results.Fig. 2 Octahydridosilasesquioxane molecule, H8Si8O1250 Structure and reactivity of silica and zeolite catalysts Table 5 Structural parametersa of octahydridosilasesquioxane calculated by various methodsb,c method/basis set ref. r(SiO) n(SiOSi) n(OSiO) obsd. 32 161.7»162.0 147.5»147.6 109.4»109.7 MP2/T(O)DZP 34 164.5(2.8»2.5) 148.8(1.3»1.2) 109.3([0.1 to[0.4) DFT-B3LYP/T(O)DZP 164.8(3.1»2.8) 149.1(1.6»1.5) 119.1([0.3 to[0.6) HF/T(O)DZP 34 162.6(0.9»0.6) 150.5(3.0»2.9) 108.4([1.0 to[1.3) a Units: bond distances in pm, bond angles in degrees.symmetry assumed. c Deviation b Oh from experiment given in parentheses. Fig. 3 shows a simulated IR spectrum, based on DFT B3-LYP results, in comparison with the experimental one.33 Only the six modes belonging to the symmetry group T1u are IR active, however, the small distortion of the molecule from the sym- H8Si8O12 Oh metry gives rise to the additional bands present in the experimental spectra.The predicted intensity of the d(OwSiwH) bending mode is approximately two times lower than the corresponding intensity of the stretching mode, while the experimen- las(SiwOwSi) tal spectrum shows lower intensity for the latter.The intensity of the symmetric mode is also underestimated. Table 6 lists all vibrational frequencies of ls(SiwOwSi) calculated at the DFT B3-LYP level in the harmonic approximation and H8Si8O12 compares them with the experimental IR and Raman frequencies.33 The agreement between observed and predicted frequencies is excellent.The only exceptions are the SiwH stretching modes. The DFT-calculated frequencies are shifted towards higher values by 60»84 cm~1. This is connected with the high anharmonicity of the SiwH bond. We conclude that, despite the too long SiwO bonds, we may expect that the shell-model potential parametrized on DFT data will yield good vibrational spectra.Fig. 3 IR spectrum of octahydridosilasesquioxane. Experimental (upper)33 and calculated with the DFT-B3LYP method in harmonic approximation (lower). A gaussian band shape has been assumed.M. Sierka and J. Sauer 51 Table 6 Experimentala and calculated with the DFT B3-LYP methodb IR- and Raman-active fundamentals of octahydridosilasesquioxane normal mode/type calculated of vibration experimentala (DFT B3LYP) A1g c l(SiwH) 2302 2367 ds(OwSiwO) 580 586 ls(SiwOwS) 456 439 Eg c ds(OwSiwO) 932 916 ls(SiwOwSi) 697 668 ds(OwSiwO) 423 419 d(OwSiwO) 84 61 T2g c l(SiwH) 2296 2357 las(SiwOwSi) 1117 1136 ds(OwSiwO) 890 889 ls(SiwOwSi) 610 590 das(OwSiwO) 414 400 ds(OwSiwO) 171 166 T1u d l(SiwH) 2277 2361 las(SiwOwSi) 1141 1145 das(OwSiwO) 881 879 das(OwSiwO) 566 559 ls(SiwOwSi) 465 456 ds(OwSiwO) 399 397 a Ref. 33. b Units: cm~1. c Raman active. d IR active. 3.2.2 Disiloxane. The question that remains is why including electron correlation deteriorates the agreement with observed SiwO bond distances ? For an answer we consider disiloxane. This is the simplest molecule, containing the H3SiOSiH3 , SiwOwSi bridge (Fig. 4). The SiwOwSi bending potential is very —at, which has consequences both for theoretically predicted and observed structures.For the former, it is not easy to get results which are converged with respect to basis set and method. For the latter, the results obtained will depend on the experimental conditions and the average structures observed make it difficult to infer an equilibrium structure for this molecule.The fact that the angle in disiloxane can apparently assume diÜerent values without a signi–cant increase in the total energy explains the rich variety of the structures found in silicates, silica and zeolites. Fig. 4 Disiloxane molecule, H6Si2O52 Structure and reactivity of silica and zeolite catalysts Baé r and Sauer35 made a careful study of the eÜect that basis set extension has on the disiloxane structure predicted by the HF and MP2 methods.Table 7 summarises these results. We include also DFT B3-LYP results from the present study. For the basis set used in this study (abbreviated ì1dœ in Table 7) we have the already known pattern : the DFT and MP2 bond distances are virtually identical and larger than the HF result (by ca. 2 pm in this case). If we extend the sp basis set and include a 2d1f set of polarization functions, the HF, DFT and MP2 bond distances decrease. Further extension of the basis set makes the HF result slightly smaller, but also reduces the MP2 distance. The MP2 result for the largest basis set (162.5^0.3 pm) is almost identical with the HF result for the standard basis sets [T(O)DZP]. From the parallel change in MP2 and DFT with the basis set size we expect that the same may be true for the DFT method. From these results we learn : (i) the HF/T(O)DZP method yields very good SiwO bond distances owing to a cancellation of neglected correlation eÜects and truncated basis sets ; (ii) methods including electron correlation need very large basis sets for predicting accurate bond lengths.With respect to the structure predictions using DFTparametrized shell-model potentials we conclude that the systematically too long SiwO bond lengths are due to the basis set used for generating the data base and not to representation of the PES by a simple shell-model interatomic potential. Unfortunately, extending the basis set to the extent required to yield accurate SiwO bond lengths is far beyond computational resources. 3.3 Simulations of IR spectra of silica polymorphs A reliable parametrization of the potential should be capable of modelling the dynamical properties of microporous silica. Silicalite, the aluminium-free form of zeolite ZSM-5, and silica faujasite are polymorphs with a broad range of SiwOwSi angles and SiwO bond lengths.Fig. 5 and 6 show their IR spectra calculated with the present DFTderived shell-model potential as well as with the HF-derived one and compare them with observed spectra. Both potentials reproduce the main features of the measured spectrum well. Every band is reproduced. Larger deviations are only found in the range 700»900 cm~1. Table 7 Equilibrium structures and barriers to linearization of disiloxane calculated by various methodsa method/basis set ref.n(SiOSi) r(SiO) *Elin obsd. 36, 37 144^0.8 163.4^1.0 B1.3 MP2(5d4f*)b 35 158.2 162.3 0.12 MP2(5d4f) 35 150.7 162.8 0.88 MP2(2d1f) 35 145.4 163.7 1.86 MP2(2d) 35 139.8 164.6 4.61 MP2(1d) 35 152.9 164.9 0.65 DFT-B3LYP(2d1f) 155.9 163.7 0.52 DFT-B3LYP(2d) 143.8 164.6 3.04 DFT-B3LYP(1d) 170.6 164.5 0.00 HF(4d3f) 35 180.0 160.7 0.00 HF(2d1f) 35 172.0 161.1 0.00 HF(1d) 35 180.0 162.6 0.00 a Units: distances in pm, angles in degrees, energies in kJ mol~1. b Polarisation functions are speci–ed in parentheses, for detailed description see ref. 35.M. Sierka and J. Sauer 53 Fig. 5 Experimental IR spectrum of silica-rich faujasite (a) and IR spectra calculated with the DFT shell-model potential (b) and (c), and the ab initio HF shell-model potential (d) ; (c) and (d) have been generated from the line spectra assuming a gaussian band shape 3.4 Space group change and phase transition The most intriguing diÜerence between the DFT- and HF-parametrized potentials concerns the structures they predict for silica sodalite and silicalite.For silica sodalite, equilibrium structures with diÜerent space groups are predicted by the two potentials.The structure found using the HF potential has a —at four-membered ring, i.e. the four oxygen atoms of the four-membered face of the sodalite cage lie in the same plane, and belongs to the Im3m space group. According to Baur38 such a ìregularœ or ìexpandedœ structure is observed when the sodalite cage is –lled by a template molecule.On the other hand, the structure obtained using our DFT-parametrized potential belongs to the I43m space group and has a ì tilted œ or ìcollapsedœ structure, two of the oxygens belonging to the four-membered face of the sodalite cage are inside and two outside the unit cell. This is the structure observed for the silica sodalite with no template in the sodalite cage.When the Im3m space symmetry group was imposed during optimization using the DFT-based potential the resulting structure proved to be a high-order stationary point with eight negative frequencies in the !-point phonon spectrum. On the other hand, using the HF-derived potential the initially ì tilted œ structure always converges to the higher-symmetry ìexpandedœ one. Silicalite shows a reversible phase transition39 at ca. 340 K from monoclinic symmetry 11)27 to orthorhombic symmetry (Pnma).40 The low-temperature mono- (P21/n clinic structure has been successfully reproduced by the empirical shell-model54 Structure and reactivity of silica and zeolite catalysts Fig. 6 Experimental IR spectrum of silicalite (a) and IR spectra calculated with the DFT shellmodel potential (b) and (c), and the ab initio HF shell-model potential (d) ; (c) and (d) have been generated from the discrete spectra assuming a gaussian band shape potential.41 The HF-–tted potential yields a minimum on the PES only for the hightemperature orthorhombic structure.1 We performed lattice energy minimizations for silicalite using the present DFT-based potential.We obtained several local minima.The two with lowest energy, [12 494.6 and [12 494.1 kJ mol~1 per unit, belong to the SiO2 same space groups as the experimental structures, (monoclinic) and Pnma P21/n11 (orthorhombic), respectively. That these are minima was checked by !-point phonon calculations. To estimate nuclear motion corrections we performed phonon calculations using 27 k-points sampling of the –rst Brillouin zone for the orthorhombic and 30 kpoints for the monoclinic structure. After including zero-point vibrational corrections the energy of the monoclinic structure is lower by 45 kJ mol~1 uc~1 than that of the orthorhombic one. We calculated the free energy, *F, of the reaction We may assume that the phase transition occurs at *F\0, silicalitemono ]silicaliteortho .which we –nd is the case at 865 K. This predicted value is considerably higher than the experimental value of 340 K.27 Quartz at 853 K shows a reversible transition from the a to the b form.42 The PES of quartz has two equivalent minima, which correspond to two twinned con–gurations, and related by a rotation of the tetrahedra around their axes. The a1 a2 , SiO4 C2 corresponding tilt angle, d, is negative in the phase and positive in the phase.a1 a2 b-quartz has a tilt angle of zero, which corresponds to a transition structure on the PES. This means that at high enough temperatures ([853 K) the tilt angle samples positive and negative values and its average is zero (b-quartz). At temperatures below this pointM. Sierka and J. Sauer 55 the system is trapped into one of the minima with either a positive or a negative (a1) (a2) average tilt angle.We performed lattice-energy minimizations for both the a and b form of quartz using three diÜerent parametrizations of the shell-model potential. In all three cases, the a structure has lower energy than the b-structure, in accord with the doubleminimum shape of the PES.Indeed, the !-point phonon frequencies calculated for bquartz using the three parametrizations of the shell-model potential show one signi–cant negative value (transition structure). Based on phonon calculations we can estimate the free energy of the reaction a-quartz]b-quartz. Fig. 7 shows the dependence of the *F value on temperature. Assuming that the transition occurs for F(b)[F(a)\0 the empirical and our DFT-derived potential give a transition temperature of 1400 K.The HF-parametrized potential yields a much lower transition temperature of 540 K. As far as the comparison of the HF- and DFT-parametrized potentials is concerned, a uniform picture emerges for all the three systems studied : the HF gives undue favour to highly symmetric structures. For quartz, the energy diÜerence between b and a is too small ; for silicalite, the low-symmetry monoclinic structure is not a local minimum; for sodalite, an equilibrium structure is predicted that seems to be not even a local minimum on the ìtrueœ PES.It seems that electron correlation makes qualitative changes in the PES of silica polymorphs. It is encouraging that this feature is retained in the DFT-parametrized shell-model potential, even if there is room for quantitative improvement.This follows from the estimated temperatures for the a»b transition in quartz (1400 K) and for the monoclinic»orthorombic transition in silicalite (865 K) which are signi–cantly higher than the observed values (853 and 340 K, respectively). Our naive scheme to get a –rst approximation of the transition temperature may also contribute to this deviation, together with neglected anharmonicity eÜects.We would Fig. 7 Temperature dependence of the free energy of the a to b quartz transition calculated with various shell-model potential parametrizations56 Structure and reactivity of silica and zeolite catalysts like to stress that a qualitatively correct PES is obtained only after including electron correlation. 3.5 DFT shell-model potential results for protonated zeolites The modelling of active protonated forms of zeolites is especially attractive since many details of their structure are not known from experiments and reliable predictions of IR frequencies and NMR chemical shifts can settle assignment problems. In particular, the absolute acidity of the active sites is not known from experiments.Br‘nsted In the silicon-rich H-faujasite all the tetrahedral sites (T atoms) are crystallographically equivalent, but there are four diÜerent oxygen positions associated with each tetrahedron, Only two of the four possible OH groups, and have O1»O4. O1H O3H, actually been detected in the IR and 1H NMR spectra.Table 8 shows the relative energies and OH harmonic stretching frequencies calculated for isolated bridging hydroxy groups in the primitive cell of faujasite (one Al per 48 T sites). Results obtained using other parametrizations of the shell-model potential are also included for comparison. In agreement with observations, all three potentials predict that the and bridg- O1H O3H ing groups are energetically preferred compared to the remaining two.However, the present DFT potential yields the smallest energy gap between the and O1H O2H groups. This is in agreement with neutron diÜraction data which indicate a small amount of groups in protonated zeolite Y.43 Also in agreement with experimental O2H assignments, all three potentials predict the stretch vibration as the low-frequency O3H (LF) IR band and the stretch vibration as the origin of the high-frequency (HF) O1H band (Table 8).Harmonic frequencies calculated with the DFT B3-LYP derived potential underestimate the observed frequencies by ca. 65 cm~1. The experimental44 diÜerence between the HF and LF bands (fundamental, 73; harmonic, 80 cm~1) (O1H) (O3H) is reproduced by the DFT-parametrized potential (79 cm~1), while the other two potentials underestimate it slightly.The empirical and HF potentials predict 36 cm~1 and 22 cm~1, respectively. Table 9 shows the bond distances and angles predicted by diÜerent potentials for the and bridging hydroxy groups in faujasite. The SiwO(H) and AlwO(H) bond O1H O3H lengths calculated with the DFT potential are ca. 1 pm longer for the group than O1H for the group and are ca. 1 pm longer than in structures obtained using the HF O3H potential. From 1H MAS NMR sideband patterns AlwH distances of 248^4 and 240^4 pm for and groups, respectively, have been deduced.45 These O1H O3H compare well with the predictions of our shell-model potential (246 and 240 pm, respectively, Table 9). Lattice-energy minimizations employing the present potential were also performed for a model of H-faujasite with an Si/Al ratio of 2.43, typical for zeolite Y.A sample of Table 8 Relative energies and OH harmonic stretching frequencies of the bridging hydroxy groups in silicon-rich H-faujasitea *Erel lOH emp. shell HF shell DFT shell emp. shell HF shell DFT shell obsd.e obsd.f modelb modelc model model modeld model harm.anharm. O1H 5.3 8.8 11.1 3772 3752 3723 3787 3623 O2H 19.8 18.0 17.1 3702 3628 3602 » » O3H 0.0 0.0 0.0 3736 3694 3644 3707 3550 O4H 23.7 27.9 24.7 3751 3697 3673 » » a Units: energies in kJ mol~1, frequencies in cm~1. b Ref. 5. c Ref. 1. d Frequencies scaled by factor 0.9. e Ref. 44. f Ref. 48.M. Sierka and J. Sauer 57 Table 9 Structural parameters of the four bridging hydroxy groups in silicon-rich faujasitea r[AlwO(H)] r[AlwH] r[SiwO(H)] r(OH) n[SiwO(H)wAl] n(SiwOwH) SAlwOT SSiwOT O1H free cluster (DFT) 194.6 240.0 172.4 97.48 142.2 112.2 178.9 164.1 DFT shell model 189.8 246.0 170.7 97.05 129.8 116.0 176.2 163.8 M4T-1AlN embedded (DFT)b 190.4 250.0 172.0 97.54 123.3 118.5 176.9 164.6 free cluster (HF) 195.4 238.8 167.9 95.51 135.4 118.4 177.2 162.4 HF shell modeld 191.4 247.6 169.6 95.4 130.2 114.5 175.8 162.3 M4T-1AlN embedded (HF)b 190.4 248.2 170.1 95.55 125.6 117.8 175.6 162.8 empirical shell modelc 191.0 238.6 169.4 100.0 131.1 123.0 173.8 160.3 O3H DFT shell model 193.5 240.2 169.3 97.87 136.5 117.1 176.6 163.9 empirical shell modelc 193.0 233.2 169.7 100.2 138.7 120.7 174.2 160.7 HF shell modeld 194.7 239.7 169.1 96.1 140.2 113.8 176.4 162.6 a Units: distances in pm, angles in degrees.b Ref. 55. c Ref. 5. d Ref. 1.58 Structure and reactivity of silica and zeolite catalysts the same module but with a small amount of sodium cations was studied by Czjzek et al.43 in a neutron powder diÜraction experiment. The model was constructed as described by Schroé der and Sauer,1 namely the Al ions were distributed over the lattice sites according to a model suggested by Klinowski et al.,46 and the protons were distributed among the diÜerent sites as which is close to the O1 : O2 : O3 : O4\8 : 2 : 4 : 0, population deduced from the neutron diÜraction data.The space group was Fd36 m assumed. Table 10 shows the mean TwO distances. The average TwO bond length is calculated for each of the four oxygen types, regardless of whether the particular oxygen is protonated or not.This is not the same type of averaging as that involved in the re–nement of the neutron diÜraction data. Nevertheless, a comparison of the neutron diÜraction results with the predictions made by the diÜerent potentials is useful. All the potentials predict the same sequence of average TwO bond distances : O1[O3[ which is caused by the diÜerent occupancies of the oxygen sites.The oxygen O2[O4 , site with the highest occupancy shows the largest mean TwO bond length. This is reasonable, since the protonation of an oxygen lengthens its TwO bonds. The observed increase of the mean TwO bond lengths compared with silica faujasite, between 2.2 for and 7.0 pm for is correctly described by the shell-model potentials.As already O4 O1, discussed (Section 3.2), the present DFT-parametrized potential overestimates the average TwO bond lengths by 2»3 pm. The prediction of the average TwOwT bond angles for the H-Y zeolite model is much better using the present potential than using either of the other two listed in Table 10.This is probably due to a large data set of highly aluminated structures used for the potential –tting. Contrary to the other potentials, the present one predicts correctly the order of mean TwOwT angles : O2[O4[ The diÜerence between observed and calculated values of the average O3[O1. TwOwT angle is least for the present DFT-derived shell-model potential. 3.6 Embedded cluster calculations for the hydroxy group in H-faujasite O1H Table 9 shows the structural parameters of the bridging hydroxy site in faujasite O1H calculated with diÜerent methods: DFT and HF free clusters, DFT- and HF-derived shell-model potentials, and –nally the combined QM»Pot method presented in Section 2.5.The cluster used in QM»Pot calculations is a four-membered aluminosilicate ring saturated with H atoms, Al-cyclotetrasilicic acid, 4T-1Al (Fig. 8). For both the HF and the DFT methods similar trends are observed on embedding the cluster in its crystal environment. The AlwO(H) distance is shortened and the AlwH non-bonding distance Table 10 Calculated and observed mean bond lengths and angles and unit cell parameter of zeolite H-Y (Si : Al\2.43)a empirical shell HF shell DFT shell modelb modelc model obsd.d Sr(O1wT)T 168.9(7.5)e 169.4(7.2) 170.5(6.8) 167.7(7.0) Sr(O2wT)T 163.0(3.1) 165.0(4.1) 166.4(3.6) 163.2(3.5) Sr(O3wT)T 165.6(4.7) 166.7(5.0) 168.3(4.6) 165.4(5.0) Sr(O4wT)T 162.3(1.5) 164.5(2.9) 165.6(2.5) 163.6(2.2) n(TwO1wT) 136([2) 138([3) 135([3) 136([2) n(TwO2wT) 152(3) 156(2) 149(2) 145([4) n(TwO3wT) 146(0) 150([2) 142([2) 140([6) n(TwO4wT) 143(2) 144([1) 144(1) 144(3) a 2484(61) 2522(59) 2520(54) 2477(51) a Units: distances in pm, angles in degrees.b Ref. 5. c Ref. 1. d Ref. 43. e Deviation from the values for pure-silica faujasite given in parentheses.M. Sierka and J. Sauer 59 Fig. 8 The 4T-1Al cluster model embedded in the faujasite lattice. The oxygen labels correspond to the diÜerent crystallographic positions.becomes longer, –nally giving similar structural parameters: 190 pm for AlwO(H) and 249^1 pm for AlwH. Similarly, for the SiwO(H) distance the diÜerence between the DFT and HF free cluster results of 4.5 pm is reduced to only 1.9 pm in the embedded clusters. The OwH distance changes only slightly when the embedding scheme is applied. The AlwO(H)wSi angle becomes smaller by 10»20° on embedding, probably owing to removal of symmetry constraints and eliminated repulsion between terminal link atoms.In Table 11 we list the deprotonation energies, OwH stretching frequencies and 1H NMR chemical shifts of the bridging hydroxy group in an H-faujasite O1H (Si/Al\47.0, 1 Al atom in the primitive cell) calculated using the combined QM»Pot Table 11 QM»Pot calculations : Propertiesa of the O1H bridging hydroxy group in high-silica H-faujasiteb DFT HF *EDP free cluster » 1212c QM//QM»Pot 1319 1348 LR//QM»Pot [121 [98 QM»Pot//QM»Pot 1198 1250 method correction 0^3 [46^2 1198 1204 vibrational correction [35 –nal value 1166^10 lOH f potential only 3723 4216d QM»Pot//QM»Pot 3646 4016e QM//QM»Pot 3729 4041 1dH 4.2 4.0 a Units: energies in kJ mol~1, frequencies in cm~1, NMR chemical shifts in ppm.b Ref. 55. c Ref. 52. d Ref. 1. e Ref. 19. f Harmonic frequencies.60 Structure and reactivity of silica and zeolite catalysts approach. The deprotonation energies are decomposed into the contribution from the quantum chemical cluster calculation (QM//QM»Pot) and a long-range electrostatic contribution from the shell-model potential (LR//QM»Pot) which includes corrections discussed in Section 2.6. Both are calculated for the structures of the parent and deprotonated systems at their respective equilibrium structures obtained with the combined QM»Pot approach.The DFT method yields lower values than the HF methods. Obviously, inclusion of electron correlation weakens the OwH bond. Sauer has suggested53 removal of the main body of the systematic error from calculated deprotonation energies by adding a constant value for a class of systems which is typical for a given method.This constant value is obtained as follows : the approximation/basis set is used to calculate deprotonation energies for related small molecules, methanol and silanol in this case, and comparison is made with accurate results, in this case obtained by Sauer and Ahlrichs.54 For the HF method with the T(O)DZP basis sets a correction of [46 kJ mol~1 was obtained.53 The DFT B3-LYP method, in combination with this basis set, yields deprotonation energies that are virtually identical with the accurate data.Hence, by fortunate error compensation, the correction is zero in this case.After applying these corrections the results of the two methods agree within 6 kJ mol~1. For comparison with experiment, we have to add nuclear motion corrections which have been estimated previously as [35 kJ mol~1.54 This yields a –nal estimate of ca. 1166^10 kJ mol~1 for both methods, which agrees well with the ìexperimentalœ value of ca. 1190 kJ mol~1 for H-Y and 1200 kJ mol~1 for HNaY.47 The combined QM»Pot approach yields absolute acidities.In another study we use this approach to simulate the dependence of the acidity (deprotonation energy) on the Si/Al ratio55 and the framework type.56 It is well known that the HF method yields stretching frequencies that are systematically too large, while the DFT method nearly reproduces observed values.In both cases, however, comparison is made between values calculated in harmonic approximation and observed frequencies including anharmonicities. Our results show the same pattern, the OH-stretching frequency from the embedded DFT calculations (QM»Pot//QM» Pot), 3646 cm~1, agrees well with the experimental value48 of 3623 cm~1, while the value obtained at the HF level, 4016 cm~1 is too large.When applying, however, the usual scaling factor of 0.9 for the HF frequencies, one obtains a reasonable result of 3614 cm~1. The frequency predicted by the DFT parametrized potential alone is almost equal to the value calculated for the free cluster at the structure found by the embedding calculation. As for the deprotonation energies, the long-range contribution is larger for the DFT method than for the HF method.The 1H NMR chemical shift for the bridging hydroxy group calculated at the DFT level, 4.2 ppm, agrees better with the experimental49 value of 4.4 ppm than the result of the HF calculations, 4.0 ppm. Note, however, that the prediction of converged values of NMR shifts requires larger cluster models for the calculation of the wavefunction. 4 Conclusions We have shown that the ion-pair shell-model potential with parameter sets derived from the data base obtained by quantum mechanical calculations yields reliable predictions for structures and properties of silica and zeolite catalysts. Our parametrization based on DFT calculations yields results which show typical eÜects of electron correlation : longer SiwO bonds and preference for less symmetric structures.The combined quantum mechanics»interatomic potential approach proved to be a powerful tool for simulations of catalytically active sites in periodic framework structures. We thank K-P. Schroé der and M. Braé ndle for providing unpublished results and useful discussion. We are also grateful to R. Ahlrichs (Karlsruhe) and J. Gale (London) forM.Sierka and J. Sauer 61 making available recent versions of the TURBODFT and GULP codes, respectively. This work has been supported by the ìFonds der Chemischen Industrieœ. References 1 K-P. Schroé der and J. Sauer, J. Phys. Chem., 1996, 100, 11043. 2 J-R. Hill and J. Sauer, J. Phys. Chem., 1995, 99, 9536. 3 J-R. Hill and J. Sauer, J. Phys. Chem., 1994, 98, 1238. 4 J. Purton, R.Jones, C. R. A. Catlow and M. Leslie, Phys. Chem. Miner., 1993, 19, 392. 5 R. A. Jackson and C. R. A. Catlow, Mol. Simul., 1988, 1, 207. 6 F. Haase and J. Sauer, J. Am. Chem. Soc., 1995, 117, 3780. 7 F. Haase and R. Ahlrichs, J. Comput. Chem., 1993, 14, 907. 8 D. A. Becke, J. Chem. Phys., 1988, 88, 1053. 9 TURBODFT: O. Treutler and R. Ahlrichs, J. Chem. Phys., 1995, 102, 346. 10 A.Schaé fer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571. 11 M. J. Frisch, G. W. Trucks, H. B. Schlegel, P. M. W. Gill, B. G. Johnson, M. A. Robb, J. R. Cheeseman, T. Keith, G. A. Petersson, J. A. Montgomery, K. Raghavachari, M. A. Al-Laham, V. G. Zakrzewski, J. V. Ortiz, J. B. Foresman, 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, S. Binkley, D. J. Defrees, J. Baker, J. P. Stewart, M. Head-Gordon, C. Gonzalez and J. A. Pople, GAUSSIAN 94, Revision B.3, GAUSSIAN, Inc., Pittsburgh PA, 1995. 12 R. Ditch–eld, Mol. Phys., 1974, 27, 789. 13 M. Haé ser, R. Ahlrichs, H. P. Baron, P. Weis and H. Horn, T heoret. Chim. Acta, 1992, 83, 455. 14 TURBOMOLE: R. Ahlrichs, M. Baé r, M.Haé ser, H. Horn and C. Koé lmel, Chem. Phys. L ett., 1989, 162, 165. TURBOMOLE and TURBODFT are available commercially from MSI: San Diego, CA. 15 F. Haase and J. Sauer, J. Phys. Chem., 1994, 98, 3083. 16 GULP written and developed by J. D. Gale, Royal Institute and Imperial College, UK, 1992»1994. 17 METAPOCS, C. R. A. Catlow, A. N. Cormack and F. Theobald, Acta Crystallogr.B, 1984, 40, 195. 18 U. Eichler, C. K. Koé lmel and J. Sauer, J. Comput. Chem., 1996, 18, 463. 19 U. Eichler, J. Sauer and M. Braé ndle, J. Phys. Chem., submitted. 20 M. Braé ndle and J. Sauer, J. Mol. Catal., 1997, 119, 19. 21 M. Leslie and M. J. Gillian, J. Phys. C, 1985, 18, 973. 22 K. Fuchs, Proc. R. Soc. L ondon Series A, 1935, 151, 585. 23 J. A. Hriljac, M. M. Eddy, A. K. Cheetham, J.A. Donohue and G. J. Ray, J. Solid State Chem., 1993, 106, 66. 24 J. W. Richardson Jr., J. J. Pluth, J. V. Smith, W. J. Dytrych and D. M. Bibby, J. Phys. Chem., 1988, 92, 243. 25 B. Marler, Zeolites, 1987, 7, 393. 26 P. R. Rudolf and J. M. Garceç s, Zeolites, 1994, 14, 137. 27 H. van Koningsveld, J. C. Jansen and H. van Bekkum, Zeolites, 1990, 10, 235. 28 R. X. Fisher, W. H. Baur, R. D. Shannon, R. H. Staley, A. J. Vega, L. Abrams and E. Prince, J. Phys. Chem., 1986, 90, 4414. 29 L. Levien, C. T. Previtt and D. J. Weidner, Am. Miner., 1980, 65, 920. 30 H. Grimm and B. Dorner, J. Phys. Chem. Solids, 1975, 36, 407. 31 D. R. Peacor, Z. Kristallogr., 1973, 138, 274. 32 T. P. E. Auf der Heyde, H. B. Bué rgi, H. Bué rgy and K. W. Toé rnroos, Chimia, 1991, 45, 38. 33 M. Baé rtsch, P. Bornhauser, G. Calzaferri and R. Imhof, J. Phys. Chem., 1994, 98, 2817. 34 J-R. Hill and J. Sauer, unpublished results. 35 M. R. Baé r and J. Sauer, Chem. Phys. L ett., 1994, 226, 405. 36 A. Almenningen, O. Bastiansen, V. Ewing, K. Hedberg and M. Acta Chem. Scand., 1963, 17, Trêtteberg, 2455. 37 M. J. Barrow, E. A. V. Ebsworth and M. M. Hardig, Acta Crystallogr. B, 1979, 35, 2093. 38 W. H. Baur, J. Solid State Chem., 1992, 97, 243. 39 C. A. Fyfe, H. Strobl, G. T. Kokotailo, G. J. Kennedy and G. E. Barlow, J. Am. Chem. Soc., 1988, 110, 3373. 40 D. H. Olsen, G. T. Kokotailo, S. L. Lawton and W. M. Meier, J. Phys. Chem., 1981, 85, 2238. 41 R. G. Bell, R. A. Jackson and C. R. A. Catlow, J. Chem. Soc., Chem. Commun., 1990, 10, 782. 42 J. D. Axe and G. Shirane, Phys. Rev., 1970, 1, 342. 43 M. Czjzek, H. Jobic, A. N. Fitch and T. Vogt, J. Phys. Chem., 1992, 96, 1535. 44 K. Beck, H. Pfeifer and B. Staudte, Microporous Mater., 1993, 2, 1. 45 D. Fenzke, M. Hunger and H. Pfeifer, J. Magn. Reson., 1991, 95, 477. 46 J. Klinowski, S. Ramdas, J. M. Thomas, C. A. Fyfe and J. S. Hartman, J. Chem. Soc., Faraday T rans. 2, 1982, 78, 1025.62 Structure and reactivity of silica and zeolite catalysts 47 V. M. Mastikhin, I. L. Mudrakovsky and A. V. Nosov, Bruker Rep., 1989, 2, 18. 48 M. W. Anderson and J. Klinowski, Zeolites, 1986, 6, 455. 49 D. Freude, M. Hunger, H. Pfeifer, Z. Phys. Chem. Neue Folge, 1987, 152, 171. 50 J. Sauer and R. Ahlrichs, J. Chem. Phys., 1990, 93, 2575. 51 B. Winkler, M. T. Dove and M. Leslie, Am. Miner., 1991, 76, 313. 52 M. Braé ndle, unpublished results. 53 J. Sauer, in Modeling of Structure and Reactivity in Zeolites, ed. C. R. A. Catlow, Academic Press, London, 1992, pp. 183»216. 54 J. Sauer and R. Ahlrichs, J. Chem. Phys., 1990, 93, 2575. 55 M. Sierka, J. Datka and J. Sauer, in preparation. 56 M. Braé ndle and J. Sauer, J. Am. Chem. Soc., submitted. Paper 7/01492I; Received 3rd March, 1997