Chapter 10 Introduction to zeolite theory and modelling

Studies in Surface Science and Catalysis 137 H. van Bekkum, E.M. Flanigen, P.A. Jacobs and J.C. Jansen (Editors) 9 2001 Elsevier Science B.V. All rights reserved.

419

Chapter 10 Introduction to zeolite theory and modelling R.A. van Santen, ~ B. van de Graaf, b and B. Smit c ~Schuit Institute of Catalysis, Eindhoven, University of Technology Eindhoven, P.O. Box 513, 5600 MB Eindhoven, The Netherlands bLaboratory for Organic Chemistry and Catalysis, Delft University of Technology Julianalaan 136, 2628 BL Delft, The Netherlands CDepartment of Chemical Engineering, Universiteit van Amsterdam Nieuwe Achtergracht 166, 1018 WV Amsterdam, The Netherlands In this Chapter some of the recent advances in zeolite theory and modelling are present. In particular we discuss the current status of computational chemistry in Br0nsted acid zeolite catalysis, Molecular Dynamics simulations of molecules adsorbed in zeolites, and novel Monte Carlo technique to simulate the adsorption of long chain hydrocarbons in zeolites. 1. I n t r o d u c t i o n In the previous edition of these book various theoretical and computational approaches to study zeolites have been discussed [1]. The focus was on methods to compute zeolite structure and stability, computer graphics to visualise the complex zeolite structures, and Molecular Dynamics and Monte Carlo calculations to simulate the adsorption and diffusion of small molecules in these materials. Although most of these topics are still of interest, many new developments have taken place since the appearance of the first edition. Novel algorithms, such as ab-initio molecular dynamics, embedding, or novel Monte Carlo techniques have been applied to zeolites and allow us to use computational for many more problems in zeolite science. Furthermore, the enormous increase in computer power allows us to study systems that are more realistic and closer to experiments. A disadvantage of this development is that it is not possible to cover the entire field in a single chapter. The limited space forces us to focus on two aspects zeolite acidity and the adsorption and diffusion of molecules in the pores. The first topic deals with the chemistry of the protons in zeolites and in this part recent developments in quantum chemistry calculations are discussed. Since molecules can not be converted if they do not adsorb or diffuse in the pores of the zeolites, diffusion and adsorption are important aspects contributing to the activity of a zeolite. These properties have been studied using Molecular Dynamics or the Monte Carlo technique and in the second part of this chapter recent developments in this area are discussed.

420 2. Q u a n t u m - c h e m i s t r y

of zeolite a c i d i t y

2.1. I n t r o d u c t i o n Current electronic structure calculations have reached a stage that they can be usefully applied to analyze the energetics of reacting systems representing the catalytically reacting site [2,3]. Here we will discuss the current status of this approach with respect to the protonation reaction in zeolites. Using transition state reaction rate theory [4] rate constants for the elementary reaction steps can be computed. Quantum-chemical calculations enable not only the prediction of interaction energies in the ground state, but also the energies that correspond to the transition state. If in addition a normal mode analysis is done in reactant as well as transition state the activation entropy can be computed. Within the harmonic approximation the corresponding expression for the activation entropy is: /kS* -- kB In I-IN-1 [1 -- e x p ( - h ~ ' j / k B T ) ] -1 1-Iff [1 - e x p ( - h u ~ / k B T ) ] -1 "

(1)

With uj and uk the normal mode frequencies of the complex in the transition state and ground state respectively, kB is Boltzmanns constant and T the temperature. The expression for the activation energy becomes: 1

Eact - /kE* + ~h

pj -

p~

(2)

AE* is the computed electronic barrier energy. The correction to AE* is due to the zero point frequencies of the nuclei. Within the harmonic approximation the normal modes can currently be routinely computed by most of the available ab-initio quantum-chemical routines. To compute such properties is also very important in case one wishes to validify computational predictions. Measurements of rate constants of elementary reaction steps would be very useful, but are rarely available for zeolites. Spectroscopic data are widely available of adsorbates interacting with protons. This holds especially for Infrared and NMR data. It has been shown for instance that proton NMR chemical shifts can be used to probe the state of adsorbed CHaOH [5]. Also Nuclear Quadruple Coupling constants have been computed [6] to probe the local coordination of H20. Of course the results depend sensitively on the accuracy or level of the calculations. For the computation of infrared spectra an additional complexity arises from the intrinsic anharmonicity of the OH modes. This is typically of the order of a few 100 cm -1 In case one wishes to use the infrared spectrum as a probe tbr the adsorption state, this shift cannot be ignored. Methods [2,7] have been developed to not only solve the anharmonicity problem in one dimension, but in several dimensions. In essence one has to compute the potential energy surface with a certain dimensionality and than one solves the equations of motion of the atoms in this potential [7]. Such approaches have also been used to develop force fields or potentials so as to predict the structure and physical properties of zeolites. This has been reviewed elsewhere [8] and will not be the subject of this paper. However to study chemical reactivity one has to exclusively use quantum-chemical approaches. An important issue is the use of clusters to model the local environment around the zeolitic proton. Until recently no quantum-chemical computer routines were available

421 to solve the protonation problem in the three dimensional solid. The main technical limitation was the need to optimize the geometry of the lattice to identify ground states and transition states with the fully relaxed structures. Molecular mechanics studies based on classical potentials have provided ample proof of the flexibility of the zeolite lattice and local relaxation upon protonation [9]. Quantum-chemical cluster calculations have only given agreement with experiment when the geometry had been allowed to relax. In fact the now popular Density Functional Theory methods became only useful, once programs as DGAUSS, DMOL, DMON or AOF were provided with routines that enables geometry optimization. This situation has only quite recently changed for the solid state problem. The breakthrough has been provided for by the implementation of the Car-Parrinello method [10-13]. It is not our intention to provide an introduction to quantum-chemical methods, but here we will shortly explain the essential electronic structural principles that are basic to each technique. This will later help to appreciate and to compare the different techniques. 2.2. C o m p u t a t i o n a l m e t h o d s To study chemical reactivity quantitatively one has to use ab-initio quantum-chemical techniques. Nonetheless the use of semi-empirical methods can be quite useful. For instance to obtain a starting structure of a transition state to be completed by a first principle computational method, transition states obtained by the much less computational demanding semi-empirical methods can be quite useful. A technique that has been used in our work is MNDO [14,15]. In the transition state a computation of the normal mode frequencies should give one imaginary mode. Most ab-initio programs have now the option to not only find the stationary point and local energy minima, but also to compute the saddle points that correspond to the transition states. The electronic structure codes can be defined in two categories. One class is based on the use of the Hartree-Fock method [16] and subsequent improvements. Gaussian [17], Gamess [18] and Crystal [19,20] belong to this class, Gaussion and Gamess being codes suitable for the study of clusters. Crystal can be used to compute properties of the three dimensional solid. Important issues here are choice of basis set and with which correction the Hartree-Fock method has to be used. Crystal does not yet have the option for geometry optimization. Recently Gaussian as well as Crystal have obtained options to include Density Functional Theory options. The essence of the Hartree-Fock method is that one-electron wave functions, corresponding to molecular orbitals, are computed using a potential where the electron-electron interaction is deduced from the electrostatic interaction with other electrons distributed in the other one-electron wave functions. It implies that correlation between the motion of the individual electron is ignored. It appears that to compute hydrogen bonds and protonation this is not an acceptable approximation and correlation corrections are important. A reasonable approach appears the use of the M011er-Plesset 2 (MP2) correction [21,22], implemented in most commercial packages. This is an important issue, because the correlation correction to the Hartree-Fock energy also contains the van der Waals interactions. When a molecule adsorbs into a zeolite it not only experiences hydrogen bonding with

422 the protons but also an important interaction with the highly polarizable oxygen atoms in the wall of the zeolite. Per oxygen atom this interaction is usually quite small, typically of the order of 10 k J/mole. Therefore, it is called physical adsorption. However, a larger molecule has contact with many channel oxygen atoms, hence the total interaction can be quite large. It gives arise to the substantial adsorption energy of molecules adsorbed in zeolites. This physical adsorption is mainly due to the dispersive van der Waals interactions and therefore very important to compute. Configuration interaction methods as the MOllerPlesset method are quite important for this purpose. Many technical aspects determine the quality of a calculation. Most important are choice of basis set and geometry optimization. Often geometrys are optimized in the Hartree-Fock approximation and energies are improved by using better approximations for chosen geometry's. For the subtle problem of protonation of polar molecules, this leads to erroneous results as we will discuss in the next section. In the Density Functional Theory methods [11,12,23-31], to be discussed now, the van der Waals problem has in essence not been solved. Density Functional Theory has initially been invented to provide a shortcut to the computational demanding Hartree-Fock theory. In Local Density Theory, the initial version of Density Functional Theory, the electron exchange term in Hartree-Fock theory is replaced by a function that depends on the electron-density similarly as a free-electron gas. Surprisingly to compute structures this appears to be a very useful approach. Not however to compute bond energies. They appear often to be overestimated by a factor of two. A dramatic improvement is found when the local density approximation is corrected for by so called density gradient correction terms to correct for electron-correlations as well as electron density dependence. Different versions of these gradient correction terms are available all having a particular advantage, dependent on systems to be studied. Cluster versions as well as solid-state versions of the DFT programs are available. For the zeolite protonation DFT methods have demonstrated themselves to be very useflll and to produce ground state interaction energies often correct with an accuracy of ~5 kJ/mole. Transition state energies may be underestimated by ~30 kJ/mole. To compute infrared frequencies the situation is not that satisfactory. To produce accurate values application of an adequate electron correlation density function is critical. The B3LYP functional appears to be best, but anharmonicity corrections have to be included. It appears that weakening of vibrational bonds by hydrogen bonding is significantly overestimated [32]. As mentioned earlier vail der Waals interactions cannot be reliably computed, they appear to be underestimated (~10 kJ/mole). Whereas geometry optimization methods have been widely implemented for clusters, this is not the case for calculations of the three dimensional solid. Car-Pai'rinello methods [10-13] intrinsically provide a solution for this, because their purpose is to study the molecular motion in the field of the computed electronic structure. Again it is essential to use a version that contains the gradient-corrected functional for the electronic correlation energy. Other issues affecting accuracy relate to the k-cut-off in the number of electronic plane-wave functions to be used in the Car-Parrinello method and the choice of pseudo potentials. So far in zeolites no transition states have been studied using Car-Parrinello techniques. They have great promise because application of

423 Car-Parrinello method do not require the cluster approximation and enable determination of optimum geometries. 2.3. T h e c l u s t e r a p p r o x i m a t i o n a n d e m b e d d i n g m e t h o d s To judge which aspects of the zeolite catalysis problem can be treated using clusters to model the Br0nsted acidic site, one has to distinguish which kinetic parameters depend on structural differences between zeolites from the question whether the chemical bonding aspects of the interaction with a zeolitic proton can be usefully studied with the cluster approximation [33]. Properties that depend strongly on zeolite structure are heats of adsorption and diffusion rates. Clearly clusters cannot be used to study these aspects of the kinetics problem. Also zeolites may control the overall kinetics of a reaction by preventing access of molecules that are too large to enter a micropore or by preventing the formation of product molecules that are too large for the zeolite micropore. Product formation may be prevented if transition states are required that are too large for the micropore. These are extremely important effects. The use of zeolite Y as a cracking catalyst is based on its reduced coke formation. This is prevented because condensation of aromatic rings is suppressed by micropore size. In the case of methanol conversion catalysis by the medium-pore size ZSM-5 leads to stable operation and formation of benzene and methylbenzenes. On wider pore zeolites methanol will also be converted to higher hydrocarbons, but zeolite operation is not stable due to formation of deactivating carboneous materials. The narrow pore chabasite type zeolite SAPO-34 [34,35] will convert methanol selectively to ethylene. So in considerations of the overall kinetics of zeolites one has to explicitly treat the structure dependent aspects, that are mainly determined by van der Waals interactions [36]. Clearly the prediction of the size of the transition state of a particular reaction implies that an estimate can be made of the minimal pore-size requirement for the reaction to proceed. Important to the choice of a cluster is how well physical properties are properly computed. Initially it has been proposed that a charged cluster should be cut from the zeolite and be embedded in an electrostatic lattice, so as to compute the long range Madelung effects correctly. It is now clear that in such an approach chemical bonding features become completely dominated by edge effects and lead to physically unrealistic results. There is now a general agreement that clusters should be chosen neutrally by terminating them with hydroxyl or hydride bonds. Studies of the deprotonation energy of such clusters show significant fluctuations in computed values [37]. Hence it is not obvious which cluster choice is best. From NMR [38] and IR [39] studies of the zeolitic proton it is well known that its chemical bond is affected by compositional changes in structure at positions that are next nearest neighbor of the tetrahedra lattice bridged by protonated oxygen atoms, that forms the Br0nsted acidic site. This implies the use of clusters at least of a size of 26 tetrahedra. These effects are not electrostatic, but covalent and can be well understood using the Bond Order Conservation principle [40]. It is for instance reflected in the dependence of the intrinsic acidity of protons on the Al-concentration on zeolite Y. As follows from calorimetric ammonia adsorption measurements [41,42] the ammonia protonation energy drops as soon as a fraction of the protons reacted with am-

424 monia. This is very similar to the allosteric effect that is found in enzymes with several sites of interaction. Clearly clusters containing only three to five tetrahedra, as often used to study transition states cannot describe these composition dependent effects. However termination of the cluster can be altered so as to influence the intrinsic acidity. In small clusters, it can be mimicked by assigning different bond lengths to terminal SiO or SiOH bonds and optimizing all other parameters [43]. For instance lengthening of terminating Sill bonds shortens neighboring Sill bonds and as a consequence weakens a neighboring proton-oxygen bond. Also the same reaction can be studied as a function of cluster size. A plot can be made of the activation energy of that reaction against deprotonation energy of the clusters [44,45] (see Figure 1). Then one can extrapolate the curve to that deprotonation- energy for which one desires to compute the activation energy that for instance can be estimated for a proton based on its infrared shift induced by a weakly adsorbing molecule [46]. One then uses the Bronsted-Polanyi principle that proposes the relation: AEac t ~- CAEreaction.

(3)

The change in activation energy is proportional to the change in reaction energy as long as there is no change in reaction mechanism. We have found [33] that once neutral clusters are large with 10 or 12 tetrahedral sites, electrostatic embedding affects the interaction energy of ammonia with a zeolitic proton with less than 10%. Many physical properties of zeolitic protons are described with high accuracy by very small clusters.

2.4. General concepts and critical evaluation of c o m p u t a t i o n a l results Calculations on the Br0nsted acidic site in a zeolite have been done with many different methods. A representative (:luster is shown in Figure 2 [6]. This cluster has five T-sites and the Br0nsted acidic site is the protonated bridging atom between Si and A1. The predicted geometry is from a DFT calculation. The calculation is obtained from a non geometry constrained relaxation. The result is very similar to that obtained with a high quality H F / M P 2 calculation. The H-A1 distance can be compared with NMR estimates and the predicted OH stretch flequency, when corrected for anharmonicity will be within 50 cm -1 of that of tile typical low alumina spectrum. Proper A1OH distances are obtained only when clusters have at least three T-sites. The deprotonation energy of such small clusters is usually a few hundred k J / g a t larger than experimental values. The H-chemical shift as well as A1-NMR Quadrapole Coupling constants can be compared well with experimental values. One deduces a charge less than 0.1 e.u. for the proton, consistent with strong covalent bonding of tile OH bond. Several analyses have been made of the dependence of the deprotonation energy on the location of the Si-OHA1 unit in the zeolite. One expects each non equivalent O positron to have a different deprotonation energy. The approaches to be taken can be quantum chemical or classical. In quantum-chemical calculations one can constrain the positions of the atoms that terminate the cluster. These terminating positions can be chosen according to crystallographically determined positions or positions obtained from structural simulations of the siliceous or aluminum phosphate topologies. For the latter accurate procedures are

425 H

H

I

I

"/O~AI/O'~" O Ot t H H AHdepr = 317.5 kcal/mol

H3Sr~ O~Al"O~'si"3 /\ H H

Extrapolation for "average" zeolite

AHdepr = 302.3 kcal/mol

AHdepr = 295.4 kcal/moi

E act 1. dehydrogenation of ethane

\Z 3\ \ \

~

-7.t

,o . . . . . .

317.5

, ...........

302.3

-7.5 9

295.4

x/

AHdepr

Figure 1. Br0nsted-Polanyi dependence of the activation energy as a function of cluster size [44].

Figure 2. A cluster representing the Br0nsted acidic site. Distances are in ppm [6].

426 available. Such cluster choices imply choices for particular sites. The implicit assumption underlying the constrained cluster model is that the changes in site structure due to the replacement of a Si-O-Si or A1-O-P unit by a Si-OH-A1 or equivalent unit are local. This has indeed been found in molecular mechanics simulations of the extended lattice. These simulations employ classical potentials deduced from quantum-chemical cluster calculations. Due to the very shallow Si-O-Si bonding potential the increase in local volume that occurs upon substitution is accommodated by small changes of the SiOSi angles. A few T-sites away from the acidic sites the structure is the same as before substitution. In such calculations a spread in deprotonation energies of 80 k J/mole is found for low A1 concentration systems [47,48]. In addition to molecular mechanics simulations of extended systems also recently quantum-mechanical Car-Parrinello calculations have been done for Offretite [49]. With in essence the same result. The largest difference in deprotonation energies oil different Offretite oxygen atoms is now 40 kJ/mole. These differences are in quite good agreement with mobility NMR measurements [50]. The conclusion of both approaches is that the main parameter that determines the difference in deprotonation energy is not the long range electrostatic interaction, as might have been expected, but differences in local geometry. In structures with short Si-O-Si distances the constraint oi1 the A1-OH-Si unit is largest resulting in slightly different A1-O-Si angles. This is the main reason for the lower acidity at sites that are most expandable. One can also use a combination of a quantum-chemical calculation and classical structure prediction, by embedding the cluster in a three dimensional lattice connected to the quantum-chemical part of the system by classical forces. In this way a smaller cluster can be used and its spatial relaxation is determined fi'om the classical part of the calculation [48]. Interaction with a weakly adsorbing molecule as CO or NO shifts the OH stretch frequency downwards with ~ll00cin-~ A strongly interacting molecule as acetonitrile shifts the OH frequency downwards by ~100cm -1 Calculations using different calculational approaches on clusters reproduce this bondweakening of the OH bond, but with different accuracy [51]. DFT seems to overestimate the bondweakening. It is important that anharmonic corrections to the frequencies are flllly included. HF overestimates frequencies by 10~). HF-MP2 calculations are to be preferred. HF-MP2 and DFT calculations agree that acetonitrile is not protonated. An important issue concerns the interpretation of the broadened downwards shifted intensity due to the weakened OH stretching frequency. Especially when strong bases as acetonitrile, H2() or methanol interact. The high intensity is due to the enhancement of the transition dipole moment. This is due to the polarization of the OH band in contact with the doubly occupied lone pair electrons of the adsorbing base. The resulting Pauli-repulsion is decreased by polarization of electrons away from the interaction region. In case of the interaction with acetonitrile the charge on the proton has increased to 0.3 e.u. [52]. Cluster calculations reproduce this quite well. One finds that the acidic nature of the zeolitic Br0nsted OH bond relates to its easiness of polarization. The OH bridging a three valent and four valent lattice cation appears highly polarizable. The broadened infrared spectrum of the perturbed OH bond in contact with a strong base is found to be split in two or three peaks. Many theoretical studies have been devoted to the interpretation of this band structure. For acetonitrile or acetone it has been proposed that the peaks describe a hydrogen

427 bonded in equilibrium with a protonated species. For methanol and H20 the bands have also been assigned to the asymmetric and symmetric OH stretching modes of the protonated methoxonium or hydronium system [53]. Now according to the most accepted interpretation it is an A, B, C band system due to Fermi resonance of hydrogen bonded base, with overtones of in-plane and out-plane OH bending modes [3]. Cluster calculations have been very useful to analyze this. They also have shown the limitations of different calculational approaches. It appears that D F T theory including gradient exchange-correlation correction terms, properly optimized, generally predicts that the hydrogen bonded species is the ground state. It is essential to analyze explicitly whether a stationary point is a local minimum or a transition state. HF-MP2 theory is essential, with a good quality basis set to describe the interaction with the proton, if one wishes to use a HF based method. The structure has to be optimized at the H F / M P 2 level. To compute frequencies full anharmonicity has to be included. For interacting acetonitrile the A, B, C band profile has been theoretically reproduced. Coupling between the OH modes and acetonitrile has been computed by solving the dynamics of the proton in the potential deduced from quantum-mechanical cluster calculations [7]. With increasing acidity of the cluster the intensity maximum shifts to lower frequencies. Also frequency shifts of molecular modes of the base molecule can be used to analyze acidity [54]. Here a scaled HF approach suffices. These are sensitive to Lewis or Bronsted acidic interactions. To summarize this discussion it may be useful to emphasize that most critical to the computation of proper OH frequency shifts is the right prediction of the structure of the interacting system. Small bond distance differences can have a large effect on the final results. As we will see it will also be strongly influenced by cluster choice. Because of the large heterolytic bond dissociation energy of the zeolitic OH bond (~1250 kJ/mol), protonation of a molecule, even a very basic one can only occur because of the presence of an additional stabilizing term: [ZH-B]I -+ [Z-. HB+]II .

(4)

Structure I represents the hydrogen bonded situation, structure II is the case where the basic molecule has become protonated. The energy change between structure I and II can be considered to consist of three terms: The deprotonation energy of the zeolite, the protonation energy of the base (the sum of these two is usually endothermic) and the electrostatic stabilization of the Zwitterionic pair of structure II. It will appear that optimization of these backbinding interactions determines the structure of groundstates and especially transition states to a considerable extent. The importance of this zeolite wall effect can be deduced from cluster calculations on the protonation of NH3. This problem has been studied for NH3 interacting with a single T-site cluster as [AI(OH)2H2] H [55-58] as well as large ring systems [33,58-62]. The problem has been studied with different qualities of basis sets, correlation energy as well as cluster relaxation. For historic reasons it may be of interest to mention that Hall et al. [63] were the first to propose that decomposition of NH + to NH3 would lead to the zeolite Br~nsted acidic

428 site we here consider. One finds that NH + becomes only formed, when two or three of its protons are close to the negatively charge (Lewis basic) oxygen atoms around aluminum. When clusters increase in size, these oxygen atoms can become non-equivalent and hence the back binding N H . . . O bonds become non-equivalent. This is observed in cluster calculations [33] as well as electrostatically embedded systems [62]. Electrostatic embedding calculations, especially of small clusters have to be considered with care. Sometimes full electrostatic potentials are used or negatively charged clusters and the embedding is used to correct for this deficiency. It appears that boundary effect corrections can be significant and can make such cluster results quite unreliable. When charges are used derived from quantum-chemical calculations and at least ten-T site clusters are used, the additional energy changes as the energy due to electrostatic embedding have so far been found to be at maximum 10% [33]. Detailed studies of cluster size are available. Brand et al. [60] used a constrained relaxation method and corrected for basis set deficiencies. Whereas the deprotonation energy only slowly converges with cluster size, it appears that the protonation energy of ammonia to ammonium is much less dependent on cluster size. This is a very general observation and derives from a compensation effect. The covalent backbinding interaction of the NH4+ protons changes in a direction opposite from that of the deprotonation energy. The same has been found in a H F / M P 2 study by Zygmunt et al. [64] on H20 interacting with clustermodels of the protonic site in ZSM-5. The most important effect of cluster size is on the geometry of the adsorbed molecule. The non-equivalence of oxygen atoms around A1 tends to make two- or three point adsorption more asymmetric. For H20 and CHaOH, Greatbanks et al. [65] find correct of 17 kJ/mol, 31 kJ/mol, to the heat of adsorption of H20 or CHaOH, respectively. They studied the effect of electrostatics on clusters of 3 T sites. 3-21 G optimized geometries were used and correlation effects were computed in this geometry. Whereas this may cause an overestimate of the embedding effect, interestingly the total adsorption energies computed of 85 kJ/mol for H20 and 110 kJ/mol for CHaOH are very similar to values found by Shah et al. [13] using Car-Parrinello Plane Wave calculations, which necessarily uses DFT. For methanol this value is ~30 kJ/mol higher than the interaction energy computed by Blaszkowski et al. [66] for methanol interacting with a 3 T acid site cluster by DFT. This amounts to a correction of 30% to the energy computed on the small cluster. For the heat of adsorption a correction of the same order of magnitude is necessary for calculations on small clusters. Compared to the single tetrahedron case and especially when clusters with tetrahedral rings are considered two new geometric options arise. Firstly for steric reasons to the backbinding molecule not always three oxygen atoms around A1 may be available, but only two. As we will later see this geometry constraint may inhibit not only formation of particular adsorption modes but also transition states. Secondly in a ring also other atoms than the one connected to A1 may interact with the protonated molecule. Such crossing-overs have been found for NH4 + adsorbed to large ring system clusters by Sauer et al. [58]. Shah et al. [67] and Nusterer et al. [68,69], using Car-Parrinello calculations on extended systems representing chabasite and sodilite found the same for methanol. When the cross connecting oxygen atoms share also a lower valent lattice ion as A1, the oxygen atoms will also be basic and lead to favorable Zwitter-

429 ionic interactions. There is ample crystallographic evidence that in many situations such situations arise [70]. However it is important that experimentally quite high loadings of molecules are used. In such a case a network of hydrogen bonded molecules will be present leading to additional stabilizing interactions. For instance cluster calculations [33,54] agree that a single H20 molecule will only hydrogen bond to a zeolitic proton. For this strong spectroscopic evidence exists, not only from H-NMR [32] but also from the analysis of Infrared spectra [54]. However when two water molecules interact the di-molecular complex will protonate [32]. The additional H20 molecule stabilizes the hydronium ion. NH + is also stabilized by interaction with a second ammonia molecule. This also accompanied by a much weaker interaction with the zeolite lattice oxygen atoms and one point adsorption becomes possible [70,71]. For hydrogen-bonded adsorption two or three point adsorption can be important. This is schematically illustrated for hydrogen bonded methanol in Figure 3a, b and c. In this structure a BrCnsted acidic proton interacts with the basic oxygen from methanol. The proton of methanol backbinds to one of the basic oxygen atoms around A1. Such bonding where a basic alcohol oxygen atom is hydrogen bonded to BrCnsted acidic proton and the alcohol proton is hydrogen bonded to a Lewis basic surface atom has been first proposed by KnSzinger et al. [72]. As mentioned earlier in Car-Parrinello Plane Wave calculations on three dimensional systems backbinding is also observed with more remote oxygen atoms further from the A1 site. This indicates that backbinding in ringsystems connected to the three-dimensional zeolitic oxygen atom is rather weak. The interaction is H(2) - 0(2) (Figure 3b is overestimated because of the equivalence of the oxygen atoms around A1. The calculations Greatbanks et al. [65] are a strong indication of this. The relative small downwards shift of the CH3OH O ( 3 ) - H(2) stretch frequency to 3540 cm-1 [73,74] provides experimental support. An analogous result has been reported for the weakening of H20 stretching frequencies [75]. It is important to realize that DFT [66] as well as HF/MP2 [5] cluster calculations overestimate this bond weakening by more than 200 cm-1. The 3250 cm-1 belongs to the asymmetric stretching frequency of proton [3] in the dimethanol complex sketched in Figure 4. Quantum-chemistry has been extensively applied to study the problem [5,13,65-67] whether CH3OH is adsorbed in the hydrogen bonded state or protonated as in the methoxonium ion (Figure 3c). Analogous studies are available for the protonation problem of H20 [27-30,64,65,70,76] acetone [77,78], acetonitrile [7,79], and dimethyl ether [80]. The earlier reported increases in computed adsorption energy due to embedding indicate that this mainly due to an enhanced H(1)- 0(3) (Figure 3a) interaction. This implies an enhanced Brcnsted acidity of the H(1)- O(1). All cluster calculations so far indicate only hydrogen bonding with this strongly basic molecules. This agrees with the current interpretation of infrared spectra and H-NMR nuclear shifts. Only the more basic pyridine and ammonia are present as protonated species. As demonstrated experimentally [78] there is a parallel between gasfase proton affinity and protonation energy of basic molecules in the zeolite. Pyridine and ammonia are molecules with the highest protonation energy.

430

c~

H

(3)

It" "" \H (2) (l) ' (2) ; A l / ; -Si "Si / /\

/\

H

\Si./~)~AI/O~-si / i N /\ I \

I

CH 3 H/ \ H \si/O~Al/O~si iN /\

/ I\

Figure 3. Adsorption modes of methanol (schematic)" a end-on adsorbed methanol; b side-on adsorbed methanol; c protonated methanol

H3C\ O ...... H_O/CH3 I I

I

"Si "/0/\

/Al

Figure 4. Di-methanol complex (schematic)

S,\

431 An important technical issue in the calculations is the calculation of local energy minimum versus local saddlepoints. The latter correspond to transition states. It has become clear that for strongly Brcnsted basic molecules the difference in energy between the protonated and hydrogen bonded molecule is of the order of 10 kJ/mol. Initially with clusters erroneous results were reported because geometries were used as obtained with the Hartree-Fock method or an artificial geometry constraint had been applied to reduce computer cost. It has become clear that optimized geometries within the H F / M P 2 approximation and a high quality basis set are essential and no symmetry constraint should be applied. Cluster calculations then show as typical result for H20 or CH3OH the double- minimum potential [5,32] as shown in Figure 5. The methoxonium of hydronium ion is a transition state structure a top of the double-minimum and the hydrogen bonded situations are local minima. Of course the local minima should become different in energy for embedded situations. Such potential energy curves are not yet available. Computed H-NMR shifts [5] and QCC constants [6] correspond to a dynamically strongly perturbed H20 or CHaOH molecule that hops between the two minima in the double-minimum potential.

\

/

k2j ,H

III

I

H

.0.,,

H...O

>-.Ar--C~,, zOnal/O,, ,,>

,~

reaction coordinate

a

reaction coordinate

b

Figure 5. Potential energies of adsorbed methanol as a function of proton positions: a. surface atoms equivalent b. surface atoms non-equivalent

A recently published neutron scattering study [76] of water adsorbed in H-ZSM-5 shows very nice agreement between computed loss frequencies and the hydrogen bonded water model. Car-Parrinello Plane Wave calculations [13,68] also indicate very small energy differences between the protonated and hydrogen bonded state. No transition state is yet available in such calculations. In a combined experimental and theoretical study [81,82] of the methane-deuterium exchange reaction, the non-equivalence of the oxygen atoms around A1 has been shown to be responsible for the difference in turnover frequency per A1 in faujasite compared

432 to ZSM-5 (MFI). Because of the larger unit cell of the MFI structure, the differences in deprotonation energy of non-equivalent O atoms is much less in the MFI structure than in the zeolite Y structure.

~AI Os, @c 00 OH

Figure 6. Transition state of hydrogen/deuterium exchange of methane [66]

The computed transition state for the deuterium exchange reaction (Figure 6), in which the carbon atom becomes five coordinated as in a carbonium-ion, with the computed reaction coordinate shows that a different oxygen atom accepts the proton from CH4D + than the oxygen atom that donates deuterium. Using the Br0nsted-Polanyi relation, equation 3 and computed deprotonation energies [44,45] the average TurnOver Frequencies for zeolite Y and ZSM-5 have been computed. The result is a significantly larger TOF for ZSM-5. Important is now also the theoretical prediction that the carbonium ion CH + is not a stable local energy minimum intermediate, but a transition state with its energy maximized as a function of the reaction coordinate. The calculations have been confirmed by Blaszkowski et al. [s31 and Evleth e.a. [84]. This is a very general result and consistent with the previous discussion. Whereas in the case of CH + backbinding to the negatively charged lattice atoms is rather strong and bond distances are such that the interaction can be considered nearly covalent, usually transition states are rather ionic. In the case of H / D exchange of methane the Bronsted-Polanyi relation between /kEact and the difference in deprotonation energies of the oxygen atoms that give H-D exchange is found to be [80]:

AEact : 0.TAEprot.

(5)

However in case one deals with a ionic transition state this dependence has become much less typically

A~act ~---0.1/kEprot.

(6)

If one now compares the difference in reactivity of different sites, the activation will be dominated by differences in deprotonation energy of the Br0nsted acidic proton. The

433 computed bond distances between transition state atom positions and lattice oxygen atoms are often a good discriminator between these situations. In calculations one can also test the Br0nsted-Polanyi relation on deprotonation energies by varying the position of the atoms that bound the cluster. Such an approach has been used for proton activation in alkanes in [83]. These transition states except for H/D exchange are usually found to be quite ionic. Another important carbocation is the carbenium ion formed in superacids by protonation of an olefin. By contrast in solid acids carbenium ions are transition states. Senchenya and Kazansky [85,86] and Pelmenschikov et al. [87] were the first to propose that carbenium ion intermediates are not stable intermediates but instead are transition states. After protonation the protonated olefin becomes converted to a stable alkoxy species (scomplex) [84-86,88-90]. One has to distinguish between the hydrogen bonded ethylene (Tr-complex) and the alkoxy complex (s-complex). Cant and Hall [91] were the first to measure for ethylene adsorption the adsorption energy of the 7r-complex. They found a value of 36 kJ/mol for zeolite Y. The generation of alkoxy species has been demonstrated by NMR spectroscopy [92-94] as well as infrared spectroscopy [95].

Potential energy, kcaYmol

isobutene

~r$

. . . .

J

/

'

__

surface ~-complex

__ __

I

\

surface

\t-buto~d,,

,~

RCH=CH2Ig)

RCH~CH2 \

H I /Ox

~

/AL ,,Six/ I Separated

Acid & Base

H I ~A[ /

/ .~ \O"si

II

+ ~

RCHCH 3

R.,C# ~

i CH3

\ /OT. / /AL\ ,,Six

\ /0. /AL\ ,,Six/

Ill

IV

H- Bonded Carbeniumion Alkox,de [omplex wffh s Base

Figure 7. Reaction energy diagram of the protonation of isobutene [44] and the corresponding steps for a therminal alkene.

Figure 7 gives the computed energy reaction diagram representative of the energy changes in the different elementary reaction steps that leads to formation of the protonated alkene. The energies given have been computed for the protonation of isobutene on a three tetrahedral cluster. It is important to realize that the energy changes are only due to the interaction of the p electronic system of the molecule upon interaction with the zeolitic proton. Hence the free state of zero interaction energy is the state of isobutene adsorbed in the micropore of a zeolite. There the interaction energy is dominated by the

434 Van der Waals interaction and is strongly micropore dimension dependent [96]. If one wishes to compare the computed protonation energy (s-complex) with the measured heat change when a molecule adsorbs from the gasfase, one has to add to the computed protonation energy the heat of adsorption of isobutene in a siliceous zeolite channel. With respect to the gas phase one now realizes that the heat of protonation is different for different zeolites because of the differences in Van der Waals interaction of the adsorbed moleclile~ with the zeolite micropore. The interaction with the zeolite proton is often nearly tile same. In Figure 7 the different interaction stages are shown. As mentioned the reference state is isobutene adsorbed in zeolitic micropore. A first step is the formation of a hydrogen bond between proton and olefin to form a 7r-bonded complex. This is a non activated process and corresponds to a relatively small interaction energy. Protonation of the olefin is an activated process. The OH bond has to be ruptured. The transition state (Figure 7) has the geometry of a carbenium ion. Protonation occurs to the primary carbon atom and most of the positive charge becomes located on the tertiary carbon atom that has become part of a planar geometry and can be considered sp2 hybridized. The positively charged side of the molecule is directed towards one of the negatively charged oxygen anions, which results in a strong electrostatic stabilization. The transition state for ethylene protonation is rather covalent and becomes more ionic for isobutene. The ground state that corresponds to the protonated ethylene molecule is the s-bounded complex, the alkoxy species. The difference in energy between this state and the transition state energy is of course dominated by the covalent O-C bond, where strength only slightly varies between the different oxygen atoms. The difference in relative stability of the carbenium ion transition states controls the activation energy differences. NMR experiments have shown that a stabilized carbenium-ionic state can be converted to stable carbenium ion. An example is the cyclopentenyl cation [94]. 3. M o l e c u l a r d y n a m i c s 3.1. I n t r o d u c t i o n In recent years, molecular dynamics (MD) has been used in an increasing number of studies on zeolite systems. These studies include both the adsorption and diffusion of guest molecules and the structure and dynamics of zeolite lattices. Some of these studies will be discussed below to give the reader an idea of the scope and limitations of MD simulations but first a short introduction to this computational tool is given. For a more thorough discussion of MD the reader is referred to the books by Allen and Tildesley [97] or Frenkel and Smit [98]. 3.2. M D as c o m p u t a t i o n a l t o o l In MD classical mechanics is used to study the evolution of a system in time. First the trajectory of the system is generated. In the often used Verlet algorithm this is done by a Taylor expansion (until third order) of the positions (r) for each of the particles (atoms) around time t: 1 d2r(t) 1 d3r(t) At + - ~ A t 2 + At 3, (7) r(t+At) - r(t)+ dr(t) dt 2 dt 2 6 dt 3 1 d2r(t) 1 d3r(t) dr(t) ~Xt + ~xt 2 At 3 (8) r(t-At) -- r(t)- dt 2 dt 2 6 dt 3 "

435 Combination of these equations, in which At is the time step of the simulation, gives: + At) = 2

(t) -

A t ) + d2r(t) At 2

dt 2

(9)

.

which can be solved in combination with Newton's equation of motion: d2r(t) --- FNewton -dt 2

7rt~

db/pot dr

(10)

In the latter equation rn is the mass of the particle (atom) and b/pot is the potential energy. The velocity of the particle follows from:

v(t)-

dr(t) r(t + At) + r ( t - At) dt = 2At

(11)

From the velocities the kinetic energy of the system can be calculated: 1 N /a/kin - ~ . = ~

(12)

In equation (12) N is the number of particles (atoms) in the system. The temperature of the system follows from: 2 b/kin

T(t) = 3 Nk.

(13)

in which kB is the Boltzmann constant. Generally, MD simulations are carried out on rather small systems (N < 104). Surface or wall effects can be avoided and the link with real systems can be preserved by implementing periodic boundary conditions. In this technique the simulation box is replicated in space to form an infinite lattice. The trick is that, when during the simulation a molecule (atom) leaves the central box, one of its images will enter the box at the opposite side. The effect of periodic boundary conditions on the thermodynamic properties is generally small but certain properties cannot be studied in this way [97]. With simulations in the N V E ensemble (constant number of particles, constant volume, constant total energy) the instantaneous kinetic energy fluctuates around an average but the total energy b/tot = b/kin(t) + b/pot(t) should be constant (conserved) during the simulation. This is only achieved after the system has been equilibrated long enough to erase all start-up effects and when the time step t is sufficiently short (order of 10 -15 s (fs)). During equilibration the temperature of the system can be set approximately by coupling it temporarily to a temperature bath. Such a coupling, in which the velocities and the positions of the particles are corrected when the temperature calculated with equation (13) deviates from the simulation temperature that is desired, is central in simulations in the N V T ensemble (constant number of particles, constant volume, constant temperature). The NVT algorithm is less expensive in computer time because it largely avoids the long equilibration period. It has been argued [1] that the N V T algorithm reflects physical reality better in simulations in which the zeolite lattice is kept fixed. To simulate the trajectory one has to calculate in each time step (see equation (10)) the potential energy and its derivatives (the forces on the particles (atoms)). In most

436 simulations on zeolite systems a force field has been used for this purpose. Each force field, and there are many, has to be validated against experimental data or against results from ab initio calculations. It is clear that a quantum chemical calculation of the potential energy and the forces would be preferable. However, at the moment the computational costs of such a simulation are prohibitive. After the trajectory of the system has been generated one can find the thermodynamic properties < A > by averaging over time: 1 n < A > - - ~ A(ti). V-] -__

(14)

where n is the number of sampling points and A(t) is the value of the particular property at a specific time. Examples are the temperature in an N V E simulation and the heat of adsorption for a guest molecule in a zeolite. When the simulation is performed with a rigid lattice the latter can be calculated from: A Hads < A~ad s >

=

RT-

=

< A/~/nb > -nt- < A/~guest > ,

< Ab/ads >,

(15) (16)

in which < m/~nb > is the average value of the non-bonded interactions and < Ab/guest > is the average strain in the guest molecule as a result of adsorption. The calculation of some other properties from the simulated trajectory will be discussed below. In general, for good statistics the number of sampling points n should be large, which makes MD simulations expensive both in computer time and in storage requirements because the trajectory has to be saved for later analysis. Of course the possibilities in MD simulations depend strongly on the length of the time step, the number of the degrees of freedom and the available computer facilities. At the moment the practical limit seems to be a total simulation time of a few nanoseconds. 3.3. A d s o r p t i o n a n d d i f f u s i o n Most MD studies in zeolite science were carried out on the diffusion of relatively small molecules in zeolite channels. The self-diffusion coefficients can be obtained from the simulated trajectories with the Einstein relationship: D - lim t--,o~

< zx~(t) 2 > 6t

(17) '

in which < Ar(t) 2 > is the averaged mean square displacement of the centre of mass of the molecules: < Ar(t) 2 > = < Jr(t0 + t) - r(t0)] 2 > .

(18)

The averaging is generally carried out over all molecules and over independent time origins to. The latter implies that different points along the trajectory are taken as origin for a part of the trajectory. Of course, the limit t --+ oe (equation (17)) cannot be reached but total simulation time should be sutficiently long to obtain accurate values for D; e.g. Bull et al. [99] estimated that a simulation time of around 200 ps is required for an accurate determination of D

437 for benzene in NaY at room temperature. In practice [100] this limits the application of MD in diffusion studies to molecules with D > 10 -1~ m2s -1 These are rather small molecules; in the n-alkane series n-hexane is the largest member for which diffusion has been studied by MD [101]. For large molecules other modelling techniques have to be used. One possibility is to use transition state theory with a hopping-model in which molecules pass an energy barrier when they migrate from one adsorption site to the other. Mosell et al. [102] have used MD simulations to obtain transmission coefficients for cage-cage migration of xenon in zeolite NaY. These transmission coefficients are then used [103] as corrections in a simulation based on transition state theory. The results are in accordance with those of conventional MD simulations. Both the value of the diffusion coefficient and that of the activation energy are in excellent agreement. This validates the use of transition state theory in simulations of molecules with diffusion coefficients too small on the time scale of conventional MD simulations. An interesting approach is also that of Maginn et al. [104] who use concepts from Brownian motion theory and transition state theory in a study of long n-alkanes, up to C20, in silicalite. Most MD simulations on diffusion of molecules in zeolites have been carried out with a rigid framework and often also with a rigid molecule. The rigid framework approximation neglects the influence of the framework vibrations on the diffusion motion. It appears that for small molecules this does not make a significant difference in the calculated values of D and AHads. The framework vibrations, however, do influence the pair distribution functions [100,105]. The calculated diffusion coefficients generally agree well with experimentally obtained microscopic diffusion coefficients (see Figure 8)). The purpose of most MD studies on diffusion in zeolites, however, was not the numerical value of the self-diffusion coefficients but to get microscopic picture of the diffusion mechanism. One of such studies [102,103] has been discussed above; a few more examples follow. The diffusion of methane in silicalite-1 has been studied by several groups [100,106-111]. Figure 9 taken from the work of Goodbody et al. [108] shows the expected anisotropy of the diffusion and its dependence on loading. Shen and Rees [112] studied the diffusion of carbon dioxide in the pure siliceous frameworks of theta-1 and silicalite-1. Their results (see Figure 10) show that the number of CO2 molecules in the simulation box does not affect significantly the diffusion in the three-dimensional channel system of silicalite-1 but it does in the one-dimensional channel system theta-1. This is ascribed to single-file diffusion of CO2 in theta-1. The activation energy for diffusion can be obtained from the temperature dependence of the diffusivity (see Figure 11). Again there is reasonable agreement between the experimental observations and the data obtained by simulation. Single-file diffusion can be inferred from MD simulations by an increase of the mean square displacement < At(t) 2 > that is proportional to the square root of the time [113,114]. In the study of Shen and Rees [112], however, the total simulation time was too short for this observation. Smirnov and Van de Graaf [111] studied the diffusion of methane in silicalite-1 and 2 with a CIEEM MD simulation. The CIEEM force field [115] has geometry dependent atomic charges and allows polarization to be taken into account in the MD simulation. It was found that the polarization of the methane molecules reaches a maximum at the edges of the channel intersections. Also there is a significant difference between silicalite-1

438 ~-. - 1 8

J

Y -19

9

D

9

rn

o~m

-20 9

C~

9 9

9

9

rn

o l

I t -21

..

2

:

~

!

I

3

4

5

6

7

1/Txl0 ~

Figure 8. Log plot of diffusion coefficients for 12 molecules of methane per unit cell: computed from MD simulations, (open triangles) vibrating framework, (solid triangles) fixed framework; experimental values, (solid squares) ZSM-5 with various aluminium content and crystal size and (open squares) silicalite versus 1/T. (reproduced with permission from ref. [99]).

439 30 25

,,=,,,.

2C t~

E e

~ 5r

O q.-

05 tl

'3O 0

5

10 N=

15

Figure 9. Diffusion tensor elements (in increasing order: c-axis, a-axis , b-axis) for methane in silicalite at 298 K as a function of loading. (reproduced with permission from ref. [108].

3~K

400

400

(a)

~o

/

(a)

373 K

t/"

~K

2O0

100 ~ ,,-,,,

o

195N

/

195 K

[

v

" " 100 r (b)

373K ,

373 K 200

323 K

100

273 K

50

220 K

~

so

~Oo time/ps

150

0

50

100

150

time/ps

Figure 10. Mean-square displacements of carbon dioxide in silicalite-1 (on the left) and theta-1 (on the right) (a) one and (b) two molecules per simulation box. (reproduced with permission from ref. [112]).

440 10 .7

10 -e,

2 r

10-~.

"-4. 10 -lo

io

is

3[o

3'5 4'o IOOWT

,'s

io

s'5

Figure 11. Temperature dependence of carbon dioxide diffusivity in silicalite: (open squares) frequence response measurements; (open circles) microscopic pulsed-field gradient NMR; (solid circles, solid squares) MD with one and two molecules per box, respectively. (reproduced with permission from ref. [112]).

and 2 with respect to the polarization of the methane molecules (see Figure 12). It was expected that this difference will be expressed in a difference in catalytic activity. Yashonath and co-workers have performed extensive MD simulations to investigate the effect of the diameter of the adsorbed molecules on the diffusion coefficient [116-122]. Figure 13 shows the diffusion coefficient of Lennard-Jones molecules adsorbed in NaCaA as a function of the diameter of these adsorbed molecules. The simulations indicate that as the diameter of the molecules is increased the diffusion coefficient first decreases. However, if this diameter has approximately the same size as the window between the cages, then the diffusion coefficient increases and a maximum is observed. A further increase of the diameter of the adsorbed molecules gives again a decrease of the diffusion coefficient. Such a maximum or "levitation effect" has been observed in simulations of various other zeolites (Y [117,120], silicalite [118,122], and VPI-5 [122]) and is both observed in simulations using a rigid [117] and a flexible [121] zeolite framework. Experimental evidence of this effect has not yet been reported. The diffusion coefficient in zeolite A is determined by the hopping of the adsorbed molecules from one cage to another through a window with a given diameter. If one inserts molecules with a larger diameter, the diffusion coefficient is expected to decrease since it will be more difficult for those bigger molecules to go through this window. Yashonath and co-workers discovered that indeed if the molecules are small compared to the window diameter the window is an energetic barrier. However, for molecules that have approximately the same diameter as the window this energy barrier disappears and hence the diffusion is larger that one would expect. Of course, if the molecules are too big the diffusion drops down to zero.

441

8.5 "T

8.4

~.~

~.0

~ f. ~-

-~ 7.6

1 '~~

~

7.0 f

E

~ 0.6 E

"~

~o. 0.4 "0 /

,

0

I

2

~

I

,

J

4 6 b -axis coordinate/,~

,

I

8

,

~

10

. 02

""~1~

0

2

4 a-axis

6

8

coordinate/A

Figure 12. Dependence of the molecular electronegativity and the molecular dipole moment (D) of methane molecules on the position of the molecules in the channel between the intersections in MFI (left; straight channel) and MEL (right). Each point corresponds to the average over all segments. The minimum and maximum values of the horizontal axis correspond to the centres of the intersections. (reproduced with permission from ref. [111]).

:,

10

442

'

'

I

'

I

i

0

I

'

I

,

I

o

oo 0 . 4 04

E o

V

o

0.2

x

a

0.0

i

0.3

I

0.5

,

0.7

0.9

I

1.1

7

Figure 13. Diffusion coefficient of Lennard-Jones particles as a function of their reduced size in the zeolite NaCaA. The size of the molecule is written as V = 21/6ass/aw, where ass is the diameter of the Lennard-Jones molecules and aw the diameter of the window. The factor 21/6 is introduced to take into account the complete repulsive part of the potential. On the x-axis the Lennard-Jones diameters of the nobles gasses are indicated. The line is a guide to the eye. The simulation results are taken from ref.[ll9].

Another interesting phenomena that has been speculated upon theoretically (see references in [123]) is the so-called resonant diffusion. These theories predict the diffusion coefficient as a function of chain length to be periodic. This periodicity is related to a matching of the size of the molecule with the periodicity of the lattice. Runnebaum and Maginn [124] have performed extensive Molecular Dynamics simulations of the linear alkanes ranging form C4 to C20 in silicalite and found evidence for this phenomena. 3.4.

Vibrational

characteristics

of adsorbed

molecules

When the molecules are treated as flexible in the MD simulation it is relatively simple to calculate the vibrational characteristics from the trajectory. The power spectrum is obtained by Fourier transform of the velocity auto-correlation function < v ( t ) 2 >. A power spectrum results from all motion because in its calculation no selection rules are applied. The IR spectrum can be obtained by Fourier transform of the dipole moment auto-correlation function < #(t)2 >. An example of a calculation [111] is given in Figure 14. It shows the power spectrum and IR spectra computed for methane in silicalite-1. In the power spectrum the band at 100 cm -1 relates to translational and rotational degrees of freedom. The two IR spectra were computed with different force fields: one with fixed atomic charges (FCA) and one with geometry dependent charges (CIEEM force field). In the latter spectrum a low-intensity band is present at 2826 cm -1 that is absent in the FCA spectrum. The band relates to the normally forbidden symmetric stretch and is observed in the CIEEM spectrum as a result of polarization of the methane molecules by the zeolite lattice.

443

C

i' J

. . . . . . . . .

0

500

i . . . .

1000

I . . . . . . . . . . . . . . . . . .

1500

2000

2500

3000

.1 3500

w a v e n u m b e r / c m -~

Figure 14. Calculated vibrational spectra of methane molecules adsorbed in MFI: power spectrum of the hydrogen atoms (a) and IR spectra from a MD run using a fixed charge approximation (b) and a CIEEM MD run (c). (reproduced with permission from ref. [111]).

444

The applications of Raman scattering in zeolite science and the calculation of Raman spectra from MD simulations for molecules adsorbed in zeolites have been reviewed recently by Br~mard and Bougeard [125]. 3.5. Zeolite s t r u c t u r e and lattice d y n a m i c s MD is also used increasingly in studies on the structure of zeolites and on lattice dynamics. Of course, such studies depend on an adequate force field being available. The advantage of MD studies over energy minimization studies is that MD is closer to physical reality: entropic effects are automatically taken into account and in the calculation of the lattice vibrations it is not necessary to invoke the harmonic approximation. The calculation of the lattice vibrations follow the same procedure as the calculation of the vibrational characteristics of adsorbed molecules. Again the results of a few studies will be discussed to illustrate the scope of MD simulations for these subjects. Yamahara et al. [126] have studied the thermal behavior of silicalite-1 with MD. In agreement with the experiment they find that the monoclinic structure is the stable one at low temperature and that the orthorhombic structure is the stable one at high temperature (see Figure 15). They note that the observation of this particular phase transition in MD simulations depends critically on the balance of forces in the applied force field. Several published force fields give only one of the two structures.

908

w

!

o

= 906 o~

@

....

"o ""

o

-

Alpha

e--

w

-- -" - Gamma

904

CrJ C

e~ 9 0 . 2 r~ 90.0 I

200

|

I

1

I

i

i

1

1

400 6oo eoo Iooo 12oo T e m p e r a t u r e (K)

Figure 15. Temperature dependence of unit cell parameters (angles) of MFI obtained from MD runs starting with the orthorhombic structure. (reproduced with permission from ref. [126]).

MD simulation has also be applied to isomorphic substitution. Oumi et al. [127] have investigated the titanium substitution in TS-1. They studied the influence of this substitution on the lattice parameters and by comparing their results with experimental data (X-ray diffraction) they concluded that T8 was the most probable site for substitution. However, in a more recent study [111] it was shown that the experimental trends in the lattice parameters can be reproduced satisfactorily with titanium atoms placed at random

445 T positions (see Figure 16). Nevertheless, the approach by Oumi et al. is quite original as it does not a priori assumes thermodynamic equilibrium as has been done in all other theoretical studies on isomorphic substitution.

20.2 a - axis

1;I 9

b - axis o 9 0 mmi n

g

9

C~- axis9 9 [;I II

[]

9

g

20.1

==~=

'<~ 20.0 d m

> 19.9 E 9 -~ t~

19,8

Q" 13.5I 13,4 l-

[

13.3t

I

I

0,00

0.01

=

I

~

0,02

I

0,03

,

'

0,04

['ri]l([Ti]+[Si])

Figure 16. Dependence of the lattice parameters on the [Ti]/([Ti] + [Si]) ratio. Open squares denote calculated (MD) values, solid squares denote experimental values. (reproduced with permission from ref. [128]).

Smirnov and Bougeard [129] have studied the lattice dynamics of some siliceous zeolites with MD. Their results show that for zeolite A the presence of sodium or potassium ions leads to drastic changes in the distribution of the window diameter (see Figure 17). In view of these results one may well question simulations on the adsorption of larger molecules in which the zeolite lattice is kept fixed. MD simulations can provide a microscopic explanation of features observed in experimental spectra. An example [128] is given in Figure 18 which provides an answer to the origin of the band at 960 cm -1 in the vibrational spectra of titanium substituted zeolites. From the power spectra calculated for the TiO4 tetrahedra it is concluded that the band relates to the localized Si-O mode in the Si-O..-Ti bridges.

446

y--A

5.5

6

6.5

7

7.5

8

Diameter (s

Figure 17. The distribution function of the window diameter of zeolite A without cations (empty-A), Na zeolite A (Na-A) and K zeolite A (K-A). (reproduced with permission from ref. [129]).

447

b

a . ,

,

i

,

200

,

i

i

,

400

,

i

!

9

9

600

W a v e n u m b e r ,

0

.

200

400

600

W a v e n u m b e r ,

=

I

.

800 r

.

i

.

.

1000

.

I

.

1200

"1

800

1000

1200

c m "~

Figure 18. At the top, calculated IR spectra for silicalite (a) and Ti containing MFI (b); at the bottom, spectra of oxygen vibrations in TiO4 tetrahedra: (a) SiO4 in silicalite, (b) in SiO4 and (c) in TiO4 tetrahedra in the Ti-containing MFI, and (d) in TiO4 tetrahedra in the model in which the same values of force constants for Ti-0 and Si-0 bonds are used. (reproduced with permission from ref. [128]).

448 4. M o n t e C a r l o s i m u l a t i o n s 4.1. I n t r o d u c t i o n In the previous edition of this book, the basic Molecular Dynamics and Monte Carlo techniques have been discussed [1]. At that time, standard simulation techniques were used for the first time to simulate the adsorption [130-132] and diffusion [133,134] in zeolites. As shown in the previous section, these simulations gave important insights in the behavior of small molecules (noble gasses, methane, or ethane) adsorbed in the pores of the zeolites. A review of these simulation results can also be found in ref. [135]. These simulations also made clear that simulating long chain molecules would require much more cpu time than one could afford. Moreover, one would have to wait many decades before computers would be sufficiently fast such that one would be able to perform such simulations. Next we discuss some of the recent attempts to develop schemes to make the simulations of such systems more efficient. Before we start with the technical aspects of the Monte Carlo technique, let us discuss the previous statements of the limitations of conventional simulation techniques in more detail. Since the conventional simulation techniques such as Molecular Dynamics and the Monte Carlo technique are equally valid for simple and complex fluid, one may wonder whether a simulation of a complex system is simply a question of changing the force fields. In some cases it may indeed be as simple as that, in particular if the computers are sufficiently powerful to deal with these more complex systems in a reasonable amount of time. For some problems, such as the sorption of long-chain hydrocarbons, however, the increase in computer time can be prohibitively large; it may take many years of super-computer time before a calculations is finished. Let us illustrate this with a specific example. Over the last decade many simulation studies on the sorption of molecules in zeolites have been published. A careful look at these studies reveals that most simulations concern the adsorption of noble gases or methane, only a few studies of ethane or propane have been published. In petro-chemical applications of zeolites, however, we are interested in the behavior of much longer alkanes such as octane and decane. The reason why only small molecules have been studied becomes clear from, for example, the work of June et al. [136] and Herngndez and Catlow [101], in which Molecular Dynamics simulations were used to investigate the diffusion of butane and hexane in the zeolite silicalite. ,June et al. showed that the diffusion of butane from one channel of the zeolite into another channel is very slow compared to diffusion of bulk butane. As a consequence many hours of super-computer time were required to obtain reliable results. In addition, these results show that the diffusion coefficient decreases significantly with increasing chain length. Hence, extrapolation of these results suggests that many years of super-computer time would be required to obtain comparable results for the longer alkanes. The above example illustrates one of the main limitations of Molecular Dynamics, in such a simulation, the approach is to mimic the behavior of the molecules as realisticly as possible. If successful, all properties will be like in nature, including the diffusion. If the molecules diffuse slowly this will be reflected in very long simulation times and in the case of long chain alkanes these simulation times turned out to be much longer than we can currently afford. In principle, one can circumvent this intrinsicly slow dynamics by

449 using a Monte Carlo technique. In a Monte Carlo simulation one does not have to follow the "natural path" and one can, for example, perform a move in which it is a t t e m p t e d to displace a molecule to a random position in the zeolite. If such a move is accepted, it corresponds to a very large jump in phase space. Again, utilization of these type of un-natural Monte Carlo moves turned out to be limited to small molecules. For example, Goodbody et al. [108] have used this Monte Carlo trick to determine the adsorption isotherms of methane in a zeolite. In such a simulation one can observe that out of the 1000 attempts to move a methane molecule to a random position in the zeolite 999 a t t e m p t s will be rejected because the methane molecule overlaps with a zeolite atom. If we were to perform a similar move with an ethane molecule, we would need 1000 x 1000 a t t e m p t s to have one that was successful. Clearly, this random insertion scheme will break down for any but the smallest alkanes. The above example, the adsorption of chain molecules in the pores of a zeolite, is used to illustrate the problems that may occur if one uses a conventional simulation techniques to more complex systems. It is interesting to note that similar problems may occur in the simulation of phase equilibria of chain molecules, simulations of polymers, or liquid crystals. For many of these systems it is relatively straightforward to implement the force fields to simulate these systems, however, the required simulation times to determine reliable equilibrium properties may be prohibitively long. These simulation times may even be so extreme that it can not be expected that increasing computer power will be of any help. To be able to perform simulations on complex systems it is therefore important to develop novel algorithms that are orders of magnitude more efficient than the conventional algorithms. In this chapter such algorithms are discussed. 4.2. S i m u l a t i o n t e c h n i q u e s We start the discussion on simulation techniques with a short review on the conventional (Metropolis) Monte Carlo technique. For a more extensive discussion the reader is referred to the book by Allen and Tildesley [97] or the more recent book by Frenkel and Smit [98]. 4.2.1. B a s i c M o n t e C a r l o s i m u l a t i o n s The prime purpose of the kind of Monte Carlo program that we shall be discussing is to compute equilibrium properties of classical many-body systems. The Monte Carlo method is a technique to generate a sequence of configurations, each of these configurations occurs with a certain probability and the Monte Carlo scheme should be developed in such a way that the desired probability distribution will be sampled. In the canonical ensemble (constant temperature, volume, and number of particles) we have to sample a Boltzmann distribution:

dV'(n) (x exp[-flb/(n)]

(19)

which states that the probability of finding configuration n is proportional to the exponent of the energy of this configuration/d(n) times the Boltzmann factor fl = 1 / k B T . One way of generating these configurations is, the basic Monte Carlo algorithm, as introduced by Metropolis et al. [137], proposed: 1. Select a particle at random, and calculate its energy Lt(rN).

450 2. Give the particle a random displacement; r ' = r + A, and calculate its new energy b/(r'N). 3. Accept the move from r N to r 'N with probability acc(o -+ n) - rain (1, exp{-fl[b/(r 'N) - / ~ ( r N ) ] } ) .

(20)

The above algorithms assume we have a simple atomic system. If we are simulating molecules rather than atoms we must also generate trial moves that change the molecular orientation. A possible approach is, for example, to consider a system consisting of N linear molecules. We specify the orientation of the ith molecule by a unit vector fii. One possible procedure to change fii by a small, random amount is the following. First, we generate a unit vector ~ with a random orientation. Next we multiply this random unit vector ~ by a scale factor 7. The magnitude of 7 determines the magnitude of the trial rotation. We now add 7~ to fii. Let us denote the resulting sum vector by t: t = 7"~ + fii. Note that t is not a unit vector. Finally, we normalize t, and the result is our trial orientation vector fi~. We still have to fix 7, which determines the acceptance probability for the orientational trial move. We have not yet indicated whether or not the translational and orientational trial moves should be performed simultaneously. Both procedures are acceptable. However, if rotation and translation correspond to separate moves, then the selection of the type of move should be probabilistic rather than deterministic. Only slightly more complex is the case of a nonlinear rigid molecule. It is conventional to describe the orientation of nonlinear molecules in terms of the Eulerian angles (~b, 0, ~). However, for most simulations, use of these angles is less convenient because all rotation operations should then be expressed in terms of trigonometric functions, and these are computationally expensive. It is usually better to express the orientation of such a molecule in terms of quaternion parameters (for a discussion of quaternions in the context of computer simulation, see [97] or [98]). If the molecules under consideration are not rigid then we must also consider Monte Carlo trial moves that change the internal degrees of freedom of a molecule. In practice, it makes an important difference whether or not we have frozen out some of the internal degrees of freedom of a molecule by imposing rigid constraints on, say, bond lengths and possibly even some bond angles. If not, the situation is relatively simple: we can carry out normal trial moves on the Cartesian coordinates of the individual atoms in the molecule (in addition to center of mass moves). If some of the atoms are strongly bound, it is advisable to carry out small trial moves on those particles (no rule forbids the use of trial moves of different size for different atoms, as long as the moves for one particular atom are always sampled from the same distribution). However, when the bonds between different atoms become very stiff, this procedure does not sample conformational changes of the molecule efficiently. In Molecular Dynamics simulations it is common practice to replace very stiff intramolecular interactions by rigid constraints. For Monte Carlo simulations this is also possible. In fact, elegant techniques have been developed for this purpose [138]. However, the corresponding MD techniques [139] are so much easier to use, in particular for large molecules, that we cannot recommend the use of the Monte Carlo technique for any but the smallest flexible molecules with internal constraints.

451 4.2.2. M o n t e Carlo s i m u l a t i o n s of chain m o l e c u l e s It is possible to perform simulations beyond the Metropolis technique. These type of simulations have been developed by Siepmann and Frenkel [140] to make Monte Carlo moves of long chain molecules on a lattice possible. This so-called configurational-bias Monte Carlo (CBMC) technique for is based on the early work of Rosenbluth and Rosenbluth [141] and Harris and Rice [142]. This technique has since been extended to continuum models by Frenkel et al. [143] and de Pablo et al. [144].

4.2.3. B e y o n d M e t r o p o l i s The general idea of biased sampling is best explained by considering a simple example. Let us assume that we have developed a Monte Carlo scheme that allows us to generate trial configurations with a probability that depends on the potential energy of that configuration: c~(o-+ n ) = flU(n)]. For the reverse move, we have a(n -+ o) = f[b/(o)]. Suppose we want to sample the N , V , T ensemble, which implies that we have to generate configurations with a Boltzmann distribution (19): To proof that the correct distribution is sampled, we have to demonstrate that detailed balance is obeyed

K ( o ~ n) = K ( n --~ o),

(21)

where K ( o --+ n) is the flow of configuration o to n. This flow is given by the product of the probability of being in configuration o, the probability of generating configuration n, and the probability of accepting this move, I<(o

-+

= H(o) •

-+

•

-+

(22)

If we imposing detailed balance, we get as a condition for the acceptance rule, acc(o-+ n) = f[b/(n)] e x p { - f i [ b / ( n ) - b/(o)]}. acc(n --+ o) /[b/(o)] A possible acceptance rule that obeys this condition is acc(o --+ n) : min (1, f[/g(o)]f[/g(n)]exp{-/5[b/(n) - / g ( o ) ] } ) .

(23)

It is instructive to apply the above to the Metropolis scheme. For the Metropolis scheme, there is no difference in the probability to generate the old or new configuration: -+

=

-+

o) =

If we substitute this equation in equation (23), we recover the Metropolis acceptance rule (20). This derivation shows that we can introduce an arbitrary biasing function f (/4) in the sampling scheme and generate a Boltzmann distribution of configurations, provided that the acceptance rule is modified in such a way that the bias is removed from the sampling scheme. We now illustrate the use of non-Metropolis sampling techniques to demonstrate how they can be used to enhance the efficiency of a simulation of a system containing chain molecules.

452

4.2.4. Configurational-Bias Monte Carlo The starting point for the configurational-bias Monte Carlo technique is the scheme introduced by Rosenbluth and Rosenbluth in 1955 [140-144]. The Rosenbluth scheme itself also was designed as a method to sample polymer conformations. A drawback of the Rosenbluth scheme is, however, that it generates an unrepresentative sample of all polymer conformations; that is, the probability of generating a particular conformation using this scheme is not proportional to its Boltzmann weight. Rosenbluth and Rosenbluth corrected for this bias in the sampling of polymer conformations by introducing a conformationdependent weight factor W. However, as was shown in detail by Batoulis and Kremer [145], this correction procedure, although correct in principle, in practice works only for relatively short chains. The solution of this problem is to bias the Rosenbluth sampling in such a way that the correct (Boltzmann) distribution of chain conformations is recovered in a Monte Carlo sequence. In the configurational-bias scheme to be discussed next, the Rosenbluth weight is used to bias the acceptance of trial conformations generated by the Rosenbluth procedure. As we shall show, this guarantees that all chain conformations are generated with

the correct Boltzrnann weight. The idea of the CBMC scheme is to grow a molecule atom by atom. During the growing, a set of trial orientations is generated and out of these trial orientations the optimum orientation is selected with the highest probability. The bias which is introduced by this growing scheme is subsequently removed by adjusting the acceptance rules. The configurational-bias Monte Carlo algorithm consists of the following steps: 1. Generate a trial conformation using the Rosenbluth scheme to grow the entire molecule, or part thereof, and compute its Rosenbluth weight W(n). 2. "Retrace" the old conformation and determine its Rosenbluth factor, W(o). 3. Accept the trial move with a probability acc(o ~ n) = mini1, W(n)/W(o)].

(24)

The generation of a trial conformation n of a chain consisting of f atoms is generated using an algorithm based on the method of Rosenbluth and Rosenbluth. An important point that we have to consider is the way in which trial conformations of a chain molecule are generated. For example, one could generate trial segments with orientations distributed uniformly on a unit sphere. However, for many models of interest this procedure is not very ef-ficient, in particular when there are strong intramolecular interactions (e.g., bending and torsion potentials). The efficiency of a configurational-bias Monte Carlo algorithm depends to a large extent on the method used of generating the trial orientations. For example, an isotropic distribution of trial directions is well suited for completely flexible chains. In contrast, for a stiff chain (e.g., liquid-crystal forming polymer), such a trial position will almost always be rejected because of the intramolecular interactions. From the preceding discussion, it follows that the intramolecular interactions should be taken into account in generating the set of trial conformations Here, we consider the case of a flexible molecule with contributions to the internal energy due to bond bending

453 and torsion. The fully flexible case then follows trivially. Consider a chain of g linear segments, the potential energy of a given conformation/4 has two contributions: 1. The bonded potential energy b/b~ is equal to the sum of the contributions of the individual joints. A joint between segments i and i + 1 (say) has a potential energy t h a t depends on the angle 0 between the successive segments. For instance, u bond i ibond (0) could be of the form ai ~ bond ( O ) = k o ( O - Oo)2 . For realistic models for polyatomic molecules, ~ bond includes all local bonded potential energy changes due to the bending and torsion of the bond from atom i - 1 to atom i. 2. The ezternal potential energy U ext accounts for all interactions with other molecules and for all the nonbonded intramolecular interactions. In addition, interactions with any external field that may be present are also included i n / 4 ext . In what follows we shall denote a chain in the absence of the external interactions as the ideal chain. Note that this is a purely fictitious concept, as real chains always have nonbonded intramolecular interactions. To perform a configurational-bias Monte Carlo move, we apply the following "recipe" to construct a conformation of a chain of g segments. The construction of chain conformations proceeds segment by segment. Let us consider the addition of one such segment. To be specific, let us assume that we have already grown i - 1 segments and are trying to add segment i. This is done in two steps. First we generate a trial conformation n, next we consider the old conformation o. A trial conformation is generated as follows: 1. Generate a fixed number (say k) trial segments. The orientations of the trial segments are distributed according to the Boltzmann weight associated with the bonded interactions of monomer i (~i ~ bond~j. We denote this set of k different trial segments by {b}k = { b l , " . , b k } , where the probability of generating a trial segment b is given by

_bond Pi

exp[--flub~

(b)] db

( b ) d b - f db exp[--/3ub~

i (b)]db. = C e x p [ - "p u bond

(25)

2. For all k trial segments, we compute the external Boltzmann factors exp[-fiu~ xt (bi)], and out of these, we select one, denoted by n, with a probability

peXt

i (bn)-

exp[--fluext (bn)] toext i (n)

(26)

where we have defined

k wext(?'t) : E exp[-/buext(bj)]" j=l

(27)

454 3. The selected segment n becomes the ith segment of the trial conformation of the chain. 4. W h e n the entire chain is grown, we calculate the Rosenbluth factor of the chain: wext(n)

g -- l-I-ext

(2s)

i=1

where Rosenbluth factor of the first monomer is defined by weXt

1 (n)-

k e x p [ - p~u ext i (rl)]

(29)

where r l is the position of the first monomer. For the old configuration, a similar procedure to calculate its Rosenbluth factor is used. 1. One of the chains is selected at random. This chain is denoted o. 2. The external energy of the first monomer is calculated. This energy involves only the external interactions. The Rosenbluth weight of this first monomer is given by ~

weXt(o) - k

ext

(o)].

(30)

3. The Rosenbluth factors of the other g - 1 segments are calculated as follows. We consider the calculation of the Rosenbluth factor of segment i. We generate a set of k - 1 orientations with a distribution prescribed by the bonded interactions (25). These orientations, together with the actual bond between segment - 1 and i, form the set of k orientations (bo, b'*). These orientations are used to calculate the external Rosenbluth factor: k

exp[-flu~Xt(bj)] 9

wext(O) -- E

(31)

j=l

4. For the entire chain the Rosenbluth factor of the old conformation is defined by g

w~

- II w Xt(o) 9

(32)

'/-1

After the new configuration has been generated and the Rosenbluth factor of the old configuration has been calculated the move is accepted with a probability acc(o --+ n) = rain[l,

WeXt(n)/WeXt(o)].

In [98] it is demonstrated that this algorithm samples the desired distribution. 4.3. A p p l i c a t i o n s Here the use of the CBMC technique is illustrated with some examples.

(33)

455

4.3.1. Simulating the energetics and siting In the introduction we have used the adsorption of molecules as an example to illustrate the type of problem one can encounter in simulating systems that exhibit slow diffusion. Smit and Siepmann have used the CBMC technique to study the energetics and siting of alkanes in silicalite [146,147]. In these simulations alkane molecules are modelled with a united atom model, i.e., a CH3 and a CH2 group are considered as a single interaction centre and the zeolite is modelled as a rigid crystal. The zeolite-alkane interactions are assumed to be dominated by dispersive interaction with the oxygen atoms, which are described with a Lennard-Jones potential. A closely related technique was used by Maginn et al. [148]. Figure 19 shows that the simulations of Smit and Siepmann and Maginn et al. predict heats of adsorption of the longer chain alkanes in silicalite that are in good agreement with the experimental data and the simulations of Bigot and Peuch [149].

,, '

i

,

i

,

i

,

i

,

i

0

i

O

lOO

O

0

E O

~ 50 !

9 Experiments 0 sim. (Smit and Siepmann) [] sim. (Maginn et al.) ,

0

i

2

i

,

4

I

,

6

I

8

,

I

10

,

I

12

No

Figure 19. The heats of adsorption of the alkanes in silicalite as a function of the total number of carbon atoms. The closed symbols are experimental and the open symbols are simulation data from Maginn et al. [148] and Smit and Siepmann [147,146].

These simulations give some interesting insights in the distribution of the alkanes over the channels. In Fig. 20 the distribution of butane and dodecane is shown graphically. It is interesting to compare the probability distribution of butane with the one of dodecane. Whereas the plots of butane show an equal density of points in the straight and zig-zag channels, dodecane has a significantly lower probability to be in the zig-zag channels. In a sense, silicalite becomes more and more "unidimensionar' with increasing length of the hydrocarbon. Bates et al. [96] performed CBMC simulations to investigate the energetics of n-alkanes in various zeolites. In particular, Bates et a l . . studied the heat of adsorption of n-alkanes as a function of the pore diameter. The results of these simulations (Figure 21) indicate an optimum sorption of n-alkanes in these all-silica zeolites occurs for zeolites with an average

456 ,~'.';.,,;,~%..

~'~ ,~,~.

,~..~: ~,,~'~

~..:..~...,.-~!

9,,~..

~.~0~.,

-a~,

e~

Figure 20. Distribution of butane (top) and dodecane (bottom) over the various channels of silicalite. At regular intervals during the simulations, a dot representing the position of the alkane molecule is plotted. The density of the dots is a measure of the probability of finding an alkane in a particular section of the zeolite. The left figures are projections on the side plane and the right figures projections on the top plane.

457 pore diameter between 4 and 5 A. Details on how the zeolite influences the conformation of these hydrocarbons in various zeolites can be found in [150]. Other applications of the CBMC technique can be found in [151].

180

160

140

120

~

1O0

~"

~'

C10

8o 60

C8 C6

40

C4

20

0

, 3

i

,

4

= 5

average

,

i

,

i

6

pore diameter

7

, 8

[A]

Figure 21. Heats of adsorption (Q) of n-alkanes in zeolites (RHO, LTA, FER, MFI, MOR, FAU) as a function of mean pore diameter (A spline is fitted to the data as a guide to the eye).

4.3.2. S i m u l a t i n g a d s o r p t i o n i s o t h e r m s Most Molecular Simulations are performed either in the micro-canonical ensemble or in the canonical ensemble. In the micro-canonical ensemble the thermodynamic variables that are kept constant are the temperature, volume, and total energy, whereas in the canonical ensemble the temperature, volume and number of particles are fixed [98,97]. To determine an adsorption isotherm we have to know the number of particles inside the zeolite as a function of the pressure and temperature of the reservoir which is in contact with the zeolite. A naive, but theoretically valid approach would be to use the Molecular Dynamics technique (micro-canonical ensemble) and simulate the experimental situation: an adsorbent in contact with a gas. Such a simulation is only possible for very simple systems. In real life experiments equilibration may take minutes or even several hours,

458 depending on the type of gas molecules. These equilibration times would be reflected in a Molecular Dynamics simulation, the difference being that a minute of "experimental" time takes of the order of 10 9 seconds on a computer. In most cases we are not interested in the properties of the gas phase, yet a significant amount of cpu-time will be spent on the simulation of this phase. Furthermore, in such a simulation there is an interface between the gas phase and the zeolite. In this interfacial region the properties of the system are different from the bulk properties in which we are interested. Since in a simulation the systems are small compared to experiments --hence the interfacial region is a relatively large part of such a s y s t e m - - we have to simulate an unnecessary large system to minimize the influence of this interracial region. Most of the above-mentioned problems can be solved by a careful choice of ensembles. For adsorption studies a natural ensemble to use is the grand-canonical ensemble (or #, V, T-ensemble). In this ensemble the temperature, volume, and chemical potential are fixed. In the experimental setup the adsorbed gas is in equilibrium with the gas in the reservoir. The equilibrium conditions are that the temperature and chemical potential of the gas inside and outside the zeolite must be equal. The gas that is in contact with the adsorbent can be considered as a reservoir that imposes a temperature and chemical potential on the adsorbed gas. We therefore have to know only the temperature and chemical potential of this reservoir to determine the equilibrium concentration inside the adsorbent. This is exactly what is mimicked in the grand-canonical ensemble; the temperature and chemical potential are imposed and the number of particles is allowed to fluctuate during the simulation. This makes these simulations different from the conventional ensembles where the number of molecules is fixed. The grand-canonical Monte Carlo method works best if the acceptance of trial moves by which particles are added or removed is not too low. For atomic fluids this condition effectively limits the maximum'density at which the method can be used to about twice the critical density. Special tricks are needed to extend the grand-canonical Monte Carlo method to somewhat higher densities [98]. The grand-canonical Monte Carlo technique is easily implemented for mixtures and inhomogeneous systems, such as fluids near interfaces. In fact, some of the most useful applications of the grand-canonical Monte Carlo method are precisely in these areas of research. Although the grand-canonical Monte Carlo method can be applied to simple models of non-spherical molecules, special techniques are required since the method converges very poorly for all but the smallest poly-atomic molecules. For chain molecules the grand-canonical ensemble can be combined with the CBMC technique. Techniqual details are discussed in ref. [152], here we discuss an application of this combination. Adsorption isotherms are of practical importance since they give information on the number of molecules adsorbed in the pores of a zeolite as a function of the pressure of the reservoir. Adsorption isotherms are also of flmdamental interest since they may signal phase transitions, such as capillary condensation or wetting, of the fluid inside the pores [153]. For example, if a system exhibits capillary condensation, one would measure a stepped adsorption isotherm with hysteresis. Steps or kinks without hysteresis are occasionally observed on flat substrates [154]. Since the pores of most zeolites are of molecular dimensions, adsorbed alkane molecules behave like a (pseudo) one-dimensional fluid. In a one-dimensional system phase transitions do not occur and therefore one would t

459 expect that for alkanes the adsorption isotherms are of type I, i.e., do not show kinks or steps. If steps occur they are usually attributed to capillary condensation in the exterior secondary pore system formed by the space between different crystals [155]. For silicalite, adsorption isotherms have been determined for various n-alkanes: for the short chain alkanes (methane--pentane) the isotherms are indeed of type I, also for decane a type I isotherm is observed. For hexane and heptane, however, a kink or step is observed.

. 5

,

,

,,i

,

.

,,i

,

,

,,i

. . . .

i

. . . .

i

. . . .

i

,

n ,

,,

'

'"1

. . . .

!

,

'"1

,

,,,i

. . . .

I

'

,,,I

. . . .

I

'

,,,

, do ~

,.o ,.-.-,

o

Z

9

E E

oA."

0. 5

0 t v 0 ~ , ,

0.010 .5

nO D Oo OOO2~

O

n

9 Richard and Rees

9 Oubinin et al.

[] [] i

9 Stach et al. 0 si-m-ulati-~ .....

,,, I

,

1 0 .3

A~

o.5

,

,,I

. . . .

I

10 1

P [kPa]

,

,

,ll

I

I

III

1 01

I

I

"1

I

~176

,

9 Rakhmatkariev et al. simulations , ,,I

,

.... lo ~''

,I

,

,,I

,

,,i

io ~'

"1'o ~ .... '1'o~

P [kPa]

Figure 22. Adsorption isotherms of hexane (left), and heptane (right) in silicalite, the closed symbols are experimental data and the open symbols the results from simulations at T = 298 [K].

Adsorption isotherms are conveniently determined from simulations in the grand-canonical ensemble. In this ensemble the temperature and chemical potential are imposed, but the number of particles is allowed to fluctuate. Adsorption isotherms of alkanes in silicalite have been simulated by Smit and Maesen [156]. The simulated adsorption isotherms for hexane and heptane are shown in Figure 22. The agreement of the simulation results with the experimental data is good at high pressures, but at low pressures deviations exist which indicate that the zeolite-alkane model may be further improved. It is interesting to note that for heptane both the experiments and the simulations show a step at approximately half the loading. Also for hexane detailed inspection of the calculated adsorption isotherm shows a kink at this loading. Since the simulations are performed on a perfect single crystal, these deviations from the type I isotherm must be due to a transition of the fluid inside the pores and can not be attributed to the secondary pore system. Smit and Maesen attribute this transition to a commensurate "freezing" in the channels of a zeolite. The length of a hexane molecule is of the order of the length of the period of the zig-zag channel. At low chemical potential, the hexane molecules move "freely" in these channels and the molecules will be part of their time in the intersections. If part of the intersection is occupied, other molecules can not reside in the straight channels at

460

the same time. At high pressures, almost all hexane molecules are exactly fitting into the zig-zag channel, they do no longer move freely and keep their nose and tail out of the intersection. In such a configuration the entire straight channel can now be tightly packed with hexane molecules (see Figure 23). This may explain the plateau in the adsorption isotherm; in order to fill the entire zeolite structure neatly, the hexane molecules located in zig-zag channels have first to be "frozen" in these channels. This "freezing" of the positions of the hexane molecules implies a loss of entropy and will therefore only occur if the pressure (or chemical potential) is sufficiently high to compensate for this loss. Further experimental evidence for this unexpected behavior of hexane and heptane can been found in refs. [157-159].

Figure 23. Hexane in silicalite: the left figure is at approximately half the maximum loading and the right figure is at almost maximum loading. The figures are projections, the zig-zag channels are in the plane of the figure and straight channels are perpendicular to the plane of the figure. The hexane molecules are represented by gray spheres and red/orange lines represent the zeolite framework.

5. C o n c l u d i n g r e m a r k s Current Density Functional Theory based cluster calculations enable predictions of transition-states and their energies of quite complicated molecules. This is a useful tool to probe acid catalyzed reaction-mechanism in zeolites. There is a need to extend this approach with embedding techniques. One expects this to become feasible in the near future. In studies on absorption and diffusion Molecular Dynamics simulations are limited to rather small systems. However, MD simulations are valuable to test other computational approaches and is the method of choice to get a microscopic picture of diffusion mechanisms. One may expect that MD simulations will replace energy minimization more and more as the tool in studies on zeolite structures and dynamics. The increase in computer power will certainly enhance the number of studies in which MD is combined with ab-initio methods and embedding.

461 We have discussed novel algorithms to make Monte Carlo simulations of chain molecules adsorbed in the pores of zeolites more efficient. In particular, it is demonstrated that the configurational-bias Monte Carlo can be used to simulate the adsorption of long chain hydrocarbons in zeolites. For these type of applications, this technique is many orders of magnitude more efficient than the conventional Monte Carlo techniques or Molecular Dynamics. An important limitation of the Monte Carlo technique is that information on transport properties can not be obtained directly. However, in some cases one can get some information on diffusion if this is an activated process, where one needs to know the free energy barrier. This free energy barrier can be computed with Monte Carlo techniques. Some preliminary work in this direction can be found in [160]. REFERENCES

.

3. 4.

10. 11. 12. 13. 14.

15. 16. 17. 18. 19. 20. 21.

R.A. van Santen, D.P. de Bruyn, C.J.J. den Ouden, and B. Smit, in Introduction to Zeolite Science and Practice, Studies in Surface Science and Catalysis, edited by H. van Bekkum, E.M. Flanigen, and J.C. Jansen (Elsevier, Amsterdam, 1991). J. Sauer, P. Ugliengo, E. Garrone, and V.R. Saunders, Chem. Rev. 94, 2095 (1994). R.A. van Santen and G.J. Kramer, Chem. Rev. 95, 637 (1995). R.A. van Santen and J.W. Niemantsverdriet, Chemical Kinetics and Catalysis (Plenum, New York, 1995). F. Haasse and J. Sauer, J. Am. Chem. Soc. 117, 3780 (1995). H. Koller, E.L. Meijer, and R.A. van Santen, 1997. E.L. Meijer, R.A. van Santen, and A.P.J. Jansen, J. Phys. Chem. 100, 9282 (1996). R.A. van Santen and A.J.M. de Man, in Molecular Modeling, edited by M.A.C. Nascimento (World Scientific, Singapore, 1994), pp. 1-43. R.A. van Santen, A.J.M. de Man, W.P.J.H. Jacobs, E.H. Teunissen, and G.J. Kramer, Catal. Lett. 9, 273 (1991). R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, and J.D. Joannapoulos, Rev. Mod. Phys. 64, 1045 (1992). L.J. Clarke, I. Stich, and M.C. Payne, Comp. Phys. Chem. 72, 14 (1992). R. Shah, J.D. Gale, and M.C. McPayne, J. Phys. Chem. 100, 11688 (1996). M.J.S. Dewar, E.G. Zoebisch, E.F. Healy, and J.J.P. Stewart, J. Am. Chem. Soc. 107, 3902 (1985). M.J.S. Dewar and E.G. Zoebisch, THEOCHEM 49, 1 (1988). W.J. Hehre, L. Radom, P.V.R. Schleyer, and Pople. J.A., Ab-initio molecular orbital therory (Wiley, New York, 1986). Gaussian 94, Gaussian Inc. Pittsburg PA, 1993. M.F. Guest and J. Kendrick, 1986, gamess Users manual, CCP1/86/1; SERC Daresbury Laboratory. R. Doven and V.J. Saunders, 1992, cRYSTAL Q2 Users manual, Gruppo di Chimica Tecnica. Universita di Torino and Daresbury, SERC Laboratory. C. Pisani, R. Dovesi, and C. Rotti, Ab-initio treatment of cystalline systems (Springer, Berlin, 1988). A. Szabo and Oslund, Modern Quantum Chemistry (MacMillan, , 1982).

462 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40.

41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54.

C. Moller and M.S. Plesset, Phys. Rev. 46, 618 (1934). DGauss, Unichem Chemistry Codes, Cray Research Inc. Eagan, MN. DMOL User Guide, version 3.0.0, Molecular Simulations, Inc. San Diego, 1995. Plane Wave User Guide, version 3.0.0, Molecular Simulations, Inc. San Diego, 1995. G. te Velde, 1993, aDF User's Guide, Vrije Universiteit Amsterdam. E.J. Baerends, D.E. Ellis, and P. Ros, Chem. Phys. 2, 41 (1973). E.J. Baerends and P. Ros, Chem. Phys. 2, 52 (1973). E.J. Baerends and P. Ros, Chem. Phys. 8, 412 (1975). T. Ziegler, Chem. Rev. 91,659 (1991). P.E. BlSchl, Phys. Rev. B 50, 17953 (1994). M. Kossner and J. Sauer, J. Phys. Chem. I00, 6199 (1996). E.H. Teunissen, A.P.J. Jansen, R.A. van Santen, and R. Orlando, J. Phys. Chem. 101, 5865 (1996). S. Nawaz, S. Kolboe, K.-P. Lillerud, M. StScker, and H.M. Oren, Stud. Surf. Sci. Catal. 61,421 (1991). C.T.-W. Chu and C.D. Chang, J. Catal. 86, 297 (1984). R.A. van Santen, J. Mol. Catal. A 107, 5 (1996). H.V. Brand, L.A. Curtiss, and L.E. Iton, J. Phys. Chem. 97, 12773 (1993). M. Hunger, M.W. Anderson, A. Ojo, and H. Pfeifer, Microporous Mat. 1, 17 (1993). D. Bartomeuf, Catalysis 1987 (Elsevier, Amsterdam, 1988), p. 177. R.A. van Santen, B.W.H. van Beest, and A.J.M. de Man, in Guidelines for mastering the properties of molecular sieves, Vol. 221 of NATO ASI Series B, edited by D. Barthomeuf, E.G. Derouane, and W. HSlderich (Plenum, New York, 1990), pp. 201-224. N. Cardona-Martinez and J.A. Dumesic, Adv. Catal. 38, 149 (1992). U. Lohse and J~inchen, in Elementary reaction steps in heterogeneous catalysis, edited by R.W. Joyner and R.A. van Santen (Kluwer, Amsterdam, 1993), p. 113. G.J. Kramer, A.J.M. de Man, and R.A. van Santen, J. Am. Chem. Soc. 113, 6435 (1991). V.B. Kazansky, M.V. Frash, and R.A. van Santen, Stud. Surf. Sci. and Catal. 105, 2283 (1997). V.B. Kazansky and I.N. Senchenya, J. Catal. 109, 108 (1989). E.A. Paukshtis and E.N. Yurchenko, Russ. Chem. Rev. 52, 426 (1983). G.J. Kramer and R.A. van Santen, J. Am. Chem. Soc. 115, 2887 (1993). J.-R. Hill and J. Sauer, J. Phys. Chem. 97, 6579 (1993). L. Campana, A. Selloni, J. Weber, A. Pasquarello, I. Papai, and A. Goursot, Chem. Phys. Lett. 226, 245 (1994). P. Sarv, T. Tuherm, E. Lippmaa, K. Keskinen, and A. Root, J. Phys. Chem. 99, 13763 (1995). K.M. Neyman, P. Strodel, S.P. Ruzankin, N. Schlensog, H. KnSzinger, and W. RSsch, Catal. Lett. 31,273 (1995). A.G. Pelmenschikov, R.A. van Santen, J. J~inchen, and E.L. Meijer, J. Phys. Chem. 97, 11071 (1993). J.D. Gale, C.R.A. Carlow, and J.R. Carruthers, Chem. Phys. Lett. 216, 155 (1993). A.G. Pelmenschikov and R.A. van Santen, J. Phys. Chem. 97, 10678 (1993).

463 55. E.H. Teunissen, F.B. van Duyneveldt, and R.A. van Santen, J. Phys. Chem. 96, 366 (1992). 56. E. Kassab, E. Gerti, and M. Allavena, J. Phys. Chem. 95, 9425 (1991). 57. E. Kassab, J. Fongreet, M. Allavena, and E.M. Evleth, J. Phys. Chem. 97, 9034 (1993). 58. J. Sauer and F. Haase, J. Phys. Chem. 98, 3083 (1994). 59. A.T. Bel| and J.N. Theodorou, J. Phys. Chem. 99, 1505 (1995). 60. H.V. Brand, L.A. Curtiss, and L.E. Iton, J. Phys. Chem. 96, 7725 (1992). 61. S.J. Cook, A.K. Chakraborty, A.T. Bell, and D.N. Theodorou, J. Phys. Chem. 97, 6679 (1993). 62. S.P. Greatbanks, P. Sherwood, I.H. Hillier, R.J. Hall, N.A. Burton, and I.R. Gould, Chem. Phys. Lett. 234, 367 (1995). 63. J.B. Uytterhoeven, L.G. Christner, and W.K. Hall, J. Phys. Chem. 69, 2117 (1965). 64. S.A. Zygmunt, L.A. Curtiss, L.E. Iton, and M.K. Erhardt, J. Phys. Chem. 100, 6663 (1996). 65. S.P. Greatbanks, I.H. Hillier, and P. Sherwood, J. Phys. Chem. 105, 3770 (1996). 66. S.R. Blaszkowski and R.A. van Santen, J. Phys. Chem. 101, 2292 (1997). 67. R. Shah, M.C. Payne, M.-H. Lee, and J.D. Gale, Science 271, 1395 (1996). 68. E. Nusterer, P.E. B15chl, and N. Schwarz, Ang. Chem. 108, 187 (1996). 69. E. Nusterer, P.E. B15chl, and N. Schwarz, Chem. Phys. Lett. 253, 448 (1996). 70. L. Smith, A.K. Cheetham, R.E. Monis, J.M. Thomas, and P.A. Wright, J. Chem. Science 27, 799 (1988). 71. E.H. Teunissen, R.A. van Santen, A.P.J. Jansen, and F.B. van Duyneveldt, J. Phys. Chem. 97, 1993 (1993). 72. H. KnSzinger, H. Biihl, and K. Kochloefl, J. Cata]. 24, 57 (1972). 73. G. Mirth and J. Lercher, in Natural Gas Conversion, edited by A. Holmen e.a. (Elsevier, Amsterdam, 1991), p. 437. 74. G. Mirth, J. Lercher, M.W. Anderson, and J. Klinowski, J. Chem. Soc. Faraday Trans. 86, 3039 (1990). 75. F. Wakabayaski, J. Kondo, K. Domen, and C. Hirose, Catal. Lett. 21,257 (1993). 76. H. Jobic, A. Tuel, M. Krossner, and J. Sauer, J. Phys. Chem. 100, 19545 (1996). 77. J. Florian and L Kubelkova, J. Phys. Chem. 98, 8734 (1994). 78. D.J. Parillo, R.J. Gorte, and W.E. Farneth, J. Am. Chem. Soc. 115, 12441 (1993). 79. L Kubelkova, J. Kortla, and J. Florian, J. Phys. Chem. 99, 10285 (1995). 80. T. Fujino, M. Kashitani, J.N. Kondo, K. Domen, C. Hirose, F. Ishida, F. Goto, and F. Wakabayaski, J. Phys. Chem. 100, 11649 (1996). 81. G.J. Kramer, R.A. van Santen, C.A. Emeis, and A.K Nowak, Nature 363,529 (1993). 82. G.J. Kramer and R.A. van Santen, J. Am. Chem. Soc. 117, 1766 (1995). 83. S.R. Blaszkowski, M.C.A. Nascimento, and R.A. van Santen, J. Phys. Chem. 100, 3463 (1996). 84. E.M. Evleth, E. Kassab, and L.R. Sierra, J. Phys. Chem. 98, 1421 (1994). 85. V.B. Kazansky, Kinet. and Cata]. 21, 159 (1980). 86. I.N. Senchenya and V.B. Kazansky, Catal. Lett. 8, 317 (1991). 87. A.G. Pelmenschikov, W.U. Zhanpeisov, E.A. Paukshtis, L.V. Malyssheva, G.M. Zhidomirov, and K.I. Zamaraev, Dokl. Akad. Nauk. SSSR 293, 915 (1987).

464 88. P.M. Viruella, C.M. Zicovich-Wilson, and A. Corma, J. Phys. Chem. 97, 13713 (1993). 89. A. Corma, P.M. Viruella, and C.M. Zicovich-Wilson, J. Phys. Org. Chem. 7, 364 (1994). 90. V.B. Kazansky, Acc. Chem. Res. 24, 379 (1991). 91. N.W. Cant and W.K. Hall, J. Catal. 52, 161 (1972). 92. J.F. Haw, B.R. Richardson, I.S. Ohiro, W.O. Lazo, and A.J. Spled, J. Am. Chem. Soc. 111,2052 (1989). 93. J.F. Haw, J.B. Nicholas, T. Xu, L.W. Beck, and D.F. Ferguson, Acc. Chem. Res. 29, 259 (1996). 94. T. Xu and J.F. Haw, J. Am. Chem. Soc. 116, 7753 (1994). 95. I. R. Forester and R.F. Howe, J. Am. Chem. Soc. 109, 5076 (1987). 96. S.P. Bates, W.J.M. van Well, R.A. van Santen, and B. Smit, J. Am. Chem. Soc. 118, 6753 (199S). 97. M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987). 98. D. Frenkel and B. Smit, Understanding Molecular Simulations: from Algorithms to Applications (Academic Press, San Diego, 1996). 99. L.M. Bull, N.J. Henson, A.K. Cheetham, and S.J. Heyes, J. Phys. Chem. 97, 11776 (1993). 100.P. Demontis and G.B. Suffritti, in Proceedings of the 9th International Zeolite Conferenence, Montreal, edited by R. von Ballmoos, J.B. Higgins, and M.M.J. Treacy (Butterworth-Heinemann, Boston, 1992), p. 345. 101.E. Hern&ndez and C.R.A. Catlow, Proc. R. Soc. Lond. A 448, 143 (1995). 102.T. Mosell, G. Schrimpf, C. Hahn, and J. Brickmann, J. Phys. Chem. 100, 4571 (1996). 103.T. Mosell, G. Schrimpf, and J. Brickmann, J. Phys. Chem. 100, 4582 (1996). 104.E.J. Maginn, A.T. Bell, and D.N. Theodorou, J. Phys. Chem. 100, 7155 (1996). 105.G. Schrimpf and J. Brickmann, J. Comput.-Aided Mater. Des. 2, 49 (1995). 106.R.L. June, A.T. Bell, and D.N. Theodorou, J. Phys. Chem. 94, 8232 (1990). 107.J.R. Hufton, J. Phys. Chem. 95, 8836 (1991). 108.S.J. Goodbody, K. Watanabe, D. MacGowan, J.P.R.B. Walton, and N. Quirke, J. Chem. Soc., Faraday Trans. 87, 1951 (1991). 109.K.S. Smirnov, Chem. Phys. Lett. 229, 250 (1994). 110.D. Dumont and D. Bougeard, Zeolites 15, 650 (1995). 111.K.S. Smirnov and B. van de Graaf, J. Chem. Soc. Faraday Trans. 92, 2475 (1996). 112.D. Shen and L.V. Rees, J. Chem. Soc. Faraday Trans. 92, 487 (1996). 113.D. Keffer, A.V. McCormick, and H.T. Gavis, Mol. Phys. 87, 367 (1996). 114.K. Hahn, J. Kaerger, and V. Kukla, Phys. Rev. Lett. 76, 2762 (1996). 115.K.S. Smirnov and B. van de Graaf, J. Chem. Soc. Faraday Trans. 92, 2469 (1996). 116.S. Yashonath and P. Santikary, J. Phys. Chem. 97, 13778 (1993). 117.S. Yashonath and P. Santikary, J. Phys. Chem. 98, 6368 (1994). 118.S. Yashonath and S. Bandyopadhyay, Chem. Phys. Lett. 228, 284 (1994). 119.S. Yashonath and P. Santikary, J. Chem. Phys. 100, 4013 (1994). 120.S. Bandyopadhyay and S. Yashonath, Chem. Phys. Lett. 223, 363 (1994). 121.P. Santikary and S. Yashonath, J. Phys. Chem. 98, 9252 (1994). 122.S. Bandyopadhyay and S. Yashonath, J. Phys. Chem. 99, 4286 (1995).

465 123.J. Wei, Ind. Eng. Chem. Res. 33, 2467 (1994). 124P~.C. Runnebaum and E.J. Maginn, J. Phys. Chem. B 101, 6394 (1997). 125.13. B%mard and D. Bougeard, Adv. Mater. 1, 10 (1995). 126.YK. Yamahara, K. Okazaki, and K. Kawamura, Catal. Today 23, 397 (1995). 127.Y. Oumi, K. Matsuba, M. Kubo, T. Inui, and Miyamoto A., Microporous Mater. 4, 53 (1995). 128~(.S. Smirnov and B. van de Graaf, Microporous Mater 7, 133 (1996). 129K.S. Smirnov and D. Bougeard, J. Mol. Struc. 348, 155 (1995). 130.H.J.F. Stroud, E. Richards, P. Limcharoen, and N.G. Parsonage, J. Chem. Soc., Faraday Trans. I 72, 942 (1976). 131.S. Yashonath, J.M. Thomas, A.K. Nowak, and A.K. Cheetham, Nature 331, 601 (1988). 132.B. Smit and C.J.J. den Ouden, J. Phys. Chem. 92, 7169 (1988). 133.S.D. Pickett, A.K. Nowak, J.M. Thomas, B.K. Peterson, J.F. Swift, A.K. Cheetham, C.J.J. den Ouden, B. Smit, and M. Post, J. Phys. Chem. 94, 1233 (1990). 134A.K. Nowak, C.J.J. den Ouden, S.D. Pickett, B. Smit, A.K. Cheetham, M.F.M. Post, and J.M. Thomas, J. Phys. Chem. 95, 848 (1991). 135.C.R.A Catlow editor, Modelling of Structure and Reactivity in Zeolites (Academic Press, London, 1992). 136P~.L. June, A.T. Bell, and D.N. Theodorou, J. Phys. Chem. 96, 1051 (1992). 137.N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.N. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953). 138.M. Fixman, Proc. Nat. Acad. Sci. USA 71, 3050 (1974). 139.G. Ciccotti and J.P. Ryckaert, Comp. Phys. Rep. 4, 345 (1986). 140.J.I. Siepmann and D. Frenkel, Mol. Phys. 75, 59 (1992). 141.M.N. Rosenbluth and A.W. Rosenbluth, J. Chem. Phys. 23, 356 (1955). 142.J. Harris and S.A. Rice, J. Chem. Phys. 88, 1298 (1988). 143.D. Frenkel, G.C.A.M. Mooij, and B. Smit, J. Phys.: Condens. Matter 4, 3053 (1992). 144.J.J. de Pablo, M. Laso, and U.W. Suter, J. Chem. Phys. 96, 2395 (1992). 145.J. Batoulis and K. Kremer, J. Phys. A: Math. Gen. 21, 127 (1988). 146.13. Smit and J.I. Siepmann, J. Phys. Chem. 98, 8442 (1994). 147.13. Smit and J.I. Siepmann, Science 264, 1118 (1994). 148~E.J. Maginn, A.T. Bell, and D.N. Theodorou, J. Phys. Chem. 99, 2057 (1995). 149.B. Bigot and V.-H. Peuch, J. Phys. Chem. 99, 8206 (1995). 150.S.P. Bates, W.J.M. van Well, R.A. van Santen, and B. Smit, J. Phys. Chem. 100, 17573 (1996). 151.S. Bandyopadhyay and S. Yashonath, J. Phys. Chem. B 101, 5675 (1997). 152.13. Smit, Mol. Phys. 85, 153 (1995). 153.R. Evans, in Liquides aux Interfaces/Liquids at interfaces, Les Houches, Session XLVIII, 1988, edited by J. Charvolin, J.F. Joanny, and J. Zinn-Justin (North Holland, Amsterdam, 1990), pp. 1-98. 154.S.J. Gregg and K.S.W. Sing, Adsorption, Surface Area and Porosity (Academic Press, London, 1982). 155.H. Stach, U. Lohse, H. Thamm, and W. Schirmer, Zeolites 6, 74 (1986). 156.B. Smit and T.L.M. Maesen, Nature 374, 42 (1995).

466

157.W.J.M. van Well, J.P. Wolthuizen, B. Smit, J.H.C. van Hooff, and R.A. van Santen, Angew. Chem. (Int. Ed.) 34, 2543 (1995). 158.D.H. Olsen and P.T. Reischmann, Zeolites 17, 434 (1996). 159.M.S. Sun, O. Talu, and D.B. Shah, J. Phys. Chem. 100, 17276 (1996). 160.B. Smit, L.D.J.C. Loyens, and G.L.M.M. Verbist, Faraday Discuss. 106, 93 (1997).