# Molecular dynamics simulations in zeolites: from deterministic to random motion

S. El Amrani, M. Kolb

<figure><img src="image_1.png" /></figure>

## To cite this version:

S. El Amrani, M. Kolb. Molecular dynamics simulations in zeolites: from deterministic to random motion. The Journal of Chemical Physics, 1992, 96, pp.4857-4869. 10.1063/1.464315. hal-00006246

HAL Id: hal-00006246

https://hal.science/hal-00006246v1

Submitted on 25 May 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Molecular dynamics simulations in zeolites: From deterministic to random motion

S. El Amrani and M. Kolb

Institut de Recherches sur la Catalyse, 2 Avenue Albert Einstein, 69626 Villeurbanne Cedex, France and
Laboratoire de Chimie Théorique, Ecole Normale Supérieure, 46 allée d'Italie, 69364 Lyon Cedex
07, France

Molecular dynamics (MD) simulations of argon in silicalite were performed in order to study the transition from the deterministic motion at short times to the random motion at long times. A characteristic length $\lambda$ (time $\tau$) is associated with this change in dynamics. At a temperature $T=229(14)$ K, we find $\lambda=0.45(8)$ Å [$\tau=0.21(6)$ ps]. A detailed analysis of the trajectories shows that the optimal choice for the elementary time step $\Delta t$ in MD is slightly below $\tau$. Little is gained by decreasing $\Delta t$ as two trajectories with almost identical initial conditions diverge exponentially fast with time. A MD algorithm with variable $\Delta t$ is proposed and used to show that the size of the zeolite channels and cavities influences how the asymptotic Gaussian random process is reached.

## I. INTRODUCTION

Zeolites are crystalline microporous structures with important applications in the petrochemical industry (Ref. 1). Their usefulness is related to the large number of catalytically active sites and to their selectivity. These properties are a direct consequence of the internal structure of the crystal in the form of channels and cavities on a molecular scale.

During the past years, an increasing number of investigations, at first experimental and more recently also theoretical, have been performed to improve the understanding of the complex catalytic processes in zeolites. Two different aspects influence the properties of zeolites as catalysts: (i) the chemical reaction on a molecular scale; (ii) the transport properties on the scale of the channel width, and upward. Adsorption/desorption on the zeolite crystal, surface diffusion, penetration into the crystal, and the internal structure determines the global diffusion properties.

Different experimental methods have been used to determine diffusion coefficients for a variety of molecular diffusers in zeolites. Macroscopic sorption-uptake measurements have determined effective diffusion coefficients of hydrocarbons in A-type zeolites many years ago (Ref. 2). Gas chromatography (Ref. 3) and a permeation method using a zeolite membrane (Ref. 4) also yield global diffusion properties. Microscopic self-diffusion coefficients have been determined with the pulsed field-gradient technique, e.g., for methane, ethane, and propane in ZSM-5 type zeolites (Ref. 5). The diffusion coefficients from this method are confirmed by quasielastic neutron scattering experiments (they probe on a shorter time scale) (Ref. 6). There may be substantial differences between microscopic and macroscopic measurements, the latter being influenced by effects such as surface diffusion and by the difficulties to enter the crystal (Ref. 7).

Theoretical investigations have calculated both thermodynamic and transport properties of small molecules in zeolite frameworks. Early studies used atom-atom potentials

to determine heats of adsorption and adsorption isotherms (Refs. 8 and 9). Potential energy maps were obtained for xenon in zeolite rho (Ref. 10) and for benzene in silicalite (Ref. 11), both by Monte Carlo simulations (MC). MC simulations were also used to calculate the water distribution in a ferrierite-type zeolite (Ref. 12). Model potentials for methane, water, and methanol in silicalite were obtained using atom-atom potentials (Ref. 13).

Molecular dynamics (MD) simulations have been used successfully to calculate the diffusion coefficients for a variety of molecules in zeolites. Simulations of xenon in silicalite were performed in Refs. 14 and 15. Methane trajectories in a cavity of NaA zeolite were compared with experimental observations in Ref. 16. Methane MD trajectories in faujasite show that diffusion is confined to the surface of the cages (Ref. 17). Global diffusion properties based on lattice models for a two-dimensional network were obtained from Monte Carlo calculations (Ref. 18). The different theoretical approaches are reviewed in Refs. 19 and 20. The influence of the external surface barriers on the sorption has been studied in Ref. 21.

The purpose of the present work is to study the changeover from the necessarily deterministic trajectory at short distances to the apparently diffusive motion over long distances. Formally, the motion is determined for all times once the initial conditions are specified because the trajectory is obtained from solving Newton's equations. However, on a sufficiently large scale, the motion looks as if the particle jumps randomly in the different directions given by the zeolite structure. Qualitatively, one expects that the smallest (largest) crystalline parameter of the zeolite structure sets a crude upper (lower) bound for the range of the deterministic (random) motion. In MD, there are two conflicting constraints—one needs to use a small elementary time step $\Delta t$ to ensure accurate numerical integration of the equations of motion and one needs a large $\Delta t$ to obtain a trajectory that diffuses through many channels and cavities of the zeolite in order to obtain a reliable
estimate of the asymptotic diffusion coefficient. A careful
analysis of the trajectories leads us to propose an optimal
choice for $\Delta t$.

## II. THE MODEL

In order to have a well-characterized system, our calculation was performed for the diffusion of a single argon atom in silicalite at a temperature of ~230 K. There are no conceptual difficulties to study several interacting molecules, other zeolites, or molecules or to work at different temperatures (Ref. 22).

The MD algorithm used here is a standard finite time step solution of Newton's equations. The interaction between the silicalite and the argon is represented by a Lennard-Jones potential $V = -A/r^6 + B/r^{12}$ between the oxygen atoms of the silicalite and the argon atom in the atom-atom approximation. The parameters $A=3347.41$ kJ (Å)$^6$ mol$^{-1}$ and $B=2585547.9$ kJ (Å)$^{12}$ mol$^{-1}$ used are taken from Ref. 7. The silicalite framework is considered to be rigid with the position of the atoms taken from Ref. 23. The interaction range between argon and the zeolite was limited to 14 Å (corresponding to ~800 oxygen atoms); the truncated potential energy then deviates about 1% from the exact value. To obtain the position of the diffuser at time $t+\Delta t$ from the one at $t-\Delta t$ and at $t$, a constant NVE Verlet algorithm was employed (Ref. 24). Theoretically, the total energy of a single diffusing particle is rigorously constant. To assure that it does not deviate appreciably from its initial value, the trajectory was adjusted at each step (Ref. 25). The temperature was calculated from the measured average kinetic energy via the relation $\langle E_c \rangle = 3/2k_B T$. $\Delta t$ is kept fixed in the first part of our work; a variable time step version of the algorithm will be used in the remainder.

The initial position and velocity of the diffuser was chosen at random, the potential energy equal to -10.69 kJ/mol, and the kinetic energy equal to 1.83 kJ/mol or equivalently 146 K. The actual temperature measured in the simulations is $T=229(14)$ K. It is barely sensitive to the time step chosen. To test the stability of the algorithm, several elementary time steps $\Delta t$ were used in the range 0.000 05–0.05 ps; for the same computing effort, the total elapsed time of the simulation varies accordingly from 150 to 150 000 ps (the actual run for $\Delta t=0.05$ ps was only up to 40 000 ps).

## III. DIFFUSION COEFFICIENTS

The diffusion coefficient was calculated from the mean square displacement as a function of time, separately along each spatial direction. The average displacement then was obtained from $R^2(t)=[X^2(t)+Y^2(t)+Z^2(t)]/3$ and the diffusion coefficient from $D=R^2/2t$ (Ref. 22). In Fig. 1(a), the data for $R^2$ is plotted against $t$ for $\Delta t=0.0005$ ps [the data are averaged along the whole trajectory $R^2(t) = \langle R^2(t_0,t_0+t) \rangle_{t_0}$, $0<t_0<t_{max}-t$]. The total length of the trajectory $t_{max}$ is 1500 ps, which requires $3\times10^6$ iteration steps. The slope of the straight line segment determines $D=4.8\times10^{-9}$ m$^2$/s, using the Einstein relation; for com-

<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) The mean square displacement $R^2$ (in Ångstroms squared) for one atom of argon in silicalite at 229(14) K is plotted as a function of time $t$ (in picoseconds). The diffusion coefficient is calculated from the slope of the straight line. The elementary time step is $\Delta t=0.0005$ ps. (b) The same data are plotted on logarithmic ($\log_{10}$) scales. There is a clear crossover from the deterministic motion at small times to the random motion at large times.</figcaption></figure>

parison, the experimental value for methane at 200 K is $D=2.5 \times 10^{-9}$ m$^2$/s (Ref. 26). The figure shows clearly that the data gets unreliable beyond 300 ps, due to poor statistics. In order to extract more information from the data, it is plotted on logarithmic scales, as shown in Fig. 1(b). There are two major advantages in using logarithmic scales: (1) any straight line portion in such a plot indicates a power law behavior $y=Ax^\alpha$. The slope of the straight line fitting the data determines the power $\alpha$ and its height determines the prefactor $A$; (2) as a log-log plot expands the data range for small $x$ and compresses it for large $x$, the statistics in the less reliable large $x$ domain improves. For the MD trajectories, the fluctuations increase with increasing $R$ (a source of additional noise is the decreasing number of independent samples for a trajectory of fixed total length). Plotting the data logarithmically counteracts this trend.

The data of Fig. 1(b) shows very clearly two scaling regimes. Well below $t\sim1$ ps, the straight line has slope 2, corresponding to Newtonian dynamics $R=vt$. Well above $t\sim1$ ps, the slope is 1 corresponding to a diffusive motion $R^2=Dt$. The intersection of the two straight line fits defines the characteristic distance $\lambda$ (time $\tau$) that separates the two regimes. For argon in silicalite at 229(14) K, $\lambda=0.45(8)$ Å and $\tau=0.21(6)$ ps. The transition region is
<figure><img src="image_4.png" /></figure>

FIG. 2. Log₁₀-log₁₀ plot of the mean square displacement R² (in Ångströms squared) as a function of time *t* (in picoseconds) for different elementary time steps Δ*t*=0.000 05, 0.0005, 0.005, and 0.05 ps. The overlap of the data in the common time range indicates that the diffusion process is not affected by the variations of Δ*t*. The last series of points has been obtained by the variable Δ*t* algorithm.

given by the range over which the slope changes, in the
present case between 0.4 and 4 Å, respectively, between 0.2
and 30 ps. The straight lines determine *v*=3.2 m/s and the
same *D* as from Fig. 1(a).

IV. TRAJECTORIES

A frequently used criterion for determining the time
step Δt in MD is the fluctuation of the total energy around
its theoretical value—a variation below 1% is usually con-
sidered acceptable (Ref. 22). Here we will use as an alter-
native criterion the stability of the velocity at short times
and the diffusion coefficient at longer times. The observa-
tion in Fig. 1(b), of a clear separation between the New-
tonian and the diffusive regimes suggests that no substan-
tial error is introduced as long as Δt lies in the Newtonian
range. In order to test this idea, we performed a series of
simulations with identical parameters and initial condi-
tions, but with Δt varying from 0.000 05 to 0.05 ps. Figure
2 shows the curves of R² vs t. The data for each Δt covers
a different time range, but in the common time regions,
they coincide to within run to run fluctuations. In partic-
ular, the crossover from Newtonian to random motion oc-
curs for the same values of R² with respect to Fig. 2 also
shows clearly that the 0.000 05 ps data yield much too
short a trajectory to reliably determine D.

In order to study the trajectories in more detail, a com-
parison between trajectories with identical parameters, but
with slightly different initial conditions has been made. To
measure the stability of a trajectory, we determine the
mean square distance ΔR²=(R₁−R₂)² as a function of
time between two trajectories with position vectors R₁ and
R₂, respectively, which initially differ by a very small
amount ΔR(t=0)=ΔR₀. Note that ΔR² compares two
trajectories at the same time t, but with different initial
conditions. Figure 3(a) shows the logarithm (log₁₀) of
ΔR(t) the distance between the trajectories as a function of
time t. The straight line behavior indicates an exponential
divergence with a characteristic time t₀=1.1 ps. A change

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

FIG. 3. (a) Logarithm (log₁₀) of squared distance ΔR² (in Ångströms squared) as a function of time *t* (in picoseconds) between two trajectories which are initially separated by a very small amount ΔR₀<1 Å. The straight line behavior indicates an exponential divergence. The different curves (from left to right) correspond to different initial separations ΔR₀ =ΔR(t=0)=10⁻ⁿ Å, *n*=3, 6, 9, 12. (b) Evolution of the logarithm (log₁₀) of the squared distance ΔR² (in Ångströms squared) as a function of time *t* (in picoseconds) between trajectories with identical initial conditions, but using different time steps Δ*t*. The curves (from left to right) correspond to the pairs (Δ*t*₁, Δ*t*₂)=(5×10⁻³, 5×10⁻⁴), (5×10⁻⁴, 5×10⁻⁵), and (5×10⁻⁵, 5×10⁻⁶) ps.

in the initial distance ΔR₀ simply shifts the curves in the
figure, suggesting a form ΔR(t) = ΔR₀e^t/t₀. The occa-
sional drop in the distance occurs when—accidentally—
the trajectories almost cross; presumably this could be
avoided if the velocity difference were included in measur-
ing the “distance” between the trajectories. A wide range
of ΔR₀ was covered in the simulations by setting ΔR₀
= 10⁻ⁿ Å with n=3, 6, 9, 12, and 15. In all cases, the
exponential behavior continues up to a separation of
ΔR₀ ~ 1 Å. This value corresponds to the onset of the dif-
fusive regime. This correspondence is also supported by the
fact that t₀ and τ have comparable values. The observed
divergence of the trajectories not only correlates with the
onset of the diffusive regime, but may be considered to be
at its origin. An analogous behavior of the trajectories has
been observed for MD of liquids (Ref. 27).

The integration time step Δt in a MD simulation in-
troduces an error of O(Δt⁴). Therefore, we expect to find
similar behavior as in Fig. 3(a) when comparing trajecto-
ries from identical initial conditions, but with a different
integration step Δt. Figure 3(b) shows the logarithm
(log₁₀) of the squared distance ΔR² between pairs of tra-
trajectories using different time steps $\Delta t_1$ and $\Delta t_2$ as a function of time $t$. The figure shows the same exponential divergence as Fig. 3(a), even though the data are much noisier. The trajectories also become randomized on a scale of the order of $\Delta R \sim 1$ Å. The observed exponential behavior implies that for an equally accurate trajectory of twice the original length, one needs twice as many digits in the calculation. In practice, this means that it is impossible to determine a MD trajectory accurately beyond a few picoseconds.

## V. VARIABLE TIME STEP MD

From the above calculations, one concludes that (i) it is not possible to calculate long trajectories accurately, but (ii) that it is not necessary either to do so to obtain a good estimate of the diffusion coefficient. How far can one increase $\Delta t$ before an appreciable change of the diffusion properties is observed? An absolute limitation for the modified Verlet algorithm is the criterion that the kinetic energy (before the correction) be positive (Ref. 25). For the set of parameters used here, $\Delta t$ has to stay below ~0.005 ps to avoid negative kinetic energies. Not surprisingly, the method fails as one approaches $\tau$, the upper limit of the Newtonian regime.

If some exceptions were acceptable, the MD could be accelerated considerably. Therefore, we implemented an algorithm that determines $\Delta t$ dynamically. One starts with a large basic time step $\Delta t_0$ (~0.05 ps), which would sometimes yield negative kinetic energies for a small fraction of the iterations (~1%). Whenever for a step from $t$ to $t+\Delta t_0$ the computed kinetic energy is negative, the diffuser returns to its position at time $t-2\Delta t_0$ and three steps of length $\Delta t_0$ are replaced by 3$n$ steps of length $\Delta t_1=\Delta t_0/n$ ($n$ is kept fixed, typically $n=10$). The following step uses the original $\Delta t_0$ again, if possible: if the refined trajectory still leads to negative energies, the procedure is refined further—the original $\Delta t_0$ is then replaced by 2$n$ steps of length $\Delta t_2=\Delta t_0/2n$. Continuing this process allows for arbitrary refinements ($\Delta t_k=\Delta t_0/kn$, $k=3,4,5,...$). The refined positions at $t-\Delta t_k$ needed to calculate the positions at $t+\Delta t_k$ are obtained by linear interpolation.

In the actual calculations, a more restrictive criterion than the positiveness of kinetic energy is used to initiate the refinement of the step size—the total energy must not deviate more than 10% from its theoretical value. This avoids the kinetic energy getting too small. With this method, the basic time step could be raised by more than one order of magnitude from 0.005 (the upper bound for the fixed $\Delta t$ procedure) to 0.05 ps, effectively gaining a factor of 5 in speed. The simulations show that when the energy is outside the tolerated range, there are often many (up to 30 for $n=10$) subdivisions necessary. The mean square displacement does not change measurably using the variable $\Delta t$ method. In Fig. 2, the mean square displacement for the modified algorithm is compared with the standard algorithm. There is clear agreement for the Newtonian as well as the diffusive range. Beyond ~100 ps, the statistics of the variable $\Delta t$ for the standard simulations with fixed small time steps is so poor that no diffusion

<figure><img src="image_7.png" /><figcaption>FIG. 4. Rescaled (dimensionless) distribution function $\Pi=\mathbf{R}(t) \cdot \mathcal{P}(\mathbf{R},t)$ plotted against the reduced mean square displacement $\rho=\mathbf{R}/\mathbf{R}(t)$. $\mathbf{R}(t)$ is the average displacement $\mathbf{R}(t) = \sqrt{\langle [\mathbf{\Sigma}, \mathbf{R}_l(t)]^2 \rangle / \langle \mathbf{\Sigma} \rangle}$ [1]. For comparison, the corresponding Gaussian distribution for a random process is superposed. The curves correspond to 1, 2, 10, 100, and 300 ps.</figcaption></figure>

coefficient can be extracted. On the other hand, the variable time step simulations yield reliable data well beyond the lattice periodicity. The curves obtained in this way show that the asymptotic diffusion regime is attained only well beyond the periodicity of the zeolite. The anisotropy explains the long transitory regime; the separate diffusion coefficients for the different lattice directions vary considerably (Ref. 22). In an intermediate stage, the average root mean square displacement is then a mixture of Newtonian, transitory, and diffusive motions in the different directions.

It is not surprising that the optimal value of $\Delta t_0$ is below, but not far below the channel width as estimated from purely energetic considerations (2–3 Å). Increasing $\Delta t_0$ further does not speed up the calculation as the number of exceptions increases very fast.

## VI. DISTRIBUTION FUNCTIONS

Calculating diffusion coefficients for zeolites requires long simulation times to assure that many cavities and channels are visited by the diffuser. The “true” diffusion coefficient is only attained in this limit. In the same limit, the distribution functions of the mean square displacement $\mathbf{R}^2(t)$ and of the length of the trajectories $s(t)$ do approach the expected Gaussian form. Using the large time step simulations, it is possible to probe a large number of different channels and to determine how the asymptotic limit is approached. Histograms for $\mathbf{R}^2(t)$ and for $s(t)$ were generated for different times $t$ varying from 0.1 to 300 ps. The results are shown in Figs. 4 and 5. Plotting the data in normalized form shows clearly that for the largest times, the Gaussian form is approached. As an interesting detail, one notices that the profile becomes almost Gaussian for $t \approx 1$ ps then again deviates from it for intermediate times before returning to the asymptotic limit for the largest times used ~300 ps. This may be interpreted as follows: the trajectory gets randomized on the scale of the effective channel width (~2 Å) of the zeolite structure, but the
<figure><img src="image_8.png" /></figure>

FIG. 5. Rescaled (dimensionless) distribution function $\phi=\delta s(t) \cdot P[s-s(t),t]$ plotted against the reduced displacement $\xi=[s-s(t)]/\delta s(t) \cdot s(t)$ is the average length of the trajectory $s(t)=\Sigma_{i}(t)/\Sigma_{i}1$, and $\delta s(t)$ is the average length of its variance $\delta s(t) = \sqrt{\{\Sigma_{i}[s_{i}(t)-\bar{s}(t)]^{2}/\Sigma_{i}1\}}$. For comparison, the Gaussian distribution for a random process is superposed. The curves correspond to 1, 2, 10, 100, and 300 ps.

asymptotic randomization only sets in on scales larger than the periodicity of the crystal (~10 Å).

## VII. DISCUSSION

In order to optimize MD calculations for diffusion in zeolites, an argon particle in a silicalite structure has been simulated using different methods. The calculation shows that reducing the elementary time step much below the effective channel width cannot improve the precision of the diffusion coefficient. In a rather broad intermediate range (roughly between the channel width and its length), the effectively measured diffusion coefficient varies continuously. The truly diffusive regime can only be reached by very long trajectories. This is related to the channel/cavity structure and the anisotropy of the zeolite (Ref. 22). For the specific example studied, the variations of the diffusion coefficients are quite significant—they depend on the total length of the trajectories used and vary by a factor of 3 in the range between 1 and 10 Å [Fig. 1(b)].

Nothing in the proposed method is specific to the example studied. It can therefore be applied to larger molecules, including rotational degrees of freedoms, without any modifications.

## ACKNOWLEDGMENT

We acknowledge many stimulating discussions with and a critical reading of the manuscript by F. Vigné-Maeder.

<footnote n="1">¹Studies in Surface Science and Catalysis: Catalysis by Zeolites, edited by B. Imelik, C. Naccache, Y. Ben Taarit, J. C. Vedrine, G. Coudurier, and H. Praliaud (Elsevier, Amsterdam, 1980), Vol. 5.</footnote>

<footnote n="2">²J. Karger and J. Caro, J. Chem. Soc., Faraday Trans. 1 73, 1363 (1977).</footnote>

<footnote n="3">³A. S. Chiang, A. G. Dixon, and H. Y. Ma, Chem. Eng. Sci. 39, 1461 (1984).</footnote>

<footnote n="4">⁴D. T. Hayhurst and A. R. Paravar, Zeolites 8, 27 (1988).</footnote>

<footnote n="5">⁵J. Caro, M. Bülow, W. Schirmer, J. Karger, W. Heink, and H. Pfeifer, J. Chem. Soc., Faraday Trans. 1 81, 2541 (1985).</footnote>

<footnote n="6">⁶H. Jobic, M. Bée, J. Caro, M. Bülow, and J. Kärger, J. Chem. Soc., Faraday Trans. 1 85, 4201 (1989).</footnote>

<footnote n="7">⁷J. Kärger and D. Ruthven, Zeolites 9, 267 (1989).</footnote>

<footnote n="8">⁸A. V. Kiselev, Adv. Chem. Ser. 102, 37 (1971).</footnote>

<footnote n="9">⁹A. V. Kiselev, A. A. Lopatkin, and A. A. Shulga, Zeolites 5, 261 (1985).</footnote>

<footnote n="10">¹⁰P. A. Wright, J. M. Thomas, A. K. Cheetham, and A. K. Nowak, Nature 318, 611 (1985).</footnote>

<footnote n="11">¹¹A. K. Nowak, A. K. Cheetham, S. D. Pickett, and S. Ramdas, Mol. Simul. 1, 67 (1987).</footnote>

<footnote n="12">¹²L. Leherte, D. P. Vercauteren, E. G. Derouane, and J. M. André, in *Innovation in Zeolite Materials Science*, edited by P. J. Grobet, W. J. Mortier, E. F. Vansant, and G. Schulz-Ekloff (Elsevier, Amsterdam, 1987), p. 293.</footnote>

<footnote n="13">¹³F. Vigné-Maeder and A. Auroux, J. Chem. Phys. 94, 316 (1990).</footnote>

<footnote n="14">¹⁴S. D. Pickett, A. K. Nowak, J. M. Thomas, B. K. Peterson, J. F. P. Swift, A. K. Cheetham, C. J. J. den Ouden, B. Smit, and M. F. M. Post, J. Phys. Chem. 94, 1233 (1990).</footnote>

<footnote n="15">¹⁵R. L. June, A. T. Bell, and D. N. Theodorou, J. Phys. Chem. 94, 8232 (1990).</footnote>

<footnote n="16">¹⁶E. Cohen De Lara, R. Kahn, and A. M. Goulay, J. Chem. Phys. 90, 7482 (1989).</footnote>

<footnote n="17">¹⁷S. Yashonath, P. Demontis, and M. L. Klein, Chem. Phys. Lett. 153, 551 (1988).</footnote>

<footnote n="18">¹⁸E. Aust, K. Dahlke, and G. Emig, J. Catal. 115, 86 (1989).</footnote>

<footnote n="19">¹⁹C. J. J. den Ouden, B. Smit, A. F. H. Wielers, R. A. Jackson, and A. K. Nowak, Mol. Simul. 4, 121 (1989).</footnote>

<footnote n="20">²⁰G. B. Suffriti and A. Gamba, Int. Rev. Phys. Chem. 6, 299 (1987).</footnote>

<footnote n="21">²¹F. Vigné-Maeder, J. Catal. 117, 566 (1989).</footnote>

<footnote n="22">²²S. El Amrani, thesis, Villeurbanne, 1990.</footnote>

<footnote n="23">²³D. H. Olson, G. T. Kokotailo, and S. L. Lawton, J. Phys. Chem. 85, 2238 (1981).</footnote>

<footnote n="24">²⁴L. Verlet, Phys. Rev. 159, 98 (1967).</footnote>

<footnote n="25">²⁵D. M. Heyses, Chem. Phys. 82, 285 (1983).</footnote>

<footnote n="26">²⁶P. Demontis, E. S. Fois, G. Suffriti, and S. Quartieri, J. Chem. Inf. Cap. Sci. 4, 30 (1990).</footnote>

<footnote n="27">²⁷P. Allen and D. Tildesley, *Computer Simulation of Liquids* (Oxford Science, Oxford, 1987).</footnote>