Heitjans P., Karger J. (Eds.). Diffusion in Condensed Matter: Methods, Materials, Models

Подождите немного. Документ загружается.

924 Reinhold Haberlandt

create

start situation

calculation:

new sites and

velocities

calculation:

fo rces

evaluation

end of the run ?



yes

end

Fig. 23.4. Schematic representation of the procedure of MD simulations (left) and

of the periodic boundary conditions. (right)

algorithms to solve the equations of motion (e. g. [4–7,10–13]). Fig. 23.4 (left)

shows the principal organization of an MD program (Verlet algorithm).

To calculate the velocities in each step by the often used velocity-Verlet-

Algorithm one needs the knowledge of the already calculated forces and also

of those corresponding to the new sites. Therefore it is suitable to calculate

the new velocities in two steps: considering ﬁrst that part, which needs only

the former forces and second the ﬁnal new velocities using the new calculated

forces.

Here – as an example only – this often used Verlet algorithm (cf. (23.20b))

and its velocity version (cf. (23.20c,23.20d)) are outlined. Adding to the Tay-

lor series

(t + h)=r

(t)+v

h +

(t)+... (23.20a)

its analogue for r

(t −h), all terms with odd powers of the time increment h

will be cancelled and one gets up to terms of the fourth order

(t + h)=2r

(t) − r

(t − h)+h

(t). (23.20b)

As the velocity version of the Verlet algorithm one obtains in a similar way

(t + h)=r

(t)+hv

(t)+

(t), (23.20c)

23 Statistical Theory and Molecular Dynamics of Diﬀusion in Zeolites 925

(t + h)=v

(t)+

(t)+F

(t + h)]. (23.20d)

23.3.3 Methodical Hints

For the complete wealth of tricks in handling MD simulations one clearly has

to consult the relevant literature (see, e. g., [5, 6]). However, the following

methodical hints aim to summarize a few of the most basic procedures:

– Periodic boundary conditions – replication of the MD box by periodic

continuation (see Fig. 23.4 (right)) – can help to overcome artiﬁcial surface

eﬀects due to spatially restricted systems.

– To avoid unnecessary calculations the interaction is taken into account

only in appropriate short distances (potential cutoﬀ and shifted forces).

To take into account long-range interactions, special techniques (e. g. the

Ewald summation) must be used.

– To save computer time, the time steps should be as large as possible (but

the energy must be maintained) and as short as necessary (fast processes).

In some cases multiple step lengths are recommended (multiple time step

method).

– In order to avoid further time consuming calculations the number of pair

separations explicitly considered is reduced (Verlet neighbour list, up-

dated from time to time).

23.4 Simulation of Diﬀusion in Zeolites

23.4.1 Introduction

Diﬀusion processes in gases, liquids, solids and at interfaces have been studied

with diﬀerent theoretical and experimental methods for a long time. The

present textbook provides a broad survey of this subject.

Bulk measurements of (self-) diﬀusion coeﬃcients D

over wide ranges of

temperature T are often described generally by the Arrhenius relation (cf

Chap. 1)

= D

exp



−



(23.21)

with the pre-exponential factor D

, an activation energy E

(in a certain

sense). This relation can be explained using the transition theory of Eyring,

but its range of validity is not quite clear.

In this section, MD simulations will be applied in helping to understand

the features of molecular diﬀusion in more microscopic detail, namely under

the inﬂuence of the conﬁned geometries in microporous materials. These mi-

croporous materials – especially zeolites – have attained a steadily increasing

926 Reinhold Haberlandt

role for both applications and fundamental research [14,20]. The study of dif-

fusion processes in zeolites is of great interest because these crystals contain

very regular internal surfaces and they have been used for many industrial

purposes [14].

Since the paper of Yashonath et al. [21], MD simulations (see Sect. 23.3,

[5, 6]) are applied to diﬀusion in zeolites to discuss the dynamics of kinetic

processes in zeolites in more detail. In the recent literature one ﬁnds some

reviews (e. g. [22–30]) and papers of several groups [39–81]).

MD simulations can use diﬀerent models for guest molecules, lattices and

their interaction and allow, especially, variations in the system parameters

that are not possible in experiments. So, interrelations and dependences can

be examined.

In this section it will be discussed how MD simulations can answer ques-

tions like, e. g.:

1. Does there exist a macroscopic behaviour (i. e. thermalization, Maxwell-

Boltzmann distribution) of the few particles diﬀusing in the narrow holes

of zeolites?

2. What is the range of validity of the Arrhenius relation (23.21)?

3. How are the evaluated diﬀusion data inﬂuenced, e. g., by

a) small changes of parameter sets (structural, potential)?

b) idealization of the models used (neglecting of cations, ﬁxed lattice)?

c) taking into account concentration gradients?

d) technical points with respect to the MD simulation?

4. How one can use additional theoretical tools – e. g. propagators – to

improve the understanding of diﬀusion processes in zeolites?

5. Can MD simulations at all describe experimental NMR data and/or

QENS data of diﬀusion?

6. What is the result for mixtures?

23.4.2 Simulations

To answer the above mentioned questions the calculations should at ﬁrst be

restricted to simple model systems. Thus, in the beginning of this section,

the diﬀusion of (simple spherosymmetric) methane guest molecules through

idealized zeolites of LTA-type neglecting the exchangeable cations (for short:

ZK4) is examined. Then we will proceed with more complex guest molecules,

mixtures and take into account the cations of the zeolite lattice additionally.

Possible improvements of the model and application to mixtures of guest

molecules will be discussed later on.

Fig. 23.5 (left) shows the general structure of zeolites of type LTA used

for the calculations. The sodalite units form a cubic lattice with large cavities

connected by so-called windows consisting of eight oxygen atoms.

23 Statistical Theory and Molecular Dynamics of Diﬀusion in Zeolites 927

Fig. 23.5. Structure of zeolites of type LTA (left: general view; right: interior view),

see text.

In Fig. 23.5 (right) the distribution of the lattice atoms around a large

cavity in the NaCaA zeolite is to be seen. Lattice atoms in front have been

removed in order to open the view into the interior of the cavity (windows

are marked by w).

One of the crucial points in determining diﬀusion coeﬃcients in zeolites is

the knowledge of the inter- and intramolecular potentials. The forces acting

between the guest molecules and the zeolite lattice as well as among the guest

molecules must be known as good as possible in order to solve Newton’s

equations for the guest molecules in zeolites.

These forces can be derived from interaction potentials (for more details

see, e. g., [1, 22–24] and references therein). Potentials can be obtained from

quantum chemical calculations [24] and/or from ﬁt procedures using experi-

mental results (for an illustrative example see [31] and [32–34]).

Alternatively – at least in the future –, the forces and potential energies

can be obtained by including quantum mechanical density functional calcu-

lations into an otherwise classical MD run [35–38]. However, such methods

are extremely computer-time consuming and can yield the trajectories only

for short times and/or small systems until now. So, they cannot be used for

long-time diﬀusion processes in large systems.

Here, for the interaction between the guest molecules CH

and the lattice

atoms of a cation-free analogue of LTA zeolites (‘ZK4’), we again use (23.18),

the well-known Lennard-Jones (LJ) (12,6) potential.

For NaCaA zeolites, additional polarization energy terms due to the

cations have been included. Using the conventional microcanonical MD en-

semble the self-diﬀusivity is calculated, while a generalized non-equilibrium

(NEMD) ensemble is used to determine the transport-diﬀusivities [48].

In order to get appropriate diﬀusion coeﬃcients, long-run simulations were

carried out with up to 6 000 000 time steps with short step lengths of 5 and 10

fs, respectively. The basic MD box contains between 8 and 343 large cavities

928 Reinhold Haberlandt

with occupation numbers varying between 1 and 11 per cavity, respectively.

The velocity-Verlet-algorithm (see (23.20c), (23.20d)) is used [5, 6].

23.4.3 Results

Computer simulations not only yield the diﬀusion coeﬃcients as the ﬁnal

results, but, additionally, residence times, velocity autocorrelation functions

(vacf), potential surfaces, density distributions, and – last but not least –

propagators and Van Hove correlation functions. Thus a better understanding

of the diﬀusion processes is possible.

To demonstrate the inﬂuence of the potential the results of two sets of

potential parameters A and B (see Table 23.2, p. 932) are compared.

Applicability of the Diﬀusion Equation

First of all one must answer the non-trivial question whether or not systems

of a few particles in narrow-sized pores obey the diﬀusion equation. In the

kinetic theory of ﬂuids this equation has been proven rigorously only in the

case of dilute gases. But, it is possible to check it here in the considered case

by computer experiments [42]. The validity of the diﬀusion equation (see,

e. g., Chaps. 1, 10, 18-20, and [2]) for an inﬁnite system can be examined by

the propagator

P (r, r

,t,0) = (4πDt)

−3/2

exp



−(r − r

)

4Dt



. (23.22)

P (r, r

,t,0) is the (conditional) probability density to ﬁnd a particle at a

site r at time t, which started at t =0atthesiter

. It is the solution of the

diﬀusion equation

∂P/∂t = D∆P ∆ : Laplacian (23.23)

with the initial condition P (r, r

,t=0)=δ(r − r

To examine the form of the solution of the diﬀusion equation, it is neces-

sary and suﬃcient that all moments f(r)

 of f(r)=(r − r

) (cf. (23.22,

23.24a)) will give the same diﬀusion coeﬃcient. Applying the deﬁnition

f(r)

 =



P (r, r

,t)f(r)

dr (23.24a)

we ﬁnd for the ﬁrst four moments ν =1, 2, 3, 4 of this distribution



| r − r



, (23.24b)



(r − r

)



=6Dt , (23.24c)

23 Statistical Theory and Molecular Dynamics of Diﬀusion in Zeolites 929

Fig. 23.6. D

from the ﬁrst four moments of the displacement; left hand side:

methane in a cation-free LTA zeolite T = 300 K, parameter set A (see Table 23.2), 3

methane per cavity, right hand side: methane in a NaCaA zeolite, set B, T = 173 K,

7 methane per cavity.



| (r − r

)



32(Dt)

3/2

√

, (23.24d)



(r − r

)



= 60(Dt)

. (23.24e)

From each one of these equations an expression for Dt can be obtained, e. g.

for the ﬁrst moment

Dt = π | r − r

|

/16 (23.24f)

D =

| r − r

|

. (23.24g)

The diﬀerentiation (cf. (23.24g)) was used instead of simple division by t

(cf. (23.24f)) because one gets a faster convergence of the results. Corre-

spondingly one gets for the second moment – the mean square displacement

(MSD) –

D =



(r − r

)



. (23.24h)

This equation is usually used to obtain the diﬀusion coeﬃcient from MD

simulations. If the diﬀusion equation is valid then the D values derived from

the diﬀerent moments should agree [42]. In Fig. 23.6 (left) an example of the

diﬀusivities in the cation-free analogue of NaCaA (‘ZK4’) is given. It can be

seen that there is good agreement for times larger than about 50 ps. That

means that in this case the hydrodynamical state is reached and we conclude

that (23.22) is then valid. The diﬀusion equation can be simply obtained by

comparing space and time derivatives of (23.22) and, hence, it must be valid,

too. In [42] this was also checked for other occupation numbers.

If the diﬀusion coeﬃcients are smaller – as in the case of NaCaA (see

Fig. 23.6 (right)) – the hydrodynamical state is reached later (for NaCaA in

about 6 000 ps) [44]. This is probably also due to the additional cations.

930 Reinhold Haberlandt

Residence Times and Velocity Autocorrelation Functions

In trying to understand the diﬀusion behaviour in more detail we deﬁne

the residence time for a given guest molecule as the time diﬀerence between

successive passages of limiting planes situated in the centre of each window

perpendicular to the window axis.

Fig. 23.7 shows the probability density of the thus deﬁned residence times

within the individual cavities for three diﬀerent sorbate concentrations, i. e.

loadings I (given in molecules per supercage), and temperatures, respectively

(parameter set A, see Table 23.2). The ﬁrst maximum at ∼ 0.3 ps corresponds

to times which are too short to allow a passage through the cavity to one

of the ﬁve other windows [43]. Therefore, this maximum must be attributed

to trajectories, which are reversed immediately after the molecule has passed

the window.

It is interesting to note that the intensity of this ﬁrst maximum increases

with increasing sorbate concentration. It may be concluded, therefore, that

the reversal in the trajectory is mainly caused by the inﬂuence of the other ad-

sorbate molecules. The second maximum corresponds to the passage through

the adjacent windows and is therefore responsible for diﬀusion.

These conclusions are conﬁrmed by minima in the velocity autocorrelation

function (vacf) (see [43]) where the dependence of these minima on temper-

ature and loading are in agreement with the above reasoning (see Fig. 23.8).

Fig. 23.1 shows a typical situation in a system where methane molecules

are enclosed in cavities of the cation-free LTA zeolite. The large cavities and

the windows connecting them are indicated by grey isopotential surfaces.

To illustrate the time development all methane molecules that were initially

situated in the left hand lower cavity are marked by dark spheres while the

other ones are fair. So the exchange can be observed.

p (a.u.)

Fig. 23.7. Probability density of residence times of particles as function of time

for diﬀerent loadings I (left) and diﬀerent temperatures T (right).

23 Statistical Theory and Molecular Dynamics of Diﬀusion in Zeolites 931

Fig. 23.8. Velocity autocorrelation function (vacf) for diﬀerent loadings (left) and

temperatures (right).

Density Distributions

The density distribution in cation-free LTA zeolites shown in Fig. 23.9 may

well be understood by isopotential lines (Fig. 23.10 and [43]) for the zeolites.

The density distributions – demonstrated in a plane through the centre of the

large cavity – show a remarkable structure. Corresponding to the potential

surface these distributions are diﬀerent from zero practically only near the

cavity wall and have maxima in the window (set A, cf Table 23.2) and in

front of them (set B), respectively.

ÅÅ

Fig. 23.9. Density distribution of CH

in a cation-free LTA with a loading of

I =3atT = 300 K monitored in a plane through the cavity centres (left: set

A; right: set B, cf Table 23.2) [51].

Inﬂuence of the Potential Parameters on D

The inﬂuence of potential parameters on the diﬀusion coeﬃcients is remark-

able. Even small parameter changes may cause signiﬁcant changes in the

932 Reinhold Haberlandt

Table 23.2. Parameter sets used for the LJ potential (cf. (23.18)).

zeolite σ in

A  in kJ/mol

LTA CH

− CH

3.817 1.232

LTA CH

− Si 2.14 0.29

LTA CH

− O (set A) 3.14 1.5

LTA CH

− O (set B) 3.46 0.81

LTA C

− O 3.775 1.536

silicalite CH

− CH

3.730 1.230

silicalite CH

− O 3.214 1.108

silicalite Xe − Xe 4.064 1.870

silicalite Xe − O 3.296 1.679

diﬀusion processes (see Fig. 23.11). It has been shown that the choice of

the σ parameter of the Lennard-Jones potential, (23.18), for the methane–

oxygen interaction has a dramatic inﬂuence not only on the value of the

diﬀusion coeﬃcient but on its concentration dependence. In [43], two diﬀer-

ent Lennard-Jones potentials – based on the sets A,B of potential parameters

shown in Table 23.2 – have been compared.

Analogously to the density distribution (Fig. 23.9), Fig. 23.10 shows the

shape of the potential surface and visualizes the diﬀerent behaviour in the

vicinity of the window, resulting from the parameter set A (left) and set B

(right), respectively. It can be seen that the potential values are high in the

centre of the large cavity and, of course, at the repulsive walls. The potential

Fig. 23.10. Potential surface plotted through one quarter of the planes considered

in Fig. 23.9 (left: set A; right: set B, cf Table 23.2).

23 Statistical Theory and Molecular Dynamics of Diﬀusion in Zeolites 933

Fig. 23.11. Diﬀusion coeﬃcient D for diﬀerent loadings of methane in the cation-

free LTA zeolite (left: diﬀerent temperatures, set B; right: T = 173 K, set A,B)

has a minimum in the window in set A, while for set B the maximum is in

front of the window and a saddle point in the centre of the window. This

threshold reduces the diﬀusivity in set B. From the larger value of σ in set

B, for the dominating oxygen–guest interaction in the zeolite there results a

narrower window between adjacent cavities. Especially for methane, which

has a similar size as the window, this causes dramatic changes. This eﬀect

seems to be the consequence of the interrelation of two eﬀects (reﬂection at

the inlet of the window, coming back after passing the window) with diﬀerent

density dependence [43, 51].

Fig. 23.11 (left) shows that the diﬀusion coeﬃcient D increases with in-

creasing mean number of guest molecules per cavity and temperature for set

B while this dependence for set A – demonstrated for T = 173 K – is re-

versed (right). Both cases show an Arrhenius-like behaviour (cf. (23.21)) (for

non-Arrhenius behaviour see, e. g., [67,68]).

For higher loadings this ﬁgure shows an interesting cross-over of the two

curves. Increasing diﬀusivities with increasing concentration is in fact the

behaviour found for paraﬃn in A type zeolites experimentally by PFG NMR

(see Chap. 10).

Inﬂuence of Lattice Vibrations on the Diﬀusion Coeﬃcients

The inﬂuence of lattice vibrations on the diﬀusion coeﬃcient in the cation-free

LTA zeolite is not very large for both of the parameter sets under consider-

ation as can be seen in Fig. 23.12 (left) [22, 58, 60, 66, 67]. Thus, the result,

initially given by Suﬀritti and Demontis [47] was corrected, who found a much

stronger eﬀect using a parameter set beyond the region of A and B. (Since

their parameters lead to a practically closed window, it is most likely that an

occasional opening of the windows by vibrations may drastically increase D

in such cases.) Moreover, Fig. 23.12 again shows how the potential parameter

σ may aﬀect the inﬂuence of another feature of the system (viz. the lattice

ﬂexibility) on the diﬀusion coeﬃcients.