Flexibility and swing effect on the adsorption of energy-related gases on ZIF-8: combined experimental and simulation study†

David Fairen-Jimenez‡ *^a, Raimondas Galvelis ^b, Antonio Torrisi ^b, Alistair D. Gellan ^a, Michael T. Wharmby ^c, Paul A. Wright ^c, Caroline Mellot-Draznieks ^b and Tina Düren ^a
aInstitute for Materials and Processes, School of Engineering, The University of Edinburgh, United Kingdom. E-mail: David.Fairen-Jimenez@Northwestern.edu
^bDepartment of Chemistry, University College London, United Kingdom
^cEaStCHEM School of Chemistry, University of St. Andrews, United Kingdom

First published on 27th June 2012

---

---

1. Introduction

Zeolitic imidazolate frameworks (ZIFs) are a very promising family of microporous metal–organic frameworks (MOFs) where the self-assembly of tetrahedral atoms (e.g. Zn and Co) and organic imidazolate linkers results in a large diversity of structures with similar topologies and energy landscapes to those found in zeolites.^1–3 In many cases, ZIFs exhibit large pore volumes and surface areas as well as good thermal stability (up to 550 °C). Furthermore, the possibility of tuning the host–guest interactions by choosing and functionalizing the organic building blocks,⁴ together with pore networks consisting of large cavities interconnected by small windows that allow for shape- and size-selective separations, make ZIFs very promising porous materials for gas separation and storage.

The adsorption of energy-related gases is of great importance in many industrial and automotive applications. This includes, for example, the storage of natural gas, which mainly consists of CH₄, as an attractive alternative to petroleum-based fuels. Moreover, alkane purification and CO₂ separation using porous materials are key processes in the petrochemical industry and in carbon capture and sequestration (CCS) technologies.

Grand canonical Monte Carlo (GCMC) simulation is a powerful technique that allows us to quantitatively predict the uptake of fluids and to screen existing or hypothetical MOF structures for specific applications. For the proposed applications, the prerequisite to precisely predict the adsorption behavior of a fluid in a porous material is the accurate description of the interaction between the fluid and the solid. Adsorption in ZIFs has been widely studied by molecular simulation with a particular focus on ZIFs with the gmelinite (GME) topology, such as ZIF-68 and ZIF-69,^5–7 but also in many other structures.^8–11 In most of these studies, where generic force fields such as UFF¹² and Dreiding¹³ were used, the simulated amount adsorbed was significantly overpredicted without significant differences between force fields.^14,15 However, Pérez-Pellitero et al. reported a notable improvement when the force field parameters (i.e. σ and ε of the Lennard-Jones potential) were adjusted to reproduce the experimental isotherms.¹⁶ On the other hand, a recent study by Babarao et al. revealed that when pore accessibility for ZIF-68 and ZIF-69 was correctly taken into account by partial pore blocking, simulations with generic force fields, such as UFF without adjusting the potential parameters, were able to match the experimental results.¹⁵

In this work, we focus on the study of ZIF-8 as a prototypical structure with sodalite (SOD) topology. ZIF-8 is comprised of large cavities (∼11.7 Å of diameter) interconnected by narrow windows (∼3.40 Å). Therefore, in the absence of structural defects and in contrast to ZIF-68 and ZIF-69, the whole porosity of the framework will be either accessible to gas molecules or not, depending on the relation between the kinetic diameter of the gas molecule and the 3.40 Å window size of ZIF-8. However, it has been experimentally observed that molecules which are theoretically too large to pass through the windows of the ZIF framework, such as CH₄ and N₂, are adsorbed into the material.^1,17 In a previous paper, we observed using in situ XRD the existence of a structural change in ZIF-8 during the adsorption of N₂ at 77 K, caused by a swing of the imidazolate linkers.¹⁸ In this case, the step in the N₂ adsorption isotherm for ZIF-8 was caused by a structural change, which was also responsible for the unexpected adsorption of large molecules in ZIF-8 due to the enlargement of the window size. Thus, a correct description of the flexibility using two different structural configurations (i.e. ZIF-8AP obtained for the empty structure at ambient pressure and ZIF-8HP obtained at high N₂ uptake), is essential to precisely describe the adsorption performance using GCMC simulations. More recently, Ania et al. combined gas adsorption experiments and GCMC simulations to study the structural change of ZIF-8 during the adsorption of N₂, O₂, Ar and CO at cryogenic temperatures, where their isotherms showed pronounced steps during the adsorption processes.¹⁹ However, it is unclear if this phenomenon also arises for other gases where no steps exist in the adsorption isotherms and if it is also extended to ambient, and therefore industrially applicable, temperatures. In the present paper, we assess the importance of the framework flexibility and the use of different generic force fields in ZIF-8 for the prediction of the energy-related gases methane, ethane, propane, butane and CO₂ at cryogenic and ambient temperatures. First, we adjust the generic force fields commonly used for the molecular simulation of adsorption in MOFs to correctly describe the solid–fluid interaction for CH₄, taking into account the influence of the flexibility on ZIF-8 during the adsorption process. Then, we study the transferability of the proposed force fields by comparing the simulated isotherms of different alkanes and CO₂ to the experimental isotherms of ZIF-8 and methane in ZIF-20.²⁰ Finally, we complement our results with dispersion-corrected DFT-D calculations to provide further insights into the different adsorption sites in the two different ZIF-8 structural configurations.

2. Experimental

2.1. Synthesis and structure of ZIF-8

ZIF-8 was prepared following a modified version of the method described by Huang et al.¹ A solution of 2-methylimidazole (6.442 g, 80.5 mmol) in methanol (100 mL) was added via a dropwise addition to a stirred solution of Zn(OH)₂ (3.904 g, 39.0 mmol) in aqueous ammonia (25%, 500 mL). Within 20 minutes a white precipitate was observed. The reaction mixture was stirred at room temperature for 48 hours to ensure the reaction was complete. The fine white crystalline powder was collected, washed (methanol–water 1:1) and dried in air (yield 75%).

ZIF-8 is formed via the self-assembly of Zn and 2-methylimidazolate linkers (MeIM), giving rise to four- and six-ring ZnN₄ clusters. In ZIF-8, the small windows (4-ring) are closed whereas the larger ones (6-ring), where the three –CH₃ groups from the MeIM are pointing, remain open giving access to the pore cavities. Fig. 1 represents the two structural configurations, ZIF-8AP and ZIF-8HP, found for this material during N₂ adsorption.¹⁸ The geometrical pore size distribution (PSD) of the structure was calculated with the method developed by Gelb and Gubbins which determines the diameter of the largest sphere that can fit into the cavities without overlapping with any of the framework atoms.²¹ Fig. S1† in the ESI shows a single peak with a pore diameter of 10.7 Å related to the single cavities that shape the structure.

2.2. Gas adsorption

CH₄ adsorption isotherms for ZIF-8 at different temperatures from 125 to 300 K were taken from the work by Yildirim and co-workers.¹⁷ CO₂ adsorption isotherms for ZIF-8 at 196 K (ethanol–dry ice mixture) (Grant GR150 thermostatic refrigerated bath) were obtained using a Hiden IGA automatic gravimetric porosimeter. Alkane (ethane, propane and butane) adsorption isotherms for ZIF-8 at 273 K and CO₂ at 300 K were determined using a bench-scale open-flow adsorption–desorption apparatus previously described in detail elsewhere.²² Prior to the measurements, the sample was degassed at 425 K using a heating rate of 5 K min⁻¹ for 4 h. The CH₄ isotherm for ZIF-20 at 273 K was obtained from the work by Yaghi and co-workers.²⁰

2.3. Grand canonical Monte Carlo simulations

Gas adsorption was simulated using grand canonical Monte Carlo (GCMC) simulations²³ implemented in the multipurpose simulation code Music.²⁴ We used an atomistic model for the two different ZIF-8 structural configurations previously reported, ZIF-8AP and ZIF-8HP, in which the framework atoms were kept fixed at their crystallographic positions.¹⁸ At each pressure, 5 × 10⁷ Monte Carlo steps were performed, the first 40% of which were used for equilibration and the remaining steps were used to calculate the ensemble averages. We used the standard Lennard-Jones (LJ) 12–6 potential to model the interatomic interactions between the framework and gases. The parameters for the framework atoms were obtained from the UFF force field¹² as well as the adjusted force fields, UFF(*) and UFF(+), developed in this work (Table 1). Mülliken partial charges for the framework atoms were taken from periodic DFT calculations. The fitting of an electronic distribution to obtain point charges, in contrast to using an electrostatic potential energy surface, has been shown to give satisfactory results for ZIF-8, a highly hydrophobic structure.²⁵ For the alkane molecules, the TraPPE force field was used to describe the interactions.^26,27 This force field uses a united-atom approach for these molecules and considers each carbon atom and its surrounding hydrogen atoms as one LJ sphere with one set of parameters. CO₂ was modelled using the TraPPE potential with charges placed on each atom derived to accurately describe the CO₂ vapor–liquid phase equilibrium.²⁸ The model approximately reproduces the first-order electrostatic and second-order induction interactions and has been successfully used to simulate CO₂ adsorption in MOFs.²⁹ Coulombic interactions were included for CO₂–CO₂ and CO₂–ZIF-8 using the Ewald summation method.³⁰ Interactions beyond 18 Å were neglected. To calculate the gas-phase fugacity we used the Peng–Robinson equation of state.³¹

---

After equilibration, density distributions were obtained by storing the center of mass positions of all the adsorbed molecules at regular intervals during the simulation. These density distributions provide valuable information about the preferential adsorption sites and the local spatial disorder of the adsorbed molecules.

In order to compare the simulation results with the experimental data, we calculated the Henry's constant for all the adsorption isotherms studied here. Furthermore, we determined the isosteric heat of adsorption, Q_st, for methane and CO₂ using the Clausius–Clapeyron equation (1), which is directly related to the potential energy and the interaction energy:

	(1)

2.4. DFT-D binding energy calculations

Density Functional Theory (DFT) calculations have made significant contributions in the realm of MOFs. It is becoming apparent that the inclusion of dispersion interaction corrections is crucial in modelling this class of materials. For example, to capture the bistability of the archetypal MIL-53 framework,³² the mechanical properties of the dense zeolitic imidazolate zni,³³ and to compare the energy landscape of Zn and LiB-based ZIFs.³⁴

The DFT calculations were performed using the CP2K/Quickstep code.³⁵ The restricted Kohn–Sham formalism with PBE exchange–correlation functional³⁶ was used with the semi-empirical dispersion correction by Grimme.³⁷ The electronic energy was minimised with the orbital transformation (OT) method.³⁸ The convergence criterion for the self-consistent field (SCF) procedure was set to 1.0 × 10⁻⁷. The nuclear and core electronic densities were modelled with Goedecker–Teter–Hutter (GTH) pseudo-potentials^39,40 and the valence electronic density represented by the hybrid Gaussian and plane-wave (GPW)⁴¹ basis sets scheme. All atoms had the split-valence double-ζ (MOLOPT-DZVP)⁴² basis sets. The plane wave cut-off was set to 400 Ry. Periodic boundary conditions and Γ-point only sampling was used. Structure optimizations were done in the P1 space group.

In order to understand the thermodynamics of the adsorption of CH₄ in ZIF-8, DFT calculations of the binding energy have been performed. The binding energy (BE) of one CH₄ molecule per unit-cell of ZIF-8 was evaluated as given in eqn (2):

BE = Eoptim(CH4 + ZIF-8) − (Eoptim(ZIF-8) + Eoptim(CH4))	(2)

where E_optim represents the energy of the system after optimization.

ZIF-8HP is unstable without a sufficient number of guest CH₄ molecules inside, so the binding energy cannot be calculated using eqn (2). In order to characterize CH₄ in ZIF-8HP, a constrained optimization has been performed by fixing the carbon atom of CH₄ in the middle of the 4-ring window. This provides an estimate of the penalty energy for the insertion of one CH₄ molecule into the 4-ring window upon the swing effect of the imidazolate linkers. The penalty energy (PE) was evaluated as given in eqn (3):

PE = Econstr(CH4 + ZIF-8) − (Eoptim(ZIF-8) + Eoptim(CH4)),	(3)

where E_constr represents the energy of the system after constrained optimization.

3. Results and discussion

3.1. Methane adsorption and force field modification

Previous analysis of the N₂ adsorption isotherm for ZIF-8 at 77 K showed the absence of a significant amount of non-porous defects or solvent molecules inside the pores, indicating that the sample was successfully activated.¹⁸ Starting from this point, we first considered the use of the generic UFF force field to simulate the adsorption of CH₄ in ZIF-8 without any modification or adjustment. Fig. 2a shows the comparison between the UFF-simulated and experimental adsorption isotherms obtained for CH₄ at two different temperatures, 125 and 300 K. By comparing the experimental and simulated isotherms, we find that, whereas at 300 K the amount adsorbed is over-predicted with respect to the experimental results, the maximum capacity is accurately predicted at 125 K. We previously considered the existence of defects in the experimental sample and the presence of solvent molecules inside the pores as the cause of the differences between the simulated and experimental uptake.¹⁴ These two possibilities would cause a reduction in the effective pore volume available for adsorption and therefore the maximum amount adsorbed. As a result, the uptake in the simulated curves can be generally corrected by a scaling factor, Φ, to fit the experimental curve. However, the comparison of the experiments and simulations at six different temperatures between 150 and 300 K (Fig. 2a and Fig. S3†) shows that the situation is more complicated and that this explanation is too simple. Firstly, the scaled isotherms do not match the experimental isotherms in the low to medium pressure range and secondly, the scaling factors change between 1.01 and 0.74 with the temperature of the isotherm used (Table S2†).

The match between the maximum adsorption capacity in the experimental and UFF-simulated adsorption isotherms at 125 K suggests that the accessible pore volume is correctly described by the molecular simulation. However, for all temperatures the simulated Henry's constants are much higher than the experimental values, indicating that the solid–fluid interactions are over-predicted (Table 2). Moreover, the Q_st values, directly related to the interaction energy, are overestimated by around 20% (Fig. 3a). To take this into account, we considered two alternatives: first, we used an adjusted UFF potential, named here as UFF(*), where σ and ε were adjusted in order to achieve a good fit over the whole pressure range of the isotherm, following a similar approach used, for example, by Pérez-Pellitero et al.¹⁰ Second, we used a modified UFF potential named UFF(+), where only ε was adjusted in order to achieve a good fit in the low pressure range of the isotherm (i.e. to correctly describe the Henry's law region, which is a direct measure of the strength of the solid–fluid interactions) rather than the adsorption capacity. Moreover, for the second case we conducted the simulations in two different structural configurations (i.e. ZIF-8AP and ZIF-8HP) to study the role of the framework flexibility (i.e. the swing of the imidazolate ring) previously observed for the adsorption of N₂, following a similar approach that we used in the study of the adsorption processes in flexible ZIF and MOF structures.^18,43,44

---

Using UFF(*), for which both Lennard-Jones parameters, σ and ε, were adjusted, leads to an excellent agreement with the experimental adsorption isotherms for all temperatures studied (Fig. 2b and S4†) and a significant reduction of the Henry's constants and isosteric heat of adsorption (Table 2 and Fig. 3a). On the other hand, when using UFF(+), for which only the Lennard-Jones parameter ε was adjusted, the ZIF-8AP configuration results were in excellent agreement with the experimental isotherms at loadings below 30 molecules per unit cell (Fig. 2c and S5†) as well as the Henry's constants and isosteric heat of adsorption (Table 2 and Fig. 3a). This indicates that UFF(+) provides a good description of the fluid–solid interactions. However, for higher loadings (i.e. above 30 molec uc⁻¹), the simulations start diverging from the experimental behavior. This is particularly pronounced for 125 K, where the simulated isotherm reaches a plateau around 33 molec uc⁻¹ whereas the experimental isotherm increases up to 37 molec uc⁻¹. Previously, we have shown for the N₂ adsorption in ZIF-8 that the structure switches from the ZIF-8AP to the ZIF-8HP configuration at higher loadings.¹⁸ Similarly, excellent agreement between the simulated and experimental CH₄ isotherms at higher loadings is achieved when using the ZIF-8HP configuration. The same is observed for all the temperatures studied although at higher temperatures there is very little difference in the uptake in the two configurations. This observation is a strong indication that the structural change is not only induced by the N₂ molecules but also by other fluids.

The need to reduce the Lennard-Jones parameter, ε, to correctly describe the experimental isotherms for ZIFs is not unique to UFF but also applies in the use of the Dreiding force field¹³ and a combination of the Dreiding and OPLS force fields proposed by Zhou et al.^8,45 (Fig. S6, ESI†). These force fields also show significant deviations to the experimental values when used at lower temperatures, where the adsorption isotherm is more sensitive to the force field parameters.

3.2. Force field transferability

We first studied the performance and transferability of the UFF, the UFF(*) and UFF(+)-modified force fields together with the effect of the structural flexibility during the adsorption of CO₂. Using the non-modified UFF force field to predict the CO₂ adsorption isotherms resulted in an over-prediction of the Henry's constants for both 196 K and 273 K (Table 2) and the Q_st values (Fig. 3b). This indicates that the CO₂⋯ZIF-8 interactions are over-predicted. Additionally, the maximum amount adsorbed is over-predicted at 273 K (Fig. 4a). Whereas the adapted force field UFF(*), for which both σ and ε were modified, presented the best performance for the prediction of the CH₄ adsorption isotherms at different temperatures, we found important discrepancies between the experimental and simulated CO₂ isotherms at both 196 and 273 K (Fig. 4b). For CO₂, the simulated isotherms under-predict the experimental values over the whole pressure range, which is more pronounced at low temperature where the isotherm is more sensitive to the CO₂⋯ZIF-8 interactions, resulting in a poor transferability of the UFF(*) force field. On the other hand, the use of UFF(+), where only ε was modified, revealed an excellent agreement in the calculated Henry's constants when using the ZIF-8AP configuration, which means that UFF(+) correctly describes the solid–fluid interactions not only for CH₄ but also for CO₂ (Fig. 4c). Small deviations in the amount adsorbed are observed when close to the saturation pressures, where they are slightly under-predicted. As in the case of N₂ and CH₄ adsorption, the use of UFF(+) on ZIF-8HP solves this deficiency, nicely predicting the maximum capacity of the experimental isotherm. In contrast to N₂ and CH₄, the shape of the CO₂ adsorption isotherm is not very sensitive to structural changes in ZIF-8.

The transferability of the force fields was also studied for the adsorption of ethane, propane and butane. Fig. 5 shows the comparison between the experimental and simulated isotherms, using UFF(*) and UFF(+) for these gases at 273 K. For UFF, the Henry's constant and amount adsorbed is over-predicted for the whole adsorption isotherm, as expected from the results for CH₄ and CO₂ (Table 2 and Fig. S7†). The differences between the UFF(*)-simulated isotherms and the experimental ones increase in the order ethane, propane and butane (i.e. they are more significant for the longer alkanes which differ more from CH₄, the fluid used to adapt the force field). In contrast, the UFF(+) simulations on ZIF-8AP fit both the experimental isotherm in the low pressure range and the Henry's constant very well. However, the prediction of the amount adsorbed deviates from the experimental data when the pressure (and therefore the loading) increases around 60, 5 and 1.5 kPa for ethane, propane and butane, respectively. When running the simulations using UFF(+) on the ZIF-8HP configuration, the simulated isotherms do not fit the experimental ones at low pressure but the maximum capacity and the higher pressure range starting at 100, 20 and 5 kPa for ethane, propane and butane, respectively, are predicted accurately. The good description of the experimental behavior for every fluid when using a combination of simulated isotherms on ZIF-8AP and ZIF-8HP indicates that the material is experiencing the same structural change as was observed during the adsorption of N₂ at 77 K.¹⁸

This structural change is in agreement with similar changes reported during the adsorption of alkanes in ZIF-7, where the gate-opening pressure evolves in the order propane < ethane < methane.^46–48 The driving force that causes the change upon gas adsorption is most likely the additional interaction the structure receives when adsorbing extra molecules. Temperature, fluid–solid and fluid–fluid interactions also play an important role in the structural change and in the shape of the adsorption isotherm.^49–51 Analogous structural and energetic transitions during gas adsorption have been analyzed by Boutin and co-workers^52,53 and others.⁵⁴

In a further step, we studied the transferability of the force field to methane adsorption in ZIF-20.²⁰ ZIF-20, with an LTA topology, presents three different cavities: α-cages, β-cages and small cube pockets with pore apertures of 2.8, 2.0 and 1.5 Å, respectively. As a result, methane molecules could be adsorbed only in the α-cages and the other cavities need to be blocked during the GCMC simulation. Fig. 6 shows the results of the comparison of the experimental isotherms with the simulated isotherms, using UFF and modified UFF(*) and UFF(+) force fields for methane at 273 K. Firstly, using UFF resulted in a significant over-prediction of the Henry's constant and experimental isotherm, whereas UFF(+) was able to accurately predict the experimental results, giving better results than UFF(*) (Table 2).

3.3. Methane adsorption mechanism

Having validated that UFF(+) provides a good description of the adsorption isotherms of different fluids at different temperatures, we validated our force field in more detail by studying the adsorption mechanism of CH₄ and the influence of the swing effect of the imidazolate ring. Thus, we compared the different adsorption sites obtained in the simulation at 125 K (i.e. the lowest temperature where an experimental isotherm of CH₄ in ZIF-8 was measured)¹⁷ with the neutron diffraction experiments from Yildirim and co-workers,⁵⁵ following a similar approach that we previously used in the study of CH₄ and H₂ adsorption.^56,57 We compared the adsorption sites by using density distributions of CH₄ adsorption (125 K) at three different uptakes in ZIF-8AP and ZIF-8HP: 15, 36 and 40 molec uc⁻¹ (Fig. 7). At low loading (i.e. 15 molec uc⁻¹), two positions (I and II, previously named IM and Channel Site, respectively in the neutron diffraction experiments) are occupied by the CH₄ molecules, correctly describing the experimental results.⁵⁵ Site I is localized on the top of the CC bond of the organic linkers whereas Site II is localized at the center of the 6-ring channel, where the three –CH₃ groups from the MeIM are pointing. These sites form a hypothetical cube inside the main cavity, previously referred to as a “nanocage”.⁵⁵ Here, four of the eight corners of the nanocage are defined as Site II (eight positions per unit cell) whereas the other four corners remain empty during the whole adsorption process. Site I is localized at the center of the nanocage edge (twelve positions, twenty four per unit cell). At 36 molec uc⁻¹, the first two adsorption sites continue to be filled but the adsorption of extra molecules and the existence of packing effects and CH₄–CH₄ interactions confine the molecules into much more localized positions. Besides this, up to four extra molecules are localized in the center of the nanocage, here called Site III, giving a total amount of 36 positions (i.e. molecules adsorbed) per unit cell in ZIF-8AP, saturating the unit cell and again matching the neutron diffraction experimental results. Neutron diffraction experiments were obtained for up to 36 molec uc⁻¹. However, the experimental adsorption of CH₄ continues up to 39 molec uc⁻¹ at 125 K,¹⁷ a higher amount adsorbed than that obtained in the simulation using the ZIF-8AP configuration. In contrast, when the simulations were performed in the ZIF-8HP configuration, we were able to increase the maximum amount adsorbed up to 40 molec uc⁻¹. In this case, CH₄ molecules start filling positions similar to the ones observed in ZIF-8AP. For the maximum uptake at 40 molec uc⁻¹, the changes in the pore shape of ZIF-8HP compared to ZIF-8AP allow the adsorption of additional molecules (Site IV) in the middle of the 4-ring windows placed at the center of the cube face, a situation that did not arise in the case of ZIF-8AP and was not observed experimentally.

3.4. DFT-D binding energy calculations

In order to evaluate the existence of a new adsorption Site IV and to quantify the solid–CH₄ interaction energy during the adsorption process, we calculated the potential energy of a CH₄ molecule in different positions of the solid using DFT calculations with the Grimme dispersion correction. Up to six different energy minima (i.e. configurations) of a CH₄ molecule in ZIF-8AP have been found, plus one extra configuration in ZIF-8HP. Table 3 summarizes the most important adsorption sites (full details are given in the ESI†). The most stable configuration of ZIF-8AP (configuration α₁: Fig. 8) with a binding energy (BE) of −21.0 kJ mol⁻¹, corresponds to a CH₄ located on the top of the CC double bond of the imidazolate linker, with one CH₄ hydrogen pointing directly towards the bond. The shortest H_(CH4)⋯CC distance is 2.7 Å, which is in very good agreement with the experimental value of 2.98 Å.⁵⁵ The strong H_(CH4)⋯CC interaction is reflected by a slight polarization of the charge density on CH₄, as demonstrated by the more positive Mülliken charge on this hydrogen in CH₄ when compared to the isolated molecule (Table S3 in the ESI†). Very similar to α₁, we found two additional configurations, α₂ and α₃, with slightly less negative BEs (i.e. a weaker interaction) of −17.2 and −16.4 kJ mol⁻¹, respectively (see ESI† for more details). The positions and values of the BE for the α₁, α₂ and α₃ adsorption sites are similar to Site I and Site II. This allows us to postulate that at T > 0 K, the experimentally observed adsorption site results in an average of these configurations, corresponding to Site I or II predicted by the GCMC simulations (Fig. 7) and observed in the experiments.⁵⁵

---

An additional stable configuration is β₁, with a BE equal to −17.6 kJ mol⁻¹, obtained for a CH₄ molecule at the center of the 6-ring pore-window (Fig. 9). It arises from the strong interactions between the charge on the aromatic CC double bond of the imidazolate ring, as indicated by the short H_(CH4)⋯CC distance of 3.2 Å, which is in very good agreement with the experimental value of 3.21 Å.⁵⁵ A further stabilising contribution comes from the interaction of CH₄ with the electronic charge on the methyl groups, as indicated by the short H_(CH4)⋯C_(CH3) distance of 3.4 Å. This configuration is related to the channel Site II, as previously described.⁵⁵ The next stable configuration (configuration γ) has a very weak BE of −2.6 kJ mol⁻¹ and corresponds to a CH₄ molecule located at the center of the unit cell, which resembles Site III observed in the previous experiments⁵⁵ and GCMC calculations. This adsorption site plays an important role only at high uptake, when the adsorption of further molecules is stabilized by CH₄–CH₄ interactions.

In the case of ZIF-8HP, the PE required for positioning one CH₄ molecule in the middle of the 4-ring window (Site IV, or Site Θ) is −1.12 kJ mol¹. The small value of this penalty energy, assuming a structural change from ZIF-8AP to ZIF-8HP, suggests that the center of the 4-ring window (Fig. 10a) is thermodynamically accessible to CH₄. This is consistent with Site IV predicted from the GCMC simulations, which becomes populated at higher loading. The main interaction stabilizing CH₄ in this position is between the hydrogens of CH₄ and the two nitrogens at opposite imidazolate linkers, with H⋯N distances of about 2.64 Å (Fig. 10b).

4. Conclusion

In this paper, we have assessed the level of agreement between the experimental and simulated isotherms of a variety of gases (CO₂, CH₄ and other alkanes) in ZIF-8 using generic and adapted force fields. We could exclude that the differences between the experimental and simulated isotherms using generic force fields result from the existence of non-porous defects. By adapting the UFF force field and reducing the Lennard-Jones parameter ε, we achieved excellent agreement between the simulated and experimental isotherms for a wide range of temperatures and a number of different fluids as we were able to reproduce the experimental adsorption mechanism of CH₄. We have shown that a comparison of the data for different temperatures and for different adsorbates is crucial when adapting a force field to produce an accurate and transferable one. Our results also show that two different structural configurations of ZIF-8 are needed to properly describe the adsorption performance of this material, even when no steps are present in the experimental isotherms, demonstrating that ZIF-8 undergoes the same structural change from ZIF-8AP to ZIF-8HP during gas adsorption that we previously observed for N₂ at 77 K. In turn, our results suggest that this flexibility will have a great impact on the diffusivity of the studied gases (i.e. CH₄ and other alkanes) at ambient temperature through the porous network and therefore on the performance of this material as a membrane. Moreover, the transferability of the adapted force field was also successfully tested on ZIF-20.

The results obtained from the DFT-D calculations were found to be consistent with GCMC and experimental observations, illustrating the nature of the CH₄ adsorption sites and their related binding energies. They show that the CC double bond of the imidazolate ligand is the most favorable adsorption site, followed by that at the center of the 6-ring pore-window. The difference in binding energy between the two is relatively small, only about 3.0 kJ mol⁻¹. On the other hand, the DFT calculations confirm the existence of a stable location of CH₄ in the middle of the 4-ring window of the ZIF-8HP structural configuration, generated by the swing effect of the imidazolate linkers. The calculations show the importance of the N atoms of the imidazolate ring and their favorable role in the interaction with CH₄. Our results indicate the role of this site in the structural change of ZIF-8 occurring with increases in the number of adsorbed CH₄ species.

Acknowledgements

We thank Prof. N. A. Seaton for fruitful discussions. D. F .-J. thanks the European Commission for a Marie Curie Intra-European Fellowship (PIEF-GA-2009-236665). C. M. D. thanks the EPSRC for an Advanced Research Fellowship. A. D. G. thanks the College of Science and Engineering, University of Edinburgh for financial support. We acknowledge the use of a HECToR HPC facility via our membership in the UK's Materials Chemistry Consortium funded by the EPSRC (EP/F067496).

References

1. X. C. Huang, Y. Y. Lin, J. P. Zhang and X. M. Chen, Angew. Chem., Int. Ed., 2006, 45, 1557–1559 CrossRef CAS .
2. K. S. Park, Z. Ni, A. P. Côté, 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 CrossRef CAS .
3. D. W. Lewis, A. R. Ruiz-Salvador, A. Gomez, L. M. Rodriguez-Albelo, F.-X. Coudert, B. Slater, A. K. Cheetham and C. Mellot-Draznieks, CrystEngComm, 2009, 11, 2272–2276 RSC .
4. A. Torrisi, C. Mellot-Draznieks and R. G. Bell, J. Chem. Phys., 2010, 132, 044705–044713 CrossRef .
5. D. Liu, C. Zheng, Q. Yang and C. Zhong, J. Phys. Chem. C, 2009, 113, 5004–5009 CAS .
6. R. B. Rankin, J. Liu, A. D. Kulkarni and J. K. Johnson, J. Phys. Chem. C, 2009, 113, 16906–16914 CAS .
7. B. Liu and B. Smit, J. Phys. Chem. C, 2010, 114, 8515–8522 CAS .
8. H. Bux, C. Chmelik, J. M. van Baten, R. Krishna and J. Caro, Adv. Mater., 2010, 22, 4741–4743 CrossRef CAS .
9. S. Keskin, J. Phys. Chem. C, 2010, 115, 800–807 Search PubMed .
10. J. Pérez-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 CrossRef .
11. W. Morris, B. Leung, H. Furukawa, O. K. Yaghi, N. He, H. Hayashi, Y. Houndonougbo, M. Asta, B. B. Laird and O. M. Yaghi, J. Am. Chem. Soc., 2010, 132, 11006–11008 CrossRef CAS .
12. A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard and W. M. Skiff, J. Am. Chem. Soc., 1992, 114, 10024–10035 CrossRef .
13. S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897–8909 CrossRef CAS .
14. D. Fairen-Jimenez, P. Lozano-Casal and T. Düren, in Characterisation of Porous Solids VIII, ed. S. Kaskel, P. Llewellyn, F. Rodriguez-Reinoso and N. A. Seaton, Royal Society of Chemistry, Cambridge, 2009, pp. 80–87 Search PubMed .
15. R. Babarao, S. Dai and D.-E. Jiang, J. Phys. Chem. C, 2011, 115, 8126–8135 CAS .
16. J. Pérez-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 CrossRef .
17. W. Zhou, H. Wu, M. R. Hartman and T. Yildirim, J. Phys. Chem. C, 2007, 111, 16131–16137 CAS .
18. D. Fairen-Jimenez, S. A. Moggach, M. T. Wharmby, P. A. Wright, S. Parsons and T. Düren, J. Am. Chem. Soc., 2011, 133, 8900–8902 CrossRef CAS .
19. C. O. Ania, E. García-Pérez, M. Haro, J. J. Gutiérrez-Sevillano, T. Valdés-Solís, J. B. Parra and S. Calero, J. Phys. Chem. Lett., 2012, 3, 1159–1164 CrossRef CAS .
20. H. Hayashi, A. P. Côté, H. Furukawa, M. O'Keeffe and O. M. Yaghi, Nat. Mater., 2007, 6, 501–506 CrossRef CAS .
21. L. D. Gelb and K. E. Gubbins, Langmuir, 1998, 15, 305–308 CrossRef .
22. J.-H. Yun, T. Düren, F. J. Keil and N. A. Seaton, Langmuir, 2002, 18, 2693–2701 CrossRef CAS .
23. D. Frenkel and B. Smit, Understanding Molecular Simulations: from Algorithms to Applications, 2nd ed. edn, Academic Press, San Diego, 2002 Search PubMed .
24. A. Gupta, S. Chempath, M. J. Sanborn, L. A. Clark and R. Q. Snurr, Mol. Simul., 2003, 29, 29–46 CrossRef CAS .
25. T. Watanabe, T. A. Manz and D. S. Sholl, J. Phys. Chem. C, 2011, 115, 4824–4836 CAS .
26. M. G. Martin and J. I. Siepmann, J. Phys. Chem. B, 1998, 102, 2569–2577 CrossRef CAS .
27. C. D. Wick, M. G. Martin and J. I. Siepmann, J. Phys. Chem. B, 2000, 104, 8008–8016 CrossRef CAS .
28. J. J. Potoff and J. I. Siepmann, AIChE J., 2001, 47, 1676–1682 CrossRef CAS .
29. Q. Yang and C. Zhong, J. Phys. Chem. B, 2006, 110, 17776–17783 CrossRef CAS .
30. P. P. Ewald, Ann. Phys., 1921, 369, 253–287 CrossRef .
31. R. C. Reid, J. M. Pausnitz and B. E. Poling, The Properties of Gases & Liquids, 4th edn, McGraw-Hill Companies, New York, 1987 Search PubMed .
32. A. M. Walker, B. Civalleri, B. Slater, C. Mellot-Draznieks, F. Corà, C. M. Zicovich-Wilson, G. Román-Pérez, J. M. Soler and J. D. Gale, Angew. Chem., Int. Ed., 2010, 49, 7501–7503 CrossRef CAS .
33. T. D. Bennett, J.-C. Tan, S. A. Moggach, R. Galvelis, C. Mellot-Draznieks, B. A. Reisner, A. Thirumurugan, D. R. Allan and A. K. Cheetham, Chem.–Eur. J., 2010, 16, 10684–10690 CrossRef CAS .
34. R. Galvelis, B. Slater, A. K. Cheetham and C. Mellot-Draznieks, CrystEngComm, 2012, 14, 374–378 RSC .
35. CP2K developers page: http://cp2k.berlios.de CrossRef CAS; J. van de Vondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103 CrossRef CAS .
36. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865 CrossRef CAS .
37. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799 CrossRef CAS .
38. J. van de Vondele and J. Hutter, J. Chem. Phys., 2003, 118, 4365 CrossRef CAS .
39. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter, 1996, 54, 1703–1710 CrossRef CAS .
40. M. Krack, Theor. Chem. Acc., 2005, 114, 145 CrossRef CAS .
41. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477 CAS .
42. J. van de Vondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105 CrossRef .
43. J. P. S. Mowat, V. R. Seymour, J. M. Griffin, S. P. Thompson, A. M. Z. Slawin, D. Fairen-Jimenez, T. Düren, S. E. Ashbrook and P. A. Wright, Dalton Trans., 2012, 41, 3937–3941 RSC .
44. J. P. S. Mowat, S. R. Miller, J. M. Griffin, V. R. Seymour, S. E. Ashbrook, S. P. Thompson, D. Fairen-Jimenez, A.-M. Banu, T. Düren and P. A. Wright, Inorg. Chem., 2011, 50, 10844–10858 CrossRef CAS .
45. M. Zhou, Q. Wang, L. Zhang, Y.-C. Liu and Y. Kang, J. Phys. Chem. B, 2009, 113, 11049–11053 CrossRef CAS .
46. C. Gücüyener, J. van den Bergh, J. Gascon and F. Kapteijn, J. Am. Chem. Soc., 2010, 132, 17704–17706 CrossRef .
47. J. van den Bergh, C. Gücüyener, E. A. Pidko, E. J. M. Hensen, J. Gascon and F. Kapteijn, Chem.–Eur. J., 2011, 17, 8832–8840 CrossRef CAS .
48. S. Aguado, G. Bergeret, M. P. Titus, V. Moizan, C. Nieto-Draghi, N. Bats and D. Farrusseng, New J. Chem., 2011, 35, 546–550 RSC .
49. T. K. Trung, P. Trens, N. Tanchoux, S. Bourrelly, P. L. Llewellyn, S. Loera-Serna, C. Serre, T. Loiseau, F. Fajula and G. Férey, J. Am. Chem. Soc., 2008, 130, 16926–16932 CrossRef CAS .
50. P. L. Llewellyn, P. Horcajada, G. Maurin, T. Devic, N. Rosenbach, S. Bourrelly, C. Serre, D. Vincent, S. Loera-Serna, Y. Filinchuk and G. Férey, J. Am. Chem. Soc., 2009, 131, 13002–13008 CrossRef CAS .
51. D. Fairen-Jimenez, N. A. Seaton and T. Düren, Langmuir, 2010, 26, 14694–14699 CrossRef CAS .
52. M. Jeffroy, A. H. Fuchs and A. Boutin, Chem. Commun., 2008, 3275–3277 RSC .
53. F.-X. Coudert, M. Jeffroy, A. H. Fuchs, A. Boutin and C. Mellot-Draznieks, J. Am. Chem. Soc., 2008, 130, 14294–14302 CrossRef CAS .
54. M. Pera-Titus and D. Farrusseng, J. Phys. Chem. C, 2011, 116, 1638–1649 Search PubMed .
55. H. Wu, W. Zhou and T. Yildirim, J. Phys. Chem. C, 2009, 113, 3029–3035 CAS .
56. J. Getzschmann, I. Senkovska, D. Wallacher, M. Tovar, D. Fairen-Jimenez, T. Düren, J. M. van Baten, R. Krishna and S. Kaskel, Microporous Mesoporous Mater., 2010, 136, 50–58 CrossRef CAS .
57. A. Rossin, D. Fairen-Jimenez, T. Düren, G. Giambastiani, M. Peruzzini and J. G. Vitillo, Langmuir, 2011, 27, 10124–10131 CrossRef CAS .

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c2dt30774j

‡ Current address: Department of Chemical and Biological Engineering, Northwestern University, USA.

This journal is © The Royal Society of Chemistry 2012