#### Abstract

The self-diffusion of hydrogen in Ca^{2+}-, Mg^{2+}- and Ba^{2+}-exchanged X zeolites (Mg_{46}X, Ca_{46}X, and Ba_{46}X) 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 m^{2}·s^{−1} to m^{2}·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 , 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).

#### 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 H_{2} 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 H_{2} 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 CH_{4} in NaX zeolites decreases with increasing loading from molecular dynamics simulation. Although several molecular dynamics simulations have been performed [11–13], simulation studies on H_{2} 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 H_{2} in Ca^{2+}, Mg^{2+}, and Ba^{2+} cations exchanged X-type zeolites (Mg_{46}X, Ca_{46}X, and Ba_{46}X, X = Si_{100}Al_{92}O_{384}). 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 (Ca^{2+}, Mg^{2+}, and Ba^{2+}) exchanged X zeolites such as Ca_{46}X, Mg_{46}X, and Ba_{46}X (X = Si_{100}Al_{92}O_{384}) 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 Ca_{46}X, Mg_{46}X, and Ba_{46}X 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.

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:

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.

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 and .

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 (Ca^{2+}, Mg^{2+}, and Ba^{2+}) 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 Ca_{46}X, Mg_{46}X, and Ba_{46}X, 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 where is the position of molecule at time and (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] where is the self-diffusion coefficient of the molecules. is the number of the diffusing molecule. Furthermore, the self-diffusion coefficient components , , and in -, -, and -direction can be calculated by replacing with , , and 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 distance from the reference atom. This concept also embraces the idea that the atom at the origin and the atom at distance may be of different chemical types, say *α* and *β*. The resulting function is then commonly given the symbol and is defined by Hansen and McDonald [26]
where and represent the numbers of *α* and *β* particles, respectively. is the model volume. represents the number of *β* particles in the range from to at the origin of th *α* particle.

#### 3. Results and Discussion

The diffusion behavior of hydrogen in Ca_{46}X, Mg_{46}X, and Ba_{46}X 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.

**(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 Ca_{46}X, Mg_{46}X, and Ba_{46}X 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.

**(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^{−8 }m^{2}/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} m^{2}/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 Figure 3). The self-diffusion coefficients increase with increase in temperature for a given loading. The activation energies for the diffusion of hydrogen in Ca_{46}X, Mg_{46}X, and Ba_{46}X 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 H_{2} 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, [28]: for , 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 , 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 . The diffusion mechanism is shown in Figure 4.

**(a)**

**(b)**

**(c)**

According to kinetic molecular theory of gas, the mean free path of gas molecules is defined as where is the Boltzmann constant, is temperature (K), is the effective diameter of the gas molecule (m), and is pressure (Pa).

With (5), the mean free path of H_{2} in Ca_{46}X, Mg_{46}X, and Ba_{46}X 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 on the loadingsfor H_{2} in Ca_{46}X, Mg_{46}X, and Ba_{46}X zeolites at 77 K as sample was shown in Figure 5. It was found that for these zeolites the values of is in range from 0.1 to 10 for full range of loadings (except the loading of 80 molecules/unit cell for Ba_{46}X 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 were maximum for Ba_{46}X zeolite and were minimum for Ca_{46}X 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 6 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 Ba_{46}X < Mg_{46}X < Ca_{46}X. 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 Ca^{2+} cations are located at site II in the supercage for Ca_{46}X; the 24 Mg^{2+} cations and 4 Mg^{2+} cations are located at site II and site II′, respectively; for Mg_{46}X and the 30 Ba^{2+} cations are located at site II in the supercage for Ba_{46}X, also as shown in Table 1. The total number of three cations at sites II and II′ is equivalent. Since the ionic radius of Mg^{2+} (0.66 Å) is smallest and the ionic radius of Ba^{2+} (1.35 Å) is largest, the windows of Ba_{46}X were seriously blocked. The diffusion of hydrogen in Ba_{46}X zeolite was the slowest of three zeolites. In X-ray crystallographic data of Mg_{46}X provided by Yeom et al. [14], Mg_{46}X zeolite was partially dehydrated. Mg_{46}(H_{2}O)_{4}-X was the final structure. The existence of partial water molecules in Mg_{46}X 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.

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 Figures 7, 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 Mg_{46}X zeolite was similar to Ca_{46}X zeolite. But there is no peak on RDF curve of hydrogen diffusion for Ba_{46}X. The reason may be that the pores in Ba_{46}X 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 Ba_{46}X 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.

#### 4. Conclusions

The diffusion behaviors of molecular hydrogen in Ca^{2+}-, Mg^{2+}-, and Ba^{2+}-exchanged X zeolites (Mg_{46}X, Ca_{46}X, and Ba_{46}X) 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 Ba_{46}X < Mg_{46}X < Ca_{46}X, 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.

#### Conflict of Interests

The author of this paper declares that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

The author gratefully acknowledges the financial support from the program for Liaoning Excellent Talents in University (LNET), China (LJQ2012016).