Force field for ZIF-8 flexible frameworks: atomistic simulation of adsorption, diffusion of pure gases as CH₄, H₂, CO₂ and N₂

Abstract

A full set of flexible force field parameters for ZIF-8 is presented, based on the AMBER, UFF parameters and the partial charges computed by the density-derived electrostatic and chemical charge method (DDEC). The parameters for the 2-methyl imidazole (MeIM) ring are adopted from the AMBER force field, while the van der Waals (VDW) parameters for organic linkers and metal centers were determined by rescaling the UFF parameters as ε = 0.635ε_UFF and σ = 1.0σ_UFF to fit the CH₄ adsorption isotherms obtained by Grand Canonical Monte Carlo (GCMC) simulations with the force field parameters to the experimental ones. The CH₄ adsorption isotherms on four different structures of ZIF-8 at 298 K obtained by GCMC simulations are compared with the experimental data. The results show that the simulated CH₄ adsorption isotherms on the ZIF-8 structure reported from the Cambridge Crystallographic Data Centre (CCDC) are closest to the ones on the ZIF-8 structure from the report of Moggach et al. To test our model, adsorption isotherms of CH₄, H₂, CO₂ and N₂ at different temperatures were computed using GCMC simulations, and the results were found to be in a good agreement with the experimental data. In the case of H₂, the equilibrium configurations obtained by GCMC simulations were statistically analyzed with ad hoc code to get probability density distribution profiles. These profiles were transformed to visual slice images, which indicate that the preferential adsorption sites of H₂ molecules in ZIF-8 are located close to the MeIM rings, where the host–guest VDW or electrostatic interactions are maximal, as revealed by the potential energy surfaces (PES). In addition, these force field parameters were confirmed to well reproduce the ZIF-8 structural properties including lattice constants, bond lengths and angles over a wide range of temperatures. The self-diffusivities at the specific loadings of adsorbed gases (CH₄, H₂ and CO₂) in ZIF-8 were calculated by the mean squared displacement (MSD) method. It was found that our self-diffusivities of H₂ are slightly higher than the ones in the literature, and our self-diffusivity of CO₂ is as about three times as the one in the literature, due to the different partial charges and the effect of different force field parameters on framework shape and flexibility in our simulations.

---

1. Introduction

In recent years, metal–organic frameworks (MOFs), a class of crystalline nanoporous materials composed of flexible organic linkers and bridging tetrahedral metal centers, were extensively investigated.^1–3 The possibility of tailoring these materials in terms of pore size or group functionality makes them extremely attractive for a wide variety of applications such as gas storage, gas purification and separation, and alkylation catalysis.^3–5 The so-called zeolitic imidazolate frameworks (ZIFs) are a particularly interesting subcategory of MOFs, based on molecular topologies that resemble zeolites, wherein transition metals (Co, Cu, Zn, etc.) replace Si-atoms and imidazolate-type linkers replace oxygen bridges in the conventional aluminosilicate structures.^6,7 Because of the strong bonding between the imidazolate linker and the metal center, many ZIFs have exceptional thermal and chemical stability.⁶ Moreover, it is possible to tune their pore size and affinity for adsorbate molecules by modifying organic linkers with different substituents or functional groups. The aim of this modification is to enhance adsorption capacity and selectivity of ZIFs for adsorbate molecules.⁸

Among large number of synthesized ZIFs to date, the prototypical ZIF-8 is the most extensively investigated.^9–11 In contrast to zeolites with relatively rigid frameworks, ZIF-8 shows an interesting structural flexibility as evidenced by several experimental and theoretical studies.^12–14 Bux et al.¹⁵ observed considerable permeability of methane in ZIF-8 using infrared (IR) microscopy despite methane possessing a larger kinetic diameter (3.8 Å) than the pore limiting diameter (PLD) of the rigid framework. The other gas molecules (N₂, C₂H₆, and C₃H₈) with kinetic diameters larger than the window size of ZIF-8 were also observed to freely diffuse through ZIF-8, due to the swing (reorientation) of imidazolate linkers that enlarges the window size (from 3.4 Å to 4.2 Å).^16–18 Zhang et al.¹⁹ recently estimated an effective aperture size for ZIF-8 between 4.0 and 4.2 Å based on the kinetic uptake measurements of short alkanes. Moggach et al.¹² reported experimentally the flexible structures of ZIF-8 over a range of pressures from 1 atm to 1.47 GPa. In addition, some simulation studies have been performed to compute gas adsorption isotherms and self-diffusion coefficients for ZIF-8 with both rigid and flexible frameworks.^20–23 Numerous force-field-based Grand Canonical Monte Carlo (GCMC) simulations of adsorption have been already carried out for MOFs including ZIF-8 with rigid framework,^20,22,24 but huge differences in the diffusion coefficients of guest molecules within rigid or flexible ZIF structures were observed, even for small guest molecules.^25,26 Hertäg et al.²¹ employed the DREIDING force field to investigate CH₄ diffusion in ZIF-8 and found the diffusivity of CH₄ in the flexible framework was smaller by a factor of ∼4 than that predicted for the rigid structure. At the same time, these authors also found that the crystal structure of ZIF-8 was not accurately predicted with this force field. As a result, the framework flexibility should be taken into consideration for an accurate study of the dynamics of guest molecules within ZIFs, and this is particularly true for predicting gas diffusivities by using molecular dynamics (MD) simulations.^23,27 The flexible framework model of Battisti et al.²⁸ matches closely the geometrical parameters of the ZIFs unit cells, but the lack of partial charges lowers the accuracy of the MD prediction for polar molecules.

Most recently, Zheng et al.²³ developed a fully flexible force field for ZIF-8, for which MD simulations accurately reproduce the experimental lattice constants and predict diffusivities for CO₂ that are in better agreement with the experiment. However, we found that CH₄ and H₂ adsorption isotherms for ZIF-8 obtained by using GCMC simulations with the nonbonded parameters of the flexible force field do not agree well with the experiment.²⁵ Zhang et al.²⁹ also developed a fully flexible force field for ZIF-8 based on the AMBER force field. They found that MD simulations with their force field parameters reproduced the structural transition of ZIF-8 at high loading of N₂ molecules as well as the ZIF-8 lattice constants. Also with the aforementioned force field, hybrid simulations combining MD with Gibbs Ensemble Monte Carlo (GEMC) reproduce accurately N₂ adsorption isotherms for ZIF-8, but the non-bonded interaction parameters between atom pair separated by two other atoms (1–4) is not explicitly mentioned in the paper and thus it is unable to determine the structure of ZIF-8 using MD simulation with the incomplete force field.

Despite of the above-mentioned progress in modeling the framework flexibility of ZIF-8, to the best of our knowledge, we are not aware of any simulation study that mimics the structural parameters of ZIF-8 and predicts gas adsorption and diffusion in ZIF-8 accurately at the same time. In this work, we developed a new force field that describes well the crystal structure of ZIF-8, as well as adsorption and diffusion of various gases in ZIF-8. Following this section, the next section outlines the force field development and simulation methodology. Next, the developed force field is used to predict the structural characteristics of ZIF-8, adsorption and diffusion behavior of some gases such as CO₂, CH₄ and H₂ in ZIF-8. Finally, the simulation results are compared with available experimental data.

2. Force field and computational details

2.1. Force field

With a cubic sodalite (SOD) topology, ZIF-8 consists of Zn metals tetrahedrically coordinated by four 2-methylimidazolate (mIM) linkers. ZIF-8 consists of large cavities (11.6 Å in diameter) interconnected by narrow six-ring windows. On the basis of experimental single-crystal X-ray diffraction (XRD) data, the window size in ZIF-8 is approximately 3.4 Å.⁶ The organic linker and metal center in ZIF-8 including atom types used in our force field parameterization are shown in Fig. 1a. Fig. 1b illustrates the crystal structure of ZIF-8 determined experimentally at ambient pressure taken from the Cambridge Crystallographic Data Centre (CCDC).³⁰ One flexible crystal structure of ZIF-8 obtained in Section 3.1 (ZIF8_55CH4UC) is shown in Fig. 1c, the other one taken from the ref. 12 is shown in Fig. 1d.

Fig. 1 (a) Atom types used in the flexible force field; (b) unit cell of the rigid framework taken from CCDC³⁰(ZIF-8); (c) unit cell of the flexible framework considered in this work (ZIF8_55CH4UC); (d) unit cell of the experimentally-established flexible framework by Moggach et al.¹² (ZIF8HL). The big yellow balls refer to the free volume of ZIF-8.

To develop the force field for ZIF-8, the parameters for bonded and nonbonded interactions are derived on the basis of the AMBER force field,³¹ the UFF force field,³² and experimental data. The “bonded” terms include bond stretching and bending, and proper and improper torsional potentials:

E_total = E_{bonded terms} + E_{nonbonded terms} (1)

E_{bonded terms} = E_bond + E_angle + E_dihedral + E_improper (2)

E_bond = K_b(b − b₀)² (3)

E_angle = K_θ(θ − θ₀)² (4)

E_dihedral = K_φ[1 + cos(nφ − φ₀)] (5)

E_improper = K_ψ[1 + cos(nψ − ψ₀)] (6)

where K_b, K_θ, K_φ and K_ψ are the force constants, b, θ, φ, and ψ are bond lengths and angles, proper and improper dihedrals, respectively, n is the multiplicity and was set to two for most dihedrals and three for N–Zn–N–C1 and N–Zn–N–C2, and b₀, θ₀, φ₀, and ψ₀ are the equilibrium values, adopted from the averaged bond lengths and angles based on experimental crystallographic data. The parameters for the organic linkers were adapted from the AMBER force field, which has been proved correct in describing the structure of 2-methyl imidazole (MeIM) ring.³¹ At the same time, the parameters for the interactions of tetrahedral ZnN₄ were adopted from the ref. 29. Table S1 and S2 (ESI†) list the optimized parameters for the bond stretching and bending and the torsional potentials.

The nonbonded interactions include Lennard-Jones (LJ) and Coulombic potentials:

E_{nonbonded terms} = E_LJ + E_COUL (7)

where σ_ij is the collision diameter and ε_ij is the energy well depth. The Lorentz–Berthelot combining rules are used for parameters of cross interactions.where q_i and q_j are the partial charges on the ith and jth atoms, respectively; ε₀ = 8.8542 × 10⁻¹² C² N⁻¹ m⁻², is the vacuum permittivity, and r is the distance between the ith and jth atoms.

In most GCMC simulation studies for MOFs, the LJ parameters were usually adopted from common force fields such as universal force field (UFF)³² and DREIDING³³ including the nonbonded interactions. However, these force fields led to overestimation of gas adsorption on ZIFs. In order to match experimental adsorption isotherms, Perez-Pellitero et al.²⁰ rescaled the UFF parameters as ε = 0.69ε_UFF and σ = 0.95σ_UFF to reproduce the CH₄ and CO₂ adsorption on ZIFs. In the case of CH₄, only van der Waals (VDW) interactions have to be taken into account.^22,25 In this work we adjusted the original UFF parameters rescaling as ε = 0.635ε_UFF and σ = 1.0σ_UFF in order to achieve better agreement between GCMC simulated adsorption isotherms and experimental data. The new VDW parameters of our force field, listed in Table S1,† were applied in the GCMC and MD simulations performed for the ZIF-8 framework. However, in the case of polar molecules, electrostatic interactions besides VDW interactions have to be taken into account in these simulations.²³

The partial charges of our model, listed in Table S1,† were computed by using the density-derived electrostatic and chemical charge method (DDEC).³⁴ Other charge models available in the literature, namely, the connectivity based atom contribution charge method (CBAC),³⁵ the repeating electrostatic potential extracted atomic charge method (REPEAT),³⁶ and ab initio calculation method³⁷ were also used in MD simulations to compute the lattice constants of ZIF-8 at 258 K. We found that the results predicted with DDEC method are in the best agreement with the experimental structure of ZIF-8 among the aforementioned charge models (see Fig. 2). Therefore, the remaining GCMC and MD simulations were performed using the atomic charges obtained by DDEC method.

Fig. 2 Lattice parameters of ZIF-8 crystal structure over a wide range of temperatures. The open circles refer to the simulation data from ref. 23. The filled squares refer to the experimental data from ref. 6 and the upper triangles refer to the experimental data from ref. 45.

H₂ and CH₄ were represented by a united-atom model using the LJ potential, while N₂ and CO₂ were mimicked as a three-site rigid model with parameters fitted to the experimental properties of bulk N₂ and CO₂, respectively.^25,38–40 The N–N bond length was 1.10 Å, and a charge of −0.482e (e = 1.6022 × 10⁻¹⁹ C is the elementary charge) was assigned on the N atom, as well as a charge of +0.964e at the center-of-mass (COM). The charges on C and O atoms were +0.6512e and −0.3256e, respectively. The C–O bond length was 1.149 Å, and the bond angle ∠OCO was 180°. Table S2 (ESI†) lists the LJ potential parameters and the atomic charges for N₂, CO₂, H₂ and CH₄. The cross interaction parameters between sorbates and ZIF-8 were obtained using the Lorentz–Berthelot combining rules.

2.2. GCMC simulation methods

The adsorption isotherms of various guest molecules on ZIF-8 were predicted using grand canonical Monte Carlo (GCMC) simulations, based on the guest–host and guest–guest force fields. All GCMC simulations were performed using the Multipurpose Simulation Code (MuSiC-4.0) developed by the Snurr group.⁴¹ Chemical potentials of gas sorbates at different temperature and pressure were converted to fugacity with the Peng–Robinson equation of state (PREOS).⁴² In many studies, the absolute adsorbed amounts obtained by GCMC simulation are converted to the excess adsorbed amounts to compare with the experimental data. In this paper, the excess adsorption N_ex is calculated using eqn (10):

N_ex = N_abs − ρ_bulkV_free (10)

where N_abs is the absolute amount, V_free represents the free pore volume of adsorbent, and ρ_bulk is the density of the sorbate calculated using PREOS at a given temperature and pressure.⁴³ The numbers of unit cells of ZIF-8 frameworks adopted in this simulation were 2 × 2 × 2, and periodic boundary conditions were applied in all three dimensions to eliminate boundary effects. The ZIF-8 frameworks were treated as rigid with frozen atoms during simulation. Cutoff radius of LJ interactions is 1.4 nm, and long-range electrostatic interactions are treated via the Ewald method. In this work, a total of 10⁷ steps were used; the first half of these moves was used for equilibration, and the remaining steps were used for calculating the ensemble averages. For pure CH₄ and H₂, three types of moves (translation, random insertion, and deletion) were used. For pure CO₂ and N₂, an additional move (rotation) was included. Every possible move was given equal probability.

2.3. Molecular dynamics simulation methods

The self-diffusion characteristics of various guest molecules in ZIF-8 were calculated using molecular dynamics (MD) simulations in the canonical ensemble (NVT) and microcanonical ensemble (NVE) with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS),⁴⁴ while lattice parameters of ZIF-8 crystal structure at different fixed temperature were calculated using MD in the isothermal–isobaric ensemble (NPT) with the same code. The numbers of unit cells of ZIF-8 frameworks adopted in this simulation were 2 × 2 × 2, and periodic boundary conditions were applied in all three dimensions. Cutoff radius of LJ interactions is 1.4 nm, and long-range electrostatic interactions are treated via the Ewald method. The simulation time step equals 1.0 fs. Note that, for each framework atom, a “scaled 1–4” policy is applied; i.e., both the VDW and electrostatic interactions between couples of bonded atoms (1–2) or between atoms bonded to a common atom (1–3) are excluded, while for the interaction between atoms separated by two other atoms (1–4) both the VDW “ε parameter” and electrostatic interaction for the given couple are divided by 2.

3. Results and discussion

3.1. Force field validation

The validation of our force field was initially focused on the unit-cell length of the simulated flexible ZIF-8 framework. MD simulations in NPT ensemble at different fixed temperature and pressure of 1.0 × 10⁵ Pa were performed using the LAMMPS package with the barostat fluctuations controlled via Nosé–Hoover method.⁴⁴ After equilibrating the system with running a 1 ns NPT simulation, the trajectories of all the framework atoms for another 1 ns were collected to compute the average unit-cell lengths of ZIF-8 crystal at the specific temperature, see Fig. 2. The results indicate that our force field parameters can better reproduce the experimental unit-cell lengths reported by Zhou et al.⁴⁵ and Park et al.⁶ than those parameters of Zheng et al.²³ The average bond length and the average bond angle values of ZIF-8 framework at 258 K and 10⁵ Pa were calculated, finding that these are in good agreement with the data obtained from X-ray diffraction experiments, see Table 1. The maximum percentile difference between experimental (16.991 Å) and simulated (16.982 Å) lattice parameter is about 0.05% (in our simulation, the fluctuation of the lattice parameter is about 0.03 Å, corresponding to 0.2% of its mean value).

Table 1 Comparison of the simulated bond lengths and angles with experimental data of ZIF-8 at 258 K and 10⁵ Pa

Average bond length (Å)
i–j	Simulation	Experiment
Zn–N	1.996	1.987
N–C1	1.343	1.340
N–C2	1.371	1.371
C1–C3	1.492	1.493
C2–C2	1.345	1.346
C2–H2	0.931	0.929
C3–H3	0.959	0.960

---

Average bond angle (degree)
i–j–k	Simulation	Experiment
N–Zn–N	109.34	109.47
Zn–N–C1	128.25	128.35
Zn–N–C2	125.30	126.40
C1–N–C2	105.11	105.24
N–C1–N	111.94	112.17
N–C1–C3	123.78	123.89
N–C2–C2	108.65	108.67
N–C2–H2	125.36	125.66
C2–C2–H2	125.59	125.67
C1–C3–H3	109.68	109.44
H3–C3–H3	109.11	109.50

---

The continuous structural transition of ZIF-8 upon gas sorption is mimicked using our force field as well as the force field of Zhang et al.²⁹ After loading 55 molecules per uc (unit cell) of CH₄ into the crystal, we run a 2 ns simulation in the NVT ensemble to equilibrate the system at 298 K. Orientational change of the imidazolate rings in the framework conformation (ZIF8_55CH4UC) is shown in Fig. 3 as compared to the original conformation of ZIF-8 framework. Apparently, the imidazolate rings in the high gas loading structure change orientation being perpendicular to the four-ring window, which is consistent with the experimental XRD data measured at high pressures.¹²

Fig. 3 (a) 2 × 2 × 2 supercell of the ZIF-8 framework reported from CCDC; (b) 2 × 2 × 2 supercell of the ZIF8-55CH4UC framework after running a 2 ns NVT ensemble MD simulation at 298 K.

3.2. Adsorption isotherms

Adsorption isotherms of four gases (CH₄, H₂, N₂ and CO₂) on ZIF-8 at 77 K and 298 K up to 100 bar were simulated and compared with the experimental ones (see Fig. 4–9). An excellent agreement was obtained between our simulation and experiment for adsorption isotherms of both polar gases and non-polar gases at 77 K and 298 K, respectively.

Fig. 4 Adsorption isotherms of CH₄ on ZIF-8 at 298 K. The isotherms in the left panel were calculated using GCMC simulations with the force field of Zhang et al.²⁹ The circles refer to the experimental data from ref. 46.

The excess CH₄ uptake isotherms obtained for four different initial structures (ZIF-8, ZIF8HL, ZIF8_ja3, ZIF8_55CH4UC, see Table S3 for details of their lattice constants and their six-ring windows see Fig. S1†) of ZIF-8 framework at 298 K generated by GCMC simulations using the force fields respectively developed by Zhang et al.²⁹ and this work, and the corresponding experimental adsorption isotherms reported by Zhou et al.⁴⁶ are compared in Fig. 4. Our simulation isotherms for ZIF-8 with the new force field are in a good agreement with the experimental data, but the simulations with the force field of Zhang et al.²⁹ underestimate the experimental values. At the same time, the effect of structure on the excess CH₄ uptake follows the order ZIF8HL > ZIF-8 > ZIF8_55CH4UC > ZIF8_ja3, and the simulation results of our structure (ZIF-8) are closest to that for CH₄ adsorption on ZIF8HL in the whole range of pressures. The crystal structure of ZIF8HL is considered as that with containing the structural transition at high gas loading and high pressure.¹² As a result, our GCMC simulation results for CH₄ adsorption on ZIF-8 might be trusted in spite of neglecting the framework flexibility and the important effect of the crystal structure on CH₄ adsorption. The CH₄ adsorption isotherms for ZIF-8 at 240 K were also computed by GCMC simulations and found to be in a moderate agreement with the experimental data reported by Zhou et al. (see Fig. S2†).⁴⁶

Fig. 5 and S3† show the excess N₂ uptake isotherms on ZIF-8 framework at 298 K and 77 K obtained by GCMC simulations using our force field, together with the experimental data reported by Perez-Pellitero et al.²⁰ and Park et al.,⁶ respectively. The agreement between the simulation isotherms and the experimental data is valid in the whole range of pressures at 298 K for the adsorption of N₂. There is also a fairly good agreement between the simulation results and the experimental adsorption isotherms of N₂ at 77 K at low pressures below 10⁻³ kPa, but the simulations performed by using the force field proposed in this work as well as those performed by Zhang et al.²⁹ underestimate the experimental adsorption of N₂ at 77 K and pressure ranging from 10⁻³ kPa to 100 kPa. This disagreement may be due to neglecting the flexibility of ZIF-8 framework in our GCMC simulations, causing that second step is not predicted by these simulations.

Fig. 5 Adsorption isotherms of N₂ on ZIF-8 at 298 K. The circles refer to the experimental data from ref. 20. The open circles refer to the simulated data from ref. 28.

The transferability of our force field has been tested by comparing the simulation results with the experimental adsorption of H₂ and CO₂ adsorption on ZIF-8 at 77 K and 298 K, respectively (see Fig. 6 and 7). In general, the simulation results obtained by our force fields are in a good agreement with the experimental data reported by Zhou et al.⁴⁶ and Park et al.⁶ for H₂ adsorption on ZIF-8 at 77 K and 298 K. Other researchers used the effective potential of Feynman and Hibbs to treat the quantum effects for H₂ at 77 K, but this aspect is neglected as the simpler potential in this work. At the same time, the simulations with the non-bonded parameters of the force field proposed in this work and the force field developed by Zhang et al.²⁹ both reproduce well the experimental adsorption of CO₂ adsorption in ZIF-8 at 298 K at pressures below 30 bar. But as the pressure increases to more than 30 bar, there is a maximum of CO₂ adsorption in the former simulated isotherms while there is not in the latter simulated ones. This difference may come from the different techniques employed to obtain the electrostatic charge distribution for the ZIF-8 framework. In general, there are higher charges on C3 atoms (−0.4526e) and lower charges on N atoms (−0.3879) in the former than those in the latter.

Fig. 6 Adsorption isotherms of H₂ on ZIF-8 at 77 K and 298 K. The upper triangles and diamonds refer to the experimental data from ref. 46 and the circles refer to the experimental data from ref. 6.

Fig. 7 Adsorption isotherms of CO₂ on ZIF-8 at 298 K. The open upper triangles refer to the data obtained by GCMC simulations with the force field of Zhang et al.²⁹ and the squares refer the experimental data from ref. 20. The open circles refer the simulated data from ref. 28.

3.3. Adsorption sites

To obtain the sorbate probability profile, the sorbate trajectory was obtained by rescaling all coordinates of sorbate molecules to the primitive unit cell via ad-hoc PBC.^23,25,47 Then, discretization of the simulation box in cubic bins (each of side 0.34 Å), we computed the probability (p) to find a sorbate molecule in the given bin, satisfying the condition where V is the volume of the primitive unit cell. The isoprobability contours are plotted using the Visual Molecular Dynamics (VMD) software.⁴⁸

H₂ probability density distribution profiles for ZIF-8 framework at 77 K at various pressures, showing the largest H₂ density (and thus the most probable adsorption sites) obtained from the Monte Carlo simulations are presented in Fig. 8. At 1 bar and 77 K, a little amount of H₂ molecules are adsorbed on ZIF-8 framework and their preferential adsorption sites are close to the imidazole rings. With increasing of gas pressure, more H₂ molecules are adsorbed in the middle of the faces of the accessible surface area of ZIF-8 framework. According to the aforementioned adsorption isotherms of H₂ on ZIF-8 at 77 K, the excess H₂ adsorption on ZIF-8 demonstrates little change in spite of pressure increasing from 30 bar to 80 bar. As shown in Fig. 8, H₂ probability density distribution profile at 30 bar is similar to that at 80 bar, indicating that the more H₂ molecules enter into the middle of the free cavity of ZIF-8 framework at high pressures.

Fig. 8 Probability density distribution (plane of XOY) for H₂ molecules adsorbed on ZIF-8 at 77 K at different pressures: (a) 1 bar (b) 30 bar and (c) 80 bar.

In addition, to understand better the main adsorption sites present in the ZIF-8 framework, potential energy surfaces (PES) representing the different energy contributions (VDW and electrostatics) exerted by the host framework on the adsorbed molecules were calculated. The potential-energy surfaces for VDW and electrostatic interactions were computed by numerical evaluation of the energy grids that were used in adsorption calculations (see Fig. 9). In the first case (VDW), the PES was obtained for the interaction of a hypothetical H₂ molecule with the solid structure by using the optimized parameters of the proposed force field. For the electrostatic surfaces we have also used the electrostatic grid computed for the adsorption process and evaluated with a hypothetical molecule with a positive charge of +1. This evaluation does not preserve the electro-neutrality of the system, but the objective of this plot is not to perform simulations of molecular adsorption but only to qualitatively reveal the surfaces of the solid which are most strongly charged. This kind of plot can help in the recognition of the regions over the surface of the solid where polar molecules could be preferentially adsorbed.

Fig. 9 Potential energy surfaces for H₂ molecules interacting with ZIF-8 at 77 K (a) VDW interaction between H₂ molecules and ZIF-8 framework. (b) Electrostatic interactions of a positive charge (+1) over the surface of ZIF-8 framework.

As shown in Fig. 9a, the main adsorption sites match the regions of maximal VDW interactions, which is consistent with the report of Perez-Pellitero et al.²⁰ It is also interesting to note that these prominent interaction sites are the same as those found for methane adsorption on ZIF-8 using neutron diffraction.⁴⁹ Since the electrostatic interactions between H₂ and the ZIF-8 framework is weak, the preferential adsorption sites for H₂ in ZIF-8 have nothing to do with electrostatic interactions, as shown in Fig. 9b.

3.4. Self-diffusivity

Before loading gas molecules into the crystal to perform the self-diffusivity computations, the empty relaxed structure of ZIF-8 was obtained through running 1 ns simulation in the NVT ensemble to equilibrate the system at 300 K and followed by a 10 ns simulation in the NVE ensemble according to the method of Zheng et al.²³ Subsequently, we loaded this structure with the specified number of guest molecules and performed a 1 ns NVT simulation to equilibrate the system. After that a 2 ns NVE simulation was carried out. The self-diffusivity (D_s) values of various guest gases within flexible ZIF-8 were calculated by using the mean squared displacement (MSD) method⁵⁰ according to the Einstein relation (eqn (11)):where r_k(t) is the position vector of guest molecule k at time t. In particular, the D_s value is averaged over all the N sorbed molecules and over multiple time origins t₀ (as symbolized by the brackets).

Fig. 10 shows the effect of gaseous sorbate loading on the self-diffusivity values for H₂ within flexible ZIF-8 at 300 K computed by MSD method. This figure reveals that H₂ self-diffusivity reaches a maximum (2.538 × 10⁻⁸ m² s⁻¹) at loading 20 H₂ molecules per unit cell (20 molecules per uc) of ZIF-8, while it is lower, 1.704 × 10⁻⁸ and 1.549 × 10⁻⁸ m² s⁻¹ at loading 3 and 55 molecules per uc, respectively. Pantatosaki et al.⁵¹ reported the simulated self-diffusivity values for H₂ within ZIF-8 at 298 K ranging from 1.1 × 10⁻⁸ to 1.8 × 10⁻⁸ m² s⁻¹ with the loading of H₂, which are slightly lower than our simulated results. This difference may be due to free orientation of the imidazole rings in the simulations with the proposed force field, resulting in the bigger size of four-ring windows in ZIF-8 framework. A comparison of self-diffusion coefficients for various gaseous sorbates (H₂, CH₄ and CO₂) in ZIF-8 at 300 K is shown in Table 2. The imidazolate rings in ZIF-8 can be classified into two types. Type I is located at the four-ring window perpendicular to the c axis, while type II is located at the six-ring window. Pantatosaki et al.⁵¹ reported that CH₄ is entrapped in the rigid ZIF-8 framework, whereas CO₂ with the less kinetic diameter is able to pass the four-ring windows presenting a tiny diffusivity value even in the rigid structure. But for similar flexible structures, their diffusivity values come up to three times higher than the ones obtained by the rigid frameworks. At the same time, the framework charge distribution also affects primarily the diffusivity behavior of sorbate molecules bearing charged centers like CO₂. According to Zheng et al.,²³ it is very important to add the proper partial charges, correct framework shape and flexibility during calculating self-diffusivity for polar molecules such as CO₂ in ZIF-8. In comparison with the simulated result in the charged-flexible ZIF-8 framework reported from Zheng et al.,²³ our simulated self-diffusion coefficient for CO₂ in ZIF-8 at 300 K is about two times higher, due to the different partial charges and the effect of different force field parameters on framework shape and flexibility in our simulations.

Fig. 10 Self-diffusion coefficients for H₂ in ZIF-8 at 300 K. The square refers to the simulated datum from ref. 28.

Table 2 Self-diffusivity of various gases in ZIF-8 at 300 K

Guest type	CO2	CH4	H2
a uc: unit cell.
Loading (molecules per uc^a)	6	10	10
D s (10^−10 m^2 s^−1)	5.463	0.399	228.0

---

4. Conclusions

In this work, a full set of flexible force field parameters for ZIF-8 have been presented, based on the AMBER, UFF parameters and the partial charges computed by DDEC method. The CH₄ adsorption isotherms for four different initial structures of ZIF-8 at 298 K have been computed using GCMC simulations with the proposed force field as well as those available in literature. The results show that the CH₄ adsorption isotherms on ZIF-8 and ZIF8HL obtained by using the proposed force field are in better agreement with the experimental data as compared to those obtained for the available force fields. To test our model, adsorption isotherms for four gases such as CH₄, H₂, CO₂ and N₂ at different temperatures were computed using GCMC simulation and were found to be in a good agreement with the experimental data. In the case of H₂, the probability density distribution profiles indicate that the preferential adsorption sites of H₂ molecules in ZIF-8 are located close to the MeIM rings, where the host–guest VDW interactions are maximal, as revealed by the potential energy surfaces (PES). In addition, it is shown that the force field parameters reproduce well the ZIF-8 structural properties including lattice constants, bond lengths and angles over a wide range of temperatures. The self-diffusivities at specified loadings of sorbed gases (CH₄, H₂ and CO₂) in ZIF-8 have been calculated. It was found that our self-diffusivities of H₂ are slightly higher than the ones in the literature, and our self-diffusivity of CO₂ is as about three times as the one in the literature, due to the different partial charges and the effect of different force field parameters on framework shape and flexibility in our simulations. In conclusion, the proposed force field is suitable to study gas adsorption as well as diffusion in the flexible structure of ZIF-8. The possibility of using this force field to other ZIF structures will be investigated in future.

Acknowledgements

This work was supported by the Natural Science Foundation of Hubei Province of China (2013CFB344), the National Natural Science Foundation of China (51272201), the program for New Century Excellent Talents in University of the Ministry of Education (NCET-13-0942) and the Fundamental Research Funds for the Central Universities, China (WUT: 2013-II-014).

References

1. M. Eddaoudi, J. Kim, N. Rosi, D. Vodak, J. Wachter, M. O'Keeffe and O. M. Yaghi, Science, 2002, 295, 469–472.
2. N. L. Rosi, J. Eckert, M. Eddaoudi, D. T. Vodak, J. Kim, M. O'Keeffe and O. M. Yaghi, Science, 2003, 300, 1127–1129.
3. J. L. C. Rowsell, E. C. Spencer, J. Eckert, J. A. K. Howard and O. M. Yaghi, Science, 2005, 309, 1350–1354.
4. T. Watanabe, S. Keskin, S. Nair and D. S. Sholl, Phys. Chem. Chem. Phys., 2009, 11, 11389–11394.
5. D. A. Gomez and G. Sastre, Phys. Chem. Chem. Phys., 2011, 13, 16558–16568.
6. K. S. Park, Z. Ni, A. P. Cote, J. Y. Choi, R. Huang, F. J. Uribe-Romo, H. K. Chae, M. O'Keeffe and O. M. Yaghi, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 10186–10191.
7. R. Banerjee, A. Phan, B. Wang, C. Knobler, H. Furukawa, M. O'Keeffe and O. M. Yaghi, Science, 2008, 319, 939–943.
8. B. Wang, A. P. Cote, H. Furukawa, M. O'Keeffe and O. M. Yaghi, Nature, 2008, 453, 207–211.
9. H. Amrouche, S. Aguado, J. Perez-Pellitero, C. Chizallet, F. Siperstein, D. Farrusseng, N. Bats and C. Nieto-Draghi, J. Phys. Chem. C, 2011, 115, 16425–16432.
10. Y. Pan, Y. Liu, G. Zeng, L. Zhao and Z. Lai, Chem. Commun., 2011, 47, 2071–2073.
11. A. Centrone, E. E. Santiso and T. A. Hatton, Small, 2011, 7, 2356–2364.
12. S. A. Moggach, T. D. Bennett and A. K. Cheetham, Angew. Chem., Int. Ed., 2009, 48, 7087–7089.
13. D. Fairen-Jimenez, S. A. Moggach, M. T. Wharmby, P. A. Wright, S. Parsons and T. Dueren, J. Am. Chem. Soc., 2011, 133, 8900–8902.
14. D. Fairen-Jimenez, R. Galvelis, A. Torrisi, A. D. Gellan, M. T. Wharmby, P. A. Wright, C. Mellot-Draznieks and T. Dueren, Dalton Trans., 2012, 41, 10752–10762.
15. H. Bux, A. Feldhoff, J. Cravillon, M. Wiebcke, Y.-S. Li and J. Caro, Chem. Mater., 2011, 23, 2262–2269.
16. H. Bux, F. Liang, Y. Li, J. Cravillon, M. Wiebcke and J. Caro, J. Am. Chem. Soc., 2009, 131, 16000–16001.
17. H. Bux, C. Chmelik, J. M. van Baten, R. Krishna and J. Caro, Adv. Mater., 2010, 22, 4741–4743.
18. H. Bux, C. Chmelik, R. Krishna and J. Caro, J. Membr. Sci., 2011, 369, 284–289.
19. C. Zhang, R. P. Lively, K. Zhang, J. R. Johnson, O. Karvan and W. J. Koros, J. Phys. Chem. Lett., 2012, 3, 2130–2134.
20. J. Perez-Pellitero, H. Amrouche, F. R. Siperstein, G. Pirngruber, C. Nieto-Draghi, G. Chaplais, A. Simon-Masseron, D. Bazer-Bachi, D. Peralta and N. Bats, Chem. - Eur. J., 2010, 16, 1560–1571.
21. L. Hertäg, H. Bux, J. Caro, C. Chmelik, T. Remsungnen, M. Knauth and S. Fritzsche, J. Membr. Sci., 2011, 377, 36–41.
22. H.-C. Guo, F. Shi, Z.-F. Ma and X.-Q. Liu, J. Phys. Chem. C, 2010, 114, 12158–12165.
23. B. Zheng, M. Sant, P. Demontis and G. B. Suffritti, J. Phys. Chem. C, 2012, 116, 933–938.
24. H. Frost, T. Duren and R. Q. Snurr, J. Phys. Chem. B, 2006, 110, 9565–9570.
25. X. Wu, X. Yang, J. Song and W. Cai, Acta Chim. Sin., 2012, 70, 2518–2524.
26. T. Chokbunpiam, R. Chanajaree, O. Saengsawang, S. Reimann, C. Chmelik, S. Fritzsche, J. Caro, T. Remsungnen and S. Hannongbua, Microporous Mesoporous Mater., 2013, 174, 126–134.
27. P. Puphasuk and T. Remsungnen, J. Comput. Theor. Nanosci., 2013, 10, 227–231.
28. A. Battisti, S. Taioli and G. Garberoglio, Microporous Mesoporous Mater., 2011, 143, 46.
29. L. Zhang, Z. Hu and J. Jiang, J. Am. Chem. Soc., 2013, 135, 3722–3728.
30. The Cambridge Crystallographic Data Centre, http://www.ccdc.cam.ac.uk Accessed March, 2013.
31. Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. M. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. M. Wang and P. Kollman, J. Comput. Chem., 2003, 24, 1999–2012.
32. A. K. Rappe, C. J. Casewit, K. S. Colwell, W. A. Goddard and W. M. Skiff, J. Am. Chem. Soc., 1992, 114, 10024–10035.
33. S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897–8909.
34. T. A. Manz and D. S. Sholl, J. Chem. Theory Comput., 2010, 6, 2455–2468.
35. Q. Xu and C. Zhong, J. Phys. Chem. C, 2010, 114, 5035–5042.
36. C. Campana, B. Mussard and T. K. Woo, J. Chem. Theory Comput., 2009, 5, 2866–2878.
37. M. K. Rana, F. G. Pazzona, G. B. Suffritti, P. Demontis and M. Masia, J. Chem. Theory Comput., 2011, 7, 1575–1582.
38. C. S. Murthy, K. Singer, M. L. Klein and I. R. McDonald, Mol. Phys., 1980, 41, 1387–1399.
39. K. M. Gupta, Y. Chen, Z. Hu and J. Jiang, Phys. Chem. Chem. Phys., 2012, 14, 5785–5794.
40. W. Zhang, H. Huang, C. Zhong and D. Liu, Phys. Chem. Chem. Phys., 2012, 14, 2317–2325.
41. A. Gupta, S. Chempath, M. J. Sanborn, L. A. Clark and R. Q. Snurr, Mol. Simul., 2003, 29, 29–46.
42. D. Y. Peng and D. B. Robinson, Ind. Eng. Chem. Fundam., 1976, 15, 59.
43. T. Duren, Y.-S. Bae and R. Q. Snurr, Chem. Soc. Rev., 2009, 38, 1237–1247.
44. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
45. W. Zhou, H. Wu, T. J. Udovic, J. J. Rush and T. Yildirim, J. Phys. Chem. A, 2008, 112, 12602–12606.
46. W. Zhou, H. Wu, M. R. Hartman and T. Yildirim, J. Phys. Chem. C, 2007, 111, 16131–16137.
47. X. Wu, J. Zheng, J. Li and W. Cai, Acta Phys.-Chim. Sin., 2013, 29, 2207–2214.
48. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33.
49. H. Wu, W. Zhou and T. Yildirim, J. Phys. Chem. C, 2009, 113, 3029–3035.
50. B. Liu, Q. Yang, C. Xue, C. Zhong and B. Smit, Phys. Chem. Chem. Phys., 2008, 10, 3244–3249.
51. E. Pantatosaki, F. G. Pazzona, G. Megariotis and G. K. Papadopoulos, J. Phys. Chem. B, 2010, 114, 2493–2503.

---

Footnotes

† Electronic supplementary information (ESI) available: Force field parameters for guest molecules and host framework of the ZIF-8 used in this work and the lattice constants of four different initial structures of ZIF-8 are listed in Table S1–S3. The N₂ adsorption isotherms at 77 K and the CH₄ adsorption isotherms at 240 K on ZIF-8 are shown in Fig. S2–S3. See DOI: 10.1039/c4ra00664j

‡ The authors declare no competing financial interest.

---