# Molecular Dynamics Study of Hydrogen on Alkali-Earth Metal Cations Exchanged X Zeolites

Authors: Du Xiaoming
Journal: International Journal of Chemical Engineering
 (2014)

---

## Abstract

The self-diffusion of hydrogen in Ca2+-, Mg2+- and Ba2+-exchanged X zeolites (Mg46X, Ca46X, and Ba46X) has been studied by molecular dynamics (MD) simulations for various temperatures and loadings. The results indicate that in the temperature range of 77–298 K and the loading range of 1–80 molecules/cell, the self-diffusion coefficients are found to range from 1.2
×
10
-
9 m2·s−1 to 2.3
×
10
-
7 m2·s−1 which are in good agreement with the experimental values from the quasielastic neutron scattering (QENS) and pulse field gradients nuclear magnetic resonance (PFG NMR) measurements. The self-diffusion coefficients decrease with loading due to packing of sorbate-sorbate molecules which causes frequent collusion among hydrogen molecules in pores and increases with increasing temperature because increasing the kinetic energy of the gas molecules enlarges the mean free path of gas molecule. The mechanism of diffusion of hydrogen molecules in these zeolites is transition diffusion. Knudsen diffusion occurs at low loading and the molecular bulk diffusion occurs at higher loading. For given temperature and loading, the self-diffusion coefficients decrease in the order Ba
46
X
<
Mg
46
X
<
Ca
46
X, due to the different sizes and locations of the divalent cations. Moreover, the effect of concentration of molecular hydrogen on self-diffusion coefficient also is analyzed using radial distribution function (RDF).


---

## Body

## 1. Introduction

Zeolites are nanoporous crystalline aluminosilicates with a rich variety of interesting properties and industrial applications [1–3] based on their structural characteristics, that is, high surface area, high thermal and hydrothermal stability, and rigid and regular systems of cavities and channels with dimensions ranging from around 3 Å to 12 Å. These characteristics give rise to their ability to transport various molecules into their channels or cavities and adsorb these molecules on the certain sites. Therefore, a proper understanding of how a molecule diffuses and adsorbs in a nanoporous medium is crucial to the successful application of the materials in adsorption and heterogeneous catalysis.In the past decades, adsorption and storage of molecular hydrogen on microporous zeolites have been investigated theoretically and experimentally [4–6]. However, the knowledge of kinetics and diffusion for hydrogen in microporous zeolites is little, and a very limited number of studies have been reported [7–12]. It is very difficult to obtain accurately the diffusion coefficient of hydrogen molecules in porous materials by using experimental methods. The quasielastic neutron scattering (QENS) and pulse field gradients nuclear magnetic resonance (PFG NMR) techniques were usually used to investigate the diffusion behavior of hydrogen in zeolites [7–9]. Fu et al. [7] have studied the diffusion of molecular hydrogen adsorbed in zeolite 13X at high coverage’s using QENS for temperatures ranging from 10 to 60 K. The diffusive motion of the adsorbed H2 could be described by a liquid like jump diffusion model above 35 K. DeWall et al. [8] also obtained the similar results to Fu et al. Bar et al. [9] reported the results that the self-diffusion coefficients of hydrogen in NaX, NaA, and CHA zeolites decreased with the decrease of the pore diameter for these zeolites using QENS and NMR techniques.However, there is a great difference for experimental results obtained from different methods. Molecular dynamics simulation is a powerful tool to provide insight into molecular-level details of the underlying adsorption and diffusion mechanisms that can guide the future rational design and synthesis of zeolite materials with significantly improved H2 adsorption capability. van den Berg et al. [10] investigated the influence of the cation distribution in sodalite zeolite on the hydrogen self-diffusion using molecular dynamics simulation; they found that the hydrogen self-diffusion coefficient was much larger for the ordered cation distribution than for the disordered distribution. Papadopoulos and Theodorou [11] concluded that the diffusivity of CH4 in NaX zeolites decreases with increasing loading from molecular dynamics simulation. Although several molecular dynamics simulations have been performed [11–13], simulation studies on H2 in ion-exchanged X-type zeolites are very scarce.In this work, we perform a systematic molecular dynamics simulation study on the diffusion properties of H2 in Ca2+, Mg2+, and Ba2+ cations exchanged X-type zeolites (Mg46X, Ca46X, and Ba46X, X = Si100Al92O384). The effect of temperature, loading, and cations on self-diffusivities is analyzed. All results are for diffusion over the temperature range 77–298 K and the loading range of 10–80 molecules/cell. Also, the predicted results are amenable to comparison with mesoscopic and microscopic experimental measurements, such as pulsed-field gradient NMR and quasielastic neutron scattering, respectively, which have been carried out on a variety of zeolite-guest systems [7–9].


## 2. Model and Methods

Structures of initial unit cell of alkali-earth metal cations (Ca2+, Mg2+, and Ba2+) exchanged X zeolites such as Ca46X, Mg46X, and Ba46X (X = Si100Al92O384) were constructed according to the results of X-ray crystallographic studies [14]. The crystals are cubic with space group Fd-3 (number 203). The unit cell parameters of Ca46X, Mg46X, and Ba46X in dehydrated form are a = 25.024, 24.671, and 22.266 Å, respectively. The pore topology of the faujasite zeolites is shown schematically in Figure 1. The structure of the faujasite system used in this investigation consists of large cavities (supercages) which have roughly spherical symmetry and a diameter around 12.5 Å with a window size of 7.4 Å (Figure 1). Each cavity is connected to four others in a tetrahedral arrangement. The structure also contains sodalite cage units linked together by double six rings. Each unit cell of faujasite is composed of 8 supercages/unit cell. The divalent extraframework cations preferentially occupy different crystallographic sites, named I, I′, II, and II′ and depicted in Figure 1 and Table 1 [14]. The simulations were performed using one unit cell with 8 supercages. Test simulations using eight unit cells gave identical results but were deemed too computational expensive to use with the Ewald summation for all the simulations.Table 1
Distribution of cations over sites [14].
Sites
Crystal
Mg46X
Ca46X
Ba46X
I
14
16
14.5
I′
4
0
1.5
II
24
30
30
II′
4
—
—Figure 1
Structure of the X zeolite framework and description of the main crystallographic sites for the extraframework cations.The parametrisation of the hydrogen-hydrogen interaction potential model is described by the 6–12 Lennard-Jones (LJ) potential. A commonly used parametric form for the hydrogen-zeolite interaction potential takes into account the short-range repulsion, dispersion, and electrostatic multipolar interactions by using a sum of Lennard-Jones and Coulombic interactions:(1)
U
s
z
=
∑
i
=
1
N
z
∑
j
=
1
N
(
q
i
q
j
r
i
j
)
+
4
ε
i
j
[
(
σ
i
j
r
i
j
)
12
-
(
σ
i
j
r
i
j
)
6
]
.The interactions between hydrogen and the oxygen, silicon, and aluminum atoms of the framework and the extraframework cations are derived and used in [15–18]. The values of the potential parameters are summarized in Table 2.Table 2
Lennard-Jones interaction parameters.
Atom type
ε (kJ·mol−1)
σ (Å)
Reference
H
0.3039
2.958
[15]
Si
0.1549
0.076
[16]
Al
0.1607
1.140
[16]
O
1.397
3.040
[16]
Ca2+
0.6492
2.98
[17]
Ba2+
1.874
3.82
[17]
Mg2+
0.0820
3.402
[18]The partial charges can be initially assigned on the basis of the electrostatic charge distribution of the bare zeolite and the isolated sorbate. The LJ interaction parameters between unlike atoms can be obtained from the parameters of the like pairs using mixing rules. The Lorenz-Berthelot mixing rules are applied to obtain the cross potentialsε
i
j
=
ε
i
·
ε
j and σ
i
j
=
(
σ
i
+
σ
j
)
/
2.The hydrogen molecule is treated as a centrosymmetrical Lennard-Jones particle, which has proven to be a valid and accurate approximation with respect to extended representation in other studies [19, 20]. To define the zeolite structure by molecular mechanics, most simulation studies are performed using the method proposed by Bezus et al. [21]. In a Kiselev-type potential the atoms from the zeolite framework are fixed at their crystallographic position. In our calculations, the zeolite frameworks were kept rigid and the extraframework cations (Ca2+, Mg2+, and Ba2+) were allowed to move freely in the zeolite. In order to get better starting structures for the barrier optimizations, we performed energy minimization calculations to determine the equilibrium structures using density function theory (DFT) methods that were successfully used to explore hypothetical FAU and LTA with dense topologies. DFT calculations were carried out with the code Quantum-ESPRESSO [22] which uses atom centered basis functions that are particularly efficient for total energy studies of very low density materials with very large unit cells. Periodic DFT calculations were carried out using the PBE exchange correlation functional within the GGA approximation. The minimized zeolites are used as adsorbent for hydrogen adsorption.All MD simulations have been performed using the simulation program code DL_POLY [23]. A simulation box with one unit cell was used. Periodic boundary conditions were applied in three dimensions in order to simulate an infinite system. The canonical (NVT) ensemble was used, with a step size of 1 fs, at 77 K, 150 K, and 298 K, using a Nosé-Hoover thermostat [24]. The initial positions for the hydrogen adsorbate were determined from grand canonical Monte Carlo (GCMC) adsorption simulations. The loading of each starting configuration was determined from adsorption isotherms at 77 K, 150 K, and 298 K to calculate the self-diffusion coefficients of hydrogen inside three ion-exchanged X zeolites. The starting configuration was equilibrated for 50 ps before collecting averages. Molecular dynamics simulations were run for 500 ps. The longer simulation runs (such as 3 ns) yielded results that were not significantly different from the reported results. The nonbond interaction cut-off distance was chosen to be 12.5 Å, 12 Å, and 11 Å which is less than half of the length of the simulation box for Ca46X, Mg46X, and Ba46X, respectively. Test runs with a larger cut-off distance did not yield significantly different results.The primary means of tracking the movement of adsorbate molecules over the course of the simulation was via the mean-squared displacement. The overall mean-squared displacement is(2)
MSD
=
〈
|
r
i
(
t
)
-
r
i
(
0
)
|
2
〉
,
where r
i
(
t
) is the position of molecule i at time t and r
i(0) is an arbitrary starting point.The mean-squared displacement results were used to calculate self-diffusion coefficients for hydrogen in each simulation using the Einstein equation. The Einstein equation for the self-diffusion coefficient in three dimensions is [25]
(3)
D
=
lim
⁡
t
→
∞
1
6
N
a
t
〈
∑
i
=
1
N
a
[
r
i
(
t
)
-
r
i
(
0
)
]
2
〉
,
where D is the self-diffusion coefficient of the molecules. N
a is the number of the diffusing molecule. Furthermore, the self-diffusion coefficient components D
x, D
y, and D
z in x-, y-, and z-direction can be calculated by replacingr with x, y, and z in (3), respectively.The radial distribution function gives a measure of the probability that, given the presence of an atom (molecule, colloid, etc.) at the origin of an arbitrary reference frame, there will be an atom with its center located in a spherical shell of infinitesimal thickness at a distancer from the reference atom. This concept also embraces the idea that the atom at the origin and the atom at distance r may be of different chemical types, say α and β. The resulting function is then commonly given the symbol g
α
β
(
r
) and is defined by Hansen and McDonald [26]
(4)
g
α
β
(
r
)
=
V
s
(
N
α
·
N
β
)
〈
∑
i
=
1
N
α
n
i
β
(
r
)
4
π
r
2
Δ
r
〉
,
where N
α and N
β represent the numbers of α and β particles, respectively. V
s is the model volume. n
i
β
(
r
) represents the number of β particles in the range from r to (
r
+
Δ
r
) at the origin of ith α particle.


## 3. Results and Discussion

The diffusion behavior of hydrogen in Ca46X, Mg46X, and Ba46X at 77 K and 150 K for loadings of 10, 20, 40, and 80 molecules per unit cell was simulated. For 298 K, loadings of 1, 5, 10, and 20 molecules per unit cell were simulated. Figure 2 displays the mean-squared displacement (MSD) of the hydrogen molecule for loading 10 molecules/cell at different temperatures which were presented as a function of time. The relationship between the observed MSD and time is linear to a very good approximation. It indicates that the normal diffusion occurs at this time scale. Also as shown in Figure 2, MSD of the hydrogen molecules increases with increasing temperature for three zeolites, illustrating that the self-diffusion of hydrogen molecules becomes quicker. For the other loadings, MSD curves exhibit the same trend as Figure 2.Mean-squared displacement (MSD) of hydrogen for loading 10 molecules/cell at different temperatures; (a) Ca46X, (b) Mg46X, and (c) Ba46X.
(a)
(b)
(c)The self-diffusion coefficients were obtained from the slope of the curves via a linear regression fit to linear portion of MSD curves to (3). Figure 3 shows the self-diffusion coefficients of hydrogen in Ca46X, Mg46X, and Ba46X zeolites at 77 K, 150 K, and 298 K. For all zeolites, the self-diffusion coefficients decrease as the loading increases. This can be explained by the decrease of mean free path of hydrogen molecules due to increasing loading. The more frequent occurrence of molecular collisions between the different molecules as higher loadings reduces the freedom of diffusion. But considering the influence of the diffusion temperature, the self-diffusion coefficients represent diverse decrease trend. For the high temperature, the self-diffusion coefficients fall drastically as the loading increases from 1 to about 20 molecules/unit cell. However for low temperature, such as 77 K, the self-diffusion coefficients fall smoothly as the loading increases. With further increase in loading, the self-diffusion coefficients reach a plateau. This can be explained by the number of adsorbed molecules in zeolite channels. The saturation of gas adsorption in micropores is rapidly reached under the conditions of the low temperature and the moderate loading. Moreover, the motion of gas molecules is strongly affected by the wall potential at low temperature. This results in the very frequent collision between gas molecules and confining wall. The diffusion of degree of freedom for gas molecules decreases minimally. Therefore, the self-diffusion coefficients are insensitive to increasing loading at the higher loading. However, for higher temperature and lower loading case, the number of gas molecules is small in micropores. And the motion of hydrogen molecules is predominately determined by the kinetic energy. The collision among molecules is predominant. The diffusion motion is quick. As further increase in loading, the collision between hydrogen molecules and pore wall is predominant. As a result, the self-diffusion coefficients rapidly fall. Schuring et al. [27] also reported that the self-diffusion coefficients of alkenes in Mordenite, ZSM-5, Ferrierite, and ZSM-22 zeolite decreased with increase of loading. An interpretation given by them was that increase in loading results in the more frequent collision among gas molecules.The self-diffusion coefficients of molecular hydrogen at different temperatures; (a) Ca46X, (b) Mg46X, and (c) Ba46X.
(a)
(b)
(c)Bar et al. [9] reported that the self-diffusion coefficient of hydrogen molecules in NaX zeolite at 110 K–140 K was in the order of magnitudes of 10−8m2/s using the quasielastic neutron scattering (QENS) and pulse field gradients nuclear magnetic resonance (PFG NMR) methods. These values are in good agreement with self-diffusion coefficients of 2.3 × 10−7–1.2 × 10−9 m2/s obtained by molecular dynamics method in Figure 3 for three zeolites. This indicated that our models and methods used in simulation were reasonable.To investigate the effect of temperature on diffusion coefficients, the self-diffusion coefficients at the same loading were compared for three zeolites (also see Figure3). The self-diffusion coefficients increase with increase in temperature for a given loading. The activation energies for the diffusion of hydrogen in Ca46X, Mg46X, and Ba46X zeolites were calculated by fitting self-diffusion coefficient only for loadings of 10 molecules/cell and 20 molecules/cell to the Arrhenius equation. The obtained values of the activation energies for three zeolites were in range from 0.35 to 2.41 kJ·mol−1. These values are extremely low, as is expected for this small molecule at relatively high temperatures, and, in fact, diffusion is probably not activated at all for H2 in these zeolites. In this case a power law dependence on temperature, as for molecular bulk diffusion, is expected. Transport is more likely closer to Knudsen diffusion at low loading and moves up to bulk diffusion at higher loading when intermolecular collisions are more important.The diffusion mechanism of gas molecules in porous media can be classified into three mechanisms according to the ratio of the pore size of porous media to the mean free path of gas molecules,d
/
λ [28]: for d
/
λ
≤
0.1, it is termed Knudsen mechanism. At low pressure, collisions are dominantly between molecules and the walls, and the free path is restricted by the geometry of the void space. For d
/
λ
≥
10, it is termed molecular diffusion (or Fickian diffusion) mechanism. When the gas pressure is high, molecule-molecule collisions dominate and the system is said to be in the normal or Fickian mechanism. The transition diffusion is defined for 0.1
<
d
/
λ
<
10. The diffusion mechanism is shown in Figure 4.Diffusion mechanisms: (a) Knudsen, (b) Fickian, and (c) Transition.
(a)
(b)
(c)According to kinetic molecular theory of gas, the mean free path of gas molecules is defined as(5)
λ
=
k
B
T
2
π
d
0
2
p
,
where k
B is the Boltzmann constant, T is temperature (K), d
0 is the effective diameter of the gas molecule (m), and p is pressure (Pa).With (5), the mean free path of H2 in Ca46X, Mg46X, and Ba46X zeolites under the conditions of temperatures and loadings can be calculated. The pore size of these zeolites is about 0.12 nm. The dependence of d
/
λ on the loadingsfor H2 in Ca46X, Mg46X, and Ba46X zeolites at 77 K as sample was shown in Figure 5. It was found that for these zeolites the values of d
/
λ is in range from 0.1 to 10 for full range of loadings (except the loading of 80 molecules/unit cell for Ba46X zeolite) in this study. This shows that the mechanism of diffusion of hydrogen molecules in these zeolites is transition diffusion. That is, the molecular bulk diffusion (or Fickian diffusion) and Knudsen diffusion are both important during hydrogen diffusion. At low loading, the mean free path of hydrogen molecules is much larger than the pore diameter of zeolites and there is less molecules in pore, which increases probability of collisions between hydrogen molecules and the walls. Therefore, collisions between hydrogen molecules and the walls are dominant. Transport is more likely closer to Knudsen diffusion. As loadings increase, the mean free path of hydrogen molecules becomes smaller and the number of hydrogen molecules in pore increases. The intermolecular collisions are more important than that between hydrogen molecules and the walls. The diffusion moves up to bulk diffusion. Also as seen from Figure 5, the values of d
/
λ were maximum for Ba46X zeolite and were minimum for Ca46X zeolite. It is indicated that the collisions between hydrogen molecules and the walls (or intermolecular collisions) are the most frequent. The diffusion of hydrogen molecules is the slowest of three zeolites.Figure 5
The dependence ofd
/
λ on the loadingsfor H2 in Ca46X, Mg46X, and Ba46X zeolites at 77 K.Figure6 shows the comparison of self-diffusion coefficients of hydrogen in three zeolites at loading of 10 molecules/cell and different temperature. It was found that self-diffusion coefficients of hydrogen on alkali-earth metal cations exchanged X zeolites decreased in the order Ba46X < Mg46X < Ca46X. This can be explained from the difference of sites and size of divalent cations in X zeolites. The number and size of the cations locating site II, site II′, and site III determine the diameter of windows of X zeolites. For alkali-earth metal cations exchanged dehydrated X zeolites in this work, the 30 Ca2+ cations are located at site II in the supercage for Ca46X; the 24 Mg2+ cations and 4 Mg2+ cations are located at site II and site II′, respectively; for Mg46X and the 30 Ba2+ cations are located at site II in the supercage for Ba46X, also as shown in Table 1. The total number of three cations at sites II and II′ is equivalent. Since the ionic radius of Mg2+ (0.66 Å) is smallest and the ionic radius of Ba2+ (1.35 Å) is largest, the windows of Ba46X were seriously blocked. The diffusion of hydrogen in Ba46X zeolite was the slowest of three zeolites. In X-ray crystallographic data of Mg46X provided by Yeom et al. [14], Mg46X zeolite was partially dehydrated. Mg46(H2O)4-X was the final structure. The existence of partial water molecules in Mg46X zeolite may decrease windows size and pores volume, which hinders the diffusion of hydrogen. Hence, the self-diffusion coefficients of hydrogen in MgX zeolite were less than that of CaX zeolite.Figure 6
Comparisons of the self-diffusion coefficients of hydrogen in three zeolites for loading 10 molecules/cell.In order to investigate the effects of loading, temperature, and cation type on the distance among hydrogen molecules, the radial distribution functions (RDF) of the centers of mass for hydrogen molecules in zeolite pores at different loading, temperatures, and zeolites were shown in Figures7, 8, and 9, respectively. As shown in Figure 7, the loading has little influence on the peak positions of RDF for CaX zeolite. Whereas the peak values of RDF curves decreased and the peak shape broadened with increasing loading. It is indicated that the concentration of hydrogen molecules in a local area decreases, and random distribution density increases, as the loadings increase. This enhances the interactions among hydrogen molecules. The collisions among hydrogen molecules become more frequent during diffusion. Hence, the self-diffusion coefficient decreases continuously. In Figure 8, with the increase of temperature, the peak position of RDF curves for CaX zeolite moves to the right. Meanwhile the peak value decreases, and the peak shape broadens. It is indicated that at high temperature hydrogen molecules cluster evolves from a more compact and ordered packing to incompact and disordered packing. The reason is that the kinetic energy of the gas molecules increases with increasing temperature. This can increase the mean free path of gas molecule, showing the very quick diffusion in the case of the higher temperature. Figure 9 shows the effect of cations in X zeolite on RDF of hydrogen diffusion at 77 K and 10 molecules/cell. RDF of hydrogen diffusion in Mg46X zeolite was similar to Ca46X zeolite. But there is no peak on RDF curve of hydrogen diffusion for Ba46X. The reason may be that the pores in Ba46X zeolite are smaller. The concentration of hydrogen molecules can be changed in such confined space, which causes the decrease of the self-diffusion coefficient for Ba46X zeolite. Yang et al. [29] thought that the self-diffusion coefficient of Lennard-Jones fluid in confined space increases with the increase of limited scale.Figure 7
Radial distribution functions of the centers of mass for hydrogen molecules for 77 K and Ca46X zeolite at different loading.Figure 8
Radial distribution functions of the centers of mass for hydrogen molecules for 10 molecules/cell and Ca46X zeolite at different temperature.Figure 9
Radial distribution functions of the centers of mass for hydrogen molecules in three zeolites at 77 K and 10 molecules/cell.


## 4. Conclusions

The diffusion behaviors of molecular hydrogen in Ca2+-, Mg2+-, and Ba2+-exchanged X zeolites (Mg46X, Ca46X, and Ba46X) were investigated by using molecular dynamics (MD) simulations for the temperatures 77 K, 150 K, and 298 K and the loading range of 10–80 molecules/cell. Our molecular dynamics results show that the self-diffusion coefficients of hydrogen in three zeolites decrease with increase in loading and increase with increase in temperature, and they decrease in the order of Ba46X < Mg46X < Ca46X, due to the different size and location of the divalent cations. The mechanism of diffusion of hydrogen molecules in these zeolites is transition diffusion. Knudsen diffusion occurs at low loading and the molecular bulk diffusion occurs at higher loading. Moreover, the effect of concentration of molecular hydrogen on self-diffusion coefficient also is analyzed using radial distribution function (RDF). It is found that as the loading and temperature increase, the concentration of hydrogen molecules in a local area decreases, and random distribution density increases. The collisions among hydrogen molecules become more frequent. Hence, the self-diffusion coefficient decreases continuously. The results and data of hydrogen diffusion properties obtained from the simulations are theoretically significant for understanding of the mechanism of hydrogen storage on microporous zeolites.



---

*Source: 701057-2014-02-23.xml*