4768
Langmuir 1996, 12, 4768-4783
Grand Canonical Monte Carlo Simulations of Adsorption of Polar and Nonpolar Molecules in NaY Zeolite Roland J-M Pellenq,†,‡ Bernard Tavitian,§ Didier Espinat,§ and Alain H. Fuchs*,† Laboratoire de Chimie-Physique des Mate´ riaux Amorphes, URA-CNRS 1104, Universite´ Paris-Sud, 91405 Orsay ce´ dex, France, and Departement de Physico-Chimie Applique´ e et Analyses, Institut Franc¸ ais du Pe´ trole, 1-4 avenue de Bois-Pre´ au, BP 311, 92506 Rueil Malmaison, ce´ dex, France Received January 10, 1996. In Final Form: June 11, 1996X In this work, we report grand canonical Monte Carlo (GCMC) simulations of adsorption of polar and nonpolar molecules such as xenon, methane, and p- and m-xylene in the NaY zeolite at different temperatures. The adsorption isotherms as well as the curves giving the evolution of the isosteric heat with coverage have been calculated and compared to their experimental counterparts and to other simulation results when available. This work is based on a new adsorbate-zeolite potential function. In the particular case of adsorption of xylene isomers, we also report a detailed structural analysis of the sorbed phases.
I. Introduction Adsorption phenomenon in zeolitic microporous materials is of great scientific interest due to applications in separation technology and catalysis, but comprehensive theory on zeolite adsorption still remains unavailable. In the last few years, a large number of computer simulation studies on the physisorption of molecules in zeolites has been reported in the literature (see section II). These works have revealed the molecular simulation technique as a powerful and complementary tool to the usual experimental characterization methods. Among all natural and synthetic zeolitic structures reported,1 few of them are used in industrial processes. In particular, zeolite Y is important for isomeric separation of xylene isomers. The molecular mechanism which underlines the separation of molecules such as m- and p-xylene in this faujasite zeolite is still not known. As a first step toward the direct simulation of mixtures (such as p- and m-xylene in NaY zeolite), it is possible to use phenomenological models such as the ideal adsorption solution theory (IAS) or the real adsorption solution theory (RAS)2 to calculate the selectivity for a given mixture composition at a given temperature. In both cases, single species adsorption isotherms are needed which can be obtained from grand canonical Monte Carlo simulation. This technique as all other simulation techniques relies on an accurate description of the intermolecular interactions between the adsorbate and the zeolitic framework and between the adsorbed molecules themselves. In this study, we have selected four different adsorbates: xenon, methane, and m- and p-xylene isomers in order to test the transferability of a new adsorbate-zeolite potential function recently introduced in GCMC simulation of adsorption of simple molecules in silicalite-13 and in AlPO4-5.4 It is indeed * To whom correspondence should be sent: tel, +33-1.69.41.75.84; fax, +33-1.60.19.33.02; e-mail,
[email protected]. † Universite ´ Paris-Sud. ‡ Present address: Centre de Recherches sur la Matie ` re Divise´e, UMR-CNRS 131, 1b rue de la Fe´rollerie, 45071 Orle´ans, ce´dex 02, France. § Institut Franc ¸ ais du Pe´trole. X Abstract published in Advance ACS Abstracts, August 15, 1996. (1) Meiers, W. M.; Olson, D. H. In Atlas of zeolite structures; Structure Commission of the International Zeolite Association, 1978. (2) Meyers, A. L.; Praunitz, J. M. AIChE J. 1978, 11, 121. (3) Pellenq, R. J.-M.; Nicholson, D. Langmuir 1995, 11, 1626. (4) Boutin, A.; Pellenq, R. J.-M.; Nicholson, D. Chem. Phys. Lett. 1994, 219, 484. Lachet, V.; Boutin, A.; Pellenq, R. J.-M.; Nicholson, D.; Fuchs, A. H. J. Phys. Chem. 1996, 100, 9006.
S0743-7463(96)00035-2 CCC: $12.00
particularly interesting to assess the relative importance of the different adsorbate/zeolite potential terms (see section IV) in the case of such a chemically complex zeolitic structure which contains four different framework species (O, Si, Al, Na+). The remainder of the paper is organized as follows. Previous simulation works on adsorption in zeolite NaY are surveyed in section 2. Section 3 is devoted to the description of the NaY zeolite structure. In particular, we focus on the aluminum distribution in the zeolite framework in relation with the countercations distribution in the zeolitic cavities. Section 4 is devoted to the different potential functions used in this work to model the adsorbent-adsorbate and adsorbate-adsorbate interactions. We outline in section 5, the grand canonical Monte Carlo simulation technique applied to adsorption. Simulated isotherms and isosteric heat curves are compared to experiment and to other simulation results for each different adsorbate in section 6. We also give a detailed analysis of the structure of the adsorbed phase in the particular case of m- and p-xylene. II. Survey of Previous Works Studies of single gas adsorption in zeolite faujasite NaY by canonical Monte Carlo (CMC), grand canonical Monte Carlo (GCMC), or molecular dynamics (MD) already exist in the literature. They mostly concentrate on adsorption of simple adsorbates such as the rare gases or small hydrocarbon molecules. Kiselev et al.5 calculated Henry’s law constants and isosteric heats of adsorption at zero coverage for nonpolar and polar small inorganic molecules and hydrocarbons adsorbed in zeolite NaX and NaY. They introduced an atom-atom approximation in order to calculate the adsorbate/zeolite potential energy. The short range interaction, i.e., the dispersion and the repulsion interaction, was modeled with a Lennard-Jones function whose parameters are obtained from the dipole polarizability, the magnetic susceptibility, and the van der Waals radius of each interacting species. It is important to note (i) that this potential model neglects interactions between the adsorbate and the so-called T atoms (silicon and aluminum) of the zeolite structure and (ii) that the electronic characteristics (polarizabilities, magnetic susceptibilities, and van der Waals radius) of zeolitic oxygen and sodium atoms are taken equal to those of correspond(5) Kiselev, A.; Pham Quang Du J. Chem. Soc., Faraday Trans. 2 1981, 77, 1.
© 1996 American Chemical Society
Adsorption on Zeolite
ing neutral and isolated species. Kiselev and co-workers also modified the Lennard-Jones function between the adsorbate and the framework oxygen and sodium atoms by introducing a multiplicative empirical parameter in order to bring the calculated heat of adsorption and the Henry constants in agreement with their experimental counterparts. For polar molecules, an allowance was also made for the electrostatic interaction between the adsorbate rigid molecular multipoles with zeolite ions. This contribution was calculated using the point moment approximation considering the adsorbate n-pole(s) at its center of mass. Yashonath et al. have reported the first (canonical) Monte Carlo study of the adsorption of one methane molecule in a NaY supercage6 using the Kiselev potential model based on the atom-atom approximation. Higher loading simulations for adsorption of methane and xenon were also carried out by Yashonath et al. using molecular dynamics.7,8 Wood and Rowlinson presented GCMC simulations of xenon and methane in zeolites NaX and NaY in order to calculate adsorption isotherms and isosteric heat curves. They used two models: a simple one based on a spherically-averaged potential for a single isolated cavity9 and a more detailed one based on the Kiselev potential;10 no long range coulombic interactions (for methane adsorption) or induction interactions were taken into account. More recently, Maddox and Rowlinson have studied the adsorption of nitrogen and methane mixtures in NaY using the GCMC technique.11 It is interesting to mention the attempt of Vigne´-Maeder in deriving NMR chemical shift for the xenon in NaY (and other zeolites) at low coverage.12 We also note the recent work of Schrimpf et al.13 on the adsorption of spherical molecules in which the zeolite framework flexibility is taken into account. Yashonath and co-workers have recently undertaken a systematic molecular dynamics study of the influence of the sorbate size on diffusivity.14 Adsorption of aromatics has also been studied by means of numerical simulations. Some years ago, Demontis et al. reported a MD simulation of benzene in NaY at medium loading.15 These authors used a more detailed adsorbate/ zeolite potential than that proposed by Kiselev and co-workers: the benzene/NaY electrostatic interaction is calculated using point charges placed on zeolite atoms and on each atomic component of the adsorbate molecules. Demontis et al. successfully localized the benzene molecules in the zeolite cavities. Henson et al. recently performed a canonical MC study of adsorption of benzene at low coverage.16 Klein et al. studied the location of (6) Yashonath, S.; Thomas, J. M.; Nowak, A. K.; Cheetham, A. K. Nature 1988, 331, 601. (7) Yashonath, S.; Demontis, P.; Klein, M. J. Phys. Chem. 1991, 95, 5881. (8) Yashonath, S.; Santikary, P. J. Phys. Chem. 1993, 97, 3849 (and references therein). (9) Wood, G. B.; Panagiotopoulos, A. Z.; Rowlinson, J. S. Mol. Phys. 1988, 63, 49. (10) Wood, G. B.; Rowlinson, J. S. J. Chem. Soc., Faraday Trans. 2 1989, 85, 765. (11) Maddox, M. W.; Rowlinson, J. S. J. Chem. Soc., Faraday. Trans. 1993, 89, 3619. (12) Vigne´-Maeder, F. J. Phys. Chem. 1994, 98, 4666. (13) Schrimpf, G.; Schlenkrich, M.; Brickmann, J.; Bopp, P. J. Phys. Chem. 1992, 96, 7404. (14) Yashonath, S.; Santikary, P. Mol. Phys. 1993, 78, 1. Yashonath, S.; Santikary, P. J. Chem. Phys. 1994, 100, 4013. Bandyopadhyay, S.; Yashonath, S. Chem. Phys. Lett. 1994, 223, 363. Yashonath, S.; Bandyopadhyay, S. Chem. Phys. Lett. 1994, 228, 284. Yashonath, S.; Santikary, P. J. Phys. Chem. 1994, 98, 6368. (15) Demontis, P.; Yashonath, S.; Klein, M. J. Chem. Phys. 1989, 93, 5016. (16) Henson, N. J.; Cheetham, A. K.; Redondo, A.; Levine, S. M.; Newsam, J. M. In Zeolites and Related Microporous Materials: State of the Art; Weitkamp, J., Karge, H. G., Pfeifer, H., Holderich, W., Eds.; Studies in Surface Science and Catalysis; Elsevier: Amsterdam, 1994; Vol. 84, p 2059.
Langmuir, Vol. 12, No. 20, 1996 4769
benzene, p-xylene, and m-xylene in NaY using a simple (0 K) energy minimization technique.17 More recently, Klein and co-workers performed low coverage molecular dynamics calculation of different aromatic molecules in NaY.18 Finally, we note the recent MD study of Schrimpf et al. on the adsorption of p-xylene in NaY.19 This work deals with a detailed analysis of the effects of loading and temperature on diffusivity. All these simulation studies reviewed here above are based on a simplistic description of the adsorbate/zeolite short range interaction and on a point charge description to account for electrostatic long range interaction. It has been recently demonstrated20 that such an approach for nonpolar adsorbates is not consistent with the nature of the bonding in zeolitic material. This work is concerned with a realistic description of the interaction between the adsorbate and the zeolitic solid in order to gain detailed insights on the adsorption phenomenon in NaY zeolite at a molecular level. III. The NaY Faujasite Structure A unit cell of faujasite has the composition NaxAlxSi192-xO384 where 0 e x e 96 (1 e Si/Al e ∞) and exhibits a cubic symmetry which has been described within the Fd3m group (a ) 24.8536 Å)21 or within the F23 group (a ) 25.0184 Å).22 This zeolite exhibits a microporous network made of large cavities (or supercages) which interconnect to give a highly microporous structure. These supercages have a spherical symmetry; their diameter is about 12.5 Å. The windows by which they interconnect are about 7.5 Å large. Each cavity is connected to four others in a tetrahedral arrangement. The structure also contains smaller sodalite (or β) cages. These sodalite cages are not accessible to xenon or methane. One unit cell of faujasite zeolite contains eight supercages. The presence of aluminum (trivalent) introduces charge defects in the zeolite framework since the T atoms (Al or Si) are exclusively in tetrahedral sites in such materials. This charge excess has to be compensated with some nonframework counterions (in this study Na+) in order to achieve global electroneutrality. Their number per zeolite unit cell depends on aluminum content and their own electric charge. It should be mentioned that all sodium atoms are fully ionized. The two structures X and Y only differ in their aluminum content; X having a lower Si/Al ratio (hence more sodium cations) than Y. For each value of the Si/Al ratio, there are several ways of distributing the aluminum atoms in the zeolite framework since diffraction studies cannot distinguish between aluminum or silicon sites. It is clear that the distribution of the nonframework sodium cations depends on the aluminum distribution in the zeolitic framework. However, X-ray diffraction studies of dehydrated X and Y zeolites have shown that there are four crystallographic sites for the nonframework cations22 defined as sites SI, SI′, SII, and SIII. Site SI and SII locations are shown in Figure 1: there are 16 SI sites at the center of the hexagonal prisms, 32 SI′ sites in the six-ring windows of the sodalite cages, 32 (17) Klein, H.; Fuess, H. In Zeolites and Related Microporous Materials: State of the Art; Weitkamp, J., Karge, H. G., Pfeifer, H., Holderich, W., Eds.; Studies in Surface Science and Catalysis; Elsevier: Amsterdam, 1994; Vol. 84, p 2067. (18) Klein, H.; Kirschhock, C.; Fuess, H. J. Phys. Chem. 1994, 98, 12345. (19) Schrimpf, G.; Tavitian, B.; Espinat, D. J. Phys. Chem. 1995, 99, 10932. (20) Pellenq, R. J.-M.; Nicholson, D. J. Phys. Chem. 1994, 98, 13339. (21) Fitch, A. N.; Jobic, H.; Renouprez, A. J. Phys. Chem. 1986, 90, 1311. (22) Uytterhoeven, L.; Dompas, D.; Mortier, W. J. J. Chem. Soc., Faraday Trans. 1992, 88, 2753.
4770 Langmuir, Vol. 12, No. 20, 1996
Pellenq et al.
Figure 2. A fragment of the zeolite structure showing the various crystal species (see text).
Figure 1. A supercage of NaY zeolite. Sites I and II for sodium counterions are shown.
SII in the six-ring windows which face into the supercages, and finally 96 SIII sites in these large cavities. As far as NaY is concerned, the experimental value of the Si/Al ratio ranges from 1.5 to 2.5. In order to avoid any problem due to partial occupancy of the sodium sites, it is often chosen in numerical simulations13,19 to consider full occupancy of sites SI and SII. Thus there are 48 sodium cations and the same number of aluminum atoms per unit cell when the Si/Al ratio is equal to 3.0. Unlike sodium cations, the choice of 48 tetrahedral sites for aluminum atoms is not univocal. The distribution adopted in our simulations was the one described by Uytterhoeven et al.22 This regular distribution is based on symmetry considerations on the secondary building units from which one can build the whole unit cell. It is characterized by two types of supercages alternated throughout the crystal; the first type has 12 aluminum atoms and the second none. This constitutes certainly a severe approximation concerning the description of the zeolitic structure since the aluminum density is certainly homogeneous. It is also assumed that the sodium ions remain fixed in the zeolite unit cell as usual in most simulations. The validity of these two assumptions will have to be tested in future works. Siliceous crystals such as quartz or the zeolite silicalite are known as iono-covalent solids: the Si-O bond in SiO2 systems exhibits a strong covalent character; the partial charges (calculated from periodic ab initio methods) are half the formal charges, i.e. -1e and +2e for oxygen and silicon atoms respectively in siliceous structures.23 Bonding in aluminosilicates such as NaY faujasite is obviously affected by the introduction of aluminum atoms in the crystal structure. Periodic ab initio studies have shown24 that the Al-O bond in Al2O3 has a stronger ionic character than the Si-O bond in SiO2 systems. For an aluminosilicate such as NaY zeolite, the partial charge for each framework species can be obtained from the electronegativity equalization method (derived from density functional theory),22 which fully takes into account the crystal periodicity. The averaged charges are -0.823e, +1.452e, and +1.228e for oxygen, silicon, and aluminum atoms, respectively. The charge excess in AlO4 induces an electrostatic field in the crystal which determines the (23) White, J. C.; Hess, A. C. J. Phys. Chem. 1993, 97, 6398. (24) Salasco, L.; Dovesi, R.; Orlando, R.; Causa, M.; Saunders, V. R. Mol. Phys. 1991, 72, 267. Borosy, A. P.; Silvi, B.; Allavena, M.; Nortier, P. J. Phys. Chem. 1994, 98, 13189.
location of the sodium cations at the various sites described earlier. Since these cations are ionized (+1e), they induce in turn an electric field on the framework atoms (Si, Al, O). Silicon and aluminum atoms are less polarizable than oxygen atoms, and one expects the oxygen atoms near to a counterion to be more polarized than other framework oxygens. Framework atoms polarizabilities can be derived from Auger electron spectroscopy (AES) data.25 From AES data, it was also possible to differentiate between different types of oxygen atoms in aluminosilicate zeolites.26 We can distinguish between three different environments for framework oxygens in structures such as NaY (see Figure 2): the types I and II oxygens have a polarizability equal to 1.2 Å3 (8.1a03) while type III oxygens have a polarizability which depends on the nature of the counterion. When the counterion is sodium, the polarizability of the type III oxygens is equal to 1.4 Å3 (9.5a03).26 This value is higher than that for types I and II and is fairly constant with the Si/Al ratio.25 It implies that the presence of the sodium cations is only felt by the nearest oxygen neighbors. Moreover the value of the polarizability of oxygen atoms near to the sodium counterions indicates that there is no bonding between these two species. The polarizabilities of aluminum and silicon are estimated to be 0.3 Å3 (2.0a03) and 0.5 Å3 (3.4a03) respectively.4,20 The polarizability of sodium cations 0.148 Å3 (1.0a03) is obtained from quantum mechanical calculations.27 The electronic properties of framework atoms (partial charges and polarizabilities) are essential quantities to derive the adsorption potential function (see section IV). Let us consider an adsorbate molecule in a supercage of the NaY zeolite. The electrostatic field created by the zeolite atoms in the supercage depends obviously on their distribution throughout the crystalline structure. It is thus essential to know the distribution of each type of oxygen framework species in the zeolite unit cell. The use of the F23 space group defines unambiguously the silicon, the aluminum, and the sodium distributions.22 It is possible to select among all oxygen locations (obtained from X-ray studies) the oxygen atoms near to the counterions (designed here as type III oxygens): a sodium cation in site I has 12 oxygen atoms as immediate neighbors while a sodium cation in site II has only six neighbor oxygen atoms. From the remainder of type III oxygen atoms, it is easy to obtain all possible locations for types I and II oxygen atoms, which have same polarizability and partial charge. (25) Pellenq, R. J.-M.; Nicholson, D. J. Chem. Soc., Faraday Trans. 1993, 89, 2499. (26) Pellenq, R. J.-M.; Nicholson, D. Stud. Surf. Sci. Catal. 1994, 87, 31. (27) Pyper, N. C. Mol. Sim. 1990, 5, 23.
Adsorption on Zeolite
Langmuir, Vol. 12, No. 20, 1996 4771 Table 1
(a) Electronic Characteristics of Zeolite Species (Partial Charges and Polarizabilities) in Atomic Units q(e) R(a03)
OI,II
OIII
Si
Al
Na+
-0.823 8.30
-0.823 9.45
1.452 3.37
1.228 2.03
1.0 1.00
(b) Electronic Characteristics of the Adsorbates Considered in this Work (Partial Charges and Polarizabilities) in Atomic Units q(e) R(a03)
Xe
H(-Csp2)
Csp2
H(-Csp3)
Csp3
CH3(-Csp2)
0.0 27.2
+0.115 4.32
-0.115 12.8
0.0 4.53
0.0 11.9
+0.115 13.5
IV. Potential Functions Used in the GCMC Simulations
(1)
The Coulombic term is calculated classically using the partial charges of each interacting species
Ucoul(rij) )
qiAqj rij
(2)
where the subscript i corresponds to an atom of an adsorbate molecule A and the subscript j corresponds to the framework atoms (O(I), O(II), O(III), Si, Al, Na+). This term is strictly pair additive. The partial charges qj for the framework atoms are those given in section III. The partial charges qi carried by the adsorbate are those reported below in the description of the adsorbateadsorbate intermolecular potential function (see Table 1). Thus, the description of the adsorbate/adsorbent Coulombic interaction is fully consistent with the adsorbateadsorbate potential function. It is obvious that the Coulombic interaction for the Xe/NaY system is zero. The Coulombic interaction between methane octopole and zeolite partial charges was neglected in this work because of the very small magnitude of the methane octopole. The adsorbate/zeolite electrostatic interaction is therefore considered only for xylene isomers. The induction energy is given by the first term of the induction multipole expansion
1 Uind(rij) ) - RiA{(E(rij)}2 2
(3)
According to the reference system given for the Coulombic interaction, RiA is the dipole polarizability of the atom i pertaining to the adsorbed molecule A and E(r) is the electrostatic field at the position occupied by the atom i due to the partial charges carried by the framework species. It can be shown that the induction energy derived from the perturbation theory for intermolecular interaction28 is not pair-additive. Even in polar liquids such as water, this term does contribute very little to the total (28) Buckingham, A. D. Adv. Chem. Phys. 1967, 12, 107.
Udisp(rij) ) -f6
C6ij rij6
- f8
C8ij rij8
- f10
C10ij
- ... + rij10 three body terms (4)
with
IV.1. Adsorbate-Zeolite Potential Function. The adsorbate-zeolite potential function used in this work is based on the usual partition of the total interaction energy derived form the perturbation theory applied to intermolecular interactions. It considers the interactions between each atom of the adsorbate molecule with each zeolite species including the sodium cations. The general expression for the adsorbate-zeolite potential function is given by28
Utotal ) Ucoul + Uind + Udisp + Urep
interaction energy29 so that the expression presented above is sufficient. Moreover, we consider that the presence of the polar adsorbate molecule in the pore does not polarize the framework atoms, so that there is no back-polarization effect. The third term of the general expression of the adsorbate-zeolite potential function corresponds to the dispersion energy. The dispersion energy is given by the following multipole expansion30
f2n ) 1 -
[∑ ] 2n
(bijrij)k
k)0
k!
exp(-bijrij)
(5)
Pellenq and Nicholson have recently shown20 that the dispersion coefficients (C6ij, C8,ij, C10ij, ...) can be all obtained from the knowledge of the dipole polarizability and an effective number of electrons which is directly linked with the partial charge of each interacting species. The subscripts i and j have the same meaning as those used for the Coulombic and the induction interactions. The f2n functions (n ) 3, 4, 5) are damping functions used to represent possible electronic exchange between the two species in interaction. Their effect is to attenuate the dispersion part of the potential function at short separations.30 These functions are parametrized with the single repulsive parameter bij corresponding to each ij pair in interaction. The quantum mechanical justifications of such damping functions can be found in ref 31. The dispersion three-body terms mentioned in the last expression involve triplets of species in which one atom pertains to the adsorbed molecules and the two others are zeolite species. It has been shown20 that the three-body dispersion coefficients involved in the three-body dispersion terms can be also obtained from the dipole polarizability and the partial charge of each interacting species. Thus these two atomic parameters are crucial for the calculations of the Coulombic, induction, and dispersion terms. The dipole polarizabilities for zeolite atoms are those given in section III. Table 1 presents the electronic characteristics (partial charge and dipole polarizability) for all species considered in this work including those for each atomic component of each adsorbate molecule. It is known that the atomic polarizability varies with the atomic net charge which reflects the atomic environment in terms of bonding. Recently, Sheraga et al.32 have put forward linear correlations between the atomic net charge and the atomic polarizability for a small number of species such as hydrogen, carbon, nitrogen, and oxygen. These authors have shown that the atomic polarizability increases as the partial charge (excess in electron) increases. We have used in this work the same linear correlation between the net atomic charge and the atomic polarizability taking the polarizability of the neutral isolated species as a reference for a partial charge equal to zero. As far as methane is concerned, the sp3 carbon atom and the hydrogen atoms have negligible atomic partial (29) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, N. A. In Intermolecular Forces; Clarendon Press: Oxford, U.K., 1981. (30) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 3726. (31) Pyper, N. C. Philos. Trans. R. Soc. 1986, A230, 107. (32) No, K. T.; Cho, K. H.; Jhon, M. S.; Sheraga, H. A. J. Am. Chem. Soc. 1993, 115, 2005.
4772 Langmuir, Vol. 12, No. 20, 1996
charges.33 Therefore their atomic polarizabilities are that of the corresponding neutral and isolated species: 11.89 a03 and 4.49 a03 for sp3 carbon and hydrogen, respectively. As far as m-xylene is concerned, atomic polarizabilities for sp2 carbon and hydrogen are calculated using the correlations given in ref 32 (see Figure 1 in ref 32) using the partial charges given in ref 34: we found 12.8a03 and 4.35a03 for sp2 carbon and hydrogen respectively. In this model for aromatic molecules,34 the methyl group of the m-xylene molecule is taken as a whole: its polarizability is 13.6 a03.35 The last term in the expression of the total adsorbateadsorbent interaction (eq 1) is the repulsive energy. This term represents the overlap between the wave functions of two interacting species at short range separations (Pauli principle). This term can be obtained in principle only from ab initio Hartree-Fock (SCF) calculations since it is not possible to derive any analytical expressions (as is the case for the Coulombic, induction, and dispersion terms using the perturbation theory). As all other short range terms, the repulsive interaction is not pairwise additive. Unfortunately, the theory on intermolecular interaction has not been able to formulate any analytical expression for the three-body exchange repulsion interactions (for a recent review on this subject see ref 36). Thus, the two body expression used for the repulsive interaction in eq 1 is an approximate expression (see eq 6), which therefore has to be considered as an “effective repulsion term”. This expression is only distance dependent. It has been shown37 that a more accurate analytical representation for the repulsive interaction should be anisotropic, i.e., should have an angle dependence in addition to the distance dependence found in the isotropic case (eq 6). From a theoretical point of view, the energy from SCF calculations strictly corresponds to the sum of the Coulombic, induction, and repulsive parts of the total potential function (see eq 1). This approach has been recently used to calculate the Coulombic and repulsive interactions for the CH4/CH4 and N2/N2 systems.38 Obviously, the SCF method requires a repeat of the calculations for a large number of molecular separations and relative orientations. However, this type of calculation faces several technical problems when applied to the gas-solid interface. The correct approach requires taking into account the crystalline periodicity.39 Less accurate treatments consider a fragment of the solid (cluster approach) whose peripheral bonds are terminated with dummy atoms (usually hydrogen atoms). The influence of capping atoms in the cluster approach has been studied recently by Hess et al.40 in the case of zeolite. These authors compare their results obtained on crystal fragments of different size with periodic calculations: they found that the partial charges converge to their value obtained with the periodic method for the largest molecular fragments. When considering solid-state periodic SCF calculations, it is clear that the number of atoms per unit cell and the size of the atomic basis set are the two limiting factors. These methods have only been applied to adsorption of small molecules adsorbed on a crystal surface whose unit cell only contains a few atoms (see, for example, (33) Ciolowski, J.; Nanayakkara, A. J. Am. Chem. Soc. 1993, 115, 11213. (34) Jorgensen, W. L.; Nguyen, T. B. J. Comput. Chem. 1993, 15, 195. (35) Bader, R. F.; Keith, T. A.; Gough, K. M.; Laidig, K. E. Mol. Phys. 1992, 75, 1167. (36) Elrod, M. J.; Saykally, R. J. Chem. Rev. 1994, 94, 1975. (37) Stone, A. J.; Tong, C. S. J. Comput. Chem. 1994 15, 1377. (38) Schindler H.; Vogelsang R.; Staemmeler V.; Siddiqi M. A.; Svejda P. Mol. Phys. 1993, 80, 1413. (39) Dovesi, R.; Roetti, C.; Freyria-Fava, C.; Apra, E.; Saunders, V. R.; Harrison, N. M. Philos. Trans. R. Soc. London 1992, A341, 203. (40) Nicholas, J. B.; Hess, A. C. J. Am. Chem. Soc. 1994, 116, 5428.
Pellenq et al.
ref 41). It is clear that the full periodic SCF treatment of adsorption of m-xylene in zeolite NaY, for example (whose unit cell contains 624 atoms), is a huge computing task out of reach of the biggest computers available. Extended Hu¨ckel theory for solid state provides an alternative solution which can be applied to large systems. This semiempirical periodic technique was used recently to calculate repulsive energy profiles of argon adsorbed in silicalite42 and small alkanes adsorbed on graphite.43 For simulation purpose, the repulsion interaction between two interacting species is usually represented with an exponential Born-Mayer term
Urep(rij) ) Aij exp[-bijrij]
(6)
Thus, in the case of adsorption in NaY zeolite, there is a pair of repulsive parameter (A and b) for each couple ij; the subscripts i and j describe the atoms of the adsorbate molecule and the zeolite species, respectively. The fact that the edge of a supercage cavity is covered with oxygen atoms (the silicon and aluminum atoms being further away in the crystal) and with sodium cations (in site II) implies that the repulsive interaction between the adsorbate and the T atoms (silicon or aluminum) can be neglected since there is no electron cloud overlap. This simplification of the adsorbate-zeolite potential function has been proven valid in the case of computer simulation of adsorption of simple nonpolar molecules in silicalite3,20 and AlPO4-5.4 Thus only the Coulombic, the induction and the dispersion terms in the adsorbate-T atoms (Si or Al) potential functions are considered in the present calculations, the adsorbate/zeolite Coulombic interaction being zero for xenon and methane. Bo¨hm and Ahlrichs44 have reported some combination rules to estimate the repulsive parameters for a pair of heteroatoms ij using those of the ii and jj pairs. We have chosen the second of the four Bo¨hm and Ahlrichs rules. This combination rule considers the following equations:
Aij ) xAiiAjj
(7)
2biibjj bii + bjj
(8)
and
bij )
Thus in the present case, it is necessary to have estimates of the repulsive parameters for the pairs (Oδ--Oδ-), (Na+Na+) as far as the zeolite atoms are concerned and the pairs (Csp2-Csp2), (Csp3-Csp3), (H-H), and (Xe-Xe) for the adsorbates considered in this study. The (Csp2-Csp2), (Csp3Csp3), (H-H) and (Xe-Xe) repulsive parameters are obtained from the fit of the repulsive part of the best shortrange (Lennard-Jones) potential function available reported in the literature with a Born-Mayer form (see section IV.2): AXe-Xe ) 1025.7Eh, bXe-Xe ) 1.728a0-1, ACsp3-Csp3 ) 322.0Eh, bCsp3-Csp3 ) 1.709a0-1, AH-H ) 1.339Eh, bH-H ) 1.995a0-1, ACsp2-Csp2 ) 35.681Eh, bCsp2-Csp2 ) 1.681a0-1, ACH3-CH3 ) 57.37Eh and bCH3-CH3 ) 1.519a0-1. The (Na+-Na+) repulsive parameters are obtained from those for the (Xe-Na+) and (Xe-Xe) pairs reversing eqs 7 and 8. The (Xe-Na+) repulsive parameters are obtained from the fit of the ab initio Xe-Na+ potential curve recently (41) McCarthy, M. I.; Hess, A. C.; Harisson, N. M.; Saunders, V. R. J. Chem. Phys. 1993, 98, 6387. (42) Pellenq, R. J.-M.; Pellegatti, A.; Nicholson, D.; Minot, C. J. Phys. Chem. 1995, 99, 10175. (43) Vidal-Madjar, C.; Minot, C. J. Phys. Chem. 1987, 91, 4004. (44) Bo¨hm, H. J.; Ahlrichs, R. J. Chem. Phys. 1982, 77, 2028.
Adsorption on Zeolite
reported by Freitag et al.:45 AXe-Na+ ) 689.0Eh and bXe-Na+ ) 2.009a0-1. We found ANa+-Na+ ) 462.83Eh and bNa+-Na+ ) 2.430a0-1. Since all types of framework oxygens have the same partial charge (δ ) -0.823e), it is reasonable to take the same (Oδ--Oδ-) repulsive parameters, whichever type of oxygens are considered. This, along with the neglect of the repulsive interaction between the adsorbates and the T atoms (Si and Al), further simplifies the calculation of the adsorbate-zeolite interaction energy. As a starting set of values, we have taken the (Oδ--Oδ-) repulsive parameters derived for oxygen atoms in silicalite zeolite (having a partial charge of -1.0e): AO-O ) 1543.5Eh and bO-O ) 2.190a0-1.20 We have chosen here to modify only the adsorbate-oxygen b parameter. Therefore, we have only one, two, and three adjustable parameter(s) for the Xe/NaY, the CH4/NaY, and m,p-xylene/NaY system(s), respectively. In this way, one keeps constant the AO-adsorbate parameters obtained with the Bo¨hm and Ahlrichs rule; all the Na+-adsorbate repulsive parameters being calculated from the Na+-Na+ repulsive parameters derived from Freitag’s work on the (Xe-Na+) pair. The modified bO-adsorbate parameters are optimized in order to obtain good estimates of the isosteric heat of adsorption at zero coverage and acceptable values for the Henry constant at different temperatures (see section VI). The optimized bO-Xe repulsive parameter is 2.0a0-1. This value is in keeping with the fact that the oxygen partial charge in NaY zeolite (-0.823e) is lower than that observed in silicalite (-1e): this value is slightly larger that the one used previously in the Xe-O potential to describe adsorption in silicalite (1.93a0-1)20 and corresponds to a weaker repulsive interaction than that for xenon in silicalite. Furthermore, we note that our value for the (Xe-O) b parameter is close to the value given by Ramdas and Thomas (1.9a0-1).46 Reversing the Bo¨hm and Ahlrichs combination rule, the corresponding oxygen-oxygen repulsive parameters are AO-O ) 1543.5Eh and bO-O ) 2.370a0-1 (the A parameter is kept constant as already mentioned). This new value for bO-O is in good agreement with that found by Boutin et al. in their adjustement of the methane/AlPO4-5 adsorption potential function4 (the oxygen partial charge in AlPO4-5 is -0.9e, close to the value found in NaY zeolite). We therefore used the set of repulsive parameters derived by Boutin et al.4 for the pair methane-oxygen (in AlPO4-5) as starting values: AC-O ) 705.0Eh, bC-O ) 1.999a0-1, AH-O ) 45.47Eh, and bH-O ) 2.176a0-1. We found that only a small adjustment of bC-O to 2.051a0-1 was needed to obtain good zero-coverage data (see section VI.2). Again this is in keeping with the lower oxygen partial charge in NaY compared to that in AlPO4-5. It is interesting to note that for the xenon/NaY and CH4/NaY systems, only a very small adjustment of the badsorbate-O repulsive parameter (from that used in previous works) is necessary in order to have a good agreement between simulation and experiment at zero coverage; this adjustment being in each case justified from a physical point of view. Thus, the adsorbate/zeolite potential model used in this work presents a good degree of transferability at least for simple adsorbate. The repulsive parameters for the m- and p-xylene/NaY systems are obtained from the same combination rule using the same set of repulsive parameters for the (O-O) and the (Na+-Na+) pairs as for xenon and methane: AO-O ) 1543.5Eh, bO-O ) 2.37a0-1, ANa+-Na+ ) 462.83Eh, and bNa+-Na+ ) 2.43a0-1. In this way, the method to derive the (45) Freitag, A.; van Wullen, C.; Staemmler, V. Chem. Phys. Lett. 1995, 192, 267. (46) Ramdas, S.; Thomas, J. M. In Specialist Periodical Reports on Chemical Physics of Solids and Their Surfaces; Roberts, M. W., Thomas, J. M., Eds.; The Chemical Society: London, 1978; p 31.
Langmuir, Vol. 12, No. 20, 1996 4773 Table 2. Dispersion and Repulsion Potential Coefficients in Atomic Units oxygen I, II Xe H(-Csp2) Csp2 H(-Csp3) Csp3 CH3(-Csp2)
C6
C8
C10
A
b
104.4 13.92 44.80 14.95 42.00 55.74
2728 251.0 1029 273.7 949.7 1215
77520 4726.0 24000 5370.4 21840 27470
1258 315.6 233.2 45.47 705.0 297.6
2.00 2.36 1.87 2.18 2.05 1.65
Table 3. Dispersion and Repulsion Potential Coefficients in Atomic Units oxygen III Xe H(-Csp2) Csp2 H(-Csp3) Csp3 CH3(-Csp2)
C6
C8
C10
A
b
119.7 16.00 51.44 17.18 48.21 63.82
3219 298.9 1219 326.8 1126 1438
93880 5886 29400 6665 26780 33690
1258 315.6 233.2 110.8 279.7 297.6
2.00 2.36 1.77 2.18 2.05 1.65
Table 4. Dispersion and Repulsion Potential Coefficients in Atomic Units sodium (Na+) Xe H(-Csp2) Csp2 H(-Csp3) Csp3 CH3(-Csp2)
C6
C8
C10
A
b
17.56 2.211 7.310 31.408 6.876 9.657
342.8 28.52 125.4 31.41 115.6 150.5
1258 406.7 2346 478.2 2122 2662
689 172.7 127.7 65.62 153.1 163.0
2.00 2.58 1.98 1.98 1.82 1.86
adsorbate/zeolite repulsive parameters is consistent whichever adsorbates (polar or nonpolar) are considered. The three badsorbate/O repulsive parameters of the p-xylene/NaY system have to be systematically decreased by 0.2a0-1 in order to have a reasonable agreement between simulation and zero coverage experiment data (see section VI). This shows the limit of the adsorbate/zeolite potential model when one has to consider complex systems. The adsorbate/zeolite potential model presented here has also been used recently to describe the adsorption of rare gases in silicalite3 and argon and methane in AlPO45.4 It has the clear advantage of providing a simple, accurate, and systematic method to calculate the coulombic, the induction, and the dispersion terms (including high order two- and three-body terms if needed). The method used to evaluate the repulsive interaction is, as already mentioned, empirical: the main uncertainty only being in the fine determination of the b adsorbate-oxygen repulsive parameter(s) estimated a priori with the Bo¨hm and Ahlrichs rule. We have shown that the adjusted repulsive parameters for the Xe/NaY and CH4/NaY systems are consistent with the value of the partial charges carried by the zeolite species. All coefficients for the dispersion and repulsion terms used for the total adsorption potential functions are given in Tables 2, 3, 4, 5, and 6. The Xe/NaY system is certainly the simplest and therefore the best system to consider in order to evaluate the contribution of the different short-range terms (induction, dispersion, and repulsion) in the total adsorbate/ zeolite potential. Preliminary calculations for the Xe/NaY system have shown that three-body dispersion terms account for less than 2% of the total energy in the supercages and it therefore has been neglected (see Figure 3). The three-body terms considered in this work are the following: XeOI,IIOI,II, XeOI,IIOIII, XeOIIIOIII, XeOI,IISi, XeOIIISi, XeOI,IIAl, and XeOIIIAl. It was not necessary to
4774 Langmuir, Vol. 12, No. 20, 1996
Pellenq et al.
Figure 3. Energy profile across a supercage (z direction). The different contributions to the total energy are shown: (circle) total energy, (square) induction energy, (triangle) three-body dispersion energy, (cross) two-body dispersion-repulsion energy. Table 5. Dispersion Potential Coefficients in Atomic Units silicon Xe H(-Csp2) Csp2 H(-Csp3) Csp3 CH3(-Csp2)
C6
C8
33.28 4.431 14.27 4.762 13.38 17.78
764.2 65.68 282.1 71.83 259.8 332.0
Table 6. Dispersion Potential Coefficients in Atomic Units aluminum Xe H(-Csp2) Csp2 H(-Csp3) Csp3 CH3(-Csp2)
C6
C8
25.28 3.405 10.90 2.402 10.22 13.43
562.3 46.83 205.7 31.41 189.1 240.4
consider the three-body terms containing one sodium cation or two T (Si or Al) atoms since these atoms have a dipole polarizability significantly smaller than that of oxygen atoms. Furthermore, since xenon is the more polarizable species among all other atomic fragment used in this work (see Table 1b), it is therefore safe to neglect all XOO interaction terms in the methane/zeolite and xylene zeolite adsorption potential (X ) Csp3, Csp2, H, CH3). The situation, as far as the magnitude of three-body xenon/ NaY terms is concerned, is very different than that found in the case of rare gases adsorbed in silicalite20 where the three-body terms were about 10-20% of the total interaction energy. The reasons for such differences are certainly difficult to assess. The induction energy appears to be 17% of the total interaction energy when the xenon atom is close to the cavity walls (see Figure 3). This term cannot be neglected in the adsorbate/solid potential as it is the case in most molecular simulations of adsorption in NaY. Again the situation with NaY is different than that encountered with the rare gases adsorbed in silicalite.20 Indeed it has been found that the induction energy was always negligible in silicalite channels. The origin of such a difference is certainly due to the (type II) sodium counterions which are on the edge of the supercages. In Figure 4, the Xe/Na+ and the Xe/O dispersion-repulsion pair potentials used in this work are compared to the
Figure 4. (a) The Xe/Na+dispersion-repulsion pair potential functions. (b) The Xe/O dispersion-repulsion pair potential functions. The dispersion-repulsion potential curves used in this work are compared to those derived within the Kiselev model: (circle) this work; (square) Kiselev model (see text).
corresponding Lennard-Jones potentials parametrized with the Kiselev method.5,13 It is clear from Figure 4a that the dispersion-repulsion energy obtained with the Xe/Na+ Kiselev potential underestimates significantly the dispersion-repulsion energy found with the Xe/Na+ potential function used in this study. However in both cases, the magnitude of the Xe/Na+ dispersion-repulsion energy is small compared to that for the Xe/O pair (see Figure 4b). It is interesting to note that both potential models are equivalent as far as the Xe/O pair is concerned as demonstrated in Figure 4b. These remarks confirm the importance of the induction energy in the description of the xenon/zeolite interaction, which depends directly upon the choice of the zeolite partial charges. IV.2. The Adsorbate-Adsorbate Interaction. An atom-atom model is used to describe adsorbateadsorbate interactions. For xenon, a usual LennardJones form is used with the standard parameters ( ) 211 K, σ ) 3.849 Å).47 This representation of the adsorbateadsorbate interaction has to be considered as a first approximation. Indeed, in the much simple case of argon adsorbed in silicalite, it has been demonstrated that highorder three-body dispersion terms involving more than one adsorbate are not negligible and have some influence on thermodynamic quantity such as the isosteric heat of (47) Barker, J. A.; Watts, R. O.; Lee, J. K.; Schafer, T. P.; Lee, Y. T. J. Chem. Phys. 1974, 61, 3081.
Adsorption on Zeolite
Langmuir, Vol. 12, No. 20, 1996 4775
Table 7. Adsorbate-Adsorbate Interaction Parameters (see text) in Atomic Units (Number in Brackets Are Power of Ten) Xe-Xe s(a0)
e(Eh)
7.273
6.6819(-4) CH4-CH4
C6C-C(Eha06)
C8C-C(Eha08)
C10C-C(Eha010)
134.9
5572
234300
AC-C(Eh)
bC-C(a0-1)
AC-H(Eh)
bC-H(a0-1)
322.4
1.71
1.334
1.994
C8H10-C8H10 sC-C(a0)
eC-C(Eh)
sC-H(a0)
eC-H(Eh)
sH-H(a0)
eH-H(Eh)
6.708
1.116(-4)
5.539
7.303(-5)
4.573
4.781(-5)
adsorption.48 Methane and m-xylene were treated atomically except for methyl groups (in m- and p-xylene), which are taken as united atoms centered on carbons. The molecules are held rigid in the calculations. The potential energy between two adsorbed molecules is the sum of Coulombic, dispersion, and repulsion terms. For methane, we have chosen the modified RMK methane potential of Meinender and co-workers49 already used in a MD and GCMC studies of methane in zeolites.4,6 This potential function was fitted to a variety of solid and gas phase properties. It consists of an exponential shortrange repulsive part between two carbons given by
Urep(rC-C) ) AC-C exp(-bC-CrC-C)
(9)
and between a carbon and a hydrogen by
Urep(rC-H) ) AC-H exp(-bC-HrC-H)
(10)
The repulsive term between two hydrogens was set to zero. The dispersion interaction acts only between carbon atoms and includes terms up to C10
{
Udisp(r) ) -f(r)
C6 r
6
C8 + r
8
}
C10 +
r10
(11)
with the damping term f(r) ) 1 for r > 5.18 Å and f(r) ) exp{-(1 - 5.18/r)2} elsewhere. This term represents the possible electronic exchange31 in a less sophisticated form than that presented in the adsorbate/zeolite potential function (see eqs 4 and 5). The octopole-octopole interactions between methane were modeled by assigning a charge of +0.143e to the hydrogen atoms and -0.572e to the carbon atom. We found that this electrostatic interaction part did not quantitatively change our results (the value of the methane octopole being small) and was therefore neglected in our simulations. This assumption has been proved valid in MC simulation of liquid methane.38 As far as p- and m-xylene are concerned, the mathematical representation was much simpler since it is based on a Lennard-Jones form (repulsive and dispersion terms) and a Coulombic term
(( ) ( ))
U(rij) ) 4ij -
σij rij
6
+
σij rij
12
+
qiqj rij
(12)
The Lennard-Jones parameters are obtained using geometric combination rules given by Jorgensen.34 These (48) Fernandez-Alonso, F.; Pellenq, R. J.-M.; Nicholson, D. Mol. Phys. 1993, 86, 1021. (49) Meinender, N.; Tabisz, G. C. J. Chem. Phys. 1983, 79, 416.
short-range parameters are in general transferable while the choice of the partial charges for a new molecule remains the main problem. The optimized potentials for liquid simulations (OPLS) derived by Jorgensen,34 reported previously for benzene and attached methyl, hydroxyl, methoxyl, and cyano groups, have been used in this work. All unsubstituted benzene carbons and hydrogens have the same partial charges ((0.115e) as in benzene itself. The partial charge in the particular case of methyl groups is also given at +0.115e (equal to that of hydrogen atoms) but differs for other substitution groups. It should be said however that these values of partial charges have been optimized to reproduce thermodynamic liquid properties at standard conditions of temperature and pressure. Further tests on supercooled different aromatic systems validate such empirical potential functions. The set of partial charges used here does not provide accurate values of the different n-pole moments carried by xylene molecules34 (the first nonzero moments are a dipole and a quadrupole for m- and p-xylene molecules, respectively). These partial charges were also used to calculate the dispersion coefficients of the adsorbate/zeolite potential function (see ref 20) so that all potential coefficients are obtained in a consistent manner. All coefficients for the dispersive and repulsive adsorbate-adsorbate potentials functions are given in Table 7. V. Grand Canonical Monte Carlo (GCMC) Simulation of Adsorption V.1. Basics. In the grand canonical ensemble, the independent variables are the chemical potential, the temperature, and the volume. At equilibrium, the chemical potential of the adsorbed phase equals the chemical potential of the bulk gas, which can be written as a function of temperature and the bulk gas pressure50
µads ) µbulk(T,P)
(13)
Consequently, the independent variables in the grand canonical Monte Carlo simulation are the temperature, the pressure of the bulk gas, and the volume of the zeolite unit cell. This set of variables is convenient because temperature and pressure are independent variables in the adsorption isotherm. Therefore, the adsorption isotherm can be obtained directly from the simulation by evaluating the ensemble average of the number of adsorbate molecules whose chemical potential equals that of the bulk gas at given temperature and pressure. Note that the bulk gas is usually assumed to obey the ideal gas law. The PVT data for the adsorbate molecules studied (50) Nicholson, D.; Parsonage, N. G. In Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: London, 1982.
4776 Langmuir, Vol. 12, No. 20, 1996
Pellenq et al.
indicate that this assumption is valid for the temperature and pressure conditions considered in this work. However, this is not a necessary assumption in the use of the grand canonical Monte Carlo method. Monte Carlo simulation in the grand canonical (µVT) ensemble is now a well established and widely used technique in the field of adsorption (see for instance ref 50). The fact that in GCMC simulations, the number of particles varies, helps to minimize ergodic difficulties. The simulation algorithm is based on the original algorithm proposed by Adams51 and Soto52 which allows some creations and destructions of adsorbed molecules in the system (see below). Each adsorbate molecule has a position vector associated with it, defined with respect to a Cartesian coordinate system. Nonspherical molecules also have in addition three orientation vectors. In the first step of the simulation cycle, a molecule (already adsorbed in the zeolite cavity) is picked at random and is moved. For xenon, this move is always a translation, while for methane and xylene isomers, moves can be translations or rotations. The decision whether to translate or reorient is made at random and with an equal probability. The potential energy of the new configuration called the trial configuration is calculated. A decision is then made whether to accept the move or to return to the old configuration according to the usual Metropolis method50 with a probability ACC PTRANS/ROT
) min(1, exp(-β
∆UTRANS/ROT ) N
(14)
where ACC ) Unew - Uold ∆UTRANS/ROT
(15)
The second step in the simulation cycle is the decision whether to attempt to add a molecule or to subtract a molecule. This decision is made at random with equal probability. If an addition is attempted, then a random position (within the volume V) and three orientation vectors are generated for the new particle. In such a creation trial, the probability of an acceptance is given by50
(
PACC CR ) min 1,
zV exp(-β∆UCR N ) N+1
)
(16)
where the initial number of molecules is N (before addition), z ) exp(βµ)/Λ and Λ is the reciprocal of the space and momenta partition function; for spherical particle Λ is given by
Λ)
h (2πmkT)1/2
(17)
If the addition is accepted, then the new molecule is kept. Otherwise, the new molecule is deleted and the system returns to the old configuration. If a subtraction is attempted, a randomly chosen particle of the system is deleted. The potential energy of the new configuration is calculated and the subtraction is accepted or rejected according to the probability (from N + 1 to N molecules)50
(
PACC DES ) min 1,
N exp(-β∆UDES N ) zV
)
(18)
Control charts in the form of plots of a number and internal energy versus the number of Monte Carlo steps were used (51) Adams, D. J. Mol. Phys. 1975, 29, 307. (52) Soto, J. L.; Meyer, A. L. Mol. Phys. 1981, 42, 971.
to monitor the approach to equilibrium. Acceptance rates for creation destruction were also followed and should be equal at equilibrium. The step length (for translation) and the step angle (for rotation) were adjusted to give an acceptance rate for translation and rotation of 50%. After equilibrium had been attained, all averages were reset and calculated over several million configurations. After a certain number of cycles, the averages tend to fluctuate about a steady value rather than converging to their final values and remaining constant. Isosteric heats of adsorption were obtained from fluctuations over the number of particles in the system N and from fluctuations of the internal energy U50
qst ) RT -
∂〈U〉 〈UN〉 - 〈U〉〈N〉 ) RT ∂〈N〉 〈N2〉 - (〈N〉)2
(19)
where the brackets denote an ensemble average. For instance, the ensemble average for any thermodynamic quantity A is
〈A〉 )
1
M
∑ Ai M i)1
(20)
where M is the number of cycles and A the value of a given quantity {U, N, UN, N2} recorded after cycle i. The isosteric heat of adsorption was split in two components: the adsorbate-adsorbent (qw) and the adsorbate-adsorbate interactions (qm). Consequently, the isosteric heat of adsorption can be calculated in two ways: directly from the above expression and from the sum of these two components. This procedure was used as a self-test for the simulation code. Because of fluctuations, it is clear that the accuracy of any quantity measured by computer simulation is accompanied of an intrinsic error. This error can be evaluated from the standard deviation technique as described by Bakaev et al.53 and is less than 5% in Xe/NaY and CH4/NaY simulations and less than 8% in xylene/NaY simulations. In addition to isotherms and adsorption heats, molecular configurations at equilibrium (snapshots), pair distribution functions, and singlet distribution functions (i.e., the density of adsorbed molecules within the volume V) were obtained. The latter were found by dividing the simulation box into a number of slices and calculating the average density in each subsystem. V.2. The Tricks of the Trade. The simulation box (of volume V) was constructed from one NaY unit cell and the periodic boundary conditions applied in three dimensions in order to simulate an infinite (macroscopic) system. The adsorbate-adsorbent potential function defined in the previous section is too complicated to be used directly in a simulation in terms of CPU time required. We therefore used a grid-interpolation procedure in which a part of the NaY zeolite unit cell is split into a collection of small cubes with side length 0.03 nm. We have calculated a grid for each atomic (or pseudoatomic) component of the adsorbate molecules considered in this study: one grid for xenon, two grids for methane (a H(Csp3) grid and a C(sp3) grid), and three grids for xylenes (a H(-Csp2) grid, a C(sp2) grid, and a CH3 grid). The total X-zeolite energy (X ) {Xe, H(-Csp3), H(-Csp2), C(sp3), C(sp2), CH3) was calculated at each corner of each cube. In order to reduce the size of the different grids, the symmetry of a unit cell was exploited. The unit cell on which the xenon grid, the two methane grids and the three (53) Bakaev, V. A.; Steele, W. A. Langmuir 1992, 8, 148.
Adsorption on Zeolite
xylene (short-range) grids were calculated, was embedded in a 27 (3 × 3 × 3) unit cell system during the grid calculation procedure prior to any GCMC simulation. The system was large enough to ensure a good convergence of the short-range (dispersion + repulsion + induction) adsorbate/zeolite interaction energy (the short-range cutoff was set at 30 Å). For xylene isomers however, a larger system (729 ) 9 × 9 × 9 unit cells) has to be considered because of the long range character of the Coulombic interaction (the long range cutoff was set at 100 Å). The grid interpolation method allows taking into account any degree of accuracy in the description of the adsorbate/ zeolite interaction energy since all the needed grids are calculated separately prior to any simulation runs. The accuracy of the interpolation based on this cube size was high except in the strongly repulsive regions of the potential where the strength of the repulsion tends to be overestimated. In the Monte Carlo calculations, however, few acceptances are to be expected in this region so that no serious error was incurred by interpolating the potential function from the grid. The minimum image convention was applied between adsorbate molecules54 during GCMC simulation runs. Therefore, according to the size of the simulation box (one single cubic NaY unit cell, a ) 24.85 Å), the adsorbate-adsorbate interaction cutoff should be less than half the shortest size of the simulation box (here 12.4 Å). This cutoff is certainly suitable for xenon and methane since there is no long range (Coulombic) interaction but is not sufficient for xylene/xylene interactions because of the multipole-multipole Coulombic interaction. The long range Coulombic interaction between m-xylene molecules was corrected using the reaction field scheme.55 However as mentioned by Jorgensen et al.,34 a generally accepted procedure to correctly represent long range interactions has not yet emerged to correct for the Coulombic interaction beyond the cutoff. These authors have shown that long range corrections modify thermodynamic properties of liquid aromatics by less than 3% at ordinary temperature and pressure conditions. The simulation program was run on in-house work stations: a HP720 computer (at CPMA, Orsay, France) and a Silicon Graphics 4400 work station (at IFP, RueilMalmaison, France), and on the C98 CRAY computer at IDRIS (CNRS Orsay, France). A typical simulation run consists on (5-10) × 106 Monte Carlo steps (the first 2 × 106 MC steps were systematically discarded before any average calculations). VI. Results and Discussion Grand canonical Monte Carlo simulations were carried out to determine the adsorption characteristics of pure gases in NaY zeolite. In all cases, we compare our results with available experimental data. This is the reason why xenon simulations were performed at 298 K only. Methane simulations were carried out at 298, 313, 323, 333, and 343 K. m- and p-xylene simulations were performed at 423, 523, and 573 K and are also compared to experiments. VI.1. Adsorption of Xenon at 298 K. Our xenonsimulated isotherm at 298 K is compared to the experimental data recently published by Lui et al.56 and Boddenberg et al.57 in Figure 5. We have also compared (54) Allen, M. P.; Tildesley, D. In Computer Simulation of Liquids: Oxford Press: Oxford, 1987. (55) Neumann, M. J. Chem. Phys. 1986, 85, 1567. (56) Shang-Bin Lui; Fung, B. M.; Tran-Chin Yang; Eng-Chin Hong; Chau-Ting Chang; Pei-Chin Shih; Fu-Hsing Tong; Tun-Li Chen J. Phys. Chem. 1994, 98, 4393. (57) Boddenberg, B.; Sprang, T. J. Chem. Soc., Faraday Trans. 1995, 91, 163.
Langmuir, Vol. 12, No. 20, 1996 4777
Figure 5. Xenon adsorption isotherms at 298 K. Comparison between simulations (this work), Rowlinson et al.,10 and experiments.57
our results to GCMC simulation data given by Rowlinson et al.10 (see Figure 5). It is clear from Figure 5 that our simulated isotherm is in good agreement with experiment in the range of pressure (2 × 104 to 105 Pa). At lower pressure (
