13554
J. Phys. Chem. B 2006, 110, 13554-13559
A Lattice-Gas Cellular Automaton to Model Diffusion in Restricted Geometries P. Demontis,* F. G. Pazzona, and G. B. Suffritti Dipartimento di Chimica, UniVersita` di Sassari, Via Vienna, 2, I-07100 Sassari, Italy ReceiVed: March 22, 2006; In Final Form: May 16, 2006
Topological constraints play a fundamental role in problems involving the transfer of heat, matter, or information in a discrete network. Flows in microporous media are the most common cases in which the thermodynamic and transport properties of adsorbate are strongly influenced by its interactions with the confining medium. Combining two local Monte Carlo operations and a Lattice-Gas Cellular Automaton, we constructed an equilibrium synchronous network of cells highly sensitive to the state of their neighborhood and able to capture the effects of confinement by means of few flexible local parameters and a parallel, local evolution rule. Results of an application of the model to the specific problem of geometrical restricted long-range diffusion are used to interpret the behaviors of particles in tight confinement.
1. Introduction Geometrical restrictions at the level of molecular length scale strongly influence the thermodynamic and transport properties of molecules subjected to the interaction with the confining medium.1 The physical phenomena occurring in these systems are multifarious including heterogeneous catalysis, biological transport, percolation, and even a dramatic change of the phase diagram.2 In this environment, the molecular mobility is ruled by the topology of its surrounding medium, which provides the energy landscape for the transport process. Increasing the number and type of molecules in the accessible volume deeply modifies this energy landscape, destroying or creating effective paths for molecular motion. Despite a great deal of effort, these phenomena, associated with diffusion in tight confinement, are still far from being understood.3 This gives rise to the quest for models, which are simple enough to be analyzed and able to capture the essential features of the real physical systems.4 In recent years, Lattice-Gas Cellular Automata (LGCA) have been used to model a wide range of physical problems, including fluid dynamics, oscillating chemical reactions, and heat flow. They consist of a network of identical cells exchanging information. At each time step, each cell is updated synchronously to all others according to an evolution rule depending only on the state of its neighborhood.5,6 In this paper, we develop a new and general LGCA consisting in a d-dimensional network of hypercubically arranged cells (applications to other spatial arrangement will follow the same lines), and then we specialize our results on the case study of diffusion in a LTA-type zeolite.7,8 We prove that our approach allows one to study the collective diffusion of species bounded to an adsorbent’s surface in large systems for large observation times avoiding the shortcomings of kinetic Monte Carlo simulations as recently pointed out.9 The key idea lies in a coarse-graining of space and time scales and, simultaneously, a reduction of the degrees of freedom. The coarse-graining of spatial scales rests on the assumption that the system can be partitioned into connected cells and each cell contains an appropriate number of sites that realize a representative number of states depending on the * To whom correspondence
[email protected].
should
be
addressed.
E-mail:
number of occupied sites. Thus, the system consists of a lattice of cells with the unit size λ. The reduction of the degrees of freedom rests on the assumption that instead of the coordinates of the particles, only a fixed number of states accessible to the particles in each cell is considered. The dynamics along the coordinates within the cell are approximated by the dynamics of the states. Each cell will be characterized by the same fixed number of exit and inner sites that sets up the allowed maximum number of occupied sites. The actual dynamic state of each cell (including its tendency to particle transfer) will be determined by the energies of its occupied sites according to a local MC sampling scheme on the allowed configuration space. In this way, one can introduce an elementary time scale larger than the time scale of the slowest processes occurring in the cell, and the equilibrium properties will depend on the value of the few parameters that completely define the thermodynamics of the model. 2. Description of the Model In this work, the concept of geometrical restriction in a transport problem has been applied to create an LGCA in which the cells exchange simultaneously a limited number of particles. Our model is constructed in such a way as to give a coarsegrained picture of a system of particles diffusing in a microporous material (such an LTA-type zeolite) in which the dimensions of the particles are such that each window connecting two adjacent cells is too narrow to allow the simultaneous intercell transfer of more than one particle through it. The main peculiarities of the model stand in the fact that the cells of the network (structured and connected, as shown in Figure 1a,b) are designed so that their exchange ability is highly sensitive to their internal state, and they all update their states simultaneously. Each cell contains two kinds of sites: exit sites available to particle exchange with neighboring cells and inner sites not available. They are different not only in the exchange ability but also in the statistical weight because energies -ex, -in (ex, in > 0) to exit and inner sites are assigned, respectively. The energy parameters according to temperature characterize the tendency of the single cell to retain its content or to transfer all or part of it to neighboring cells. Additional constraints can be added to the exchange process by means of a transmission
10.1021/jp061783z CCC: $33.50 © 2006 American Chemical Society Published on Web 06/17/2006
LGCA to Model Diffusion in Restricted Geometries
J. Phys. Chem. B, Vol. 110, No. 27, 2006 13555 The system evolves in discrete time steps of arbitrary duration τ. Each cell r ∈ L contains K g ν distinguishable sites labeled as (r, i), i ) 1,..., K; we assume that an exclusion principle applies so that the site (r, i) has only two states coded by the Boolean variable si(r, t) ) 1, 0 if the site is occupied or not, respectively, at time t. Therefore, the occupancy of cell r at time t is defined as K
n(r, t) ) Figure 1. (a) Two representations of a unit cell with ν ) 6, K ) 16; three-dimensional sketch (up) and string representation (down). Grey numbered cubes represent exit sites, whereas the sphere represents the ensemble of inner sites. Assigned numbers estabilish the correspondence between exit sites in the 3D sketch (grey cubes) and in the string representation (grey squares). The topology of exit sites in one cell is the same as the topology of the cells in the entire system; therefore, it must be specified. The inner sites are distinguishable, but their spatial arrangement is ignored (this is why we employ a sphere for their 3D representation and white numbered squares in the string). (b) Spatial arrangement of cells in a 3D system portion.
probability κ between each pair of adjacent cells, which can affect the frequencies of transfers. Several levels of dependence of the parameters ex, in, κ on local observables are also possible, to enhance or depress the cell sensitivity. Thus, the cells exchange their contents (particles, heat, or information in general) by means of a dynamics strictly dependent on the way we model the small set of parameters defined above. To demonstrate the basic ideas, we applied the model to the problem of diffusion of particles in tight confinement, and we found that even maintaining a relatively low level of cell sensitivity (i.e., all parameters constant), the model shows several diffusion profiles depending on temperature, whereas a completely different behavior is observed by simply inverting the assigned values of ex, in. At low temperatures, the model predicts a condensation of particles in the cells that gives rise to a phase transition between a low- and high-mobility regime: this is mirrored by the appearance of a singularity on the reduced variance of the occupancy number. The critical loading (or concentration) at which the transition occurs is strictly connected to the architecture and energetics of the single cell: its value depends on the average ability of the cell to make accessible an exit site if ex < in (i.e., when the inner sites bind more strongly the guest particles than the exit sites), whereas its value depends on the ability to vacate an exit site if ex > in (i.e., when all of the particles tend to occupy exit sites, thus making transfer a rare event). 2.1. Structure of the Network. Particle diffusion in a microporous material is a suitable framework in which our model can be illustrated. The model system consists of N identical, distinguishable particles diffusing in a d-dimensional hypercubic lattice L of Ld cells with periodic boundary conditions. Each cell is labeled by the coordinates of its center r ∈ L. The lattice has connectivity ν ) 2d: that is, each cell in L is connected to, and can exchange particles with, its ν nearest-neighboring cells labeled as r + λej, j ) 1,...,ν where e1,...,ed is the basis set of d orthogonal unit vectors of the lattice, ei+d ) -ei, i ) 1,..., d, and λ is the fixed distance between centers of adjacent cells. Therefore, we define the neighborhood I(r) of cell r as the set constituted by the cell r plus the cells with which r can exchange particles: I(r) ) {r + λej}j ) 0,...,ν, i.e., the same cell and its nearest-neighbor cells (e0 is the null vector).
si(r, t) ∑ i)1
(1)
ranging from 0 to K, whereas the configuration of sites is stored in the vector s(r, t) ) {si(r, t)}i)1,...,K. The geometrical restriction mentioned in Sec. 2 is represented by imposing a constraint on the particle traffic: each pair of adjacent cells, say r and r + λej with j ∈ [1,ν], can exchange at most one particle at each time step. Let us explain how to implement this condition. Let r be a reference cell. In the lattice, the cell r is surrounded by the first-neighboring cells r + λe1, r + λe2, ..., r + λeν. The cell r has K sites; ν of them are exit sites. Each exit site connects r to a particular neighboring cell: for each j ∈[1,ν], (r, j) and (r + λej, j + d) are the only (adjacent) sites connecting the adjacent cells r and r + λej they belong to, respectively. Then the constraint on particle traffic arises naturally in the following form: at each time step the cell r: • can employ only its (r, 1) site to exchange one particle with only the site (r + λe1, 1 + d) belonging to the r + λe1 cell and vice versa; • can employ only its (r, 2) site to exchange one particle with only the site (r + λe2, 2 + d) belonging to the r + λe2 cell and vice versa; • ... • can employ only its (r,ν) site to exchange one particle with only the site (r + λeν,ν + d) belonging to the r + λeν cell and vice versa. Analogous conditions are imposed in the same way to all cells in L . In this notation, the sums j + d are intended as sums modulo ν (see Figure 1a). As a consequence of this constraint, the change in the occupation of a cell between two consecutive time steps is always an integer in the interval [-ν,ν]. 2.2. Structure of the Cell. Given that a cell can host a maximum of K particles, and ν sites among K are used to allow particle transfers to neighboring cells, a cell in L consists of ν exit sites available to particle exchange with neighboring cells, and K - ν inner sites not available. In a cell containing n particles, a number of K!/n!(K - n)! accessible configurations of n particles in K site is allowed. To make each cell sensitive to the amount of particles it contains, we differentiate the statistical weight of the two kinds of sites by assigning them a different energy (lattice-gases with two nonequivalent sites were extensively studied by Chvoj et al. and Tarasenko et al. in refs 10 and 11 and by Yashonath et al. in ref 16). We denote -ex, -in (where ex, in > 0) the energy of exit and inner sites, respectively. Therefore, the allowed configurations are differently weighted depending on the temperature T. To measure the tendency of the cell to transmit particles outside, we define in a cell of occupancy n the accessibility Rn as the equilibrium probability that a specific exit site is occupied, divided by the occupancy n. If ex ) in, exit and inner sites are equivalent and a direct calculation gives Rn ) 1/K independent of n and T (compare this value with Fickian diffusivity for the diffusive model in ref 6). Setting ex * in, the accessibility Rn varies with n. This
13556 J. Phys. Chem. B, Vol. 110, No. 27, 2006
Demontis et al.
is enough to generate a wide range of behaviors of the system. Although inside a cell the K sites are distinguishable, the cell interior is not as detailed as a lattice except in the sites that join two cells. As can be seen in Figure 1a, only the topology of its exit sites is specified (because they connect the cell with its topologically well defined neighborhood), whereas we choose to ignore all about spatial arrangement of its inner sites. This is equivalent to neglecting the time-scale relative to intracell motion, because we are interested only in intercell (long-range) diffusion. In this way, at each time step we simulate an intracell configuration change by mixing the guest particles according to a MC sampling scheme13 in which each particle can try to jump in a site chosen among all K sites inside the cell. This approach is less time-expensive and allows the system to relax toward equilibrium faster than the finer procedure in which the particles can jump only between neighboring sites in a cell interior detailed as a lattice. In doing so, ex, in represent the parameters that model the mean interactions of the n(r, t) guest particles among themselves and with the host cell. The model is flexible enough to permit a setting of the energy parameters ex, in as functions of properly chosen local observables to represent averaged local multiparticle interactions. As long as we are interested in the thermodynamics of equilibrium systems, all parameters (the energies ex, in and the transmission probability κ) may not exhibit any explicit time-space dependence but must be chosen in such a way that the system preserves ergodicity. For the sake of simplicity, in this outline of the model we will refer to ex, in, and κ as constants. Further generalizations are also possible. An important consequence of neglecting part of the internal cell structure is that the location of each particle is defined by the position of its containing cell, and inside of each cell, the adsorbed particles are in fact undistinguishable. 2.3. Time Evolution. We suppose that the system is placed in a large heat bath at temperature T, so that the model generates a random walk in the canonical ensemble. In a cell the occupied site i has potential energy - i ) -ex, -in if it is an exit or inner site, respectively. Following Chvoj et al.,10 to study the system thermodynamics, we do not consider (i) the kinetic energy of diffusing particles (constant at given T), (ii) the proper energy of substrate, and (iii) the entropy of individual particles connected with their internal degrees of freedom. Thus, the Hamiltonian of a cell with J occupied exit sites and n - J occupied inner sites results
HJ,n ) -Jex - (n - J)in
(2)
and the Hamiltonian of the entire system at time t reads H(t) ) ∑r∈LHJ,n(r, t). The evolution rule E follows the typical scheme of LGCA models,6 i.e., at each time step it consists of two moves: • a local mixing rule to change the internal configuration of a cell according to its present state only and acting simultaneously and independently on each cell by means of a randomization operator R • a propagation rule to generate transfers among adjacent cells, where the occupancies of all cells are updated synchronously, by means of a propagation operator P . Therefore, the evolution operator reads E ) R o P, and the system evolves according to a RP-dynamics as intended in ref 6. The mixing consists of the application of a randomization operator R ) R pos o R id that first changes the configuration of sites (R pos) and then mixes (R id) the identities of particles inside of each cell.
Figure 2. (a) Example of randomization in the cell of Figure 1a with n ) 3. Inside the cell and independently of other cells, the R pos operator performs a Monte Carlo procedure on the input configuration of occupied sites to generate a new configuration accordingly to the cell Hamiltonian defined in eq 2. Next, the identities of the particles in this new configuration are mixed by the R id operator. (b) Scheme of the propagation. In this cross section of the system, the cell represented in Figure 2a is highlighted at the center and surrounded by other cells. The squares are the exit sites numbered as in Figure 1a. Big numbers inside the white circles indicate the number of particles residing in the inner sites of each cell. Black circles indicate particles. Each particle in an exit site will attempt a jump along the direction of the arrow with probability κ (see text for details) if the destination site is empty. Otherwise, no intercell migration can occur (as indicated by the cross).
At each cell r ∈ L, the R pos operator changes the particle positions by means of site-to-site jumps according to an Arrhenius-type dynamics12 where the activation energy of a jump is the energy barrier i that the particle residing in the (r, i) site has to overcome in jumping to a different site, so that the jumps are treated as activated processes. This operation can be performed by the following Monte Carlo algorithm:13 (i) select the K sites in r at random order; (ii) for each selection, let the selected site be the leaVing site i; a destination site f among all K sites belonging to r is selected, then a particle jump from i to f occurs with probability pRif ) si(r, t)(1 - sf(r, t))φi. φi is the probability that the particle has enough energy to overcome the energy barrier i of the leaving site i. We set φi as11
( )
φi ) A exp -
i kBT
(3)
where A is a normalization constant and kB is the Boltzmann constant. The optimal choice for A11 is A ) exp(min/kBT), where min ) min(ex, in). This algorithm produces an output configuration that satisfies conservation and distinguishability of particles. We recall that in step (ii) of randomization, each particle can choose a destination site among all remaining sites in the cell, because we did not specify any intracell structure. As inside of the unit cell the n adsorbed particles must be undistinguishable, the R id operator permutes their identifications in a way randomly chosen on the n! possible (and equiprobable) ones. Therefore the effect of the operator R : s f sR is to transform an input cell configuration s to an output one sR among the K!/n!(K - n)! possible ones according to the Hamiltonian (2). The randomization applies simultaneously and independently to all cells in L and satisfies detailed balance:10,11 that is, for each pair of sites i, f belonging to the same cell, 〈pRif - pRfi 〉 ) 0. Figure 2a summarizes the randomization operation. After randomization, the propagation operator P allows the intercell particle transfers by acting simultaneously on all pairs
LGCA to Model Diffusion in Restricted Geometries
J. Phys. Chem. B, Vol. 110, No. 27, 2006 13557
of adjacent exit sites connecting adjacent cells. Let a, b be such a pair of sites. Let their occupancies after randomization be sRa , sRb , respectively. Then the operator P : sRa , sRb f sPa , sPb swaps the occupancies sRa , sRb with probability pPab + pPba where
pPab ) sRa (1 - sRb )κ
(4)
The parameter κ acts as a transmission probability between the leaving cell (to which the site a belongs), and the destination cell (to which the site b belongs). We can formulate κ as
( )
κ ) κ0A exp -
ex kBT
(5)
where A is the same as in eq 3, and ex in the Boltzmann factor is the proper choice because the propagation operation involves always an exit site. Here κ0 ∈ (0;1] ∈ R is an additional barrier to intercell migration. If we set κ0 ) 1, then no additional barrier acts on diffusion. If ex, in, κ are constants then it is easy to compute the Grand-Canonical partition function of the single cell, thus obtaining predictions on thermodynamic properties and (using a probabilistic approach under some simplifying assumptions about the correlations in the random walk motion) self-diffusivity. Under all these conditions, the propagation satisfies detailed balance: if a and b are two adjacent exit sites connecting two adjacent cells in L , then 〈pPab - pPba〉 ) 0. The propagation operation is summarized in Figure 2b. 3. Results and Discussion The most simple system to study is one in which ex, in, κ are constant parameters. To test our LGCA on a model of a real system, we map on the lattice an LTA-type zeolite, the ZK4, a framework that satisfies the topology requirements of the model illustrated here.14 This system consists of a simple cubic array of nearly spherical cavities with an internal radius of about 5.7 Å connected to six neighboring cavities by nearly circular windows of about 4.2 Å in diameter. We assume that each cavity can host a maximum of K ) 16 particles. In our previous MD simulations15 on this system, we observed that very high local density fluctuations are one of the ways in which simple fluids confined in porous media manifest their altered dynamic behavior. We evidenced also that the molecules adsorbed in the R-cages of ZK4 spend most of their time oscillating around the adsorption sites. We found four different regimes at different time scales: quasi-free motion on a small time scale (up to 0.2 ps), oscillations about the preferred positions in the cavity (0.2-3 ps), intra-cavity diffusion (some tens of picoseconds), and long-range diffusion. In this last condition, the long-time motion of the adsorbed molecules can be treated as stochastic, as the time spent in each R-cage is long enough to cancel, to a great extent, any correlation with the past history of the migration dynamics. Thus, adsorbed molecules make a random walk between cages whose properties may be very much complicated by the presence of other molecules and the interactions with the surrounding walls of the cage. This sharp separation of time scales between the elementary transition processes (atomic motion between adsorption sites in the same cage) and the time lapse between transitions to different cages, whose sequence leads to the diffusion phenomena in which we are interested, enables us to represent diffusion on our lattice of cells.
On the basis of these considerations, we will treat the migration process by invoking the Markovian hypothesis that the current state of the system is fully determined by its state at one past time rather than by its entire past history. We stress that in our model, we neglect all the details of the intracage motion. We reduce all of the physical processes of diffusion in a microporous material into a sequence of particle jumps from one cell to a neighboring cell. Therefore, to represent the intercage migration, the only quantity we have to determine at each time step is the number of particles that are leaving each cell, then which particles are actually entering to a neighboring cell. We simulate this process by means of the twostep rule illustrated before; focusing on a cell at a particular time step, (i) we decide how many particles will be able to try a jump to neighboring cells, and (ii) we perform the allowed jumps. A particle will be able to try a jump to a neighboring cell if it reaches an exit site, and the jump from this exit site to the exit site of the destination cell is allowed with the probability κ defined in eq 5 if and only if the destination (exit) site is empty. In our model, we distinguish between the energies -ex, -in of the exit and inner sites inside of a cell to have a simple criterion to mimic a behavior reminiscent of that in the real systems. In a “real” ZK4,15 the complex dynamics determined by the interactions between molecules and between molecules with the zeolite’s inner surface will have as a final result that some molecules will be in the condition to attempt to leave the R-cage whereas some others certainly will not. Although in systems of noninteracting hard spheres all of the configurations of the guest particles in a cell of occupancy n have the same statistical weight, the presence of different energies between the exit and inner sites inside of a cell causes different configurations to be differently weighted. In their work, Yashonath et al.16 made use of a differentiation in binding energies of some adsorbing sites in a Metropolis Monte Carlo model of particles diffusing in a framework of two-dimensional confining multiparticle cells; Chvoj et al. and Tarasenko et al.10,11 studied an Arrhenius Monte Carlo model of particles diffusing on a surface with two types of sites nonequivalent in adsorption energies. In both models, the differentiation among sites allows to sample differently weighted configurations of particles in the adsorption sites. In this way, depending on the difference in energy, different kinds of diffusion profiles with loading are generated. In our model, we apply this differentiation between exit and inner sites to see how different accessibilities of exit sites influence the dependence of the diffusivity on loading. To test our model, the size of the system and observation time-scale were chosen to get a high statistical accuracy in the computation of equilibrium averages and diffusivities. A total of 128 numerical simulations were carried out on a grid of Ld ) 323 cells with ν ) 6, K ) 16, κ0 ) 1, each one spanning a time scale of 105 time steps. We found nonnegligible finite size effects for grid sizes less than L ) 8. Following Yashonath et al.,16 the energy parameters were taken to be ex ) 10 kJ/mol and in ) 20 kJ/mol for exit and inner sites, respectively. In Figure 3, curves of diffusivity vs loading 〈n〉 at several temperatures are shown. Similarly to ref 16, changes in temperature led to different types of diffusional behavior, which correspond to the I, II, IV, and V types observed by Ka¨rger and Pfeifer17 in the PFG-NMR measurements of intracristalline self-diffusion coefficient depending on sorbate concentration. Let us consider first the trend of D vs 〈n〉 with ∆ ) in ex ) 10 kJ/mol. At T ) 300 K, the curve for 0 < 〈n〉 j 11
13558 J. Phys. Chem. B, Vol. 110, No. 27, 2006
Figure 3. Self-diffusivity D vs loading 〈n〉 for a system of 32 × 32 × 32 cells with ν ) 6, K ) 16, and κ0 ) 1. The energy parameters are ex ) 10 kJ/mol and in ) 20 kJ/mol, independent of n. The simulation results refer to the following temperatures: T ) 300 K (9), T ) 460 K (b), T ) 600 K (2), T ) 2500 K (1). (Inset) Effect of inverting the energy parameters (i.e., setting ex ) 20 kJ/mol, in ) 10 kJ/mol) on the diffusivity. For a direct comparison, in the upper part of the figure, the five different diffusivities observed by Ka¨rger & Pfeifer17 are reported.
shows a behavior of type V. For higher loadings because the system approaches saturation by increasing the loading, D reaches a maximum around some njmax and then decreases to zero, so that in the 〈n〉 interval, 8 < 〈n〉 e 16 the trend is of type IV. At 〈n〉 ) njmax, the average accessibility of an exit site and the probability of having two adjacent exit sites simultaneously occupied are balanced so as to give the greatest probability of migration. By removing a small number of particles from the system, the average accessibility decreases and the diffusivity decreases too. By adding a small number of particles to the system, the average accessibility crosses the threshold above which the probability of two simultaneously occupied adjacent exit sites results enhanced, therefore lowering the probability of migration. The trend of this curve is qualitatively analogous to the behavior of diffusivity of methane in ZK4 at T ) 300 K.18 The type III diffusivity trend is proper to a situation in which the balance persists for some interval njmax < 〈n〉 < nj′max, where nj′max is a loading above which the diffusivity must decrease to zero as the system saturates. The case study of Figure 3 does not cover this behavior because, to obtain such a curve, the choice of constant values of ex, in is not sufficient, and a finer tuning of the accessibility Rn is needed (e.g., modeling the dependence of in on the number of the occupied inner sites) to depress the average accessibility for njmax < 〈n〉 < nj′max. If we look at the behavior of the diffusivity at T ) 460, 600 K, we note that for low loadings, the curve raises with the temperature, whereas in the end of the curve, the diffusivity slightly decreases. The reason is the following: for low loadings, most of the cells have low occupancy; given that in > ex, an increase in temperature implies a gain in accessibility of exit sites because the term ∆/kBT decreases and the statistical weights of exit and inner sites become more similar. However, as the loading is low, the probability of having two adjacent exit sites simultaneously occupied (such an event allows no transfer) is also low. Therefore, the gain in accessibility of exit sites produces a large gain in diffusivity for low loadings. For high loadings most of the cells have high occupancy; a temperature increase raises the accessibility of exit sites, but at the same time, as the loading is high, the probability of having two adjacent exit sites simultaneously occupied is also high, and this causes a decrease of the diffusivity with the temperature for high loadings.
Demontis et al.
Figure 4. Reduced variance for the same system as in Figure 3. Dashed line: T f ∞, which corresponds to in ) ex. Solid lines: from up to down, T ) 522, 261, 174, 130, 105 K. Around 〈n〉 ) 10 (which corresponds to K - ν), if the temperature is sufficiently low, a phase transition occurs. At T ) 174 K, the transition is diffuse but becomes more neat by further temperature lowering. (Inset) Reduced variance for inverted energies (see Figure 3). In this case, the critical loading is 〈n〉 ) 6, corresponding to ν.
For T ) 2500 K, the behavior is of type I: this is because the high temperatures tend to minimize the difference in the statistical weight of exit and inner sites and the system approaches the ideal condition represented by T f ∞, where the particles diffuse as noninteracting mutually exclusive particles; under such conditions, a linearly decreasing trend of D vs 〈n〉 is expected19 because the accessibility Rn is constant. In conclusion, for ∆ > 0, by raising the temperature, the D vs 〈n〉 curve seems to vary continuosly from the trend described above in the case of T ) 300 K to the trend of type I at T ) 2500 K, passing through the type II at some intermediate temperature. The inset of Figure 3 shows how inverting the energy parameters (i.e., setting ∆ ) -10 kJ/mol) produces diffusivity trends substantially different from the previous case. From 〈n〉 ) ν to K, the diffusivity can only decrease because each cell has all or almost all of the exit sites occupied, so the particles hinder each other during the propagation operation. As can be seen, at the temperatures T ) 300, 460, 600 K, the trends are essentially of type IV, whereas by increasing the temperature, the curves seem to get near the type II. According to the theoretical model developed by Chvoj et al.,10 we found that, depending on the Boltzmann factor exp(-|∆|/kBT), at a certain value of 〈n〉, a phase transition between a low and a high diffusivity regime occurs. We describe this phase transition by characterizing first the ideal phases between which the phase transition proceeds for the system under study when ∆ > 0 and κ ) 1. According to ref 10, the whole system can be divided in two subsystems: the particles occupying the exit sites form the first phase, and the particles occupying the inner sites form the second phase. Then we study both phases as separate systems. The exit sites belong to the phase characterized by the energy barrier ex whereas the inner sites belong to the phase characterized by the energy barrier in. These characteristics do not depend on the number of particles in the respective subsystems. At low loadings, most particles are in the inner sites and the system behaves predominantly as the isolated second phase. For 〈n〉 > K - ν, most of the inner sites will be occupied, and diffusion proceeds only via the exit sites: the system behaves like the isolated first phase. Therefore, close to 〈n〉 ) K - ν, we observe the phase transition between the first and the second phase. For T f 0, the two phases are completely separated, so that the phase transition is a point transition; this reflects on a singularity (a cusp) in the reduced variance σ2/〈n〉 ) (〈n2〉 〈n〉2)/〈n〉 (see Figure 4), which is a measure of the thermody-
LGCA to Model Diffusion in Restricted Geometries namic tendency to accept a new particle in the cell.15 When the temperature is increased, we observe a diffuse phase transition as intended in ref 10; because the equilibrium state of the system consists of the first and second ideal phases coexisting in their respective equilibrium concentrations, the transformation of the second into the first phase proceeds not exactly at loading 〈n〉 ) K - ν, but more properly within a finite 〈n〉 interval around 〈n〉 ) K - ν. The curve of the reduced variance around this loading becomes smooth. Because we can no longer properly assign a transition point, we shall speak about a diffuse transition around 〈n〉 ) K - ν. Moreover, we found that diffusion profiles and thermodynamic properties change completely with the sign of ∆. For low T, if ∆ > 0, the particles will tend to saturate the inner sites. Then, increasing 〈n〉, the diffusivity regime changes from low to high diffusivity around 〈n〉 ) K - ν, i.e., the loading at which the probability of finding a cell with occupancy n ) K - ν + 1 is not negligible (in such a cell, it is impossible to avoid that a particle reaches an exit site). If instead ∆ < 0, the particles will preferably tend to occupy exit sites (see insets of Figures 3 and 4). Then, by increasing 〈n〉, the diffusivity collapses to very low values around 〈n〉 ) ν, where the probability that a cell will have all ν exit sites occupied is large. Therefore, the probability of having two adjacent exit sites simultaneously occupied will increase, and particle migration will be a rarer event. Thus, if ∆ < 0, similarly to the previous case, at low loadings, most particles will occupy the exit sites and will diffuse through them: the system behaves predominantly as the isolated first phase; for 〈n〉 > ν, most of the exit sites will be occupied, and the system behaves as the isolated second phase. Therefore, the diffuse phase transition takes place around 〈n〉 ) ν. 4. Conclusions In conclusion, our LGCA is able to capture the most essential features of confinement by means of a simple probabilistic scheme that satisfies detailed balance by means of its evolution rule consisting of two consecutive Monte Carlo moves locally applied to each cell of the network at each time step. Although we have presented calculations only for a specialized case, the LGCA introduced here may have a scope for wider applications. Indeed, our work introduces a stochastic modeling and computational framework, easy to code on a parallel machine, capable of effectively sampling larger length and time scales
J. Phys. Chem. B, Vol. 110, No. 27, 2006 13559 than conventional simulations while still retaining the essential microscopic details. Acknowledgment. This work has been carried out with financial support provided by Italian Ministero dell’Istruzione, dell’Universita` e della Ricerca, by Universita` degli Studi di Sassari, and by INSTM. References and Notes (1) Klafter, J.; Drake, J. M. Molecular Dynamics in Restricted Geometries; John Wiley and Sons: New York, 1989. (2) Huwe, A.; Kremer, F.; Behrens, P.; Schwieger, W. Phys. ReV. Lett. 1999, 82, 2338. (3) Conner, W., Fraissard, J., Eds.; Fluid Transport in Nanoporous Materials: Proceedings of the NATO AdVanced Study Institute, held in La Colle sur Loup, France, 16-28 June 2003; NATO Science Series II: Mathematics, Physics and Chemistry - Vol. 219; Springer: Berlin, Germany, 2006. (4) Coppens, M.-O.; Bell, A. T.; Chakraborty, A. K. Chem. Eng. Sci. 1998, 53, 2053. (5) Chopard, B.; Droz, M. Cellular Automata Modeling of Physical Systems; Cambridge University Press: Cambridge, England, 1998. (6) Boon, J.-P.; Dab, D.; Kapral, R.; Lawniczak, A. T. Phys. Rep. 1996, 273, 55; Dab, D.; Lawniczak, A. T.; Boon, J.-P.; Kapral, R. Simul. Pract. Theory 2000, 7, 657. (7) Beerdsen, E.; Smit, B.; Dubbeldam, D. Phys. ReV. Lett. 2004, 93, 248301. (8) Tunca, C.; Ford, D. Chem. Eng. Sci. 2003, 58, 3373. (9) Beerdsen, E.; Dubbeldam, D.; Smit, B. Phys. ReV. Lett. 2005, 95, 164505. (10) Chvoj, Z.; Conrad, H.; Cha´b, V.; Ondrejcek, M.; Bradshaw, A. M. Surf. Sci. 1995, 329, 121; Chvoj, Z.; Conrad, H.; Cha´b, V. Surf. Sci. 1997, 376, 205. (11) Tarasenko, A. A.; Jastrabı´k, L. Surf. Sci. 2002, 507-510, 108; Tarasenko, N. A.; Tarasenko, A. A.; Bryknar, Z.; Jastrabı´k, L. Surf. Sci. 2004, 562, 22. (12) Katsoulakis, M. A.; Majda, A. J.; Vlachos, D. G. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 782. (13) Novotny, N. A. In Annual ReViews of Computational Physics IX; Stauffer, D., Ed.; World Scientific: Singapore, 2001; pp 153-210. (14) Fritzsche, S.; Haberlandt, R.; Ka¨rger, J.; Pfeifer, H.; Wolfsberg, M. Chem. Phys. Lett. 1992, 198, 283. (15) Demontis, P.; Fenu, L.; Suffritti, G. B. J. Phys. Chem. B 2005, 109, 18081. (16) Bhide, S. Y.; Yashonath, S. J. Chem. Phys. 1999, 111, 1658; Bhide, S. Y.; Yashonath, S. J. Phys. Chem. B 2000, 104, 2607. (17) Ka¨rger, J.; Ruthven, D. M. Diffusion in zeolites and other microporous materials; John Wiley and Sons: New York, 1992. (18) Dubbeldam, D.; Beerdsen, E.; Vlugt T. J. H.; Smit, B. J. Chem. Phys. 2005, 122, 224712. (19) Auerbach, S. M.; Jousse, F.; Vercauteren, D. P. in Computer Modelling of Microporous and Mesoporous Materials; Catlow, C. R. A., van Santen, R. A., Smit, B., Eds.; Elsevier: Amsterdam, 2004; pp 49108.
