Kinetic effects in predicting adsorption using the GCMC method – using CO₂ adsorption on ZIFs as an example

Abstract

Using force field parameters developed and validated for zeolitic imidazolate frameworks (ZIFs) and carbon dioxide (CO₂) independently from adsorption data, we predicted CO₂/ZIFs adsorption isotherms using the Grand Canonic Monte Carlo (GCMC) method. The results are in sharp contradiction: the calculated adsorption data agree well with the experimental data for SOD-type ZIF-8, but are more than 100% higher than the experimental data for GME-type ZIFs. Using non-equilibrium molecular dynamics simulations and potential of mean force (PMF) calculations, we reveal that the discrepancies are due to the kinetic blockage which is significant for ZIF-68 but negligible for ZIF-8. This study demonstrates that a force field developed independently from the adsorption data can be used to predict the adsorptions accurately; and the kinetic factor must be considered if the bottlenecks exist in the adsorption paths due to geometric and energetic features of adsorbate and adsorbent. It could be very misleading if the force field parameters are adjusted by fitting the GCMC simulation data to experimental data without considering the kinetic factors.

---

1. Introduction

Zeolitic imidazolate frameworks (ZIFs), a subclass of metal–organic frameworks (MOFs), have been proposed as a potential carbon capture and storage (CCS) material. ZIFs exhibit high CO₂ sorption capacities with high selectivity for CO₂ over CH₄, O₂, H₂, and CO.^1–8 Computational simulations have been applied to study the CO₂ adsorption mechanisms of several types (such as GEM and SOD type) of ZIFs using force field based methods.^8–17

However, contrasting arguments on the accuracy of the force field methods have been reported. Some research groups^10,11 claimed that results obtained for gmelinite (GME)⁶ type ZIFs (ZIF-68 and ZIF-69) using the universal force field (UFF) are in good agreement with experimental data, whereas other groups^12–14 demonstrated that adsorption isotherms predicted using the UFF and DREIDING force fields are significantly overestimated. Liu and Smit¹³ adjusted force field parameters to reproduce experimental isotherms for ZIF-68 and ZIF-69. However, Babarao et al.¹⁴ argued that the overestimates are due to the inaccessibility of small channels in the GME-type of ZIF material. Han et al.¹⁵ developed a new force field for metal–organic frameworks (MOFs) and ZIFs, and obtained excellent CO₂ adsorption on sodalite (SOD) type (ZIF-8) and GEM type ZIFs. McDaniel et al.^16,17 recently reported a new force field that satisfactorily reproduces CO₂ isotherms on different types of ZIFs.

The predictions were generally conducted by using the Grand Canonic Monte Carlo (GCMC) method, which works based on equal chemical potentials between gas phase and adsorbed phase regardless of kinetic factor in the adsorption process. Although the predicted adsorption isomers are in good agreement with experimental data in many cases,^20–22 it does not warrant that the kinetic factor in the adsorption can be always neglected. This issue is often blurred by the force field quality underlying all atomistic simulations. Because UFF and DREIDING force fields are not parameterized for the adsorption purposes, it is tempting to adjust the force field parameters to reach good agreements^18,19 with experimental data. However, the fact that good prediction can also be obtained by modifying the simulation model¹⁴ indicates that the force field quality is not the only factor that impacts the predictive power of GCMC method.

In this work, we examine the kinetic aspect in GCMC simulations by using CO₂ adsorption on ZIFs as an example. We decouple the coupling between the force field quality and kinetic factor by using independently developed and validated force fields. The Murthy–Singer–McDonald (MSM)^23,24 force field for CO₂ was rigorously validated by calculating the vapor–liquid equilibrium (VLE) coexistence curves. By fixing the parameters for CO₂, we derived the force field parameters that describe the interactions between CO₂ and the ZIF frameworks based on high-level quantum mechanics (QM) ab initio data. The QM calculations were conducted at the second-order Moller–Plesset perturbation method (MP2)²⁵ level, but calibrated using the double excitations and perturbative treatment of the triple excitation method CCSD(T)²⁶ with the complete basis set (CBS), which have been known accurate for representing intermolecular interactions based on previous studies.^15,27 Using the validated force fields, we carried out GCMC simulations for two types ZIFs: the GME-type ZIF68, ZIF-69, ZIF-78 and ZIF-79 and the SOD-type ZIF-8. The results were dramatically different: the predicted adsorption isotherms are significantly overestimated for the GME-type ZIFs but fairly accurate for the SOD-type ZIF-8. We examined the kinetic factors in these two different types of ZIFs by using non-equilibrium molecular dynamics (NEMD) simulations and potential of mean force (PMF) free energy calculations. The calculations reveal that the GME-type ZIF channels with small diameters are kinetically blocked by adsorbed CO₂ molecules, while in SOD-type ZIFs such kinetic blockage does not exist.

In the following sections, we first explain how the calculations were done, then present and discuss the computational results, and finally make a summery of this work.

2. Method and model

2.1 ZIF Structures

GEM-type ZIF-68, ZIF-69, ZIF-78, ZIF-79 and SOD-type ZIF-8 were studied in this work. ZIF-68 is formed with a zinc cation (Zn²⁺) coordinated by two 2-nitro-imidazolate (nIM) and two benzimidazole (bIM). Three GME-type ZIFs, ZIF-69, ZIF-78 and ZIF-79, are derived from ZIF-68 by replacing benzimidazole by chlorobenzimidazole (cbIM), nitrobenzimidazole (nbIM) and methylbenzimidazole (mbIM) respectively. ZIF-8 has the Zn²⁺ coordinated by four methylimidazolate (mIM). Table 1 lists the compositions and important physical properties of these materials.

Table 1 Composition and porous characteristics of the ZIFs studied in this work

	Composition	Density^a g cm^−3	dpore^a Å	Surface area^b m^2 g^−1	Free volume^b cm^3 g^−1
a Obtained from ref. 4 and 6.b Calculated using the Connolly volume and surface method as implemented in Materials Studio with 0.75 Å grid intervals and a Connolly radius of 2.25 Å (N2). The data in parentheses are calculated without the small channels.
ZIF-68		1.033	10.3	1972 (975)	0.560 (0.339)
ZIF-69		1.149	7.8	1938 (942)	0.471 (0.282)
ZIF-78		1.175	7.1	1914 (949)	0.487 (0.292)
ZIF-79		1.073	7.5	1879 (927)	0.500 (0.289)
ZIF-8		1.141	11.6	2444	0.510

---

The simulations were conducted on 2 × 2 × 2 super cells constructed using the experimental XRD data^2,4,6 with Material Studio.²⁸ Fig. 1 shows the projections on the X–Y plane and along the Z-axis for (a) ZIF-68 and (b) ZIF-8 models. The GME-type ZIF-68 comprises two one-dimensional channels (aligned in the Z-direction). The small channel comprises small HPR and large GME cages alternatively. The large channel comprises KNO cages. Five adsorption sites are defined in this paper. The center of the HPR cage is denoted as HPR_C. The center and bridge of a GME cage are respectively denoted as GME_C and GME_B. The center and wall of the KNO cage are denoted as KNO_C and KNO_W respectively. The SOD-type ZIF-8 has only one type of cages, each is a polyhedron consisting of 8 faces of six-member rings and 6 faces of four-member rings. The cages are connected with each other through these faces to form three-dimensional channels. Two adsorption sites are defined as the center SOD_C and wall SOD_W of the SOD cage for analysis.

Fig. 1 The topologic features of (a) GME-type ZIF-68 and (b) SOD-type ZIF-8. ZIF-68 consists of small and large one-dimensional channels, three types of cages, and five denoted adsorption sites. ZIF-8 consists of extended three-dimensional channels, one type of cage, and two denoted adsorption sites.

2.2 Ab initio calculations

The van der Waals (VDW) dimers consisting of CO₂ and benzene (B), chlorobenzene (CB), nitrobenzene (NB), methylbenzene (MB), methylimidazole (mIM), and zinc–ammonia complex (ZN) respectively were used as model molecules for the ab initio calculations. The approximate resolution of the identity MP2 (RI-MP2)²⁵ method with the def2-TZVPP basis set were used for geometry optimization. Analytical frequency calculations were performed at the same level to verify the optimized structures were in energy minimums. Based on the optimized structures, RI-MP2 energies were calculated with various basis sets: def2-QZVPP, aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ.^29,30 The RI-MP2 energies were extrapolated to the CBS limit by applying the two-point Helgaker extrapolation scheme³¹ as follows:where X and Y denote the cardinal numbers of the two basis sets used, X = 3 and Y = 4 were used in this work. The data for CCSD(T)/CBS were obtained from the RI-MP2/CBS results²⁷ as follows:

In this work, the “small basis” is aug-cc-pVDZ. The basis set superposition errors (BSSEs) were corrected by the counterpoise method³² for all energetic data. The ab initio calculations were performed using the TURBOMOLE program.³³

2.3 Force field

The rigid model was used for both CO₂ and ZIFs. The intermolecular interactions are represented by pair-wise Coulomb and Lennard-Jones (LJ) 12-6 terms:

The Lorentz–Berthelot combination rules were used to obtain the parameters between unlike atoms.

The partial charges and LJ parameters for CO₂ were taken from the MSM model,²⁴ while the C–O bond length fixed at 1.18 Å and the C–O–C bond angle fixed at 180°. The geometries of ZIFs were fixed at the experimental data. The atomic ESP charges³⁴ on ZIFs were computed at the B3LYP/6-31G* level³⁵ using fragments taken from the ZIFs. Each fragment contains 4 zinc centers and 12 ligands. With the parameters for CO₂ and partial charges of ZIFs fixed, the LJ parameters for ZIFs were obtained by fitting the calculated ab initio potential energy data.

2.4 Monte Carlo simulation

The Gibbs Ensemble Monte Carlo (GEMC)³⁶ simulations at constant volume and temperature were conducted to calculate the vapor–liquid–equilibrium (VLE) coexistence curves of CO₂. Two boxes representing the vapor and liquid phases with a total of 250 CO₂ molecules were used in the simulations. The GEMC moves included swapping molecules between the vapor and liquid boxes, volume exchanges, and translations and rotations of molecules in each of the two boxes. Each of the GEMC simulations included 2 × 10⁶ steps for equilibration and 2 × 10⁶ steps for data collection.

The GCMC simulations^37,38 were used to predict the adsorption isotherms by specifying chemical potentials of CO₂ at specific pressures. The chemical potentials were calculated by using the Widom insertion method³⁸ and the values are listed in the ESI (Table S1†). The GCMC moves included insertions, deletions, translations and rotations of the adsorbate molecules. Each of the GCMC simulations included 1 × 10⁷ steps for equilibration and 1 × 10⁷ steps for data collection. The isosteric heat of adsorption was calculated using the fluctuations in number of adsorbate molecules N and potential energy U:³⁹

In the MC simulations, the LJ interactions were evaluated using a 12.8 Å cutoff with tail corrections. The electrostatic energies were calculated by particle-mesh Ewald (PME) summation with a 12.8 Å real-space cutoff. The block-average method³⁸ was used to estimate uncertainties. The MC simulations were carried out using the Towhee 4.16.8 program.⁴⁰

2.5 Non-equilibrium molecular dynamics (NEMD) and potentials of mean force (PMF)

Non-equilibrium molecular dynamics (NEMD) simulations were conducted on ZIF-68 to investigate the kinetic features of adsorption and desorption. In these simulations, a series of NVT MD simulations were performed on a model comprising vapor–solid–vapor phases along the Z-direction (the direction along the channels). The solid phase was represented by a slab model with 2 × 2 repeat units along the X and Y directions, as well as two and a half repeat units along the Z-direction, as shown in Fig. S1.† The vapor phase was represented by two slabs on both sides of the solid phase and each is 10 Å thick. The process of desorption was simulated with the solid slab filled with CO₂ molecules based on the equilibrated configuration of the GCMC simulations. In these simulations, the vapor phase was emptied every 100 ps to maintain a density gradient in a step-wise manner. The decay of CO₂ density in the solid slab was measured as a function of simulation time. During adsorption simulation, the slab model was initially empty and the vapor phase was filled with CO₂ according to its equation of state at 100 kPa. The vapor phase was refilled to the same state every 100 ps. The increase in adsorbed amount in the slab was measured as a function of simulation time.

The potential of mean force (PMF)⁴¹ was calculated by placing a molecule along a path that connects two cages in the small channel of ZIF-68 or two cages of ZIF-8. A series of umbrella samplings with harmonic force constant of K_i = 5000 kJ mol⁻¹ nm⁻² was conducted to sample the free energy curves. The strong harmonic force was necessary due to strong interactions between CO₂ and the ZIF surfaces, accordingly the umbrella window spacing was set to 0.1 Å to ensure sufficient overlaps between any two windows. In each window, 2 ns of NVT molecular dynamics (MD) simulation was performed to evaluate the potential energies. The weighted histogram analysis method (WHAM)^42,43 was used for evaluating the averaged PMF values.

The MD simulations were performed using the GROMACS software package (version 4.0.3).^44,45 The time step used in the MD simulations was 1.0 fs. The LJ potential was evaluated using a 12.0 Å cutoff with tail corrections. The electrostatic interactions were calculated by PME⁴⁶ summation with a real-space cutoff of 12.0 Å. The temperature was controlled using a Berendsen thermostat⁴⁷ with a coupling constant of 0.1 ps.

3. Results and discussions

3.1 Ab initio data

The optimized structures of the VDW dimers are shown in Fig. 2. The structures can be categorized into five types: a CO₂ molecule positioned at the top and side of the aromatic rings are labeled as “RT” and “RS” respectively; a CO₂ molecule positioned at the head, top, and side of substitution functional groups are labeled “SH,” “ST,” and “SS” respectively. These structures represent three types of interactions:²⁷ electron donor and acceptor, hydrogen-bond and stack π–π interaction.

Fig. 2 Optimized structures of C₆H₆⋯CO₂, C₆H₅NO₂⋯CO₂, C₆H₅Cl⋯CO₂, C₆H₅CH₃⋯CO₂, C₄H₆N₂⋯CO₂, and Zn(NH₂)₄⋯CO₂ complexes obtained at the RI-MP2/def2-TZVPP level of theory.

Five dimers (B-RT, B-RS, NB-SH, CB-ST, and MB-SH) were selected to scan the computational methods. The binding energies of these dimers calculated at various levels of theory are summarized in Table 2. Comparison of the binding energies calculated using different level of theory with those obtained using the high-end CCSD(T)/CBS method indicates that the RI-MP2/def2-QZVPP level of theory yields results close to that obtained at the CCSD(T)/CBS level of theory, with much less computational expenses. The RI-MP2/def2-QZVPP method was used for sampling of potential energy surfaces to derive force field parameters. However, the calculated energy data were corrected using a set of scaling factors. Using the CCSD(T)/CBS values as reference data, we obtained the scaling factors as follows: 0.91 for B-RT, 1.14 for B-RS, 0.95 for CB-ST, 1.23 for NB-SH and 1.19 for MB-SH.

Table 2 Binding energies (in kJ mol⁻¹) calculated at different levels of theory for five representative dimers

	B-RT	B-RS	CB-ST	NB-SH	MB-SH
RI-MP2/def2-TZVPP	−9.87	−3.16	−7.32	−9.84	−3.13
RI-MP2/def2-QZVPP	−11.88	−4.02	−8.34	−11.65	−3.48
RI-MP2/aug-cc-pVDZ	−10.53	−3.62	−7.58	−11.10	−3.03
RI-MP2/aug-cc-pVTZ	−12.36	−4.31	−8.85	−11.95	−3.76
RI-MP2/aug-cc-pVQZ	−12.94	−4.52	−9.35	−12.29	−3.94
RI-MP2/CBS	−13.36	−4.68	−9.71	−12.53	−4.08
CCSD(T)/aug-cc-pVDZ	−7.97	−3.52	−5.76	−10.69	−3.10
CCSD(T)/CBS	−10.80	−4.58	−7.90	−12.12	−4.15

---

The binding energies obtained at the RI-MP2/def2-QZVPP//RI-MP2/def2-TZVPP level are summarized in Table 3. Among different dimer configurations, RT shows the most strong binding energies, which follow the order of C₆H₅NO₂ < C₆H₆ < C₆H₅Cl < C₆H₅CH₃ < C₄H₆N₂, indicating the interaction strength is correlated with the electron-pushing power of the substitution groups. The binding energies of the RS type of dimers are significantly weaker than the RT type of dimers.

Table 3 Binding energies (in kJ mol⁻¹) between CO₂ and benzene (B), chlorobenzene (CB), nitrobenzene (NB), methylbenzene (MB), methylimidazole (mIM), at different configurations (see text). The energies are calculated at the RI-MP2/def2-QZVPP level for the optimized structures obtained at RI-MP2/def2-TZVPP

	B	CB	NB	MB	mIM
RT	−11.88	−11.90	−11.30	−14.45	−17.01
RS(1)	−4.02	−4.71	−5.72	−4.39	−5.34
RS(2)	—	−4.33	−4.97	−4.31	—
SH	—	—	−11.65	−3.48	−2.84
ST(1)	—	−8.34	−12.91	—	—
ST(2)	—		−12.50	—	—
SS	—	−8.71	−9.82	—	—

---

3.2 Force field parameterization

The VLE coexistence curves calculated using the MSM force field²⁴ are compared with the experimental curves in Fig. 3. From 220 K to 290 K, the predicted VLE curves agree very well with the experimental data, which demonstrates that force fields accurately represent CO₂ molecular interactions in both the vapor and liquid phases.

Fig. 3 Comparison of the experimental and calculated VLE curves of CO_2. The calculations are based on the MSM force field.²⁴

Twelve (12) atom types are defined for ZIF atoms, these atom types together with the optimized parameters are listed in Table 4. The charge parameters were expressed in bond-charge increment,⁴⁸ determined from the ab initio ESP charges on fragments taken from the ZIFs. The fit quality is satisfactory, as shown in Fig. 4.

Table 4 Atom types, LJ 12-6 parameters, and bond-charge increment parameters for ZIFs

LJ 12-6 parameters			Bond-charge-increments
Atom type	R^0 (Å)	ε (kJ mol^−1)	Atom pair	ΔQ
Zn	2.80	0.5191	N_Ar–Zn	−0.1865
N_Ar	3.75	0.3139	C_Ar–N_Ar	−0.0283
C_Ar	3.75	0.3897	C_Ar–C_Ar	0.0000
H_Ar	2.73	0.1461	C_Ar–H_Ar	−0.0969
C_Ar_Cl	3.60	0.4567	C_Ar–C_Ar_Cl	0.0461
C_Ar_Ni	3.50	0.5442	C_Ar_Cl–Cl	0.1428
N_Ni	3.00	0.4182	C_Ar–C_Ar_Ni	0.0660
O_Ni	3.25	0.7497	C_Ar_Ni–N_Ni	0.1847
Cl	3.80	1.1260	N_Ni–O_Ni	0.4169
C_Ar_Me	3.75	0.3897	C_Ar–C_Ar_Me	−0.0952
C_Me	3.98	0.3349	C_Ar_Me–C_Me	−0.0070
H_Me	2.73	0.0892	C_Me–H_Me	−0.0510

---

Fig. 4 Comparison of force field (FF) and quantum mechanical (QM) atomic charges calculated for the fragments of ZIFs. Each fragment consists of 4 zinc centers and 12 ligands.

With the CO₂ parameters and the charge parameters of ZIF atoms fixed, the LJ parameters for ZIF atom types were optimized by fitting the ab initio binding energy data for dimers. The 24 LJ parameters were obtained by fitting the energy curves calculated for different configurations of the dimers along the probing paths (Fig. 5). For each probing path, two orientations (parallel and perpendicular to the probe path) of CO₂ were calculated. A total of 321 energy data points were used to fit the LJ parameters. The fit quality is reasonably satisfactory as shown in Fig. 5. The unsigned differences between the ab initio values and the force field values for all structures are less than 2 kJ mol⁻¹ (as Fig. S2† shown).

Fig. 5 Comparison of the ab initio (solid lines) and FF (open symbols) potential energy curves for CO₂ interacting with model molecules along probing paths. CO₂ is oriented parallel (blue) and perpendicular (red) to the probing path.

3.3 Adsorption isotherms

The experimental and calculated adsorption isotherms of CO₂ on the ZIFs at 298 K and from 0 to 100 kPa are present in Fig. 6. The calculated data are based on the optimized force field parameters and GCMC simulations. Sharp differences in the predicted adsorption isotherms are obtained. The predicted curves are reasonably close to the experimental data⁹ for SOD-type ZIF-8, but about 100% overestimated for GME-type ZIF-68, ZIF-69, ZIF-78, and ZIF-79. However, good agreements with the experimental data can be obtained if the small channels in the GEM-type ZIFs are excluded.

Fig. 6 Comparison of the experimental and simulated adsorption isotherms for CO₂ on (a)ZIF-8, (b) ZIF-68, (c) ZIF-69, (d) ZIF-78 and (e) ZIF-79 at pressures ranging from 0 kPa to 100 kPa.

The calculated isosteric heats (Q_st) as function of pressure from 0 kPa to 100 kPa are shown in Fig. 7. For ZIF-8, the initial Q_st value is about 33 kJ mol⁻¹ and it converges to about 18 kJ mol⁻¹. For the GME-type ZIFs, the initial Q_st value is approximately 40 kJ mol⁻¹ to 45 kJ mol⁻¹. The curve drops quickly as the pressure increases to less than 10 kPa and the value converges to about 26 kJ mol⁻¹. Overall, the Q_st values of GME-type ZIFs are about 7–10 kJ mol⁻¹ higher than that of the SOD-type ZIF-8, consistent with the adsorption amount.

Fig. 7 Calculated isosteric heats of adsorption for CO₂ molecules are pressures ranging from 0 kPa to 150 kPa.

Fig. 8 shows the populations of molecules at different adsorption sites as functions of pressure on (a) ZIF-8 and (b) ZIF-68. There are two pieces of information to be noted for this figure. First, there is a clear sequence of the sites to be taken as the pressure increases. For ZIF-8, small amount of CO₂ molecules are adsorbed on SOD_W sites first, and then most molecules are adsorbed on SOD_C sites as the pressure increases. For ZIF-68, CO₂ molecules are adsorbed on HPR_C sites (in small channel) at very low pressure, then on GME_B (in small channel) and KNO_L (in large channel) sites as pressure increases to about 1 kPa, GME_N (in small channel) sites start to populate around 10 kPa, and finally on KNO_C sites (large channel). Secondly, the amounts of molecules on different sites are different at different pressure. It is interesting to note that the amounts at high pressures are affected by the free volumes and surface areas available. On ZIF-8, SOD_C sites six-member-rings interconnecting large cages with limited surface area and free volume, therefore, majority molecules are adsorbed on SOD_W sites. Approximately half of the adsorbed CO₂ molecules are on KNO_L and KNO_W sites of the large channels in ZIF-68, and another half distributed on the HPR and GME sites of small channels.

Fig. 8 Numbers of adsorbed CO₂ molecules at different adsorption sites for (a) ZIF-8, (b) ZIF-68.

There is a correlation between the order of adsorption sites and the binding energies. Fig. 9 are snapshots of the adsorbed molecules on different sites in ZIF-68 and ZIF-8. On ZIF-68, the initial adsorption site is HPR, which corresponds to very large adsorption heat, indicating very strong binding energy at this site. Only one CO₂ molecule can be accommodated and the molecule is positioned parallel to the channel. The molecule interacts with six NO₂ groups. According to ab initio data (Table 3), the binding energy of CO₂ with the NO₂ group is about 12 kJ mol⁻¹, indicating strong binding energy on HPR-C site. The CO₂ molecule interacts with one NO₂ group and with the aromatic rings on the GME_B and KNO_W sites respectively. According to the ab initio data, the CO₂ binding energy with the NO₂ group is similar to that with the benzene ring. At the GME_C and KNO_C sites, the CO₂ molecules interact mainly with the hydrogen atoms. The ab initio data show that binding is significantly weaker than that with the NO₂ group and aromatic ring. Therefore, the CO₂ molecules adsorbed on the GME_C and KNO_C sites are weak. Although the mIM-RT dimer exhibits large binding energy, the top of imidazole rings do not expose to adsorbate well so that the CO₂ molecules mainly interact with the methyl groups and aromatic hydrogen atoms. This explains why the heat of adsorption on ZIF-8 is generally weaker than that on ZIF-68.

Fig. 9 Snapshots of the GCMC simulations on different adsorption sites for (a) ZIF-68 and (b) ZIF-8.

3.4 Kinetic factors

Fig. 10 illustrates the calculated desorption curves as functions of simulation time up to 4 ns for ZIF-68. The simulation data are obtained with initial loads corresponding to equilibrium loads at 3000 kPa. In this case, there are approximately 330 and 180 CO₂ molecules in the large and small channels respectively. The population curve of the large channel decreases much faster than that of the small channels. Decomposition of the data of small channel to GME and HPR cages indicates that desorption occurs only in the GME cages and that the population of molecules in the HPR cages remains constant. At the end of the 4 ns simulation, the number of CO₂ molecules in the large channels is reduced from 330 to 54, in the small channels the number of molecules is reduced from 162 to 110 in the GME cages, and the number of molecules in HPR cages is barely reduced from 18 to 16. A close analysis indicates that the molecules in the small channel are blocked by the molecules in the HPR cages.

Fig. 10 (a) Numbers of CO₂ molecules in different channels as a function of the NEMD simulation time. The curves are obtained from desorption simulations for ZIF-68, starting from the equilibrium configurations of 3000 kPa. Curves for the small channels are decomposed to curves for the GME and HPR cages. (b) The snapshot at 4 ns that show how the molecules are locked in the GME cages because both-ends are blocked by HPR cages.

Fig. 11 shows the uploaded CO₂ as a function of simulation time at 100 kPa. The number of CO₂ molecules in the large channels increases rapidly. At 4 ns, 44 CO₂ molecules are found in the large channels but only 6 CO₂ molecules are found in the small channels, and the molecules are all populated in the HPR cages that are exposed to the vapor phases and none of them is in the inside of the small channels. Extension of the adsorption simulation to 20 ns shows that the total numbers of molecules in the large and small channels are 51 and 7, respectively. Again, the molecules in the HPR cages block other molecules enter the small channels.

Fig. 11 (a) Numbers of CO₂ molecules as a function of the NEMD simulation time for the adsorption process under a vapor phase pressure of 100 kPa. The curve for the small channels is decomposed to curves for the GME and HPR cages. (b) The snapshot at 4 ns that show how the GME cages are inaccessible because of the blocked of HPR cages.

The PMF curves of moving one CO₂ molecule from HPR cage to GME cage in the small channel of ZIF-68 and from one SOD cage to another crossing a six-member ring of ZIF-8 are given in Fig. 12. The cages were filled with CO₂ molecules according to 100 kPa adsorption data for these calculations. The energy barrier height of moving one molecule from HPR cage to GME cage is about 27 kJ mol⁻¹ and the reverse energy barrier height is about 10 kJ mol⁻¹. It is relatively easy to move a CO₂ molecule from GME cage to HPR cage, but it is much more difficult to take the molecule out of the HPR cage. As we have seen from the binding energy analysis, the interaction between CO₂ and HPR cage is very strong, and the HPR cage can only accommodate one molecule. Once a molecule enters the HPR cage, the molecule is trapped and it blocks the entire channel, as indicated by the NEMD simulations discussed above. On the SOD type ZIF-8, the PMF curve shows only an 11 kJ mol⁻¹ free energy barrier for moving one molecule from one SOD cage to another.

Fig. 12 Calculated potential of mean forces (PMFs) of moving a CO₂ molecule (a) from HPR cage to GME cage in the small channel of ZIF-68 and (b) from one SOD cage to another crossing a six-member ring of ZIF-8. The dots indicate entropy-corrected free energy barriers for one molecule escaping from current cage to any adjacent cage.

Because the PMF calculation was restricted along a path connecting two cages, the calculated PMF measures the free energy curve for the sampled path only. To estimate the free energy barriers for one molecule escapes from one cage to any adjacent cage, the probability of escaping along any paths must be counted. Assuming the number of configurations is W at the distance corresponding to the energy barrier, the total number of configurations would be N × W, where N is the number of paths. Using Boltzmann equation, the entropy contribution is:⁴⁹

S = k_B ln(NW) = k_B ln(N) + k_B ln(W) (6)

Note that the last term is included in the PMF calculations. The first term on the right side contributes −k_BT ln(N) to the free energy barrier height at temperature T. The entropy-corrected energy barriers are indicated by dots in Fig. 12(a) and (b) for ZIF-68 and ZIF-8 respectively. For ZIF-68, the change is only −1.7 kJ mol⁻¹, because the path is one dimensional (N = 2). This would reduce the free energy barrier heights to ca. 8 and 25 kJ mol⁻¹ for moving one CO₂ molecule in and out of the HPR cage. For ZIF-8, the change is significant, ca. −5.2 kJ mol⁻¹, due to three-dimensional paths (N = 8). The actual free energy barrier would be even lower due to the fact that another type of paths connected by the 4-member-ring (N = 6) have not been included in the analysis. Since the energy barrier of crossing the 4-memebr-ring would be higher than that for the 6-memebr-ring due to steric effects, we estimated the correction to the total free energy barrier would not be very significant. Nevertheless, the low free energy barrier of ca. 5 kJ mol⁻¹ or less explains why ZIF-8 does not exhibit any kinetic blockages for CO₂ adsorption.

4. Conclusion

By developing and validating force fields independently from experimental adsorption data, we are able to evaluate the GCMC simulation protocol independently from the force field quality that has strong impact to the prediction of adsorption curves. The MSM²¹ force field is validated by simulating VLE data CO₂. The interactions between CO₂ and ZIFs are parameterized using ab initio CCSD(T)/CBS energy data. Using the result force fields, the predicted adsorption isotherms are in good agreement with experimental data for SOD-type ZIF-8, but significantly overestimated for GME-type ZIFs. This sharp discrepancy cannot be attributed to the force field quality.

Using NEMD simulations we found that the adsorption and desorption in the one-dimensional small channels are significantly slower than that in the one-dimensional large channels in GME-type ZIF-68. Furthermore, the calculated PMFs indicate that the small channels of ZIF-68 are blocked by adsorbed CO₂ molecules in the HRP cages. Quantitatively, the free energy barrier is about 8 kJ mol⁻¹ for loading a CO₂ molecule into an HRP cage but it is about 25 kJ mol⁻¹ for removing the molecule from the HPR cage. In ZIF-8, the free energy barriers crossing the 6-member ring connector is only about 5 kJ mol⁻¹.

These findings explain the origin of the discrepancy. The small channels of GME-type ZIFs are completely blocked by the adsorbed CO₂ in the out-most HRP cage. By excluding the small channels, one obtains isotherms in good agreement with the experimental data for all GME-type ZIPs. This work demonstrated that a force field developed independently from the adsorption data can be used to predict the adsorptions accurately; and the kinetic factor must be considered if some kind of bottleneck exists in the adsorptions due to geometric and energetic features of the adsorbate and adsorbent. Another conclusion can be drawn is that it could be very misleading if the force field parameters are adjusted by fitting the GCMC simulated data to experimental data without considering the kinetic factors.

Acknowledgements

This work was partially funded by the National Science Foundation of China (no. 21073119 and 21173146), and National Basic Research Program (no. 2014CB239702).

References

1. J.-R. Li, R. J. Kuppler and H.-C. Zhou, Chem. Soc. Rev., 2009, 38, 1477.
2. A. Phan, C. J. Doonan, F. J. Uribe-Romo, C. B. Knobler, M. O'Keeffe and O. M. Yaghi, Acc. Chem. Res., 2009, 43, 58.
3. H. Hayashi, A. P. Cote, H. Furukawa, M. O'Keeffe and O. M. Yaghi, Nat. Mater., 2007, 6, 501.
4. R. Banerjee, A. Phan, B. Wang, C. Knobler, H. Furukawa, M. O'Keeffe and O. M. Yaghi, Science, 2008, 319, 939.
5. B. Wang, A. P. Cote, H. Furukawa, M. O'Keeffe and O. M. Yaghi, Nature, 2008, 453, 207.
6. R. Banerjee, H. Furukawa, D. Britt, C. Knobler, M. O'Keeffe and O. M. Yaghi, J. Am. Chem. Soc., 2009, 131, 3875.
7. 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.
8. H. Fang, H. Demir, P. Kamakoti and D. S. Sholl, J. Mater. Chem. A, 2014, 2, 274.
9. A. Yazaydin, R. Q. Snurr, T.-H. Park, K. Koh, J. Liu, M. D. LeVan, A. I. Benin, P. Jakubczak, M. Lanuza, D. B. Galloway, J. J. Low and R. R. Wills, J. Am. Chem. Soc., 2009, 131, 18198.
10. D. Liu, C. Zheng, Q. Yang and C. Zhong, J. Phys. Chem. C, 2009, 113, 5004.
11. A. Sirjoosingh, S. Alavi and T. K. Woo, J. Phys. Chem. C, 2010, 114, 2171.
12. R. B. Rankin, J. Liu, A. D. Kulkarni and J. K. Johnson, J. Phys. Chem. C, 2009, 113, 16906.
13. B. Liu and B. Smit, J. Phys. Chem. C, 2010, 114, 8515.
14. R. Babarao, S. Dai and D.-E. Jiang, J. Phys. Chem. C, 2011, 115, 8126.
15. S. S. Han, D. Kim, D. H. Jung, S. Cho, S.-H. Choi and Y. Jung, J. Phys. Chem. C, 2012, 116, 20254.
16. J. G. McDaniel, K. Yu and J. R. Schmidt, J. Phys. Chem. C, 2012, 116, 1892.
17. J. G. McDaniel and J. R. Schmidt, J. Phys. Chem. C, 2012, 116, 14031.
18. 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.
19. D. Liu, Y. Wu, Q. Xia, Z. Li and H. Xi, Adsorption, 2013, 19, 25.
20. J. Fu and H. Sun, J. Phys. Chem. C, 2009, 113, 21815.
21. Y. Sun, T. Ben, L. Wang, S. Qiu and H. Sun, J. Phys. Chem. Lett., 2010, 1, 2753.
22. L. Wang, Y. Sun and H. Sun, Faraday Discuss., 2010, 151, 143.
23. A. Hirotani, K. Mizukami, R. Miura, H. Takaba, T. Miya, A. Fahmi, A. Stirling, M. Kubo and A. Miyamoto, Appl. Surf. Sci., 1997, 120, 81.
24. C. S. Murthy, K. Singer and I. R. McDonald, Mol. Phys., 1981, 44, 135.
25. F. Weigend and M. Haser, Theor. Chem. Acc., 1997, 97, 331.
26. K. Raghavachari, G. W. Trucks, J. A. Pople and M. Head-Gordon, Chem. Phys. Lett., 1989, 157, 479.
27. K. D. Vogiatzis, A. Mavrandonakis, W. Klopper and G. E. Froudakis, ChemPhysChem, 2009, 10, 374.
28. Materials Studio, 4.0V, Accelrys Inc., San Diego, CA, 2005.
29. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796.
30. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007.
31. A. Halkier, T. Helgaker, P. Jorgensen, W. Klopper, H. Koch, J. Olsen and A. K. Wilson, Chem. Phys. Lett., 1998, 286, 243.
32. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553.
33. R. Ahlrichs, M. Bar, M. Haser, H. Horn and C. Kolmel, Chem. Phys. Lett., 1989, 162, 165.
34. M. C. Breneman and K. B. Wiberg, J. Comput. Chem., 1990, 11, 361.
35. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. J. A. Montgomery, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challa-combe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian 03, Revision D.01, Gaussian, Inc., Wallingford CT, 2004.
36. A. Z. Panagiotopoulos, Mol. Phys., 1987, 61, 813.
37. D. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Cambridge University Press, Cambridge, U.K., 2000.
38. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, New York, 1987.
39. R. Q. Snurr, A. T. Bell and D. N. Theodorou, J. Phys. Chem., 1993, 97, 13724.
40. Martin, M. G. MCCCS, Towhee, http://towhee.sourceforge.net/, 2006.
41. G. M. Torrie and J. P. Valleau, Chem. Phys. Lett., 1974, 28, 578.
42. S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman and J. M. Rosenberg, J. Comput. Chem., 1992, 13, 1011.
43. J. S. Hub, B. L. de Groot and D. van der Spoel, J. Chem. Theory Comput., 2010, 6, 3713.
44. D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. Berendsen, J. Comput. Chem., 2005, 26, 1701.
45. D. van der Spoel, A. R. van Buuren, E. Apol, P. J. Meulenhoff, D. P. Tieleman, A. Sijbers and K. Feenstra, Gromacs User Manual, version 4.0.1, Gromacs, Groningen, The Netherlands, 2009, http://www.Gromacs.org.
46. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089.
47. H. J. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684.
48. H. Sun, J. Phys. Chem. B, 1998, 102, 7338.
49. R. M. Neumann, Am. J. Phys., 1980, 48, 354.

---

Footnotes

† Electronic supplementary information (ESI) available: Chemical potential of CO₂, MSM force field parameters for CO₂, the models used in MD simulation, and comparison of energy differences between the FF and ab initio results. See DOI: 10.1039/c4ra03768e

‡ Present Address: School of Chemical and Environmental Engineering, Shanghai Institute of Technology, Shanghai 200240, China.

§ Present Address: Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Kyoto 615-8520, Japan.

---