# Molecular-Dynamics Simulation of Self-Diffusion of Molecular Hydrogen in X-Type Zeolite

Authors: Xiaoming Du
Journal: Journal of Chemistry
 (2013)

---

## Abstract

The self-diffusion of hydrogen in NaX zeolite has been studied by molecular-dynamics simulations for various temperatures and pressures. The results indicate that in the temperature range of 77–293 K and the pressure range of 10–2700 kPa, the self-diffusion coefficients are found to range from 1.61 × 10−9 m2·s−1 to 3.66 × 10−8 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 increasing pressure due to packing of sorbate-sorbate molecules which causes frequent collusion among hydrogen molecules in pores and increase with increasing temperature because increasing the kinetic energy of the gas molecules enlarges the mean free path of gas molecule. The activated energy for hydrogen diffusion determined from the simulation is pressure-dependent.


---

## Body

## 1. Introduction

Zeolite having the faujasite (FAU) framework structure type is a microporous material with regular crystalline structure. It is widely employed as molecular sieves in various industrial processes, such as adsorption, catalysis, purification, and separation of mixtures [1, 2]. The synthetic low silica FAU zeolites, namely, X-type zeolites, have the property that the diameter of the pores can be controlled by changing the size and charge of the exchangeable cations. This offers a variable internal surface area, which is a useful attribute for the physisorption of gases.Recently, zeolites are attractive materials for the safe and efficient storage of hydrogen [3–5]. One of the potential advantages of using zeolites as hydrogen storage materials is their thermal and chemical stability, which has been shown to be useful in catalytic applications. This is exemplified by their lack of combustibility both in air and, in particular, in a hydrogen atmosphere. The hydrogen storage capacities of X-type zeolites have been reported experimentally and theoretically [6–9]. However, little attention was previously paid to the transport properties of adsorbed hydrogen molecules which play a central role in catalytic and separation processes [10] that take place within zeolite cavities [11]. Understanding the host-guest interactions that control molecular diffusion may suggest new materials with advanced performance. We have recently calculated [12] self-diffusion coefficients for hydrogen in all-silica ZSM-5 over the temperature range of 77–293 K and the pressure range of 14–3360 kPa. The self-diffusion coefficients are found to range from 1.2 × 10−9 m2/s to 3.8 × 10−7 m2/s. The diffusion of the hydrogen in all-silica ZSM-5 zeolite channels is anisotropic. Fu et al. [13] have studied the diffusion of molecular hydrogen adsorbed in zeolite 13X at high coverage using quasielastic neutron scattering (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. [14] also obtained similar results to Fu et al. Bar et al. [15] 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. Van Den Berg et al. [16] investigated the influence of the cation distribution in sodalite zeolite on the hydrogen self-diffusion using molecular-dynamics simulation; they found that the hydrogen diffusion coefficient was much larger for the ordered cation distribution than for the disordered distribution. However, very little work on diffusion behavior of hydrogen in FAU zeolites as a function of loading was done so far.In this paper, we investigate the effect of temperature and pressure on self-diffusivities determined using molecular-dynamics method for hydrogen molecules adsorbed in X-type zeolite. All results are for diffusion over the temperature range of 77–293 K and the pressure range of 10–2700 kPa. Computational data on the diffusivities of hydrogen in zeolites are presented and compared with the experimental data.


## 2. Model and Methods

Structure of initial unit cell of NaX zeolite (Na88Si104Al88O384) was constructed according to the result of X-ray crystallographic studies [17]. The crystal is cubic with unit cell parameters a=25.099 Å and space group Fd3- (number 203). The pore topology of the NaX zeolite is shown schematically in Figure 1. NaX is a 3-dimensional zeolite and has a 12-membered ring pore opening of 7.4 Å in diameter providing access to a supercage of 12.4 Å in diameter. In dehydrated NaX zeolite, there are four Na+ cation sites. 2.7, 26.6, 29.3, and 29.3 Na+ cations are located at the site I in (the hexagonal prism), I′ (in the β-cage, close to the hexagonal window to the hexagonal prism), site II (located in the sixth ring of the β-cage), and site III (in the supercage, close to a square window between two other square windows), respectively [17].Figure 1
Structure of the NaX zeolite framework and description of the main crystallographic sites for the extraframework cations.All MD simulations have been performed using the simulation program code DL_POLY [18]. DL_POLY is a package of molecular simulation routines written by Smith and Forester, copyright CCLRC for the Central Laboratory of the Research Councils, Daresbury Laboratory at Daresbury, Nr. Warrington (1996). Lennard-Jones potentials are used for adsorbate-adsorbate and adsorbate-framework interactions. The interactions between hydrogen and the oxygen, silicon, and aluminum atoms of the framework and the extraframework cations are derived and used in [19–21]. The values of the potential parameters are summarized in Table 1. The Lorenz-Berthelot mixing rules are applied to obtain the cross potentials εij=εi·εj  and σij=(σi+σj)/2.Table 1
Lennard-Jones parameters used for hydrogen adsorption on NaX zeolite.
Atom type
ε (kcal/mol)
σ (Å)
Reference
H
0.0726
2.958
[20]
Si
0.0370
0.076
[21]
Al
0.0384
1.140
[21]
O
0.3342
3.040
[21]
Na+
0.1008
1.746
[22]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 [22, 23]. To define the zeolite structure by molecular mechanics, most simulation studies are performed using the method proposed by Bezus et al. [24]. In a Kiselev-type potential the atoms from the zeolite framework are fixed at their crystallographic position. In our calculations, the zeolite framework was kept rigid and the nonframework cations (Na+) 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 [25] 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.For the simulations, 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 0.001 ps, at 77 K, 160 K, 195 K, 230 K, and 293 K, using a Nosé-Hoover thermostat [26]. 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, 160 K, 195 K, 230 K, and 293 K and 10–3000 kPa [27] to calculate the self-diffusion coefficients of hydrogen inside NaX zeolite. The starting configuration was equilibrated for 50 ps before collecting averages. Molecular-dynamics simulations were run for 500 ps. The longer simulation runs yielded results that were not significantly different from the reported results. The nonbond interaction cut-off distance was chosen to be 13 Å which is less than half of the length of the simulation box. 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(1)MSD=〈|ri(t)-ri(0)|2〉,
where ri(t) is the position of molecule i at time t and ri(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 [28]
(2)D=limt→∞16Nat〈∑i=1Na[ri(t)-ri(0)]2〉,
where D is the diffusion coefficient of the molecules and Na is the number of the diffusing molecules. Furthermore, the diffusion coefficient components Dx, Dy, and Dz in x-, y-, and z-direction can be calculated by replacing r with x, y, and z in (2), 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 distance,r, 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 [29]:
(3)gαβ(r)=Vs(Nα·Nβ)〈∑i=1Nα(niβ(r)/4πr2Δr)〉,
where Nα and Nβ represent the numbers of α and β particles, respectively, Vs is the model volume, and niβ(r) represents the number of β particles in the range from r to (r+Δr) at the origin of ith α particle.


## 3. Results and Discussion

Figure2 displays the mean-squared displacement (MSD) of the hydrogen molecule of diffusion in NaX zeolite at 77 K for different pressures which are presented as a function of time. Figure 3 displays the same data at 10 kPa for different temperatures. 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. The self-diffusion coefficients were obtained from the slope of the curves via a linear regression fit to linear portion of MSD curves to (2). The calculated self-diffusion coefficient for hydrogen at 77 K for pressures of 10 and 100 kPa is 5.87 × 10−9 and 1.90 × 10−9 m2·s−1, respectively (Figure 2). These values are in good agreement with reported experimental diffusion coefficients obtained by Bar et al. [15] from the quasielastic neutron scattering (QENS) and pulse field gradients nuclear magnetic resonance (PFG NMR) methods. They are also in good agreement with those obtained from previous MD simulations [30]. It can be seen from Figure 2 that as hydrogen pressure increases there is an overall decrease in the diffusion coefficients for all temperatures. That is, the diffusivity decreases with increasing loading. Thus with each successive increase in pressure there is less available space for the fast moving hydrogen molecules; hence, the movement becomes hindered. At low pressure, the molecules have the highest mobility. NaX zeolite has a 12-membered ring pore opening of 7.4 Å in diameter that provides access to a supercage of 12.4 Å in diameter. Therefore, molecules have the freedom to move in all directions from one supercage to another cage (intercage). Due to their small size, hydrogen molecules have the good ability to move from one cage to another. At high pressure, we observe reduced translational mobility due to a reduced volume for the individual molecule. The diffusivity also decreases with increasing pressure. The sorbate-sorbate interaction further contributes to the decreased mobility with increasing pressure. Similar results are also found for alkane diffusion in mordenite zeolite [31].Figure 2
Mean square displacements of the center of mass of hydrogen diffusion at different pressures: (a) 10 kPa, (b) 100 kPa, (c) 600 kPa, (d) 1700 kPa, and (e) 2700 kPa at 77 K.Figure 3
Mean square displacements of the center of mass of hydrogen diffusion at different temperatures: (a) 77 K, (b) 160 K, (c) 195 K, (d) 230 K, and (e) 293 K at 10 kPa.In general as the simulation temperature increases there is an overall increase in diffusion coefficients for the hydrogen molecules. With comparison of the values of diffusion coefficients at a certain pressure and different temperatures as shown in Figure4, the diffusion coefficients increase with the increase in temperature. This shows that the diffusion of hydrogen in NaX zeolite is a thermally activated process. Figure 5 shows the relationship of diffusion coefficients changing with the inverse of the temperature at different pressure. It appears that the diffusion coefficients only for 10 and 100 kPa versus 1/T follow a linear relation. The activation energy for the diffusion of hydrogen in NaX is calculated by fitting diffusion coefficients for 10 and 100 kPa versus 1/T (Figure 5) to the Arrhenius equation [32]. The values obtained are 2.01 and 1.95 kJ/mol for 10 and 100 kPa, respectively. The activation energy of hydrogen in zeolite NaX was found experimentally to be 1–4.0 kJ/mol [15], and it was 2 kJ/mol for NaA zeolite [33]. Schuring et al. [31] considered that the activated energy is possibly dependent on pressure of gas molecules due to the increase of diffusion coefficient with the increase in temperature and the decrease in pressure in studying the diffusion of alkene in MFI zeolite. To solve the question, they proposed that the activated energy should be calculated at low pressure.Figure 4
The self-diffusion coefficients of hydrogen at different temperatures.Figure 5
Arrhenius plots of the self-diffusion coefficients of hydrogen in NaX zeolite. The corresponding equations and least squares fits are given on the graphs.In order to investigate the effects of pressure and temperature on the distance among hydrogen molecules, the radial distribution functions (RDFs) of the centers of mass for hydrogen molecules in zeolite pores at different pressures and temperatures are shown in Figures6 and 7, respectively. Only one peak is observed which is centered around 2.5 Å–3.3 Å in RDF plots of the hydrogen molecules, indicating the adsorbed hydrogen in zeolite pores is still unordered. For 77 K, at the pressures of 10 kPa and 100 kPa the distances are both about 2.9 Å which is approximate to the dynamics diameter of hydrogen molecule, 2.958 Å. The probability of the collusion among hydrogen molecules in pores is less than that of the collusion between hydrogen molecules and pore walls. The diffusion process of hydrogen molecules is controlled from the collusion between hydrogen molecules and pore walls. The distance decreases to about 2.7 Å and 2.5 Å for 600 kPa and above 1700 kPa, respectively, as shown in Figure 6. It is indicated that at high pressure hydrogen molecules are induced to change orientation or conformation to achieve a more compact packing. The collusion among hydrogen molecules in pores becomes frequent. It can be argued that this is also the reason for the very slow diffusion in the case of higher pressure. For the effect of temperature, in the case of 10 kPa, as shown in Figure 7, the distances among hydrogen molecules are about 2.9 Å at the temperatures of 77 K, 160 K, and 195 K. This increases to about 3.1 Å for 230 K and 3.3 Å for 293 K, respectively. 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 higher temperature.Figure 6
Radial distribution functions of the center of mass for hydrogen molecules at different pressures: (a) 10 kPa, (b) 100 kPa, (c) 600 kPa, (d) 1700 kPa, and (e) 2700 kPa at 77 K.Figure 7
Radial distribution functions of the center of mass for hydrogen molecules at different temperatures: (a) 77 K, (b) 160 K, (c) 195 K, (d) 230 K, and (e) 293 K at 10 kPa.To obtain a more detailed investigation of the information on the short-time dynamics of the system, the velocity autocorrelation function (VACF) of the diffusion of hydrogen in zeolite is analyzed. The velocity autocorrelation function is defined as(4)Cvv(t)=〈vi(t)·vi(0)〉〈vi2(0)〉,
where vi(t) and vi(0) denote the velocity of particle i at time t and time zero (i.e., t=0), respectively. Again, averaging is performed over all time origins in the trajectory.Figures8 and 9 show the velocity autocorrelation functions of the diffusion of hydrogen in NaX zeolite for different pressures and temperatures, respectively. In general, the structure of the velocity autocorrelation function can provide information on the short and intermediate time dynamics of the system. To quantify this short-time dynamical behavior of hydrogen molecules, we consider two quantities: (i) The time, τn, at which Cvv(t) first turns negative represents the average time at which the sorbate first encounters a repulsive barrier and is often closely related to the crossover time from ballistic to diffusion motion; (ii) τm, the position of the first minimum in Cvv(t) indicates the time scale on which a particle reverses its velocity as a result of repulsive encounters. As collisions with the wall and other sorbates increase, the correlation function decays to zero. The values of τn and τm in Figures 8 and 9 are listed in Tables 2 and 3. The magnitude of τn provides an estimate of the probable time, on average, that a molecule first encounters a repulsive barrier. From Table 2, we see that as pressure increases, τn continuously decreases. This may be attributed to the increasing frequency of intermolecular collision and collision between molecules and the walls of the pore. This trend in τn is reversed with increasing temperature (Table 3). The reason is possibly that as temperature increases, the dynamical energy of a molecule increases, which shortens the crossover time from ballistic to diffusion motion. Also seen from Tables 2 and 3, the values of τm are always larger than that of τn. However, τm values do not exhibit the good consistence with the variation of pressure and temperature.Table 2
The values ofτn and τm for the different pressures at 77 K.
10 kPa
100 kPa
600 kPa
1700 kPa
2700 kPa
τ
n (ps)
0.184
0.174
0.162
0.118
0.106
τ
m (ps)
0.220
0.358
0.300
0.403
0.261Table 3
The values ofτn and τm for the different temperature at 10 kPa.
77 K
160 K
195 K
230 K
293 K
τ
n (ps)
0.184
0.156
0.148
0.146
0.133
τ
m (ps)
0.220
0.180
0.198
0.180
0.200Figure 8
The velocity autocorrelation function of hydrogen molecules for the different pressure at 77 K.Figure 9
The velocity autocorrelation function of hydrogen molecules for the different temperature at 10 kPa.


## 4. Conclusions

The diffusion behavior of molecular hydrogen in NaX zeolite was investigated by using molecular-dynamics (MD) simulations. Our molecular-dynamics results show that the self-diffusion coefficient of hydrogen in NaX decreases as the pressure increases; this is in agreement with the earlier quasielastic neutron scattering (QENS) and pulse field gradients nuclear magnetic resonance (PFG NMR) intracrystalline self-diffusivity measurements of the other group [15]. The reason is found from RDF of hydrogen molecules that at high pressure hydrogen molecules are induced to change orientation or conformation to achieve a more compact packing. The collusion among hydrogen molecules in pores becomes frequent, showing the very slow diffusion in the case of the higher pressure. The temperature dependence of the self-diffusion coefficient shows an increase with increasing temperature between 77 K and 293 K. An explanation for the effect 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 higher temperature. The activated energy for hydrogen diffusion determined from the simulation is pressure-dependent. The value obtained below 100 kPa is about 2 kJ/mol. 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: 545367-2013-08-13.xml*