J. Phys. Chem. 1994,98, 13339-13349
13339
Intermolecular Potential Function for the Physical Adsorption of Rare Gases in Silicalite Roland J-M. Pellenqt and David Nicholson* Department of Chemistry, Imperial College of Science, Technology, and Medicine, London, SW7 2AY, United Kingdom Received: June 4, 1994; In Final Form: October 3, 1994@
We present a simple method for the derivation of two-body and three-body dispersion coefficients for incrystal atoms from the knowledge of their dipole polarizability and effective number of electrons. This method is checked by comparison with results from quantum mechanical calculations and used to derive two- and three-body coefficients for the dispersion interaction of argon, krypton, and xenon adsorbed in the siliceous zeolite silicalite- 1. Repulsive parameters for the argon/silicalite system are obtained by fitting experimental data over a wide range of temperature and the full scale potential constructed in this way (referred to as P N l ) is shown to perform well in predicting other argon data. The parameters for the repulsive energy for the kryptordsilicalite and xenonhilicalite systems obtained using combination rules and the PN1 potential for these adsorbates (without any parameter adjustment) are also found to be successful in predicting low coverage properties. We compare the performance of the PN1 function with the Kiselev adsorption potential, widely used in the field of modeling adsorption in zeolite cavities, and show that the latter tends to overestimate thermodynamic properties and also predicts a wider pore than the new potential. The energetics of adsorption are discussed in terms of site location and shape and the effective size of zeolite oxygen atoms.
1. Introduction Numerical simulation is being increasingly used in order to understand adsorption phenomena at the microscopic level. At present there are no statistical mechanical theories sufficiently well developed to be able to predict adsorption properties for heterogeneous surfaces.' Monte-Carlo and molecular dynamics simulations provide an alternative way to model physical adsorption in which ideally the statistical mechanics is not subject to any approximations. The input to a simulation is the potential energy function. To develop this function for adsorbate-adsorbent interactions, the adsorbent is described as an assembly of atoms and the interaction energy is calculated between all the adsorbent and adsorbate atoms. Although the statistical treatment can be envisaged classically, it is clear that the calculation of the interaction energy has to be considered within the framework of quantum mechanic^,^ since the intermolecular interaction potential depends on the electronic properties of each particle. A full scale quantum mechanical calculation of the interaction energy is not feasible for any real system of interest, and approximate methods need to be employed. A conventional approach to this problem is to calculate long-range interaction terms (coulombic, induced and dispersive interactions) from perturbation theory and to couple these with a variational theory such as ab initio (SCF-HartreeFock) methods to calculate the short range exchange and polarization-exchange repulsive2 interactions. Successful calculations have been achieved for the interactions between small molecules based on this type of a p p r ~ a c h . ~ However, the repulsive term is still subject to an uncertainty of about 15%. When a solid phase is involved, as is the case with adsorbate-adsorbent interaction, further difficulties arise because either a small cluster of atoms must be chosen with attendant uncertainties associated with capping, or else a method which uses Bloch functions must be adopted which generally implies a lower level of approximation. In this work we have Present address: Laboratoire de Chimie Physique des Matkriaux Amorphes, batiment 490, Universitk de Paris-Sud 2 Orsay, 91405 cedex, France. Abstract published in Advance ACS Abstracts, November 15, 1994. @
0022-365419412098-13339$04.50/0
therefore chosen a more empirical approach to the repulsive term, adopting the Mayer form for atom-atom terms and using experimental data to obtain parameters. At the same time we have relied heavily on the insights gained from ab initio work on other systems. The immediate aim of the present study was to construct potential functions for rare gas adsorbates in silicalite which represented the state of the art using the minimum number of adjustable parameters. Silicalite-1 is the pure siliceous form of the ZSM-5 class of zeolites; it possess a three-dimensional network of intersecting channels and exhibits a solid-solid phase transition at 340 K from a monoclinic to an orthorhombic form. A first objective was to gain insight into the contributions originating from higher order two-body and three-body dispersion terms, since the balance between the various contributions may be expected to be different in a nonuniform heterogeneous system from that in homogeneous uniform systems. A second objective was to construct the potential in such a way that transferability to other more complex systems would be possible. To this end we concentrate initially on the argon-silicalite system and subsequently demonstrate that our potential functions can be transferred to other rare gases using combination rules established elsewhere. The rare gas silicalite system is particularly suitable for this purpose for two reasons: (i) Since silicalite is the pure silica form of ZSM5, problems arising from charged lattice sites in uncertain locations are minimized. (ii) There is an extensive body of good quality experimental data against which theoretical predictions can be tested. A third longer term objective was to use the potential functions developed as input to simulation studies covering a range of temperatures, and adsorbate loadings. The relevant work has been reported in preliminary form elsewhere,23 and a paper giving fuller details is in press. Many previous calculations of adsorption in silicalite and similar materials have been based on the Kiselev potential or on potential functions constructed in a similar way: i.e., essentially by assuming that oxygen is the only important lattice species. We demonstrate that this method has limitations (as indeed anticipated by its originators), although capable of giving 0 1994 American Chemical Society
13340 J. Phys. Chem., Vol. 98, No. 50, 1994
Pellenq and Nicholson
a remarkably good account of adsorbate-adsorbent interactions in the case of argon adsorption. In the next section we set out the basic perturbation theory formalism for calculating dispersion interactions including twoand three-body terms to high orders of polarizability, a discussion of repulsive interaction is included since this enters into the damping terms. In the following section we describe methods for deriving all the parameters necessary to calculate the dispersion coefficients for interactions involving condensed phases. Parameters determined by these methods are compared with available data from quantum mechanical calculations. The remainder of the paper is devoted to potential energy for adsorbate interacting with silicalite. After a brief section outlining the widely used Kiselev potential we describe the development of a full scale potential (designated PN1 for convenience of reference); we use experimental data for argon in order to establish repulsive parameters. The model is compwed with the Kiselev model and its predictive performance is tested for argon against experimental data not used in the determination of repulsive parameters and for other rare gases. We conclude with a discussion of the strengths and weaknesses of the potential function developed here compared to that obtained with the Kiselev model.
first dispersion coefficients can be then rewritten:
e
= ;hm[af(iw)$(iw) 15
+ g(io)a:(io)]d o
(8)
When w = 0, a(io)becomes the static I-th pole polarizability (compare eq 6) D)
n=
1
Therefore the dynamic I-th pole polarizability can be rewritten as a Pad6 approximant:
2. Dispersion Energy Two-Body Dispersion Terms. Consider a rare gas atom in interaction with an atom of the solid surface. Time-dependent perturbation theory gives the long-range dispersion interactions in terms of correlations between charge fluctuations on the interacting species! The dispersion multipole expansion for two ground-state, spherically symmetric atoms or molecules A and B, separated by a distance R can be expressed as
where 771 is an average l-th pole transition energy. Moreover, making use of the identity 2 -
ab &=- 1 a b nJm(a2 x2)(b2 x2)
+
+
eq 11 with eq 5 gives
e+ e+ e+
AB
- - - ... R6
Ra
R'O
e
The C, coefficients are the dispersion coefficients where describes the interaction between two instantaneous dipoles, the interaction between a quadrupole and a dipole, and C$ includes the interaction between an octopole and a dipole and between two quadrupoles. These coefficients can be rewritten in the following
Therefore, the dispersion coefficients can be rewritten:
e
=P ( 1 , l )
e
(2)
+
= CAB(1,2) P ( 2 , l )
=P ( 2 , 2 )
(3)
+P(1,3) +P(3,l)
(4)
where 1, 2, and 3 denote dipole, quadrupole, and octopole respectively, and P ( 4 , 4 >=
+ '
(24 21 2)' 2Jmat(iw)aE(io) d o 4(2ll)!(2Z2)! n O
Here a/ are the 2l-pole dynamic polarizabilities, 1 {1,2,3}] at frequency w defined as
E
(5)
[I,,
Thus the C, coefficients can be obtained if the set {ali,azi, a3i,vli, v*j, ~ 3 is~ known ) for each interacting species. Three-Body Dispersion Terms. Terms involving triplets of species A, B, and C are obtained from the perturbation theory at the third and fourth order. For a triplet, the total dispersion energy is nonadditive:
12:
where fn' is the oscillator strength which is a dimensionless quantity characteristic of the n-pole transition from the ground state to an excited state with the corresponding energy dk. The
where represents a sum of several terms describing the three-body interaction. The three-body interaction is about 1020% of the total interaction energy for condensed phases.6These terms are not only distance dependent like the two-body terms but also angle dependent. Therefore, the sign of the three-body interaction depends on the geometrical configuration, each contribution being the product of a geometrical term with an
J. Phys. Chem., Vol. 98, No. 50, 1994 13341
Physical Adsorption of Rare Gases in Silicalite electronic term. The general form of each three body function is7
where Wac is the geometrical factor and Zac is the electronic function. The general expression for the electronic component for each three-body term is Zmc(11,12,13) = -&a$o)a;(io)a@o) 1
do
(19)
Using the Pad6 approximation as before gives
The calculation of the three-body dispersion terms requires the same set of parameters as for the two-body terms. A general expression for the geometrical factor Wac(11,12,13) is given in ref 5 from which the following formulas (with relevant permutations of the mixed terms) can be derived with R12, R23, R31 the sides of a triangle in which the triplet of atoms are vertices and 41, 42 ,43 are the interior angles. W(DDD) = 3R;;/Ry:R;f(
1
+ 3 cos dl cos 42cos 43) (21)
W(DDQ) = 3T;Ri;R;f[(9
cos 43- 25 cos 343) (41
The Repulsion Interaction. At the intermolecular separations that occur in condensed phases, molecular wave functions overlap. When the overlap between molecules is substantial, the long-range multipole expansion used in perturbation theory is no longer valid at short range. The perturbation theory approach is questionable at these distances. In the first place, exchange of electrons between the two molecules can no longer be ignored so it is not possible, as in the long-range treatment, to take unperturbed wave functions and antisymmetrized products must be used instead to obey Pauli’s principle. The corresponding energy, the “valence repulsion term,” is proportional to the square of the overlap integral between the orbitals of the two molecules. Since the electronic wave functions decay exponentially with distance at long range, the valence repulsion term is also quite well represented by exponentials in the intermolecular distance at the separations relevant to physical interaction. The usual treatment for the short-range interaction is to assume pairwise additivity and to fit quantum mechanical calculations with an atom-atom Born-Mayer term which is a sum of exponentials of the form
ureP=
c
A exp(-brij)
itj
where A and b are repulsive parameters. There is therefore a pair of repulsive parameters for each different pair of interacting species. Because overlap of electron clouds cannot be neglected in the region of equilibrium separation, it is necessary to correct dispersion terms. To achieve this, Toennies et aL9 derived an alternative expansion, valid for nonzero overlap of the wave functions. This can be written as a damped multipole expansion in the
COS
- 4211 [3 + 5 cos 24311 (22) where (29)
The fourth-order perturbation theory gives a dipole-dipoledipole interaction from which general expression is
The fourth-order dipole-dipole-dipole term is also a sum of products of geometrical terms with electronic factors which are given by8
where a(’) is the dipole polarizability at the vertex i.
In eq 29 b is the repulsive parameter in (27). The damping functions fzn are always positive and are unity at distances sufficiently large for overlap between wave functions to be negligible, but become less than unity at distances for which wave function overlap is appreciable. For equilibrium separations, where the overlap is appreciable, the damping functions decrease rapidly with increasing n. At large distances, where the dispersion damping functions approach unity, the l/Rn dependence ensures that terms with large n are unimportant. It should be stressed that the three-body term in the total dispersion interaction should also be damped, but no suitable expression is available.
3. Parameters for the Dispersion Interaction The Effective Number of Electrons. If qr is taken to be the lowest excitation energy for the 2l-th radiation, the approximant for al(iw) given in eq 11 is everywhere smaller than the true value.5 Conversely, if is taken to be
then al is everywhere greater than the corresponding true value. The general form for the Sl function known as the k-th-order sum rules is
Pellenq and Nicholson
13342 J. Phys. Chem., VoE. 98, No. 50, 1994
The method of Tang et uL5 provides upper and lower bounds for all the dispersion coefficients provided that the set of parameters {ali, azi,a i , Sli(0), Szi(O), Ssi(O)} is known for each interacting species i . According to the Thomas-Reiche-Kuhn sum rules, Sl(0) is equal to the total number of electrons in the system but this is found to greatly overestimate the dispersion coefficients. Rtzer suggested'O that, the number of electrons to be used to calculate c 6 should be the number of electrons in the outermost subshell on the grounds that only these electrons are polarized. Calculations" using time-dependent perturbation theory for the dipoledipole dispersion coefficient support Pitzer's assumption and the demonstration5that (30) is an upper bound. We therefore replace Sl(0) in (30) by ne^, the effective number of electrons. From eqs 2 and 13, the first dispersion coefficient for interaction between like species is
/
0
4
6
10
8
electronic structure
Substituting eq 30 and rearranging gives
Figure 1. Effective number of electronsI6 versus number of electrons
Ntff= 1 6 [ e ] 2 / 9 [ a f ( 0 ) ] 3
(33)
Recently, several calculations have appeared in which the timedependent perturbation theory and a coupled Hartree-Fock scheme have been used to obtain frequency-dependent polarizabilities and dispersion coefficients. Both neutral and charged species have been investigated.I3-l6 In Figure 1, ne^, taken from these data, for neutral isolated atoms, is plotted against the number of outer-shell electrons. In heavier atoms the innershell electrons make some contribution to Neff." For example, considering the sequence Ne, Ar,and Xe in which each atom has six p-electrons in the outer subshell, the effective number of electrons is 4.455,6.106, and 7.906, respectively.16 Each of the plots in Figure 1 can be fitted by a simple second-order polynomial:
Ntff= aN2
2
+ bN ic
(34)
+
where N is the number of electrons in the outer (s p) shells for each species. The coefficients u, b, and c are listed in Table 1. This approach can be extended to atoms in molecules or in crystals. In molecules, the polarizability resides in the bonds4 and Nee includes the valence electrons together with lone pair contribution^.^^ For purely ionic solids, Neffis close to that of the isoelectronic inert g a ~ . ' ~ For + ' ~partly covalent structures, such as silicates, there is no theoretical study of the dispersion interaction. We suggest here a general rule which can be applied to bonded atoms or crystal species whatever the degree of ionicity or covalency of the structure. We assume that the effective number of electrons in the more negative species is increased by the partial charge and that in the more positive species it is reduced by the deficit of charge, thus Neffis found from
where N 2 is the effective number of electrons of the species A in the isolated and neutral (ground) state and 141 is the magnitude of the partial charge on A due to the transfer of electrons during the bonding process. In this manner, the dispersion interaction between a third species (nonbonded) and the species A will
in the outer shell for isolated closed-shell atoms.
TABLE 1: Parameters for the Fit of N a vemu the Electronic Structure H, He Li to Ne Na to Ar K to Kr Rb to Xe
a
b
-0.109 -0.02017 -0.0293 -0.0118 -0.0237
0.9330 0.7075 0.9966 1.0059 1.1767
C
0 0.0857 0.0062 0.0011
0
depend on the nature of the bond(s) which involve(s) the species A. For example in silica, where Si and 0 cany a partial charge of +2 and -1, re~pectively,~~
@g= 3.52 - 2 = 1.52 = 3.65 + 1 = 4.65 Equation 36, as a basis for finding Neg, appears to be reasonable in the light of different calculations for isolated crystal species. Thus Neffi,s the same for an isoelectronic sequence, but varies in an homologous sequence, as found from quantum mechanical calculations for ionic crystal species.18 Quadrupole Polarizability and the Sz(0) Parameter. In order to calculate Cg one needs the sum rule S2(0), and the quadrupole static polarizability a2(0) in addition to Neff and al(0). &(O) is defined as a sum of oscillator strengths associated with an average quadrupole transition energy? Unfortunately, az(0) and Sz(0) are known experimentally for only a few isolated species such as alkali, alkaline earth metals and rare gases. KiselevZ0proposed the following approximate expression for the quadrupole-dipole dispersion coefficient between species AandB
Thus, according to the definition of 771, (eq 30) the Kiselev expression for depends only on the static dipole polarizability and on the effective number of electrons, although we note that the expression for pairs of like atoms does not depend on the effective number of electrons but only on the static dipole
e
J. Phys. Chem., Vol. 98, No. 50, 1994 13343
Physical Adsorption of Rare Gases in Silicalite 12
,
I
I
0,
I
TABLE 2: Parameters Used for Isolated and Neutral Species (All Data in atomic units) H
He Li Be C N 0
F
Ne Na Mg
Ar K
Ca Kr Xe
0
1
2
3
4
5
6
7
a
g
Number of effective electrons
Figure 2. Correction factors for the dipole-quadrupole dispersion coefficients (0) and quadrupole-quadrupole dispersion coefficients C$ (0)for isolated closed-shell atoms, plotted against Neff taken from Table 2.
e
polarizability of the species. A similar expression derived by Narayan,21 although different for a heterogeneous pair, is proportional to the Kiselev formula (by a factor of 2.4) for like atoms. Values of obtained from (37) seriously underestimate those from correlated Hartree-Fock calculations for isolated neutral atoms.22 We therefore define Kg ('1) as the ratio of obtained from quantum mechanics to those found from (37). Figure 2 shows this correction factor plotted against the effective number of electrons. The correction factor for incrystal or for in-molecule atoms can be obtained from this plot once Nee has been found from eq 34. For pairs of like atoms, eq 15 equated with eq 37 multiplied by this correction factor gives A A 112
2 4 [ 9 S f a f ] ' " = Kgaf[(%S1 )
+ (9S1Aa,A 3 1141
(38)
St(0) can be expressed in terms of the effective number of electrons (= St(0)) and the static dipole p~larizability,'~ as A A 112
S; = [9& a, 1
(39)
Substitution of (39) into (38) yields a quadratic equation in (4)li2:
Octopole Polarizability and the $(O) Parameter. In order to calculate the Clo dispersion coefficient according to (16) we need &(O) and a3(O) in addition to the quantities already determined. &(o) is defined5 as a sum of oscillator strengths associated with the octopole transition (averaged) energy. An approximate expression for Clo in terms of first-order quantities can be written20
ado)
Ka
Kio
4.500 1.385 164.5 37.84 11.88 7.423 5.412 3.759 2.676 165.0 71.32 11.075 293.1 153.9 16.75 27.32
1.635 2.066 1.162 1.944 2.423 2.552 2.653 2.750 2.830 1.414 2.236 3.289 1.602 2.129 3.512 3.602
2.974 3.918 1.248 3.287 4.868 5.282 5.724 6.21 1 6.650 1.934 4.283 9.5 11 2.327 4.189 11.054 11.553
Unlike the corresponding quadrupole equation (eq 37) the expression for pairs of like atoms now depends on the effective number of electrons. As for the estimates of obtained with this equation seriously underestimate quantum mechanical values for isolated neutral atoms.22 As before, we define K ~ as o the ratio between the values obtained from quantum calculations and those given by eq 41 and derive a quadratic equation in (a3(0))lI2:
e,
e
A A A 112 A - K!(atSf)1/2(at)1/2 (S,S,/a,) a3
K'a;(Sf)'l2
e
e
NCf 0.824 1.430 0.773 1.420 2.649 3.179 3.656 4.080 4.455 0.988 1.874 6.106 1.030 1.961 7.305 7.901
= 0 (42)
where A A 2 = 1.6406K10(af/Sf)'12- 0 . 6 2 5 (AS 2A112 / ~ ) ( % / a l ) (43)
In this equation, S;' can be expressed in terms of known quantities by the relationlg
Comparison with Previous Work. The preceding treatment for the dispersion interaction provides a simple means of calculating dispersion coefficients not only for simple and isolated systems, such as alkali, alkaline earth, or rare gas atoms, but also for other open-shell isolated systems such as carbon, nitrogen, oxygen, or fluorine. Moreover, it is readily extended to bonded or to crystal species in modeling ionic, semiionic or molecular crystal structure. The only inputs required are the dipole polarizability and the effective number of electrons for each interacting species. From these two quantities, we have shown that all the two-body and three-body coefficients can be obtained. A knowledge of the dipole polarizability is crucial: as previously shown,23 the dipole polarizability of crystal species, at least for anions, strongly depends on the ionic environment and is very different from the tabulated polarizability for corresponding isolated systems. Since the effective number of electrons is the same for a free ion or for a crystal ion, the atomic environment is mainly characterized by the dipole polarizability, which varies greatly from one environment to another. Calculations of c 6 and Cs, based on the method described here, have been made for a number of neutral and ionic species, and are compared with ab initio results24in Table 4 using the parameters in Tables 2 and 3. Calculations for c6, Cg, and c10 are usually in close agreement with ab initio results: our results for c6 are generally about 10% too high while those for c8 are
13344 J. Phys. Chem., Vol. 98, No. 50, 1994
Pellenq and Nicholson
TABLE 3: Parameters Used for In-Crystal Species (AU Data in atomic units) Lit
Na+ K+ Mgz+
Ca2+ Siz+ 0- (SiOz) Oz- (MgO)
F- (LiF) C1- (NaC1)
1.430 4.455 6.106 4.455 6.106 1.520 4.656 5.656 5.084 6.550
0.192 1.002 5.339 0.486 3.193 2.360 8.030 12.01 5.367 19.98
1.947 2.842 3.291 2.842 3.291 2.002 2.894 3.171 3.011 3.398
3.508 6.708 9.416 6.708 9.416 3.662 7.003 8.654 7.683 9.811
TABLE 4: Two-Body Dispersion Coefficients for Rare Gas-Crystal Species Interactions and for Crystal Species-Crystal S ecies Interactions. Comparison with ab Initio Calculations
E
C8
c 6
system He-Li+ He-Na+ He-F- (in LiF) He-Cl- (in NaCl) Ne-Li+ Ne-Na+ Ne-F (in LiF) Ne-Cl- (in NaCl)
ab initio
this work
ab initio
this work
covalent crystal. The model uses a Lennard-Jones form of potential function u = - - + -C
R6
B R12
(45)
where R is the distance between interacting atoms and C and B are the dispersion and repulsion constants, respectively. Thus the dispersion multipole expansion is reduced to a single term which must account for higher order two-body and three-body interactions. The dispersion constant C is identified with the f i s t term in the dispersion expansion given in eq 1. The repulsive constant B is obtained from the condition that the potential function is minimum at the equilibrium separation of the interacting pair, the equilibrium separation being calculated as the sum of the van der Waals radii:
+
B = 1/2C(rl rJ6
Therefore, the parameters in this model are the dipole polarizability, the number of effective electrons, and the van der Waals radius of each species. In order to reduce the number of parameters, only the interactions with oxygen atoms of the silicalite crystal were included. This presupposes that the influence of silicon atoms can be allowed for by fitting the parameters for the oxygen atoms in silicalite. Furthermore, the Ar-Li+ effective number of electrons was taken as the total number of Ar-Na+ electrons in the considered species. The parameters for Ar-F- (in LiF) adsorbate species were assumed to be known quantities. Thus Ar-Cl- (in NaCl) if one uses the van der Waals radius of the isolated and neutral Li+-Li' oxygen atom in its ground state,2gthe only remaining parameter Na+-Na+ F--F- (in LiF) is the oxygen polarizability. The latter can then be obtained Li+-F- (in LiF) by fitting experimental thermodynamic quantities such as the Cl--Cl- (in NaC1) Henry law constant or the isosteric heat of adsorption in the Na+-Cl- (in NaC1) limit of zero coverage, using the orthorhombic structure of silicalite. The potential parameters for some adsorbate molgenerally about 5% too low. In view of the range of variations ecules and the silicalite framework oxygen are given in which can occur in quantum mechanical calculations, the present reference.28 It is interesting to note that the oxygen dipole method can claim acceptable accuracy for the prediction of in polarizability used in the Kiselev model (0.85 A3) is significantly crystal dispersion coefficients. A more detailed paper on this method of calculating dispersion coefficients is in p r e p a r a t i ~ n . ~ ~ smaller than that derived from Auger electron spectroscopy data (1.2 A3);23the value of 0.85 A3 corresponds to that of an isolated and neutral oxygen. It has been suggested30that the neglect of the silicon is justified because the framework silicon atoms have 4. The Kiselev Intermolecular Potential for the a very small polarizability and, being further away from the Adsorption of Argon in Silicalite pore center compared to the oxygen atoms, contribute very little A complete adsorption potential function for rare gas adsorpto the total dispersion energy. From the Auger data we have tion in silicalite includes repulsion, dispersion (two and three evaluated23the dipole silicon polarizability to be 0.38 A3 which body), and induction interactions and clearly requires extensive is about one-third of the value obtained for the oxygen calculation. For this reason most calculations describing the polarizability. As discussed in reference 23, the value of the interaction between one adsorbate molecule and the zeolite oxygen polarizability used in the Kiselev model is not compatcrystal rely heavily on effective potential functions which are ible with the partially covalent-ionic nature of the Si-0 bond mathematically simpler. in silicate. The silicon dipole polarizability will be negligible only if one considers the silicalite structure as purely ionic, The effective potential introduced by Kiselev et al. (referred which contradicts the starting assumption that the silicalite to here as the Kiselev model) has been widely used in modeling crystal is a purely covalent structure. Moreover, we have seen adsorption in silicalite and related materials; it is therefore useful that the use of the total number of electrons of a given species to compare the results of the present work with this model. In this modelz8 the total interaction energy is represented by two as its effective number of electrons does not give accurate results terms: an attractive term and a repulsive term. The model for the dispersion interaction. Thus it appears that this model neglects all other types of interaction. This approximation was does not accurately describe the electronic properties of the justified by the fact that dehydroxylation of silica gel surface framework species. From eq 46, it is expected that a small sharply diminishes the adsorption of polar molecules; the change in the van der Waals radii would have a significant effect silicalite channels are considered as hydroxyl-free surfaces. on the repulsive constant and consequently on the whole Furthermore, the neglect of the induction interaction (and of adsorbate-zeolite interaction. Some authors have used the the electrostatic interaction for a molecule such as nitrogen) Kiselev model with a slightly different set of atomic parameters; implies that the framework atoms do not carry any partial for example June et aL31 fitted the silicalite oxygen van der charge: the silicalite crystal is then considered as a purely Waals radius and found a value of 1.58 A. 0.291 1.346 5.126 14.11 0.622 2.832 10.29 27.10 1.769 8.292 34.35 97.66 0.077 1.514 19.14 1.066 159.5 12.43
0.295 1.427 5.543 15.20 0.675 3.220 11.95 31.80 1.862 9.141 37.55 107.3 0.075 1.588 21.02 1.109 171.4 13.52
1.948 12.10 71.7 311.1 5.05 30.0 163.9 654.1 26.10 142.3 753.5 2680 0.273 12.32 352.8 11.41 5404 271.8
2.084 12.35 70.99 297.6 4.875 28.54 157.8 637.5 23.24 130.3 704.7 2776 0.269 10.70 325.5 10.42 5086 241.0
Physical Adsorption of Rare Gases in Silicalite
J. Phys. Chem., Vol. 98, No. 50, 1994 13345
5. The PN1 Intermolecular Potential for the Adsorption of Argon in Silicalite The Induced Interaction. This is rarely a predominant source of attraction between molecules and can be often neglected in comparison to other contributions.6 The dipole polarizability term can be written
TABLE 5: Parameters (in atomic units) Used in Calculating the Ar-Silicalite Potential Two-Body Coefficients c 6
CS
ClO
A b
Ar-0
Ar-Si
50.15 966 2.21 104 711.6 2.045
15.1 26 1
-
1422 2.090
Three-Body Coefficients Ar-0-0
where E , j is the electrostatic field due to the frameworkj-th oxygen atom at the position i occupied by the argon and Ei,k is the same quantity for silicon atoms. The framework atoms carry partial charges of - l e and f 2 e for the oxygen atoms and the silicon atoms, re~pectively,~~ as previously mentioned. The Dispersion Interaction. The dispersion interaction is considered as a multipole expansion which includes damped two-body and three-body terms. For the argon-oxygen interaction which is expected to give the larger contribution, the twobody expansion includes terms up to the Clo coefficient. In the argon-silicon expansion, terms up to the c8 coefficient are included. The three-body expansion contains all terms based on dipoles and quadrupoles (and the relevant permutations) and includes contributions from the two triplets {Ar,O,O} and {Ar,O,Si}. Table 5 gives the two-body and three-body dispersion coefficients calculated by the method described earlier. It is interesting to note that eq 40 gives a 2 O = 24.9 au in good agreement with azo= 27 au from modified electron gas calculations on quartz.33 The Repulsive Interaction. The four repulsive parameters {Ah,o,bh,oeh,Si,bh,si} required by eq 27 were obtained here by fitting to experimental data including the isosteric heat of adsorption qst in the limit of zero coverage and the Henry law constant. The isosteric heats of adsorption can be obtained experimentally from extrapolation of microcalorimetry measurem e n t ~or~from ~ adsorption isotherms.35 The first method is preferable because it is a direct measurement while the latter needs a set of adsorption isotherms over a wide range of temperature. The Henry law constant is obtained from isotherm slopes in the low-pressure region.
It is worth noting that it is often experimentally difficult to get low-temperature Henry law constant.36 These two quantities can be expressed in terms of the integrals:
I , = j u exp[-P u(r)]dr
(49)
Z2 = j e x p [ - b u(r)ldr
(50)
where u is the total potential function, and are given by
+
qst = -[Zl/Z2] RT;
kH = (Z2 - V)/Q,VRT
(51) where Vis the volume of a unit cell and eZis the zeolite number density. Calculations were performed over a unit cell embedded at the centre of 3 x 3 x 3 unit cells containing 5184 0 and 2592 Si atoms. The summations were cut off after 16.5 A. The Henry law constants and heats were found from 20 x 20 x 20 point Gaussian quadrature3' over l/4 of the central unit cell for the orthorhombic unit cell (taking advantage of the pnma symmetry) and over a full unit cell for the monoclinic structure. The validity of the integration method was checked against previously published r e s ~ l t s . 4 ~Details 3 ~ ~ are given in ref 58.
101 325 325 450 870 1.05 x 1.50 4.68 x -4.18 -3.05 -3.05
103 103 103 103 103 103
Ar-Si-0 30.2 71.8 97.2 135 232 320 434 1.03 x -1.25 x -9.09 x -2.71 x
103 103 103 103
SCF calculations for the adsorption of He on NaCl and LiF surface^^^,^^ provide a basis for imposing preliminary constraints
on {A,b} pairs such as the magnitude of the cationic repulsive parameters with respect to the anionic ones. A trial and error method controlled by four Henry law constants at four different temperatures was then used to obtain better repulsive parameters. Calculations at temperatures less than 340 K used a monoclinic structure while calculations above 340 K used an orthorhombic ~ i l i c a l i t e .It~ ~is clear from the mathematical expression used for the repulsive interaction that the total energy is very sensitive to small variations of the b parameters which then can be considered as coarse parameters. Consequently, the A parameters were adjusted in order to fine tune the repulsive energy. In Table 6 we present the Gaussian quadrature integration results for argon in silicalite at zero coverage. The repulsive parameters for the PN1 potential function found from the fit of the four Henry law constant written in bold in Table 6 are = 711.6 Eh, bh9O = 2.045~0-', and Ah,'' = 1422 Eh, b*ssi = 2.09 a0-l. The Kiselev parameters for the Ar/silicalite system are ~ h - 0= 130.3 K and oh-0 = 3.029 A. The calculated kH values obtained with the PN1 function at 195, 306, 323, and 358 K are in good agreement with the experimental data while the theoretical kH value at 423 K is underestimated. It should be mentioned that the experimental points on the isotherm, obtained by gravimetry at 423 K, are scattered because the amount of adsorbed argon at this temperature is tiny and therefore hard to detect. Henry law constants are also hard to measure at low temperature because the linear region is in a range of pressure not easily detectable by most pressure devices. The values predicted for qst are also in good agreement with experiment. It is clear from Table 6 that the PN1 model gives better results than the Kiselev model over a wide range of temperature. The Henry law constants obtained with the Kiselev model are consistently overestimated compared to the experimental data by a factor 3 or 4, while the qst values (although acceptable) are overestimated at low temperature and underestimated at high temperature. Having obtained A and b for heterogeneous interactions, the following combination rule, based on an SCF study4 of repulsive interactions between nonbonded species A and B, was used to calculate the corresponding constants for oxygen and silicon pairs. Ahqo
u;;(R) = (AAA B)112 exp[-2bAbBR/(bA
+ bB)]
(52)
13346 J. Phys. Chem., Vol. 98, No. 50, 1994
Pellenq and Nicholson
TABLE 6: Comparison of Theoretical and Experimental Results at Zero Coverage for Kiselev Potential (KI)and PN1 Potential (ThisWork) k&nmol g-I Pa-' qsJd mol-' structure temp/K KI PNl expt KI PN 1 expt * 14.8 14.3 14.4 monoclinic 77 91.4 19.1 monoclinic 87 7.10 0.42* 1.67 monoclinic 195 9.6(-5) 3.69(-5) 3.71(-5) monoclinic 306 34.0(-7) 16.3(-7) 14.9 15.4 16.0 f 0.9 15.5(-7) monoclinic 323 25.0(-7) 9.72(-7) 6.83(-7) orthorhombic 358 14.4(-7) 4.95(-7) 4.27(-7) orthorhombic 423 6.43(-7) 1.53(-7) 2.18(-7)
400
t'
I
I
I
I
I
I
I
-1
O t Y
-400
-
-800
-
-1200
-
'1
\ h
P w
400 0
53
400
F c k!+-800 -1200
-1600
I
I
I
I
I
I
1
I
-1600
-10
-5
0
5
10
y-axis I A Figure 4. Energy cross section along the straight channel (orthorhombic structure): (0)total PN1 energy, (0)2BArO, (V)2BArSi, (0) induced energy, (A) 3BArO0, (0) 3BArOSi, (W) Kiselev potential energy.
Using the A r - h repulsive parameters' Ah = 328.85 Eh and bh = 1.918 a-',this expression gives Ao = 1543.5 Eh , bo = 2.190 a-'; As' = 6163.4 Eh, bsi = 2.395 a-'.
6. Energetics and Adsorption Site Distribution for Argon in Silicalite The whole potential surface is a rather complicated object which cannot be represented in two dimensions. Its global properties are reflected in the values of qst(0)and kH already discussed. In order to bring out some of the factors responsible for the differences between the PN1 and the Kiselev potentials we illustrate some aspects of the potential surface in Figures 3-5. Figures 3 and 4 show the variation of potential energy across the centre of the straight channel and along the straight channel, respectively. Separate contributions from two-body Ar-0 Born-Mayer repulsion (2BArO), the corresponding Ar-Si term (2BArSi), three-body terms (3BArO0, 3BArOSi), and induced energy are shown. The main contribution to the total potential, as expected, comes from the 2BArO term, which gives values within 10% of the total potential for all positions of the Ar probe. The 2BArSi contribution is attractive and typically about 12% of the total but gives a similar contour pattern to that of the Ar-0 two-body potential. Each threebody terms reflects the anisotropy of the system: when the probe is close to the center of a pore, the terms involving quadrupoles almost exactly cancel the three-body fourth-order term, as is found to be the case for uniform When the probe is close to the cavity surface, however, this cancellation does not occur and the total three-body interaction can differ by as much as 20% from the triple dipole terms. In total, the 3BArOO term
+
almost cancels out with the 2BArSi term when the adsorbate is not near to the cavity walls. The sign of the three-body interaction depends on the geometrical configuration of a triplet: for example, the interaction is attractive if the atoms are on a line. The induced energy is also strongly anisotropic in this sense; however, as anticipated, its contribution is very small. There are two striking differences between the PN1 and the Kiselev potential: (i) The PN1 pore is narrower (with respect to Ar). The PN1 energy profile exhibits a single rather than a double minimum across the narrowest section of the straight pore (Figure 3) because the framework atoms are more polarizable. Thus the dispersion part of the PN1 function is more attractive than in the Kiselev function. However, the cross section at the center of the pore shows that the 2BArO term is less attractive than the Kiselev term. Clearly the repulsion interaction is much stronger in the PN1 function than in the Kiselev function in agreement with preliminary semiempirical quantum mechanical c a l c ~ l a t i o n . ~(ii) ~ . ~The ~ PN1 potential profile along the axis (Figure 4) is more complex than the Kiselev potential. Potential maps for the straight channel are shown in Figure 5 for the PN1 model and for the Kiselev model. They are in general agreement with those already p ~ b l i s h e d ~for~ argon ,~~,~~ in silicalite using the Kiselev model and similar adsorption site distributions have been obtained for methane in s i l i ~ a l i t e . ~ ~ , ~ ~ Thus it appears that adsorbates having the same size and symmetry will have similar distributions of sites in the zeolite channels; in other words, the topology of the site distribution is imposed by the zeolite geometry, while the actual energy depth and shape of each site are a consequence of the electronic
Physical Adsorption of Rare Gases in Silicalite
J. Phys. Chem., Vol. 98, No. 50, 1994 13347
Figure 5. Potential maps for argon in a straight channel of silicalite-1: (a) PN1 model, (b) Kiselev model. The length of each map is 1.6 nm and width 0.78 nm; gradations on the horizontal scale are at 0.0332 nm. Energylk is given in K, and the contours are separated by 100 K.
properties of both the adsorbate molecule and the framework atoms. The deepest adsorption sites occur at the centre of the straight channel, while the slightly more open regions near the channel intersections have a small maximum in this cross section. Again, it is clear from these potential surfaces that the PN1 model describes narrower cavities than the Kiselev model. We have seen that the Kiselev potential overestimates Henry law constants. This effect is confirmed and amplified with more polarizable adsorbates (see below). Equation 51 shows that the Henry law constant is roughly proportional to the integral of the Boltzmann probability over the volume of the unit cell. The fact that the Kiselev model overestimates the Henry law constant at a given temperature in comparison to the PN1 model can be written as
where ui is the adsorption energy at the position i of the adsorbate. This expression will be true if, for some argon positions, the adsorption energy is more negative than the energy obtained with the PN1 model. Another possibility is that the Kiselev model generates more adsorption sites than the PN1 function. The similar topology of the energy map in both cases rules out this latter possibility. Inspection of Figures 5 and 6 indicates that the Kiselev potential has lower adsorption energies than the PN1 potential. Consequently, one may conclude that the Kiselev function is too attractive. This is despite the fact that the attractive part of the Kiselev potential function is obtained from a smaller value of the oxygen polarizability. Clearly, the repulsive part of the Kiselev potential is too weak. The repulsive parameter in the Kiselev potential is given by eq 46 from which it is apparent that the defective repulsion energy would increase if the sum of van der Waals radii was increased. Since the argon van der Waals radius is known accurately the Kiselev repulsive energy can only increase if the van der Waals
13348 J. Phys. Chem., Vol. 98, No. 50, 1994 TABLE 7: Parameters (in atomic units) Used in Calculating the Kr-Silicalite Potential Two-Body Coefficients Kr-0 Kr-Si c 6 71.4 26.7 1.19 x 103 399 c.8 ClO 4.15 x 104 A 603.8 1207 b 1.907 1.98I Three-Body Coefficients Kr-0-0 Kr-Si-0 55.6 104 214 226 333 1.06 103 382 1.61 x 103 -3.44 x 103 -1.71 x 103 -653 TABLE 8: Parameters (in atomic units) Used in Calculating the Xe-Silicalite Potential Two-Body Coefficients Xe-0 Xe-Si c 6 103.7 39.0 1.94 x 103 648 c.8 ClO 7.94 x 104 A 1.26 x 103 2.51 x 103 b 1.932 2.007 Three-Body Coefficients Xe-0-0 Xe-Si-0 83.8 220 152 877 339 877 1.02 x 103 393 5.47 x 103 523 2.09 x 103 5.39 x 103 637 5.39 x 103 3.12 x 103 3.76 x 104 -8.14 x 103 -2.12 x 104 -2.61 x 103 -6.83 x 103 -1.00 x 103 -6.83 x 103 radius of the framework oxygen is increased. This will inevitably reduce the accessible space for the adsorbate in the zeolite channels in line with the result obtained with the PN1 function. June et al., using the Kiselev model for the dispersion interaction, adjusted the van der Waals radius of the zeolitic oxygen atom to a value of 1.58 A in order to obtain zero coverage results comparable with experiment for the adsorption of xenon, c&, sF6, and in silicalite. Recent periodic ab initio calculation^^^ have shown that the radius of the oxygen atom in quartz is about 2 A, much larger than the tabulated van der Waals radius (1.52 A) for the neutral and isolated oxygen atom29 used in the Kiselev model. Effective potential functions developed for alkanes in silicalite5' also have framework species similar in size to those obtained from ab initio calculations. Thus it appears that the method of calculating the repulsion interaction in the Kiselev model is unsuitable. Another consequence of increasing the repulsive part of the interaction potential is that it produces less corrugated potential surfaces. The omission of the silicon atoms, included in the full calculation, also contributes to this effect. It is clear that the correct energy surface in confined geometry such as zeolites cannot be obtained by using atomic sizes of free species. Calculations with the PN1 model show that framework symmetry is significant for the calculation of the Henry law
Pellenq and Nicholson constants but negligible for the isosteric heat. The results obtained with the monoclinic structure and those obtained with the orthorhombic structure differ by as much as about 70% for k~ but less than 3% for qst;the orthorhombic structure always gives smaller absolute values in both cases. The differences decrease as the temperature and the adsorbate size increase. This illustrates the influence of symmetry losses between the orthorhombic structure and the monoclinic crystal due to small atomic displacements. However, calculations using the Kiselev model52have shown that at room temperature small variations of some atomic positions in the orthorhombic structure do not affect the values of the Henry law constants for methane in silicalite while some large variations are obtained for large molecules such as p-xylene. The effect of the change of geometry on the calculations of the Henry law constants is due to the closer confinement in the PN1 model so that small variations in the position of framework species lead to significant modifications of the energy surface.
7. Intermolecular Potential for Adsorption of Krypton in Silicalite The dispersion coefficients for krypton were obtained by the method outlined earlier. The repulsive parameters were found from the combination rules in eq 52 with the Kr-Kr repulsive parameters determined from the repulsive part of the effective Lenuard-Jones potential6 as A&*&= 236.2 Et,,and bKr,& = 1.689 The PN1 parameters are listed in Table 7. The parameters for the Kiselev model are28EG,O = 158.62 K and OQ,O = 3.154 A. The isosteric heat of adsorption at 77 K calculated with the PN1 potential function on the monoclinic phase (17.6 kJ mol-') compares well with its experimental counterpart (17.4 kJ mol-') measured by ~ a l o r i m e t r y .Again ~ ~ the value obtained with the Kiselev potential (20.1 kJ mol-') overestimates this property. The single experimental valueJ4reported for kH = 1.86 x mmol g-' Pa-' at 250 K compares with k~ = 1.20 x mmol g-' Pa-' (PN1) and kH = 15.6 x mmol g-' (Kiselev model). The PN1 results are particularly gratifying here in view of the fact that there are no adjustable parameters. 8. Intermolecular Potential for Adsorption of Xenon in Silicalite The potential Parameters for xenon were obtained in the same manner as for krypton. The repulsive constants for the XeXe interaction are6AxeSxe= 1025.7 Eh and bxe,xe = 1.728 a-'. The PN1 parameters are listed in Table 8. The parameters for the Kiselev model are cxe,0 = 209.07 K and axe,0 = 3.278 A. The isosteric heat of adsorption at 121 K calculated with the PN1 potential function on the monoclinic phase is 27.4 kJ mol-' which compares well with its experimental counterpart (26.6 & 1. kJ mol-') obtained from the isosteric method.55 The calculated value obtained with the Kiselev potential (31.0 kJ mol-') again significantly overestimates the experimental heat. Two experimental values have been reported56 for the Henry law constant: k~ = 2.2 x low2mmol g-' Pa-' at 195 K which compares with k~ = 3.74 x lom2mmol g-' Pa-' obtained with the PN1 potential and k~ = 87.6 x mmol g-' for the Kiselev potential and k~ = 2.25 mmol g-' Pa-' at 273 K mmol g-' Pa-' and k~ = compared with k~ = 2.90 x 39.0 x mmol g-' Pa-' from the PN1 and Kiselev models, respectively. Again there are no adjustable parameters used in the PN1 function. 9. Summary and Conclusions We have presented a method to calculate two- and threebody dispersion interaction for isolated and crystal species which requires only a knowledge of dipole polarizability and effective
Physical Adsorption of Rare Gases in Silicalite number of electrons for the interacting species. Elsewhere we have shown how the former property can be found for lattice species using Auger data,23and in this paper we give a method for obtaining ne^. These methods accurately reproduce available data for isolated and crystal species. Also, we have applied the general form of intermolecular interactions to the adsorption of one molecule in silicalite. This model includes all the main features required by the general theory of intermolecular interactions. The four parameters for the Born-Mayer form of repulsive interaction between argon and the zeolitic structure have been obtained from experimental Henry law constants at four temperatures. Both the full scale potential function (PN1) for argon and the corresponding Kiselev potential function have been used to calculate zero-coverage thermodynamic quantities not used for fitting. It was found that the PN1 function performed well compared to the Kiselev function which overestimates experimental values. It is shown that the PN1 model describes a narrower and softer zeolite wall than the Kiselev model. The analysis of the contributions to the total energy has shown that the overall potential is govemed by the two-body dispersion-repulsion term for the argon-framework oxygen pair. The three-body contributions almost cancel out with the two-body dispersion-repulsion term for the argonsilicon pair in the center of the pore but are not negligible on the edge of the cavity and reflect the anisotropy of the system. The induction energy is negligible. Combination rules for the repulsive parameters were used to derive full scale potential functions for krypton and xenon in silicalite. The PN1 function is again able to predict available zero coverage data accurately in contrast to the Kiselev model which again tends to overestimate the experimental data. The calculations repeated on both silicalite structures (monoclinic below 340 K and orthorhombic above) show that this structural phase transition has a significant effect on the PN1 calculations of the Henry law constants reflecting the effect of the confinement but does not affect results for the isosteric heat of adsorption. The relatively minor contributions from three-body and induction terms suggest that an effective two-body model might be found which would provide a reasonably accurate potential for adsorbate silicalite interactions. However the Kiselev model, as anticipated by its authors,28can only be regarded as a first approximation since it is not able to reproduce the narrower pores predicted by the full potential model. The present work also highlights the fact that accurate prediction of qst(0) is a necessary but not a sufficient condition in evaluating a potential function, since this property is largely determined by the depth of the sites but less by the shape of the potential. Elsewhere the results of grand canonical Monte-Carlo simulations are reported in which the PN1 model has been used to study adsorption of rare and nitrogen58 in silicalite, argon and methane in AlP04-5,59 and methane and xenon in Acknowledgment. We thank the CEC for providing equipment and support to R.J-M.P. References and Notes (1) Bakaev, V. A. Surface Sei. 1988,198,571. Bakaev, V. A,; Steele, W. A. J. Chem. Phys. 1993, 98, 9922. (2) Schindler, H.; Vogelsang, R.; Staemmler, V. ; Siddiqi, M. A,; Svejda, P. Mol. Phys. 1993, 80, 1413. (3) Buckingham, A. D. Adv. Chem. Phys. 1967, 12, 107. (4) Claverie, P. Intermolecular Interactionsfrom Diatomics to Biopolymers; Pullman, B., Ed.; Wiley: New York , 1978; Chapter 2. ( 5 ) Tang, K. T.; Norbeck, J. M.; Certain, P. R. J. Chem. Phys. 1976, 64, 3063. (6) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, N. A. Intermolecular Forces, Clarendon Press: Oxford, U.K., 1981. (7) Bell, R. J. J. Phys. B 1970, 3, 751.
J. Phys. Chem., Vol. 98, No. 50, 1994 13349 (8) Macrury, T. B.; Linder, B. J. Chem. Phys. 1970, 54, 2056. (9) Tang, K. T.; Toennies, J. P. J. Chem. Phys. 1984, 80, 3726. (10) Pitzer, K. S . Adv. Chem. Phys. 1959. (11) Buckingham, A. D.; Fowler, P. W.; Hutson, J. M. Chem. Rev. 1988, 88, 963. (12) Thakkar, A. J.; Hettema, H. ; Wormer, P. E. S. J. Chem. Phys. 1992, 97, 3252. (13) Spackman, M. A. J. Chem. Phys. 1991, 94, 1295. (14) Halgren, T. A. J. Am. Chem. SOC. 1992, 114, 7827. (15) Cambi, R.; Cappelletti, D. ; Liuti, G. ; Pirani, F. J. Chem. Phys. 1991, 95 1852. (16) F'yper, N. C. Philos. Trans. R. SOC. London 1986, A320 107. (17) Caillet, J.; Claverie, P. Acta Crystallogr. 1975, A31, 448. (18) Fowler, P. W. Mol. Simulation, 1990, 4, 313. (19) Koutselos, A. D.; Mason, E. A. J. Chem. Phys. 1986, 85, 2154. (20) Kiselev, A. V.; Poskus, D. P. Zur Fiz. Chim. 1958, 32, 2854. (21) Narayan, R. J. Phys. Chem. Solids 1977, 38, 1097. (22) Varandas, A. J. C.; Dias da Siva, J. J. Chem. SOC.,Faraday Trans. 2 1986, 82, 593. (23) Pellenq, R. J-M.; Nicholson, D. J. Chem. SOC.,Faraday Trans. 1993, 89, 2499; Pellenq, R. J-M.; Nicholson, D. Stud. Surf: Sei. Catal. 1994, 87, 31. (24) Fowler, P. W.; F'yper, N. C. Mol. Phys. 1986, 59, 317. (25) Pellenq, R. J-M.; Nicholson, D., manuscript in preparation. (26) CRC Data Book for Chem. Phys. 1993. (27) Fowler, P. W. ; F'yper, N. C. Proc. R. SOC. London A 1985, 398, 377. (28) Kiselev, A. V.; Lopatkin, A. A.; Shulga, A. A. Zeolites 1985, 5, 261. (29) Bondi, A. J. Phys. Chem. 1964, 68, 441. (30) Pickett, S. D.; Nowak, A. V.; Thomas, I. M.; Peterson, B. K.; Swift, J. F.; Cheetham, P. A. K.; den Ouden, C. J. J.; Smit, B.; Post, M. F. M. J. Phys. Chem. 1990, 94, 1233. (31) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 1508. (32) White, J. C; Hess, A. C. J. Phys. Chem. 1993, 97, 6398, 8703. (33) Jackson, M. D.; Gibbs, G. V. J. Phys. Chem. 1988, 92, 540. (34) Rouquerol, F.; Rouquerol, J.; Everett, D. H. Thermochim. Acta 1980, 41, 311. (35) Barrer, R. M. Zeolites and Clay Minerals; Academic Press: London, 1978. (36) Grillet, Y., private communication, 1991. (37) Press, W. H.; Flanery, B. P.; Turkolsky, S . A,; Vetterling, W. T. Numerical Recipes; Cambridge University Press: Cambridge, U.K., 1989. (38) Hutson, J. M.; Fowler, P. W. Surf: Sei. 1986, 173, 337. (39) Fowler, P. W.; Hutson, J. M. Phys. Rev. B 1986, 33, 3724. (40) Reichert, H., private communication, 1991. (41) van Koningsveld, H.; Jansen, J. C.; van Bekkum, H. Zeolites 1987, 7, 564; 1990, 10, 235. (42) El Amrani, S.;Vignt-Maeder, F.; Bigot, B. J. Phys. Chem. 1992, 96, 9417. (43) Goodbody, S. J. ; Watanabe, K.; MacGowan, D.; Walton, J. P. R. B.; Quirke, N. J. Chem. Soc., Faraday Trans. 1991, 87, 1951. (44) Bohm, H. J.; Ahlrichs, R. J. Chem. Phys. 1982, 77, 2028. (45) Barker, J. A. Phys. Rev. Lett., 1986, 57, 230. (46) Pellenq, R. J-M.; Pellegatti, A,; Nicholson, D.; Minot, C., manuscript in preparation. (47) Pan, D.; Mersmann, A,, Fundamentals of adsorption; Mersmann, A., Scholl, S . , Eds.; United Engineering Trustees, 1991; p 645. (48) June, R. L; Bell, A. T.; Theodorou, D. N. J. Phys. Chem., 1990, 94, 8232. (49) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1991, 95, 8866; 1992, 96, 1051. (50) Dovesi, R.; Roetti, C.; Freyria-Fara, C.; Apro, E.; Saunders, V. R.; Harrison, N. M. Philos. Trans. R. SOC. London 1992, A341, 203. (51) Nicholas, J. B.; Trouw, F. R.; Mertz, J. E. Iton, L. E.; Hopfinger, A. J. J. Phys. Chem., 1993, 97, 4149. (52) Li, J.; Talu, 0. J. Chem. Soc., Faraday Trans. 1993, 89, 1683. (53) Llewellyn, P. L.; Coulomb, J.-P.; Grillet, Y.; Patarin, J.; Lauter, H.; Reichert, H.; Rouquerol, J. Langmuir 1993, 19, 1846. (54) Ruthven, D. M.; Tezel, F. H.; Devgun, J. S . ; Sridhar, T. S . Can. J. Chem. Eng. 1984, 62, 526. (55) Biilow, M. ; Hirtel, U.; Muller, U. ; Unger, K. K. Ber. Bunsenges. Phys. Chem. 1990, 94, 74. (56) Chen, D. Ph.D. Thesis, Imperial College, University of London, 1992. (57) Pellenq, R. J-M.; Nicholson, D. Langmuir, in press. (58) Pellenq, R. J-M., Ph.D. Thesis, Imperial College, University of London, 1994. (59) Boutin, A,; Pellenq, R. J-M.; Nicholson, D. Chem. Phys. Lett. 1994, 219, 484. (60) Pellenq, R. J-M.; Fuchs, A,; Espinat, D.; Tavitian, B., to be published.
