J . Phys. Chem. 1990, 94, 4329-4334 only the presence of parent FeS and a sextet that may be ascribed to a-iron. These data suggest that electrochemical reduction of FeS in these cells proceeds rapidly to iron and that any intermediate formed is present, at any one time, in only very small amounts: FeS
Li
[LiFeS]
Li
FeO
+ LizS
The only species observed in the cell during discharge are parent FeS, very small iron particles that exhibit superparamagnetism at room temperature, and larger iron particles that yield a magnetic sextet of intermediate field strength. To emphasize this point, the Mossbauer spectra of Li2FeSzat room temperature and 77 and 4.2 K are shown in Figure 6 . The preparation and characterization of Li2FeS2by single-crystal X-ray diffraction has been described previously.22 The room-temperature Mijssbauer spectrum shows an irresolved set of two doublets of unequal intensity. This is in agreement with the earlier work of Melandres and TanLz3 The 4.2 K spectrum is remarkably complex. The X-ray crystal structure of LizFeS2 shows the presence of four inequivalent iron sites in the ratio 3:3:1:1. The four sites are not resolved in the room-temperature Mossbauer spectrum. At 4.2 K the different sites apparently yield different hyperfine fields, leading to the observed complex spectrum. There is no evidence in the spectra of Figures 1 and 2 for the formation of LizFeSz in the FeS cathode during discharge.
J.
(22) Batchelor, R. J.; Einstein, F. W. B.; Jones, C. H. W.; Fong, R.; Dahn, 1988, 37, 3699. (23) Melandres, C. A,; Tani, B. J . Phys. Chem. 1978, 82, 2850.
R. Phys. Rev. B
4329
The results of experiments on chemically lithiated FeS are shown in Figures 4 and 5 . Here there is little evidence for the formation of superparamagnetic iron particles. For 1 equiv of lithium, the spectrum at room temperature shows unreacted FeS and the presence of a magnetic sextet corresponding to an intermediate field strength. When this sample is cooled to 4.2 K, this latter sextet shows a very significant change to a larger field strength, while the FeS peaks remain essentially unaffected. Thus, chemical reduction appears to produce particles of iron that are too large to exhibit superparamagnetism but that are small enough to yield an internal magnetic field significantly smaller than that for bulk a-iron at room temperature. When this is cooled to 4.2 K, a field close to that of a-iron is observed. It is worthy of note that the internal magnetic field for the iron particles at 4.2 K in the chemically lithiated FeS is again greater than that of bulk a-iron. This result was consistently obtained for different samples of FeS reduced with n-BuLi. Surface effects may again play a role here.
Acknowledgment. We thank Dr. R. R. Haering for suggesting this work, Dr. Ray Batchelor for his assistance in computational aspects, and Dr. Tom Birchall of the Department of Chemistry, McMaster University, Hamilton, Ont., for supplying us with a copy of the Mossbauer fitting routine GMFP. We are grateful to the Natural Sciences and Engineering Research Council of Canada for support of this work. Registry No. FeS, 1317-37-9; Li, 7439-93-2; LiAsF,, 29935-35-1; AI, 7429-90-5; Fe, 1439-89-6; Li2S, 12136-58-2; Li2FeS2,59217-78-6; LiFeS2,59217-77-5; Li,,5FeS2,79176-50-4; propylene carbonate, 108-32-7; cyclohexane, 110-82-7; n-butyllithium, 109-72-8; hexane, 110-54-3.
Molecular Dynamics Studies on Zeolites. 4. Diffusion of Methane in Silicalite Pierfranco Demontis, Ettore S. Fois, Giuseppe B. Suffritti,* Dipartimento di Chimica, Universitri di Sassari, Via Vienna 2, I-071 00 Sassari, Italy
and Simona Quartieri Istituto di Mineralogia e Petrografa. Universitci di Modena, Via S . Eufemia 19, I-41100 Modena, Italy (Receioed: June 5, 1989; In Final Form: October 23, 1989)
The diffusion of methane in silicalite was simulated by molecular dynamics using a simplified model but very long (0.2 ns) trajectories. The calculated diffusion coefficient, 6.58 X m2 s-l, resulted in good agreement with experiment [(6.5 1.0) X m2 s-l], and some details of the diffusive motion were evidenced. Structural changes induced by the sorbate and experimentally detected were also found in the simulated silicalite.
Introduction In some recent papers'-' the application of molecular dynamics (MD) to the simulation of structural and dynamic properties of zeolites was illustrated. In this laboratory, the positions and vibrations of water molecules in the cages of natrolite,'J the atomic coordinates and the crystal symmetry of dehydrated natrolite3 and Linde zeolite 4A," and their dynamical behavior were satisfactorily (1) Demontis, P.; Suffritti, G . B.; Alberti, A,; Quartieri, S.; Fois, E. S.; Gamba, A. Gazz. Chim. Ital. 1986, 116, 459. (2) Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A. In Dynamics of Molecular Crystals; Lascombe, J., Ed.; Elsevier: Amsterdam, 1987; p 699. (3) Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A. Zeolites 1987, 7 , 522. (4) Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A. J . Phys. Chem. 1988, 92, 867. (5) Shin, J. M.; No, K. T.; Jhon, M . S. J . Phys. Chem. 1988, 92, 4533. (6) Leherte, L.; Lie, G. C.; Swamy, K. N.; Clementi, E.; Derouane, E. G.; Andrt, J. M. Chem. Phys. Lett. 1988, 145, 237. (7) Yashonath, S.; Demontis, P.;Klein, M. L. Chem. Phys. Lett. 1988,153, 551.
reproduced by using simple model potentials. Fixed framework simulations were recently carried out by Leherte et aL6 on the diffusion of water in ferrierite, by Shin et al.s on the dynamics of the Na+ ions in A-type zeolite, and by Yashonath et al.' on the mobility of methane in zeolite N a y , which was studied also by Monte Carlo method.s In the present paper, MD calculations extended to the diffusive motion of methane in silicalite are described. The choice of this zeolite was suggested by the importance of ZSM-5 as molecular sieve and catalyst. Indeed, in the simulation no aluminum atom was included, and the simulated material was in fact silicalite. The same model potential (harmonic form) that was used for anhydrous natrolite and zeolite A was adopted also for the silicalite framework, while methane molecules were represented by point particles interacting with the framework and among themselves via suitable Lennard-Jones potentials. (8) Yashonath, S.;Thomas, J. M.; Nowak, A. K.; Cheetham, A. K. Nature 1988, 331, 601.
0022-3654f 90 f 2094-4329SO2.50f 0 0 1990 American Chemical Society
4330 The Journal of Physical Chemistry, Vol. 94, No. 10, 1990
In order to compare simulated and experimental results? the loading of the methane was kept relatively low (24 molecules per crystallographic cell, while the maximum loading can be estimated as about 40 molecules per cell). The motions of the methane particles were slow and highly irregular, so that statistics and averaging problems forced us to extend the length of the runs up to 200 ps, or 0.2 ns. However, a satisfactory agreement with the experimental values was finally reached. In order to get a better understanding of the influence of the framework motion on the diffusion process, a test with fixed framework was also performed. Some differences and analogies with the unconstrained system were evidenced, as will be discussed below. A simulation of the bare silicalite was also performed for reference. In spite of the simple model used, slight structural changes in the framework induced by the sorbate and detected by X-ray and spectroscopic experiments were found also in the simulated system. Method and Models The model potential proposed by the authors and used for MD calculations on anhydrous natrolite and zeolite A was adopted, in the harmonic form, also for silicalite. This model is described in detail in refs 3 and 4, and assumes that the potentials for S i 4 and 0-0 interactions are represented by quadratic functions of the displacement from a given equilibrium bond distance. No other possible contacts are included, the initial topology of the framework bonds is retained during the MD simulation, and only first neighbors are considered as interacting atoms. The parameters used for silicalite were the same as for zeolite A.
This model is certainly crude and may be interpreted as a second-order approximation of a Taylor expansion of a realistic potential which, when used in MD simulations, could be unmanageable. On the other hand, it has been shown in our previous works that it can describe reasonably the structural and dynamical features of zeolitic frameworks, so that the effects of the vibrations of the framework upon the diffusion of methane would be satisfactorily reproduced. In this view, no special physical meaning is attached to the harmonic model. As for methane molecules, they were represented by soft spherical particles, mainly because attention was focused on the general features of the diffusion, so that a model as simple as possible allowing long simulation runs was desirable. Methane-methane interactions were represented by a 20-6 Lennard-Jones potential between the centers of the molecules, derived by Matthews and Smithlo from experimental data. The minimum falls at 3.88 A and is -0.431 kcal mol-' deep. For the methane-framework potential, a slight modification of the one used by Ruthven and Derrah" for a transition-state theory study of the diffusion of CH4 in 5 A zeolite was adopted. Ruthven and Derrah assumed that the interaction between methane and the zeolite framework can be approximated by a pair potential (a Lennard-Jones function plus a spherical polarization energy term) between point particles representing methane molecules, cations, and the oxygen atoms (to which a charge of -e/4 was assigned) of the framework of the 5A zeolite. In the present work, as no cation was considered, the polarization term was small and was dropped. N o adjustment of the original parameters for both methanemethane and methaneframework potentials was attempted. Their analytical forms, giving energies in kilocalories per mole with distances in angstroms, are the following &H4-CH4 = 1.0992 X 10'lr-20- 2.09835 X 103r-6 (1) &H4+
= 2.28 X 1 O 6 F l 2 - 1.33
X IO3F6
(2)
where Of are the oxygen atoms of the framework. The value of ~~
~~
( 9 ) Caro, J.; HoEevar, S.; Kaerger, J.; Riekert, L. Zeoliles 1986, 6, 213. (10) Matthews. G. P.; Smith, E. B. Mol. Phys. 1976, 32, 1719. ( 1 1) Ruthven, D. M.; Derrah, M. I. J . Chem. SOC., Faraday Tram. 1 1972, 68, 2332.
Demontis et al. the minimum of the second potential is -0.194 kcal mo1-I for r = 3.885 A. In the MD calculations on silicalite all the atoms were left free to move without symmetry constraints under the action of the above-described potentials. The periodicity of the crystal was simulated by the usual minimum image convention and the equations of motion were integrated by means of modified Verlet's algorithm.I2 Structural analysis was performed separately for the zeolite framework and for methane particles. The coordinates of the Si and 0 atoms of the framework were referred at each step to the asymmetric unit by symmetry group operations and stored in order to derive their frequency distributions. These single-coordinate distributions represent an useful tool for structural analysis of simulated crystals. Indeed, if an atom maintains a mean position consistent with the given symmetry group, the distribution functions of its coordinates are Gaussian-like, with half-width close to the average displacements that can be derived from diffraction data. A bimodal (split) or multimodal distribution function for a coordinate means that the simulated structure is not consistent with the symmetry group operation on that coordinate and that probably it is ordered but with lower symmetry. Very large or asymmetrical distributions may arise from disorder; and diffusion is evidenced by flat distributions. Radial distribution functions (rdf) were not used for the structural analysis of the silicalite framework because in lowsymmetry crystals they can be difficult to be interpreted, but they were calculated to elucidate some features of the methane diffusion. Indeed, methane-methane rdf can yield information about dimers or clusters that can be formed in the cavities, and methane-Of rdf can explain some properties of the diffusive motion, like the average distance from the walls of the channels. A more detailed analysis of the diffusion mechanism was devised in order to elucidate the relative motion of the sorbed molecules: the distribution of dimers, trimers, ..., n-mers (up to n = 24) simultaneously present was evaluated, along with the distributions of the time of life of the n-mers from which the mean lives were derived. Moreover, the velocity autocorrelation function of the methane molecules, and the distribution of their velocities and temperatures was evaluated. The mean coordinates and temperature factors of the framework atoms were calculated following the same procedure used in refs 3 and 4. IR and power spectra were derived from the Fourier transform of the total dipole moment autocorrelation function and of the velocity autocorrelation function, r e s p e c t i ~ e l y . ~ . ~ The diffusion coefficient of methane in zeolite channels can be evaluated by using different procedures, which are in principle equivalent but numerically more or less efficient. The most used technique lies on the Einstein formula (3)
where t is time and r is the position of a methane particle. Equation 3 implies that mean square displacement in diffusive media is linear with time, so that D can be derived by the alternative expression 1 . d 6 I-- dt
D = - Iim - (Ir(t) - r(0)l2)
(4)
Finally, it is well-known that D may be calculated also by integrating the unnormalized velocity autocorrelation function
--
1
D = lim - J ' ( u ( O ) ~ ( t )dt) 3
0
(5)
but in this way no evident numerical improvement is obtained, unless the autocorrelation function comes rapidly to zero,I3 and (12) Swope, W. C.; Andersen, H. C.; Berens, C. H.; Wilson, K. R. J . Chem. Phys. 1982, 76, 637.
Diffusion of Methane in ZSM-5 in any case computing autocorrelation functions for long trajectories is often time and storage consuming. Beyond the technical problems of computing the diffusion coefficient, some attention must be paid to the physical meaning of the Einstein formula, from which eqs 3, 4 and 5 are derived, for the diffusion of molecules in zeolites. This formula is valid for Brownian motion in a tridimensional homogeneous medium, but the molecules diffusing in zeolite channels and cavities are constrained to move in a biased or hindered way, so that isotropic diffusion coefficients cannot be defined microscopically. On the other hand, in most cases experimental measurements of diffusion coefficients are performed in polycrystalline or powdered samples by means of various techniques giving an average over all directions and on a scale for which the internal structure of the zeolite is not resolved. In other words, these experimental techniques “see” the diffusion in zeolites as a continuous medium. Moreover, some of them (e.g., pulsed N M R ) evaluate directly the orientationally averaged mean square displacements of the diffusing molecules, and the diffusion coefficient is approximated with the Einstein formula, as the time intervals relative to such displacements are known. It may be concluded that the use of the above-mentioned formulas to obtain the diffusion coefficients from MD simulations could be justified if the only meaning attributed to that quantity is the possibility of comparison with experimental data. This problem will be the object of further study in future.
MD Simulations Single-crystal X-ray studies of ZSM-5 structure were carried out by Kokotailo and c o - w ~ r k e r and s ~ ~by~ Lermer ~~ et a1.16 The resulting diffraction patterns are consistent with the orthorhombic Pnma space group. However, it was suggested16that some disorder around m could be present in the structure. For MD simulations, a system corresponding to two crystallographic cells, superimposed along c, with cell parameters a = 20.076 A, b = 19.926 A, and c = 13.401 A was used. The framework atoms were 576, and no aluminum atom was included. This seems a reasonable assumption, as ZSM-5 is a high-silica zeolite with Si/AI ratios sometimes larger than 80. In this way no counterion had to be included in the structure. The spherical particles representing CH, molecules were initially located in positions occupied by CH3 and CH2 groups in tetrapropylammonium ZSM-5, whose single-crystal structure was studied by van Koningsveld et a1.I’ and by Chao et a l l s In order to compare simulated and experimentalg diffusion coefficients, three methane molecules per channel intersection were added to the MD system, i.e., 12 molecules per unit cell, resulting in a total of 600 particles. The time step used in MD runs was 1 fs. This small value was chosen in order to ensure a good energy conservation, in view of the need for very long trajectories, which had been stated in a previous work.’ Indeed, the fluctuations of total energy were less than 0.1%. Four MD simulations 200 ps long were performed: the first for the bare silicalite, with empty channels, was run at 280.5 K in order to verify the ability of the harmonic model to reproduce the most relevant structural and dynamical properties of this zeolite; the second was run for the full system, with 24 methane particles, at a temperature of 298.1 K; the third was run with a monoclinic cell (see below) at 302.6 K as a consequence of the ( 1 3 ) Zwanzig, R.; Ailawadi, N. A. Phys. Rev. 1969, 182, 280. (14) Kokotailo, G.T.; Lawton, S. L.; Olson, D. H.; Meier, W. M. Nature (London) 1978, 272, 437. ( 1 5 ) Olson, D. H.; Kokotailo, G.T.;Lawton, S. L.; Meier, M. V. J . Phys. Chem. 1981, 85, 2238. (16) Lermer, H.; Draeger, M.; Steffen, J.; Unger, K. K. Zeolites 1985,5, 131. (17) van Koningsweld, H.; van Bekkum, H.; Jansen, J. C. Acta Crystallogr. B 1987, 43, 127. (18) Chao, K-J.; Lin, J-Ch.; Wang, Y . ; Lee, G. H. Zeolites 1986, 6, 35. (19) Demontis, P.; Suffritti, G.B.; Quartieri, S.; Fois, E. S. Manuscript in preparation.
1rhe Journal of Physical Chemistry, Vol. 94, No. 10, 1990 4331
results of the second run; the last was carried out with fixed framework in order to check the differences between mobile and fixed framework models at 298.4 K. The ratio of the computing time per step for the first, second, and fourth runs was 1:3.25:2.2. It appears that, for the harmonic model, letting the framework move required only about 50% more time with respect to keeping the framework fixed.
Results and Discussion The results of the MD simulation of the bare silicalite assuming the Pnma space group yield mean values of the calculated coordinates in good agreement with the experimental ones, the standard error being 0.008 A. On the other hand, the distribution functions of the coordinates are not all satisfactory: though they are all unimodal, some of them, for they coordinate, are asymmetrical, leading to negative values of diagonal elements pZZof the anisotropic temperature factors matrixZo8. In particular, this finding is apparent for they coordinates of the four oxygen atoms (023 - 026) lying on the mirror plane normal to b at y / b = 1/4, which is characteristic of the space group Pnma. Moreover, the calculated anisotropic and isotropic temperature factors are generally slightly larger than the corresponding experimental ones, while they should be smaller, as was shown in ref 4. It may be concluded that the simulated structure belongs to the Pnma group but presents also some degree of disorder about m, as assumed by experimentalists.I6 The values of the calculated and experimental coordinates and temperature factors, along with a more detailed discussion, will be reported elsewhere. Once the methane molecules were introduced and a new MD run was started, the structure of the silicalite framework was monitored, and it seemed to remain close to the empty one for more than 100 ps; after the distribution functions of t h e y coordinate of several atoms began to split, and finally, starting from about 150 ps, all y coordinate distributions became bimodal. As these distributions were calculated by assuming that the space group was Pnma, it is evident that the simulated silicalite-methane system presented a symmetry lower than orthorhombic. In particular, the mirror plane normal to b was apparently lost, and it was argued that a space group compatible with the simulated coordinates would be P2,/n, as suggested by Lermer et aLl6 Sorbate-induced change of the crystal symmetry of ZSM-5 from orthorhombic to monoclinic space groups has been observed,2’,22 though with minor shifts in cell parameters (in particular /3 angle was reported in the range 90.4-90.6°),2’ and it was likely that the MD simulation was somewhat able to reproduce this effect. In order to verify better this finding, a monoclinic unit cell was assumed and the angle 0was optimized by minimizing the total energy of the system. The best value for /3 was 90.4O, in agreement with experiment. A new MD run 200 ps long was started at room temperature (302.6 K). Structural analysis was performed assuming P2,/n space group and, even though all coordinate distribution functions were symmetrical with respect to the mean value, some disorder was still present, and probably the true symmetry group of simulated system is only Pi. Nevertheless, the slight deformation of the unit cell does not affect the diffusion of the methane, and the study of this process was performed using the results of the last M D run, with monoclinic cell. In this paper emphasis is given to the diffusion process, and further details on MD studies of structural properties of silicalite with or without sorbed methane will be discussed in a work in preparation.ls It is to be remarked only that the simulated IR spectrum of the bare silicalite was in reasonable agreement with the experiment and that no detectable intensity was present in it below 300 cm-’. (20) Willis, B. T. M.; Pryor, A. W. Thermal Vibrations in Crystallography; Cambridge University Press: Cambridge, U.K., 1975. (21) Wu, E. L.; Lawton, S. L.; Olson, D. H.; Rohman Jr., A. C.; Kokotailo, G. T. J . Chem. Phys. 1979, 83, 2777. (22) Fyfe, C. A,; Kennedy, G.J.; De Schutter, C. T.; Kokotailo, G.T. J . Chem. SOC.,Chem. Commun. 1984, 541.
4332 The Journal of Physical Chemistry, Vol. 94, No. 10, 1990 500
Demontis et al. 5001
1
i
400-
r, A 30011
I
I
3001
’
h
I
J
/
0 1
I
l
200t
4
v
h
Y
V
I
B
root
/ ,,./”
I
X
I
~50 75 100 125 time ( p s ) Figure 1. Mean square displacements of the methane molecules (A2)vs time (ps) for simulations with mobile framework (continuous line) and with fixed framework (dashed line). The smoothness of the curves is due to the length of the MD run, entailing an average over several thousands of time origins.
“0
25
The experimental data on the diffusive process of methane in silicalite are essentially of two kinds: evaluation of the diffusion coefficients and calorimetric measurement of the sorption heat. Self-diffusion coefficients of methane in silicalite and in several samples of ZSM-5 with various aluminum content have been determined by Caro et aL9 by the N M R pulsed field gradient technique. At room temperature and at a sorbate concentration of 3 molecules per channel intersection, they were of the order of (6.5 f 1.0) X m2 S-I, depending on the different crystal size and origin of the samples. Sorption heat of methane in H-ZSM-5 was measured by Papp et aLZ3and, more recently, by VignE-Maeder and A u r ~ u for x~~ silicalite, and are quoted as 6.7 and 5.0 kcal mol-I, respectively. From the calculations reported in this work, a value of 4.5 kcal mol-l was obtained. This value includes the correction (0.1 kcal mol-]) for the truncation of the CH4-Of potential at 7 A, estimated by extending the cutoff distance to about 30 A. The calculated sorption heat compares favorably with the data of ref 24 for silicalite and can be considered reasonably good, because the data of ref 23 concern the H-ZSM-5, where extra polarization energy due to the acidic sites is present. In order to discuss the results of the MD simulation of diffusion of methane in silicalite, the trend of the mean square displacement of methane vs time, which is shown in Figure 1 (continuous line), has to be considered first. It appears that (Ir(t) - r(0)12)behaves approximately as a straight line, with slight oscillations. An estimate of the self-diffusion coefficient can be evaluated from eqs 3 and 5 , and from the oscillations of the values of the derivative present on the right-hand side of eq 4 the order of magnitude of the error can be derived. The value of D = (6.58 f 0.01) X m2 s-l results from eq 3, in good agreement with experiment [(6.5 f 1.0) X m2 s-]]. It should be stressed that the reported error of the computed D is the standard deviation from linearity of the mean square displacement vs time. The real accuracy, depending on the approximations of the adopted model, should be considerably less than that reported. Experimental data on the diffusion of CH4 in silicalite are reported in ref 9 for temperatures ranging from 160 to 300 K, so that the temperature dependence of the calculated diffusion coefficient could be checked. This test is important, because there are examples of systems where the trend of the calculated D vs the temperature was not as good as a single-temperature result could induce one to expect.25 MD simulations on this point are in progress. In Figure 1 the mean square displacements calculated by switching off the framework oscillations are also reported (dashed (23) Papp, H.; Hinsen, W.; Do, N. T.; Baerns, M. Thermochim. Acta 1984. 82, 137.
(24) VignB-Maeder, F.;Auroux, A. J . Phys. Chem. 1990, 94, 316. (25) Rahman, A . Phys. Rea. 1964, A136, 405.
50 75 100 125 time ( p s ) Figure 2. Mean square displacements of methane molecules along the x axis, y axis, and Z axis, corresponding to the a, b, and c cell parameters of the orthorhombic cell. The cell parameters of the monoclinic cell are chosen so that this correspondence is maintained ( w = x , JJ, 2). Continuous and dashed lines with the same meaning as in Figure 1. “0
25
line). This calculation was performed in order to evaluate the effect of the coupling between the framework and the sorbate molecules on the computed diffusion coefficients. The result is that, in the case of silicalite, mean square displacements are about 20% smaller for the fixed framework than for the moving one, at room temperature, resulting in a diffusion coefficient D = 5.41 X m2 s-’. This overall effect hides part of the larger discrepancies that appear when the diffusion is considered in more detail. This is evident when the motions along the three channel systems are observed. In Figure 2 the trends of mean square displacements of methane along the Cartesian coordinates vs time are reported both for mobile (continuous lines) and for fixed (dashed lines) framework simulations. It appears that the diffusion is larger along y , corresponding to the direction of the straight channels of silicalite, is smaller but relevant along x, the direction of the sinusoidal channels, and is almot negligible along z, for which no channel is present, but some diffusion is possible through the intersections of the channels along x and y . Indeed, direct inspection of the methane trajectories allowed to observe a few “jumps” along z during 200 ps. The differences among the diffusion along the Cartesian coordinates can be evidenced by computing the “monodimensional” diffusion coefficients from the slopes of the curves reported in Figure 2. For the oscillatine framework these coefficient are D- = 3.59 X 1 0 ~ 9 m 2 s - 1 , D , = ~ 4 . 190X- 9 m 2 s - 1 , a n d D Z = f.22X*10-9m2 s-1.
If this finding reflects the real values, it might be suggested that a preferential orientation of the crystallites along y or, in an case, normal to z, could enhance the efficiency of ZSM-5 as catalyst or molecular sieve. When the framework is kept fixed, the monodimensional diffusion coefficients are D, = 4.93 X 1O-’ m2 s-l; D,, = 9.66 X 10-9 m2 s-l; D, = 1.64 X m2 s-l. The differences with the values reported above deserve some comment. One of the reasons for this effect could be the negligible coupling between the diffusive motion and the vibrations of the framework. As a consequence, the slowly moving methane molecules “see” the framework as a set of thermal ellipsoids, so that their motion is not hindered along the straight channels parallel to b, where they “float” close to the channel axis (see below), while along a thermal vibrations tighten up the diameter of the sinusoidal channels, slowing down the guest molecules by collisions. Moreover, the enhanced diffusion along b can be explained by the energy exchange with the moving framework, resulting in a larger temperature distribution of the sorbate, as will shown below. From the methane-Of radial distribution function (Figure 3) it can be deduced that methane molecules move preferentially close to the axes of the channels. Indeed, the first maximum in the rdf occurs at about 4.1 A, a distance close to one-half the channel diameter (measured as distance between the wall oxygen nuclei),
The Journal of Physical Chemistry, Vol. 94, No. 10, 1990 4333
Diffusion of Methane in ZSM-5
-?
lt
J 2
: 6
4
8
10
r(A) Figure 3. Methane-oxygen atoms of the framework radial distribution function.
n’ of n-mers Figure 5. Frequency distributions of simultaneously present clusters made of n methane molecules: ( 0 )n = 1; (0)n = 2; (B) n = 3; (*) n = 4. Frequency distributions are normalized by dividing by the total number of steps.
1 5t
31
I
2
6
4
a
10
r(A) Figure 4. Methane-methane radial distribution function.
while the secohd maximum (at 5.8 A) corresponds to the distances with the other oxygen atoms of the eight-membered rings of the channels, and the third maximum, at 8.2 A, can be referred to a cylindrical shell containing, among the other ones, the surface oxygens of the adjacent channels. This situation could be related to the model by Derouane, AndrC, and L ~ c a s ~for ~ , van ~ ’ der Waals interactions between molecules and curved surfaces. Following this model, a molecule possessing a van der Waals radius close to the pore radius will appear asfloating, that is, will be subject to zero radial force, and thus it will achieve high mobility. This is approximately the case for methane in silicalite, as the minimum in the Ormethane potential function is about 3.9 A. As a consequence, the diffusion rate should be controlled by collision between methane molecules and should increase by lowering the loading. Both experimental and theoretical results are still insufficient, in our opinion, to draw final conclusions on this topics, and further investigations are going on in this laboratory. The methane-methane rdf (Figure 4) exhibits a unique maximum at r = 3.75 A, a distance slightly smaller than the minimum of the CH4-CH4 potential function. This maximum could reflect frequent collisions as well as permanent or transient methane dimers or clusters. In order to elucidate this point, a detailed analysis of the methane particles trajectories was performed, as explained in the second section of this paper. First of all, a clear-cut definition of “methane dimer” was required. From direct inspection of the trajectories, it is reasonable to state that a dimer existed when two methane molecules remained closer than 4.5 A. Many dimers oscillating without exceeding this distance were found. Moreover, clusters containing three, four, or more methane molecules, linear or branched (near (26) Derouane, E. G.; Andrt, J.-M.; Lucas, A. A. Chem. Phys. Lett. 1987, 137, 336. (27) Derouane, E. G.;Andr6, J.-M.; Lucas, A. A. J . Card. 1988, 110, 58. (28) Berne, B. J. In Physical Chemistry, an Aduanced Treatise; Eyring, H., Henderson, D., Jost, W., Eds.; Academic Press: New York, 1971; Vol. VIIIB, p 6 5 2 ff.
time ( p s ) Q.5
b
n-mers Figure 6. (a) Distributions of the lifetime for the clusters made of n methane molecules. Symbols and normalization as in Figure 5 . (b) Mean lifetimes of the clusters vs n.
channel intersections), were observed. For each cluster of n molecules ( n = 2 , ..., 24), it was assumed that all its subclusters of n‘ < n molecules did not enter in the number of the cluster of n’molecules (e.g., the two or three dimers discernible in a trimer were not enumerated as dimers). It should be reminded that in a channel cross section only one methane molecule can be accommodated, and, by setting the 24 disposable particles evenly spaced in the channels, their mean distance would be about twice a molecular diameter, so that moving molecules are forced to collide very frequently. Moreover, kinetic energy exchange with the channel wall atoms is sufficiently large and slow to give rise to more or less long-lived dimers or clusters. This picture of the diffusive process of methane in silicalite emerges from the following analysis of the methane trajectories. Henceforth, for simplicity a cluster containing n methane molecules will be called “n-mer”. The frequencies of the number of n-mers simultaneously present in the simulated system are reported
4334 The Journal of Physical Chemistry, Vol. 94, No. 10, 1990
-lo
1
2
3
4
5
6
.i
time ( p s ) Figure 7. Velocity autocorrelation function for methane molecules. .-ii
v)
e
2
4 Y
.-E frequency ( c m - 1) Figure 8. Power spectrum of methane molecules up to 200 cm-I. For higher frequency intensites are negligible, except for a weak line a t 350 cm-I. The intensity scale is arbitrary.
in Figure 5. It appears that, for instance, the most probable number of isolated molecules present in an instant is about 8, while for dimers the same number is 2, and only one for the n-mers with n > 2. For n > 4 the frequencies are smaller and they are not reported in the figure. One may conclude that about one-third of the diffusing molecules is, on the average, isolated, while the others are associated in clusters, two of which, at least, are dimers. The larger cluster observed (lasting only a few time steps) contained 15 molecules. In Figure 6, the distributions of the lifetime of the n-mers and their mean lifetimes are shown. Detectable fractions of monomers last up to about 1.5 ps before being captured to form a cluster (in one case, a life of 5.5 ps were observed), and the mean life is 0.43 ps. For dimers, the largest lifetime is about 1 ps, with a mean life of 0.2 ps. These values are sufficient to allow the dimers to oscillate a few times. The mean lives of trimers and tetramers are 0.15 and 0.12 ps, respectively, and most of them can be considered the result of multiple collisions, as they, on the average, hardly oscillate before they decompose. Further features of the dynamical behavior of the diffusing particles can be outlined by considering the velocity autocorrelation function, which is reported in Figure 7, and the power spectrum of the methane molecules (Figure 8). The power spectrum exhibits a relatively large band approximately in the range 20-1 10 cm-' with a peak of 5 3 cm-'. From direct inspection of the time dependence of the intermolecular distances, it was found that the frequencies above -30 cm-l correspond to the vibrations of the clusters (mainly dimers), while the frequencies below this value correspond more likely to the oscillations of the molecules due to the collisions with the channel walls, in particular in the sinusoidal channels.
Demontis et al. All these vibrations are reflected in the behavior of the velocity autocorrelation function (vacf), which shows a damped oscillatory trend. This vacf is more like a vacf for a one-dimensional system of Lennard-Jones particles than a vacf typical of Lennard-Jones liquid,28 as should be expected because the methane molecules are constrained to move in narrow channels. No evident resonance between methane and framework vibrations was found, except for a small peak at 350 cm-' in the methane spectrum (not shown in Figure 8), corresponding to the lowest frequency band of the framework. As was anticipated above, one can conclude that the coupling between the motion of the sorbed molecules and the framework vibrations is small. Another point of interest for the study of the dynamics of diffusing particles is the distribution of velocities and temperatures. Shin et aL5 found a large scatter among the temperatures of Naf ions diffusing in A-type zeolite, using a fixed framework model. In the case of this work, the large number of collisions allows a good equilibration of the simulated system, and the calculated distribution of the methane velocities resulted to be Maxwellian-like for the mobile framework simulation. For fixed framework the same trend is observed, but the velocity distribution is narrower. It seems that the vibrating framework acts by exchanging kinetic energy with the diffusing molecules and broadening the velocity distribution. More precisely, the methane molecules behave as a small canonical ensemble in a "thermal bath" aroused from the framework motion. Indeed, in spite of the small oscillations of the temperature of the whole system, the distribution of the methane molecules temperature was markedly large, with standard error u = 54 K, in line with the theoretical value for the temperature fluctuations in a canonical ensemble [u, = (2/3)1/2N-'/2T,for N particles at a temperature T; for the present calculation ul N 51 K].29 Concluding Remarks The results reported in this paper confirm that, among the various theoretical research method^,^^^^' molecular dynamics can be a useful and versatile technique for investigations in the field of the zeolites and, in particular, of the diffusion of sorbed molecules. The model used in this work, though very simple and approximate, is able to reproduce with good accuracy the measurable quantities related to the diffusion of methane in silicalite. Moreover, it allows one to attempt a detailed microscopic description of the diffusive process, from which some useful conclusions can be drawn. First of all, the different diffusivities along the crystallographic axes, if confirmed experimentally, could be exploited to enhance the effectiveness of ZSM-5 as catalyst or molecular sieve, once a technique able to give preferential orientation to the crystallites is devised. Second, the consideration of the details of the methane molecules motions outlined above could give suggestions about some aspects of the reactivity of the molecules sorbed in zeolites. Work is in progress along three lines: a check of the dependence of the diffusion coefficient on the temperature, an evaluation of the effects of different loadings on the diffusive process, and an attempt to generalize the results of these simulations to the study of the kinetics of diffusion and reactivity in molecular sieves. Acknowledgment. We thank Prof. A. Alberti who encouraged this work. This research is supported by Minister0 della Pubblica Istruzione (60% and 40%). Registry No. CH,, 74-82-8 (29) Lebowitz, J. L.; Perkus, J. K.; Verlet, L. Phys. Reo. 1967, 153, 250. (30) Sauer, J.; Zahradnik, R. Int. J . Quantum Chem. 1984, 26, 793. (31) Suffritti, G. B.; Gamba, A. I n f . Rea. Phys. Chem. 1987, 6, 299.
