316
J. Phys. Chem. 1996, 100, 316-326
Molecular Dynamics Simulation of Single-File Systems K. Hahn and J. Ka1 rger* UniVersita¨ t Leipzig, Fakulta¨ t fu¨ r Physik und Geowissenschaften, Linne´ strasse 5, D-04103 Leipzig, Germany ReceiVed: June 28, 1995; In Final Form: September 18, 1995X
The diffusion of particles in single-file systems is investigated by means of molecular dynamics simulation. For three different system configurations, the mean square displacements are calculated and compared with the result of analytical considerations. For a straight tube the mean square displacement is proportional to the time t. For a straight tube, where random forces act on the particles, and for a periodic tube potential it is proportional to the square root of time.
I. Introduction Molecular diffusion in microporous solids has become an attractive subject of both fundamental investigations and applied research.1-17 Their regular structure and numerous technical applications as catalysts, cation exchangers, adsorbents, and host materials for advanced technologies have made zeolites a particularly important candidate for such studies. In addition to extensive experimental studies,1-3 in the last few years a number of theoretical papers have been devoted to diffusion studies in zeolites including Monte-Carlo simulations,4-7 applications of the theory of absolute reaction rates,8,9 and molecular dynamics simulations.10-17 Molecular simulations proved to be a versatile tool for correlating the microscopic structural and energetic conditions of the adsorbent-adsorbate systems to their diffusion properties and were in particular found to confirm the order of magnitude of the results of PFG NMR self-diffusion studies.13-15,17,18 So far, molecular dynamics simulations have mainly been carried out in zeolites of type A, X, or silicalite/ZSM-5. These zeolites contain a three-dimensional pore network, and molecular transport in such systems is subject to the laws of normal diffusion, in particular to the Einstein relation
〈z2(t)〉 ) 2Dt
(1)
indicating that for sufficiently large observation times the molecular mean square displacement increases in proportion with the observation time t. Such a relation must hold for any arbitrarily chosen direction of the z-coordinate. In isotropic systems (zeolites X and A) the diffusivity D defined by this relation does not vary with the direction of the z-coordinate. In anisotropic media such as zeolite ZSM-5, the diffusivity generally changes with changing orientation. A completely different situation may be encountered with zeolites containing unconnected straight channels. Owing to recent progress in zeolite synthesis, there are already a substantial number of representatives of this special type of zeolite, including ZSM-12, -22, -23, -48; AlPO4-5, -8, -11; L; Omega; EU-1; and VPI-5.20 Practical problems associated with the transport properties of such systems have become a topic of current research interest.19,21-24 In the case where the channel diameter and the diameter of the diffusing particles are of the same order, the particles have no possibility to change their relative order. This special case of diffusion is called single-file diffusion. Due to the fact that mutual passages of the particles are excluded, a displacement X
Abstract published in AdVance ACS Abstracts, December 1, 1995.
0022-3654/96/20100-0316$12.00/0
of a given particle will necessitate the shift of many other particles in the same direction. Therefore, there is a high degree of mutual correlation between the movement of different particles. This correlation will lead to an essential change in the behavior of single-file systems compared with that of normal diffusion. A first comprehensive treatment of the influence of singlefile behavior on zeolite diffusion and catalysis is presented in ref 21. It is based on Monte-Carlo simulations where the diffusing particles are assumed to be random walkers with a well defined step length, and it includes results of previous theoretical studies and model considerations.25-27 As a main result of these considerations, in the long-time limit the mean square displacement is found to be proportional to the square root of the observation time rather than to the observation time itself, so that in the case of single-file diffusion, Einstein’s relation (eq 1) must be replaced by a relation of the type
〈z2(t)〉 ) 2Fxt
(2)
where the factor of proportionality F has been termed the mobility factor of single-file diffusion. As to our knowledge, corresponding molecular dynamics simulations have not yet been carried out. It is true that the zeolites ferrierite and mordenite, which have been considered in molecular dynamics simulations in refs 16 and 28, are good candidates for single-file studies. However, the molecular displacements considered in these studies are too small to reveal single-file behavior. There are a number of crucial questions which can only be answered by the help of molecular dynamics simulations. For example, without the rigorous treatment by molecular dynamics simulations, it must in particular remain unknown whether eq 2 holds quite generally for single-file systems or whether it is a particular consequence of the random walker model. Moreover, only in this way may the influence of the microscopic features of the system under study on the diffusion behavior be studied rigorously. For elucidating these features there is no necessity to apply particular zeolite potentials. In fact, MD simulations with single-file systems turn out to lead to a number of computational problems that have no equivalence in the case of normal diffusion. It is therefore essential to choose interaction potentials as simple as possible. It turns out that with the simple potentials used in this study most aspects of the global behavior of singlefile systems may be successfully explained. This paper is organized as follows. The next section describes the methods used for the simulation and the potentials of the © 1996 American Chemical Society
Molecular Dynamics Simulation of Single-File Systems
J. Phys. Chem., Vol. 100, No. 1, 1996 317
tube and of the diffusing particles. Also some problems related to the simulation are discussed. Section III gives a short overview of relations which may be derived analytically for single-file systems by statistical arguments. The results of the molecular dynamics simulations and their interpretation are presented in section IV. Some analytical calculations are given in the appendices. II. Model, Potentials, and Simulation Parameters The simulations are carried out using the commonly used velocity-Verlet algorithm.29-31 In most cases a time step of h ) 1.5 × 10-14 s is used. In some simulations, especially with strongly varying pore diameters, numerical instabilities occur with this time step; in such cases a smaller time step is used. For the interaction between particles a shifted force LennardJones potential
{
Vs(r) )
(( ) ( ) )
4s 0
σs r
12
-
σs r
6
- c1(r2 - r2c ) - c2
for r < rc
(3)
for r > rc
is used, where the cutoff is at r ) rc ≡ 2σs. r denotes the separation between the centers of the particles. The constants c1 and c2 are chosen in such a way that both potential and force vanish at r ) rc. The quadratic form of the additional term is used to avoid the calculation of the square root in the force loop. For the parameters σs and s we have used the values s ) 164.0 kB K and σs ) 0.383 nm (with kB denoting the Boltzmann constant), which happen to coincide with the values of krypton as given in the literature.30,32 Correspondingly, for the particle mass the value of krypton (m ) 83.8 u, with u denoting the unified atomic mass unit) has also been used. The infinite tube is simulated by a cylindrically symmetric potential V(F,z) with length L and a steep wall at the tube radius. The variable F denotes the distance of the center of the particle from the tube axis. Periodic boundary conditions are used in the z-direction. The simplest potential which could hold the particles in the tube is a hard core potential, but this type of potential leads to problems with energy and momentum conservation in molecular dynamics simulations. Therefore, a slightly more complicated potential has been chosen. An inverse Lennard-Jones potential
Vl(ρ,z) )
{
(( ) ( ) )
4l 0
ρ σl
12
-
ρ σl
6
+ l
for ρ > ρc
(4)
for ρ < ρc
was found to be appropriate for our purpose. The z-dependence of the potential, if any, is contained in σl ) σl(z); i.e., the tube radius may be considered to be a function of z. The cutoff is at F ) Fc ≡ σl/(21/6), where the potential has its minimum. The parameter l ) 150 kB K is chosen to be on the order of s. Figure 1 represents schematically the energy profiles used in this study and the resulting steric conditions within the tubes. With eq 3 and the profile represented in Figure 1a, the particle diameter is seen to be on the order of σs. Likewise, with eq 4 and Figure 1b, the part of the tube that may be occupied by the particle centers may be seen to have a diameter on the order of 2σl. Hence, the effective tube diameter σeff is on the order of 2σl + σs. Single-file behavior may be expected for tube diameters smaller than twice the particle diameters, which may be transferred to the condition
Figure 1. Schematic representation of the energy profiles used in this study (a, particle-particle interaction; b, particle-tube interaction) and steric conditions within the tube resulting from them (c). The shaded area in Figure 1c represents the range which may be occupied by the particle centers. The two indicated particles assume outmost positions within the tube. For simplicity, the tube is assumed to be unstructured; i.e., σl is independent of z.
1 σl < σs 2
(5)
The calculations are carried out in the microcanonical ensemble; i.e., the energy should be conserved. To fulfill this requirement, the average of the energy of 20 measurements during a time interval of 10-11 s is determined. If the deviation from the initial value is greater than 0.5%, the energy is reset to the initial value by rescaling the velocities. The average time between such corrections of the energy turned out to be on the order of (0.5-1.0) × 10-9 s. Due to the fact that the particles are ordered in a onedimensional tube, the use of neighborhood tables is not necessary. The positions of the particles are stored in a onedimensional array, and for a given particle only the next two or three particles in the array may be inside the cutoff radius rc. For this reason the calculation time is only proportional to the particle number N. Therefore, much larger particle numbers could be considered than in simulations in three dimensions. These large particle numbers are in fact indispensable due to the special structure of the problem. In a tube with a finite length the shift of a particle has to remain limited, no matter whether periodic boundary conditions are used or not. This is due to the strong correlations between the shifts of the particles. For example, if the leftmost particle would have a shift over say the total tube length L to the right, all other particles must also move to the right. However, this would be in contrast to the requirement that, for an unstructured tube and without external fields acting on the particles, the total momentum and hence the center of mass must be preserved. One may conclude, therefore, that, for sufficiently long times, the mean square displacement 〈z2〉 in a given system will approach a finite value 〈z2〉∞, which must be markedly smaller than the file length L. This situation is illustrated by Figure 2. It shows the mean square displacements 〈z2〉 for four systems with different small
318 J. Phys. Chem., Vol. 100, No. 1, 1996
Hahn and Ka¨rger All numerical results are obtained with a computer code written by the authors in the C-programming language. The average computation time per particle and time step on a INDIGO2 system, which we have used in our simulations, is about 5 × 10-6 s. Correspondingly, a typical run with 4000 particles over 5 × 10-8 s then requires about 20 h. Section III. Results of Analytical Treatment Single-file systems may be described theoretically in different ways. Often a random walker model is used,25-27 either directly or by considering the diffusion of the vacancies between the particles. The mean square displacement obtained from these models is
〈z2〉 ) Figure 2. Time dependence of the mean square displacement 〈z2〉 for simulations with different small numbers of particles. 〈z2〉 tends to a limiting value 〈z2〉∞. For comparison, the result for 10 000 particles is also given.
TABLE 1: Limiting Mean Square Displacements 〈z2〉∞ and the Instants of Time t∞ When the Limit Is Reached for Different Numbers of Particles N
t∞/s
50 100 200 400
2.0 × 10 3.5 × 10-10 8.0 × 10-10 1.5 × 10-9 -10
〈z2〉∞/m2 (simulation)
〈z2〉∞/m2 (eq 6)
2.1 × 4.0 × 10-17 9.0 × 10-17 1.6 × 10-16
1.96 × 10-17 3.91 × 10-17 7.82 × 10-17 1.56 × 10-16
10-17
numbers of particles (N ) 50, 100, 200, and 400) in comparison with those determined for a system with a very large number (N ) 10 000) of particles, which is used as a reference system. The simulation is done for a straight tube as in section IV.C. The relative occupancy is set to θ ) 0.2, and the temperature is set to T ) 300 K in all cases. After an initial period of time where all systems behave the same, the systems leave the reference curve and approach their limiting values 〈z2〉∞. In appendix A we calculate an analytical estimate for this quantity. As one can see from Table 1, where both the values from the simulation and from the calculation are given, the analytical result
1 Nd2 〈z2〉∞ ) (1 - θ)2 2 6 θ
(6)
is in good agreement with the MD simulation. One should note here that eq 6 applies to all cases simulated later in this paper, provided that the total momentum vanishes or is small. Moreover, eq 6 does not depend on the temperature, on the specifics of the tube potential, or on the type of particle motion (e.g. random walker or deterministic motion). In Table 1 the instants of time t∞ when the systems reach the limit are also given. The quantity t∞ is defined as the time when the first maximum in the curve is reached. For the case of a straight tube, t∞ is proportional to the particle number N. Of course, for systems with the same particle numbers but with other types of particle motion, the dependency on the time may differ. As one can see from Figure 2, the instants of time when the mean square displacements leave the reference curve are up to one order of magnitude smaller than t∞. Therefore, to avoid the influence of the finite longitudinal extension of the tube, a large number of particles has to be chosen. All calculations presented in the following sections are carried out with at least 4000 particles.
1 - θ 2 2 1/2 1/2 l t θ τπ
( )
(7)
where τ is the time step between successive jump attempts and l is the jump length. θ denotes the relative site occupancy of the single file. The essential result of these model considerations is that the mean square displacement in single-file diffusion is proportional to the square root of the observation time rather than to the time itself as for normal diffusion. On the basis of eq 7, in analogy to Einstein’s relation (eq 1), a single-file mobility factor
F)
〈z2〉 2xt
)
1 - θ l2 θ x2πτ
(8)
has been introduced. The random walker models can be simulated numerically by a Monte-Carlo simulation, where many random walkers jump in a one-dimensional array of possible positions and where a jump is forbidden if the destination is already occupied by another random walker. The results of such simulations [21] agree with the theoretical result (eq 7). In ref 33 molecular transport in a single-file system is derived from the motion of an isolated particle, i.e. of a particle without interaction with other particles, in the same environment. These considerations clearly include the case of a random walker. They may be applied, however, also to any other mode of molecular propagation including continuous movement with respect to both the space and time scale. The only condition is that the particle movement proceeds ideally in one dimension. As a result, the mean square displacement of a particle in the single-file system may be shown to be related to the absolute value 〈|s|〉 for the shift of an isolated (“noninteracting”) particle in the same tube by the relation
〈z2〉 )
(1 - θ) 〈|s|〉 c
(9)
where c ) N/L is the particle concentration. The quantity (1 - θ)/c may be visualized as the mean clearance between adjacent particles. IV. Simulation Results A. Unstructured Tube. The simplest single-file system accessible for molecular dynamics simulations is a long unstructured tube, i.e. a tube with a radius that does not depend on z. We have chosen a value of σl ) 0.12 nm, which, with the given value of σs ) 0.383 nm, clearly fulfills condition 5, which is required for ensuring single-file behavior. It could in fact be verified during the calculations that there were no mutual passages of the particles with this parameter setting. MD
Molecular Dynamics Simulation of Single-File Systems
J. Phys. Chem., Vol. 100, No. 1, 1996 319 mean absolute velocity is 〈|Vz|〉 ) (2kBT/πm)1/2 and thus the mean absolute shift is 〈|s|〉 ) (2kBT/πm)1/2t. Then, from eq 9, the mean square displacement in the single-file system follows to be
〈z2〉 )
Figure 3. Time dependence of the mean square displacement 〈z2〉 for different temperatures in a straight tube (θ ) 0.2).
simulations in unstructured tubes yield particularly large displacements. They have been, therefore, carried out with 10 000 particles to minimize the influence of the finite tube length. In one case, even a system with 40 000 particles is simulated. The consideration of such a large number of particles is probably exceptional for MD simulations in adsorbateadsorbent systems. Figure 3 shows the mean square displacement for three different temperatures (100, 300, and 500 K), as a function of the observation time. In all cases, a sorbate concentration of one particle per 1.915 nm tube length was considered. If the particles are assumed only to attain positions on the channel axis, this would correspond to a relative occupancy of θ ) 0.2. After an initial ballistic phase where 〈z2〉 ∼ t2, the mean square displacement is found to increase in proportion to the time t. The transition between the two regimes occurs at a time which is on the order of the time of free flight of the particles; i.e., the “long-time” behavior starts after the first interaction between adjacent particles. It is remarkable that there is no indication that the observed time dependence 〈z2〉 ∼ t tends to approach the proportionality 〈z2〉 ∼ xtas predicted by eqs 7 and 8. The inflections which in fact became apparent from Figure 3 are in reality caused by finite size effects. For 10 000 particles, according to eq 6, the limiting value for the mean square displacement becomes 〈z2〉∞ ) 3.91 × 10-15 m2, and extrapolating the results given in Table 1, this displacement is reached after time t ) 3.5 × 10-8 s and the influence of the finite size effects begins approximately one order of magnitude earlier, i.e. at t = 3 × 10-9 s. That the observed deviation from the proportionality 〈z2〉 ∼ t is really due to finite size effects is convincingly demonstrated by the simulation at 500 K, which has been carried out with 40 000 particles. In this case disturbing effects due to the finite size of the single-file system are only to be expected for t g1.2 × 10-8 s, which is still beyond the range of our representation. Though being in striking contrast to the random walker model and eq 7, the observed proportionality 〈z2〉 ∼ t, however, may be easily explained by using eq 9. In fact, for an unstructured tube the random walker model does not apply, since the movement of the individual particles is deterministic rather than chaotic, so that the z-component of the momentum is clearly conserved. Therefore, for a single particle in the tube the velocity Vz is constant and the shift of the particle is s ) Vzt. Having an ensemble of systems with single particles, where the velocity is Gaussian distributed as in the single-file case, the
( )
1 - θ 2kBT σ θ s πm
1/2
t
(10)
where the concentration was substituted by c ) θ/σs. Equation 10 is of the type of Einstein’s relation (eq 1) and allows the determination of a diffusivity D. Inserting the relevant data into eq 10, at T ) 500 K one obtains a diffusion coefficient D ) 1.36 × 10-7 m2 /s. On the other hand, linear regression 〈z2〉 ) 2Dt + b of the molecular dynamics simulation results in the time interval 10-10 to 10-9 s, yielding D =1.3 × 10-7 m2/s. Hence, the simulation results and the theoretical prediction on the basis of the analytical relations of section III are found to be in excellent agreement. This is particularly noteworthy, since in the simulations clearly also velocity components perpendicular to the tube axis have been considered. The simulations presented in this section have demonstrated that molecular transport in single-file systems need not necessarily lead to a proportionality between the mean square displacement and the square root of the observation times as required by eq 7. In fact, the simulations confirmed the more general eq 9, which predicts proportionality between the mean square displacement and the observation time if an isolated particle would be able to move freely within the single-file system. In real systems, in particular in zeolites with onedimensional channels, such a situation does not occur. Neither can the zeolite lattice be assumed to be ideally rigid, nor can the individual channels be considered as unstructured tubes. The influence of a deviation from these idealized conditions will be considered in the following sections. B. Unstructured Tube with Random Force. The simplest way to include the influence of lattice vibrations is to consider them as a random force acting on the diffusing particles. This random force is simulated by a random change of the particle velocities Vz after every 10-12 s. The change of the velocities is achieved in the following way: The velocity Vi of particle i is transformed to
Vˆ i ) Vi + u
(11)
where u is Gaussian distributed. The square width 〈u〉 of this distribution is set equal to a2〈V2〉, where 〈V2〉 is the width of the velocity distribution of all particles, i.e. 〈V2〉 ) kBT/m. The parameter a represents the relative strength of the random force. In our calculations a is set equal to 0.15. The mean of the random velocity 〈u〉 is chosen to be Vi((1 - a2)1/2 -1). In Appendix B it is shown that this choice ensures that the velocity distribution of the Vˆ i after the transformation is again Gaussian and has the same parameters as the original distribution of the Vi. In Figure 4 the results of simulations for particles in an unstructured tube subjected to a random force are given. Again, the mean square displacements for systems at 100, 300, and 500 K with concentrations as chosen in section IV.A (θ ) 0.2) are calculated. The time behavior of the mean square displacement is completely different from that without random forces. After the ballistic regime which is the same in both cases, there is a broad transition region at t ) 10-11 to 10-9 s where the slope of the curve decreases. For t g 10-9 s the mean square displacement finally becomes proportional to the square root of the time.
320 J. Phys. Chem., Vol. 100, No. 1, 1996
Figure 4. Time dependence of the mean square displacement 〈z2〉 for different temperatures in a straight tube with random forces (θ ) 0.2).
Hahn and Ka¨rger
Figure 6. Time dependence of the mean square displacement 〈z2〉 in a straight tube with random forces for different tube diameters (θ ) 0.2, T ) 300 K).
occupancy,
θ)
Figure 5. Time dependence of the mean square displacement 〈z2〉 for different relative occupancies in a straight tube with random forces (T ) 300 K).
Again, the knowledge of the behavior of a single noninteracting particle is the key to understanding the behavior of the single-file system. The successive transformations Vˆ i ) Vi + u lead to a stochastic motion of the single particle over long times. The correlation coefficient Kn between the initial velocity and the velocity after n transformations is found to be Kn ) (1a2)n/2 (see eq B13). Thus, the velocity of the single particle after e.g. t ) 10-9 s, i.e. after 1000 transformations, is completely uncorrelated with the initial velocity (K1000 ) 1.1 × 10-5). With respect to this time scale, therefore, the particle will follow a random walk, so that the mean absolute shift 〈|s|〉 will be proportional to the square root of the time. For smaller times, obviously the correlation has not yet completely vanished and therefore the proportionality to t1/2 is not reached immediately after the onset of the single-file behavior. In Figure 5, the results for different relative occupancies are presented. The mobilities F as obtained from the simulation are 8.3 × 10-14, 3.4 × 10-13, and 1.4 × 10-12 m2/xs for θ ) 0.8, 0.5, and 0.2, respectively. From eq 9 one would expect that the mobilities for these occupancies follow the proportion 1:4:16. The results of the simulation are in fact in reasonably good agreement with this prediction. This agreement is remarkable, given that the definition of the relative occupancy θ is not without ambiguity. In our definition of the relative
Nσs L
(12)
we assume that the diameter of a particle is equal to the parameter σs and that the particles are arranged strictly in one dimension. However, the minimal distance between two adjacent particles may be as small as 0.8σs, and due to the finite tube diameter, in reality more particles may be packed into a region of a given length than in the truly one-dimensional case. Obviously, more rigorous definitions, which take account of all these influences, would not allow such an easy estimate of the diffusivity as that provided by eq 9. The fact that the magnitude of the channel diameter does not substantially influence the simulation results (as long as the single-file condition is preserved) is demonstrated by the simulation results shown in Figure 6. In all cases, the total number of particles, and hence the relative occupancy according to our operational definition (eq 12), has been kept constant while the effective tube diameter has been enhanced from 0.503 to 0.683 nm. In all cases the simulations yielded the remarkable result that the single-file mobility remains essentially the same. Figure 7 compares the results of molecular dynamics simulation of the mean square displacement in a single-file system with those of an individual (i.e. non-interacting) particle. Via eq 9, the simulation results of the individual particles are used for a prediction of the single-file behavior, which for times t g 3 × 10-11 s (i.e. after the ballistic period) is found to be in excellent agreement with the results of the direct calculation. C. Tube Structured by a Periodic Potential. In the previous case the expected single-file behavior, i.e. 〈z2〉 ∼ t1/2, was due to an additional random force acting on the particles. Now, instead of adding a random force, the diameter of the tube is assumed to vary periodically. For this, the tube radius becomes a function of z:
(2πzλ )
σl(z) ) σl,o + ∆σl cos
(13)
Two cases with different periodicities λ are considered, and in both cases the amplitude ∆σl is chosen in such a way that eq 5 is fulfilled even at the positions with the largest tube diameter. In the first case the periodicity λ ) 6 nm is much larger than the particle diameter (0.383 nm) and the amplitude is set equal
Molecular Dynamics Simulation of Single-File Systems
Figure 7. Time dependence of the mean square displacement 〈z2〉 of the particles in a single-file system and of an isolated (non-interacting) particle and comparison with the behavior of a single-file system predicted via eq 9 from the non-interacting particle. The calculations were carried out for a straight tube with random forces (θ ) 0.2, T ) 300 K).
J. Phys. Chem., Vol. 100, No. 1, 1996 321
Figure 9. Time dependence of the mean square displacement 〈z2〉 for different relative occupancies in a periodic tube potential with λ ) 6.0 nm (T ) 300 K).
Figure 8. Time dependence of the mean square displacement 〈z2〉 for different temperatures in a system with periodic tube potential with λ ) 6.0 nm (θ ) 0.2).
Figure 10. Time dependence of the mean square displacement 〈z2〉 of the particles in a single-file system and of an isolated (non-interacting) particle and comparison with the behavior of a single-file system predicted via eq 9 from the non-interacting particle. The calculations were carried out for a periodic tube potential with λ ) 6.0 nm (θ ) 0.2, T ) 300 K).
to ∆σl ) 0.03 nm. As in the previous cases, σl,0 is chosen to be equal to 0.12 nm. The results for this case are presented in Figures 8-10. After a transition region, which is much smaller than in the case with random forces, at t g10-10 s the mean square displacement becomes proportional to xt. This means that the motion of a single non-interacting particle must imply some randomness. This is true although no explicit randomness is put into the system. The random change of the velocity is now due to the interaction with the tube potential. In contrast to the cases considered in sections IV.A and IV.B the potential depends on z and, therefore, there is a force in the z-direction. In Figure 8 the mean square displacements for different temperatures (100, 300, and 500 K) are given. The dependency of the mean square displacement on the temperature is not as pronounced as in the previous cases. Similarly, a less pronounced dependency on the relative occupancy is found. Figure 9 yields a ratio of 1:2.7:9.1 between the mobilities F, while from eq 9 a ratio 1:4:16 would have been predicted. Finally, Figure 10 demonstrates that eq 9 completely fails in predicting the single-file behavior from the simulations with an isolated (noninteracting) particle. The noninteracting particles behave
like in the case with the straight tube, and no stochastic behavior is obtained in the investigated time scale. This behavior is discussed later in this section. The second case considered in this section is a periodic tube potential with a periodicity on the order of a particle diameter. The periodicity is now set to λ ) 0.4 nm, and the amplitude is chosen to be ∆σl ) 0.02 nm. In this case, generally, the mobilities F are much smaller than in the case of the longer periodicity. This is an immediate consequence of the fact that the z-force due to the interaction with the tube potential is proportional to λ-1. The “random” force acting on the particles is, therefore, stronger for a shorter period length than for a longer period. Again, the dependency on the temperature, as presented in Figure 11, is only weak, as already found in the system with long period length. The dependency on the relative occupancy, however, results to be much more pronounced (Figure 12). The mobilities F now follow the ratio 1:5:25. In contrast to the system with long periodicity, the dependency is stronger than expected from eq 9. At present, we are not able to provide a simple explanation for these differences.
322 J. Phys. Chem., Vol. 100, No. 1, 1996
Hahn and Ka¨rger
Figure 14. Typical path of a particle within one period of the tube potential. In region A the particle becomes faster, and in region B the particle loses nearly the same velocity that it had gained in region A.
Figure 11. Time dependence of the mean square displacement 〈z2〉 for different temperatures in a system with periodic tube potential with λ ) 0.4 nm (θ ) 0.2).
Figure 12. Time dependence of the mean square displacement 〈z2〉 for different relative ocupancies in a system with periodic tube potential with λ ) 0.4 nm (T ) 300 K).
In the case of a periodic tube potential, as in the case of stochastic forces, the molecular mean square displacement in single-file systems is found to increase in proportion with the square root of the observation time. There is, however, no possibility to apply the analytical expressions given in Section III for the prediction of the single-file behavior from the simulations with non-interacting (i.e. isolated) particles. This complication may be most likely attributed to the fact that motion in the z-direction is now influenced by the transverse degrees of freedom. That is, during the interaction with the tube there is an energy transfer between longitudinal and transverse degrees of freedom. For a non-interacting particle in the system with a long periodicity, the motion does not become random in the considered time scale (see Figure 10). Test calculations showed that the velocity of a given noninteracting particle changes quasiperiodically around the initial value but does not become random. Figure 14 schematically represents a typical trajectory. While moving over a period length, the particle often interacts with the tube. Therefore, successive changes of the z-velocity are strongly correlated. The particle gains some kinetic energy while moving through the first half of the period length and loses nearly the same kinetic energy while moving through the second half. Thus, the z-velocity changes periodically as observed in the simulation. The situation changes completely if the particles interact with one another. Then, the normal series of interactions with the tube are disturbed at random positions by the particle-particle interactions, and the velocity changes after a certain number of particle-particle interactions also become random. Note that the randomness results in addition to that caused by the particleparticle interaction itself. It is evidently due to this reason that in the single-file case in contrast to eq 9 the mean square displacement is observed to be proportional to the square root of the observation time, while for the isolated particle 〈z2〉 is proportional to t. In the case of short periodicity the situation is different. Now, the z-forces which act on the particle during successive instants of interaction appear to be uncorrelated and already the non-interacting particle moves randomly. This is illustrated in Figure 15, showing that the mean square displacement of a non-interacting particle increases in proportion with the observation time. Hence, eq 9 correctly predicts the time dependency of the mean square displacement in the single-file system, though with a too high mobility. This deviation may be explained by the particle-particle interactions which lead to additional stochastic interactions with the tube and, hence, to a more random behavior, resulting in a lower mobility. Conclusions
Figure 13. Time dependence of the mean square displacement for periodic potentials with different amplitudes: ∆σ ) 0.01, 0.02, and 0.03 nm. The periodicity is λ ) 0.4 nm (θ ) 0.2, T ) 300 K).
Figure 13 shows the mean square displacements for different amplitudes ∆σ. As expected, the mobilities are found to decrease with growing amplitude of the potential.
For the first time, molecular dynamics simulations have been applied to study the phenomena of single-file diffusion. In systems where the isolated (“non-interacting”) particle moves randomly, the mean square displacement in the single-file system is found to increase in proportion with the square root of the observation time. This is the well-known behavior for singlefile systems, so far predicted by statistical arguments for a system consisting of random walkers. It could be shown, however, that also deviations from such a time dependence are possible: In systems where the noninteracting particle moves in a strictly deterministic way the mean square displacement of
Molecular Dynamics Simulation of Single-File Systems
J. Phys. Chem., Vol. 100, No. 1, 1996 323 and vanishing total momentum is calculated using some simple assumptions. The calculations are done with pointlike particles, and the result is then transformed to include also the case of particles with finite size. Due to the periodic boundaries, the particle positions are repeated periodically
xi ) xi+N - L
(A1)
Consequently, in any region of length L there are exactly N particles and the particle density is found to be
c)
N L
(A2)
In a single-file system the order of particles is the same at all instants of time Figure 15. Time dependence of the mean square displacement 〈z2〉 of the particles in a single-file system and of an isolated (noninteracting) particle and comparison with the behavior of a single-file system predicted via eq 9 from the non-interacting particle. The calculations were carried out for a periodic tube potential with λ ) 0.4 nm (θ ) 0.2, T ) 300 K).
interacting particles (i.e. in the single-file case) becomes proportional to the time. Furthermore, in the cases of a straight tube and of a straight tube with random forces acting on the particles, the simulation results were found to agree quantitatively with a simple correlation rule (eq 9), allowing the calculation of the mean square displacement of the single-file system from that of non-interacting particles. This agreement is not found for a periodic tube potential. In this case the interaction between the particles obviously leads to another sequence of interaction between the particles and the tube. Thus a “stochastic” force is generated which leads to the proportionality between mean square displacement and square root of the observation time. More generally, in the case of a tube with a periodic potential, the disagreement between the correlation rule (eq 9) and the simulations may be attributed to the fact that the simulation is done in three dimensions, whereas the theory assumes a truly one-dimensional system. Since due to the z-dependence of the tube potential the degrees of freedom in the transverse and in the longitudinal direction (i.e. in the direction of diffusion) are coupled, both methods cannot be expected to yield identical results. Since analytical methods fail to describe single-file phenomena in structured potentials quantitatively, molecular dynamics simulations appear to provide the best tool for predicting the diffusional behavior in zeolites with one-dimensional pore structure. Although the computational expense, as explained in section II, is very large, the simulation of more realistic systems (i.e. the use of realistic zeolite potentials) is possible and will be done in the near future.
xN - L < x1 < x2 ... < xN-1 < xN
Due to the vanishing total momentum, the center of mass is conserved. The origin is chosen in such a way that
1
In this appendix the limiting value 〈z2〉∞ of the mean square displacement in a single-file system with periodic boundaries
∑xi ) 0
(A4)
For the calculation of the limiting value of the mean square displacement, it is assumed that the particle distributions at t ) 0 and t f ∞ are completely uncorrelated and that all allowed particle distributions have the same probability. For simplicity, the initial positions at t ) 0 are denoted by xk and the final positions at t f ∞ are denoted by yk. Clearly, eqs A1, A3, and A4 are fulfilled also for yk. To formulate the expression for the mean square displacement, it is easier to use new variables xˆ k and yˆ k, which fulfill only eq A3 but do not fulfill eq A4. The original variables may be expressed in terms of the new ones as
xk ) xˆ k -
yk ) yˆ k -
1
N
1
N
∑xˆ i N i)1 ∑yˆ i
(A5)
N i)1
i.e. this transformation only corresponds to a shift of the origin. Using the new variables, the limiting value of the mean square displacement for the lth particle may be written as
〈z2l 〉∞ ) 1
∫Ldxˆ N∫xˆxˆ -Ldxˆ N-1...∫xˆxˆ -Ldxˆ 2∫xˆxˆ -Ldxˆ 1∫0Ldyˆ N∫yˆyˆ -L × N
2 0
3
N
dyˆ N-1...
2
N
∫yˆ -Ldyˆ 2∫yˆ -Ldyˆ 1 yˆ 3
yˆ 2
N
N
((
yˆ l -
N
1
N
N
)(
∑yˆ i
N i)1
- xˆ l -
N
1
N
∑
))
2
xˆ i N i)1 (A6)
The weight w is defined by
w) Appendix A: Limiting Value of the Molecular Mean Square Displacement in a Finite Single-File System, Subject to Periodic Boundary Conditions
N
N i)1
w Acknowledgment. Financial support by the Deutsche Forschungsgemeinschaft (SFB 294) and by the Alexander von Humboldt Foundation and the Max Planck Society is gratefully acknowledged. We are obliged to S. Fritzsche and R. Haberlandt, Leipzig, to D. M. Ruthven, Fredericton, to D. N. Theodorou, Berkeley, and to R. Q. Snurr, Berkeley/Leipzig/ Evanstone, for stimulating and helpful discussions.
(A3)
ˆ xˆ dxˆ 2∫xˆ -Ldxˆ 1 ∫0Ldxˆ N∫xˆxˆ -Ldxˆ N-1...∫xˆdx-L N
N
3
N
2
(A7)
N
Because of the periodic boundaries all particles are equivalent, and for any l ) 1, ..., N the relation
〈z2〉∞ ) 〈z2l 〉∞
(A8)
324 J. Phys. Chem., Vol. 100, No. 1, 1996
Hahn and Ka¨rger
holds true. Therefore, eq A6 gives the mean square displacement for any given particle and also the mean square displacement 〈z2〉∞. In the following we calculate 〈z2〉∞ by calculating 〈z2N〉∞. The integrals in eqs A6 and A7 can be calculated more easily by using the further substitution
rk ) xˆ k - (xˆ N - L) sk ) yˆ k - (yˆ N - L)
k ) 1, ..., N - 1
Changing the order of summation and integration, the solution of the integrals from eq A13 is straightforward. We give here only the results of the calculations:
I1 )
(A9)
I2 )
Then, the weight (eq A7) becomes
∫0 dxˆ N∫0 drN-1...∫0 dr2∫0 dr1
w)
)
L
L
r3
I3 )
r2
(A10)
〈z2〉∞ ) 〈z2N〉∞ )
∫Ldxˆ N∫0LdrN-1...∫0r dr2∫0r dr1∫0Ldyˆ N∫0d dsN-1... 3
(
1
N-1
1
N-1
3
2
i)1
i)1
)
2
(A11)
The integrand may be rewritten
(
N-1
∑ i)1
N
si -
1 N
N-1
∑ i)1
)
2
ri )
N-1 N-1
1
(
N2
∑ ∑ j)1 i)1
N-1 N-1
rirj +
∑ ∑ sisj j)1 i)1 N-1 N-1
2 N-1
1
(
) 2
N N-1
I4 )
L2N+2 (N - 1)(3N - 7N + 2) 1 (A15) 24 ((N - 1)!)2 N
Inserting these results into eq A14, one obtains
N-1 N-1
∑ ∑ risj) j)1 j)1
〈z2〉∞ )
r2i + 2 ∑ ∑ ∑ i)1 i)1 j)i+1 N-1
rirj +
N-1 N-1
(A12)
Ld 1 ) (1 - θ)2 6 θ
N-1
2
Nd2 1 ) (1 - θ)2 2 6 θ
i)1
I2 ) N-1 N-1
∫0 dxˆ N∫0 drN-1...∫0 dr1∫0 dyˆ N∫0 dsN-1...∫0 ds1 ∑ ∑ rirj i)1 j)i+1 L
r2
L
L
s2
N-1
I3 )
∫0Ldxˆ N∫0LdrN-1...∫0r dr1∫0Ldyˆ N∫0LdsN-1...∫0s ds1 ∑ risi 2
2
i)1
I4 ) N-1 N-1
∫0Ldxˆ N∫0LdrN-1...∫0r dr1∫0Ldyˆ N∫0LdsN-1...∫0s ds1 ∑ ∑ risj i)1 j)i+1 2
2
(A13)
1 1 〈z 〉∞ ) 2 2(2I1 + 4I2 - 2I3 - 4I4) w N
(A14)
(A18)
According to eq A18, the limiting value of the mean square displacement becomes proportional to both the length L and the particle number N if the relative occupancy θ is constant. On the other hand, with a constant particle number the limiting value depends strongly on the relative occupancy. Appendix B: Probability Distribution of a Random Velocity Change Which Conserves the Particle Velocity Distribution The random force acting on the particles is simulated by a random change of the velocity. The transformation of a given velocity V is done by adding a random velocity u, which is Gaussian distributed with
the mean square displacement becomes 2
(A17)
1 L2 〈z2〉∞ ) (1 - θ)2 6 N
∫0Ldxˆ N∫0LdrN-1...∫0r dr1∫0Ldyˆ N∫0LdsN-1...∫0s ds1 ∑ r2i 2
1N-1 (1 - θ)2L2 6 N2
In the limit of a large number of particles, one obtains
Introducing the integrals
L
(A16)
For N ) 1, the mean square displacement vanishes, as it should be for a single particle with vanishing total momentum. The above result was obtained for pointlike particles. For particles with nonvanishing diameter, the result must be modified. The behavior of a single-file system of N particles with diameter d in a length L is the same as the behavior of a system of N pointlike particles in a length L - Nd.33 Therefore, one has to replace the length L by (1 - θ)L, where L is now the length in the system with particles having diameter d. θ ) Nd/L is the relative occupancy. Taking into account this modification, the mean square displacement becomes
N-1 N-1
s2i + 2 ∑ ∑ sisj - 2 ∑ risi - 4 ∑ ∑ risj) ∑ i)1 i)1 j)i+1 i)1 i)1 j)i+1
I1 )
1N-1 2 L 6 N2
2
∫0s ds2∫0s ds1 N ∑ si - N ∑ ri
...
1
1 L2N+2 (N - 1)(2N - 1) 6 ((N - 1)!)2 N
〈z2〉∞ )
2 0
w
1 L2N+2 (N - 1)(N - 2) 8 ((N - 1)!)2
2
LN (N - 1)!
Using the new variables from eqs A9, one obtains from eq A6
1
1 L2N+2 (N - 1) 3 ((N - 1)!)2
g(u,V) )
1
x2πσ2u
(
exp -
)
(u - µu)2 2σ2u
(B1)
Molecular Dynamics Simulation of Single-File Systems
J. Phys. Chem., Vol. 100, No. 1, 1996 325
The mean and width of the distribution g(u,V) may depend on the initial value of the velocity V. The velocity distribution f(V) of the initial velocities itself is Gaussian. The transformation Vˆ f V + u must not change this distribution. To ensure that the Gaussian distribution of the velocities remains the same after the transformation of all particles, the width and the mean of the distribution of u must be chosen in a specific way. Let the width of the velocity distribution f(V) be equal to σ, and let the mean be µ ) 0. (Clearly, this should be fulfilled for an ensemble of particles in thermal equilibrium.) The width σu of the distribution g(u,V) is chosen to be
σu ) aσ
(B2)
where the parameter a represents the relative strength of the random force and must be in the range 0-1. Then, the mean µu must obey the equation
µu ) -V ( Vx1 - a2
(B3)
to fulfill the above condition, that the distribution of the velocities before and after the transformation is the same. With a ) 0 the distribution g(u,V) approaches a Dirac delta function and the original value of V remains unchanged. With a ) 1, independent of the original value V, the transformed value V + u is distributed according to the Gaussian distribution f(V). With the positive sign in the definition of µu the mean of the transformed velocity V + u will be near the original value V, while in the case of a negative sign the mean of V + u is near -V. In the calculations presented in section IV.B only the transformation with the positive sign is used. Now we show that the above transformation conserves the original Gaussian distribution. To conserve the distribution, the number of particles going away from a given V and the number of particles transformed to this V must be the same, i.e. the equation
∫-∞∞du f(V) g(u,V) ) ∫-∞∞du f(V-u) g(u,V-u)
inserted
1
x2πa2σ exp
(
∞ du × ∫ -∞ 2
1 )
(
)
2 2 2Vu - u2 (u - x1 - a (V - u) + (V - u)) 2σ2 2a2σ2
x2πa σ
∫∞ du exp ×
2 2 -∞
2Vua2 - u2a2 - (V - u)2(1 - a2) ( V2 + (V - u)Vx1 - a2
1 )
(
∫∞ du exp 2 -∞
x2πa2σ
and
f(V-u) )
1
x2πσ2 1
)
x2πσ
(
exp -
Kn )
2
∫
µi,j )
∫-∞∞dV0...dVn ViVjw(V0,...,Vn)
The density distribution is found to be n
G(Vi,cVi-1,aσ) ∏ i)1
)
µ0,n ) cnσ2
∫
(
(B12)
Inserting these moments into eq B9, one finally obtains for the correlation
(B6) Kn ) cn ) (1 - a2)n/2
(B13)
References and Notes
)
2
2Vu (u - µu(V-u)) ) du exp -∞ 2σ2 2σ2u x2πσ2u ∞
(B11)
where G(V,m,σ) is the Gaussian distribution of V with mean m and variance σ. c is the abbreviation for (1 - a2)1/2. The moments used in eq B9 are found to be
V2 2Vu - u2 exp 2σ2 2σ2
2Vu 1 ) -∞du exp g(u,V-u) 2σ2 1
(B10)
µn,n ) σ2
2Vu - u 2σ2
( )
(B9)
µ0,0 ) σ2
one obtains from eq B4 the condition ∞
xµ0,0µn,n
w(V0,...,Vn) ) G(V0,0,σ)
2
) f(V) exp
µ0,n
where the moments are defined as
(B4)
(B5)
)
Hence, the condition from eqs B4 and B7 is fulfilled and the transformation conserves the Gaussian distribution f(V). A measure for the randomness of the motion of a particle undergoing n successive transformations is the correlation Kn between the initial velocity V0 and the final velocity Vn. This correlation is
)
( ) ( ( )
exp -
2a2σ2
)
(B8)
(V - u)2 2σ2
(u - V(1 - x1 - a2))2
)1
must hold true. Using
∫-∞∞du g(u,V) ) 1
2a2σ2
(B7)
To show that the right-hand side of eq B7 really becomes equal to 1, the values for σu and µu according to eqs B3 and B2 are
(1) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolites and other Microporous Solids; Wiley: New York, 1992. (2) Rees, L. V. C. Proceedings of the 10th International Zeolite Conference, Garmisch-Partenkirchen, 1994; Weitkamp, J., Karge, H.-G.., Pfeifer, H., Ho¨lderich, W., Eds.; Elsevier: Amsterdam, 1994; p 1133. (3) Chen, N. Y.; Degnan, T. F.; Smith, C. M. Molecular Transport and Reactions in Zeolites; VCH: New York, 1994. (4) Ruthven, D. M. Can. J. Chem. 1974, 52, 3523. (5) Theodorou, D. N.; Wei, J. J. Catal. 1983, 83, 205. (6) Fo¨rste, C.; Germanus, A.; Ka¨rger, J.; Pfeifer, H.; Caro, J.; Pilz, W.; Zikanova, A. J. Chem. Soc., Faraday Trans. 1987, 83, 2301.
326 J. Phys. Chem., Vol. 100, No. 1, 1996 (7) Rajadhyaksha, R. A.; Pitale, K. K.; Tambe, S. S. Chem. Eng. Sci. 1990, 45, 1935. (8) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1994, 98, 11948. (9) Ruthven, D. M.; Eic, M.; Xu, Z. In Studies in Surface Science and Catalysis; Ohlmann, G., Pfeifer, H., Fricke, R., Eds.; Elsevier: Amsterdam, 1991; Vol. 65, p 233. (10) Yashonath, S.; Demontis, P.; Klein, M. L. Chem. Phys. Lett. 1988, 153, 551. (11) Yashonath, S.; Demontis, P.; Klein, M. L. J. Phys. Chem. 1991, 95, 5881. (12) Fritzsche, S.; Haberlandt, R.; Ka¨rger, J.; Pfeifer, H.; Wolfsberg, M. Chem. Phys. Lett. 1990, 171, 109. (13) Demontis, P.; Suffriti, G. B.; Quartierie, S.; Fois, E. S.; Gamba, A. J. Phys. Chem. 1988, 92, 867. (14) Goodbody, S. J.; Watanabe, K.; MacGowan, D.; Walton, J. B. P.; Quirke, N. J. Chem. Soc., Faraday Trans. 1991, 87, 1957. (15) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 8232. (16) Leherte, L.; Vercauteren, D. P.; Derouane, E. G.; Lie, G. C.; Clementi, E.; Andre´, J.-M. Zeolites: facts, figures, future, Proceedings of the 8th International Zeolite Conference, Amsterdam, 1989; Jacobs, P. A., van Santen, R. A., Eds.; Elsevier: Amsterdam, 1989; p 773. (17) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1992, 96, 1051. (18) Hong, U.; Ka¨rger, J.; Kramer, R.; Pfeifer, H.; Seiffert, G.; Mu¨ller, U.; Unger, K. K.; Lu¨ck, H.-B.; Ito, T. Zeolites 1991, 11, 816. (19) Ka¨rger, J.; Keller, W.; Pfeifer, H.; Ernst, S.; Weitkamp, J. Microporous Mater. 1995, 3, 401.
Hahn and Ka¨rger (20) Meier, W. M.; Olsen, D. H. Atlas of Zeolite Structure Types; Butterworths-Heinemann: London, 1992. (21) Ka¨rger, J.; Petzold, M.; Pfeifer, H.; Ernst, S.; Weitkamp, J. J. Catal. 1992, 136, 283. (22) Lei, G.-D.; Sachtler,W. M. H. J. Catal. 1993, 140, 601. (23) Nivarthi, S. S.; McCormick, A. V.; Davis, H. T. Chem. Phys. Lett. 1994, 229, 297. (24) Shen, D.; Rees, L. V. C. J. Chem. Soc., Faraday Trans. 1994, 19, 3017. (25) Fedders, P. A. Phys. ReV. B 1978, 17, 40. (26) van Beijeren, H.; Kehr, K. W.; Kutner, R. Phys. ReV. B 1983, 28, 5711. (27) Ka¨rger, J. Phys. ReV. A 1992, 45, 4173. (28) Nowak, A. K.; den Ouden, C. J. J.; Pickett, S. D.; Smit, B.; Cheetham, A. K.; Post, M. F. M.; Thomas, J. M. J. Phys. Chem. 1991, 95, 848. (29) Verlet, L. Phys. ReV. 1967, 165, 98. (30) Allen, M. P.; Tildesley, D. J. Computer simulations of liquids; Clarendon Press: Oxford, 1987. (31) Haberlandt, R.; Fritzsche, S.; Peinel, G.; Heinzinger, K. Molekulardynamik: Grundlagen und Anwendungen; Vieweg: Braunschweig/Wiesbaden, 1995. (32) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular forces: their origin and determination; Clarendon Press: Oxford, 1981. (33) Hahn, K.; Ka¨rger, J. J. Phys. A: Math. Gen. 1995, 28, 3061.
JP951807U
