Understanding the Adsorption and Diffusion of Carbon Dioxide in

---

5004

J. Phys. Chem. C 2009, 113, 5004–5009

Understanding the Adsorption and Diffusion of Carbon Dioxide in Zeolitic Imidazolate Frameworks: A Molecular Simulation Study Dahuan Liu, Chengcheng Zheng, Qingyuan Yang, and Chongli Zhong* Laboratory of Computational Chemistry, Department of Chemical Engineering, Beijing UniVersity of Chemical Technology, Beijing 100029, China ReceiVed: October 22, 2008; ReVised Manuscript ReceiVed: December 30, 2008

In this work, atomic partial charges in the framework atoms of two typical zeolitic imidazolate frameworks (ZIFs), ZIF-68 and ZIF-69, were calculated using density functional theory, and a suitable force field for describing CO2 adsorption in ZIFs was identified. On the basis of this, a combined grand canonical Monte Carlo (GCMC) and molecular dynamics (MD) simulation study was performed to investigate the adsorption and diffusion behaviors of CO2 in ZIFs. The results show that the small pores formed by the nIM linkers in ZIF-68 and ZIF-69 are the preferential adsorption sites for CO2 molecules, with the corners formed by the phenyl rings in the large pores being less preferential adsorption sites. This work demonstrates that the chlorine atoms in cbIM linkers in ZIF-69 lead to enhanced binding energy but reduced diffusivity for CO2, the electrostatic interactions produced by the frameworks are important that cannot be ignored, and, down to 180 K, no steps are found in isotherms. In addition, this work demonstrates that the diffusion of CO2 in ZIFs is likely to be much slower than that in other MOFs. 1. Introduction Zeolitic imidazolate frameworks (ZIFs) are a subclass of metal-organic frameworks (MOFs) with tetrahedral networks that resemble those of zeolites with transition metals (Co, Cu, Zn, etc.) linked by imidazolate ligands.1-8 They show exceptional chemical and thermal stability1,8 and exhibit great promise for a variety of applications especially for gas storage. To date, several investigations have been performed on gas adsorption in ZIFs. For example, Wu et al.6 studied the H2 adsorption sites and binding energies in ZIF-8 using a combined experimental and computational method, Yaghi and co-workers measured the adsorption isotherms of CO2 in several ZIFs,7,8 and the results show that ZIFs are promising candidate materials for CO2 capture. However, little is known about the adsorption features of CO2 molecules as well as its diffusivity in ZIFs, which is crucial for developing new ZIFs materials for practical CO2 capture applications. In contrast to the extensive investigations on other MOFs,9-11 the studies on ZIFs are farther limited to date. Therefore, this work extended our previous work on MOFs11-16 to ZIFs to explore the adsorption and diffusion characteristics of CO2 in ZIFs, which may provide useful information for designing new ZIFs with improved CO2 capture capability. 2. Models and Simulation Method 2.1. ZIFs Structure. In this computational study, two typical ZIFs, ZIF-68 and ZIF-69, were selected and their structures were constructed from the X-ray diffraction (XRD) data7 using Materials Studio Visualizer17 as shown in Figure 1. ZIF-68 and ZIF-69 have the same primitive topology with Zn linked by different imidazolate/imidazolate-type linkers: ZIF-68 contains benzimidazole (bIM) and 2-nitroimidazole (nIM), whereas ZIF69 contains 5-chlorobenzimidazole (cbIM) and 2-nitroimidazole (nIM). The structural properties for these two ZIFs are given * To whom correspondence should be addressed. E-mail: zhongcl@ mail.buct.edu.cn.

Figure 1. Crystal structures of the ZIFs used in the simulations: a) ZIF-68, b) ZIF-69 (Zn, pink; O, red; C, gray; N, blue; Cl, green; H, white).

in Table 1. The accessible surface area (Sacc) and total free volume (Vfree) of each ZIF material were estimated using the Atoms Volume & Surfaces calculation within the Materials Studio package.17 A probe molecule with the diameter equal to the kinetic diameter of CO2 (3.3 Å)18 was used to determine the accessible surface area, whereas the total free volume not

10.1021/jp809373r CCC: $40.75  2009 American Chemical Society Published on Web 02/27/2009

Zeolitic Imidazolate Frameworks

J. Phys. Chem. C, Vol. 113, No. 12, 2009 5005

TABLE 1: Structural and Limiting CO2 Adsorption Properties of the ZIFs Studied in This Work materials

unit cella, Å

cell angle (°)a

dporea Å

Fcrsya g/cm3

Saccb mb/g

Vfreeb cm3/g

qst0c kJ/mol

ZIF-68 ZIF-69

a ) 26.6407, b ) 26.6407, c ) 18.4882 a ) 26.0840, b ) 26.0840, c ) 19.4082

R ) β ) 90, γ ) 120 R ) β ) 90, γ ) 120

10.3 7.8

1.033 1.000

1283 1251

0.588 0.574

23.32 29.04

a

Obtained from the XRD crystal data.7 b Calculated with the Materials Studio package.17 c Obtained in this work with molecular simulation.

TABLE 2: LJ Potential Parameters for CO2 and the ZIFs Used in This Work Atoms

LJ parameters σ

CO2_C CO2_O ZIF_Zn ZIF_C ZIF_O ZIF_N ZIF_H ZIF_Cl

3.43 3.12 2.46 3.43 3.12 3.26 2.57 3.52

(Å)

ε/kB

(K)

52.84 30.19 62.40 52.84 30.19 34.72 22.14 114.23

occupied by the framework atoms was calculated by a probe size of 0.0 Å. 2.2. Force Field. In this work, various atomistic force fields have been tested and we found only the universal force field (UFF) of Rappe´ et al.19 can reproduce well the experimental CO2 adsorption isotherms in both ZIF-68 and ZIF-69 (section 3.1) without refining the potential parameters. Therefore, this force field, different from that used in our previous work, 11 was adopted to describe all of the interactions, including those between CO2 molecules as well as those between CO2 and ZIFs. In this case, CO2 was represented as a rigid linear molecule with three charged Lennard-Jones (LJ) sites located on each atom, and a combination of the site-site LJ potential and Colombic potentials was used to calculate all of the interactions. For consistency, the atomic partial charges in both CO2 molecules and the atoms in the frameworks of ZIFs were calculated using the density functional theory (DFT) method given later. Furthermore, in this work all of the LJ cross interaction parameters were determined by the Lorentz-Berthelot mixing rules. The LJ parameters are given in Table 2, and the atomic partial charges calculated using DFT are given in the following paragraphs. 2.3. Density Functional Theory Calculations. Atomic partial charges for the atoms of adsorbate and adsorbent are required as input parameters in the simulations. In this work, the partial charges of CO2 and the ZIFs were extracted by DFT calculations using the Accelrys Dmol3 code.17 The calculations

Figure 3. Comparison of simulated and experimental32 bulk PVT behavior of CO2 at different temperatures.

were performed on the periodic models of ZIFs using the PW91 GGA density functional20 and the double numerical basis set plus d-functions (DND), which has been considered as a viable alternative to the ab initio methods and can provide reasonable accuracy for computation.21 The partial charges for the atoms were extracted using the Mulliken charge partitioning method,22 which has been successfully employed to study the adsorption of CO2 in MOFs recently.23,24 The calculated charges for CO2 are qo ) -0.272 e and qc ) 0.544 e, and those for the ZIFs are given in Figure 2. 2.4. Grand Canonical Monte Carlo Simulation. In this study, the conventional GCMC simulation was performed to calculate the adsorption of CO2 in ZIFs. Similar to previous works,11-16,25-29 the ZIFs were treated as rigid with atoms frozen at their crystallographic positions during simulations. The simulation box contained 12 (2 × 2 × 3) unit cells with periodic boundary conditions applied in all three dimensions, and no finite-size effects existed by checking the simulations with larger boxes. A cutoff radius was set to 12.8 Å for the LJ interactions, and the long-range electrostatic interactions were handled using the Ewald summation technique with tin-foil boundary condition.30 For each state point, GCMC simulation consisted of 1.0 × 107 steps to guarantee the equilibration, followed by 1.0 × 107 steps to sample the desired thermodynamics properties. The isosteric heat of adsorption qst was calculated from

qst ) RT -

Figure 2. Atomic partial charges of a) ZIF-68 and b) ZIF-69 (Zn, pink; O, red; C, gray; N, blue; Cl, green; H, white).

〈UffN〉 - 〈Uff〉〈N〉 〈N2〉 - 〈N〉〈N〉

-

〈UsfN〉 - 〈Usf〉〈N〉 〈N2〉 - 〈N〉〈N〉

(1)

where R is the gas constant, N is the number of molecules adsorbed, and 〈 〉 indicates the ensemble average. The first and second terms are the contributions from the molecular thermal energy and adsorbate-adsorbate interaction energy Uff, respectively. The third term is the contribution from the adsorbent-adsorbate interaction energy, Usf. A detailed description of the simulation methods can be found in ref 30 and our previous work.11 2.5. Molecular Dynamics Simulation. Equilibrium molecular dynamics (MD) simulation was used to investigate the diffusion behaviors of CO2 in ZIFs at 273 K. The temperature was held constant with a Nose´-Hoover chain (NHC) thermostat as formulated by Martyna et al.31 The simulation boxes used in

5006 J. Phys. Chem. C, Vol. 113, No. 12, 2009

Liu et al.

Figure 4. Comparison of simulated and experimental7 adsorption isotherms of CO2 in ZIF-68 and -69 at 273 K.

the MD simulations are the same as those in the GCMC simulations, and periodic boundary conditions were also applied in all the three dimensions. The long-range electrostatic interactions were evaluated using the Ewald summation method and LJ interactions were calculated with a 12.8 Å cutoff radius. The velocity Verlet algorithm was used to integrate Newton’s equations of motion. The time step used in the MD simulations was taken as 1.0 fs. The adsorbed CO2 densities obtained in previous GCMC simulations were used in MD simulations. Each MD simulation was performed with the randomly insertion of the molecules into the ZIFs lattices first, and then relaxed using approximately 100 000 NVT Monte Carlo moves and 500 000 MD steps following relaxation. After equilibration, 2 500 000 MD steps were run to sample the diffusion properties of interest. At least 10 independent simulations were performed for each loading to estimate the statistical error. During each simulation, the trajectory of the system was saved every 100 steps to subsequently calculate the self-diffusion coefficient Ds by meansquare displacements (MSD) method.12,30 3. Results and Discussion 3.1. Validation of the Force Field. To study the adsorption and diffusion of CO2 in ZIFs, the force field has to be validated. First, the bulk PVT behavior of CO2 gas at different temperatures was calculated using NPT MC simulations to validate the force field for CO2. The results are shown in Figure 3. It can be found that the simulation results are in well agreement with the experimental values.32 Then GCMC simulations were performed to calculate the adsorption isotherms of CO2 in both ZIF-68 and ZIF-69 and compared with the available experimental data.7 The results in Figure 4 show that both experimental adsorption isotherms can be reproduced well by the force field adopted in this work, indicating that the combination of the atomic partial charges obtained from the DFT calculations and LJ parameters from UFF works well for describing CO2 adsorption in both ZIF-68 and ZIF-69. As a result, we can use these parameters to do further study. 3.2. CO2 Adsorption Sites in ZIF-68 and ZIF-69. The knowledge of adsorption sites of adsorbate molecules in a porous material is important in understanding the adsorption mechanism. Therefore, center of mass (COM) probability distributions of CO2 at different pressures in ZIF-68 and ZIF69 were calculated on the basis of all of the configurations recorded during the GCMC simulations, as shown in Figure 5. From parts a and b of Figure 5, it is obvious that at low pressures CO2 molecules are mainly adsorbed in the small pores formed by the nIM linkers (the circle with solid line in Figure 5), with a slight accumulation in other regions of

the frameworks. Whereas with the increase of pressure, more and more CO2 molecules begin to distribute in the large pores formed by the bIM (for ZIF-68) or cbIM (for ZIF-69) linkers (the circle with dot line in Figure 5), as shown in parts c-f of Figure 5. This indicates that the small pores formed by the nIM linkers are the preferential sites for CO2 adsorption in both ZIF-68 and ZIF-69. Furthermore, it was observed that with increasing pressure, CO2 molecules were first adsorbed in the corners formed by the phenyl rings in the large pores (parts c and d of Figure 5), and then occupied the center of the pores (parts e and f of Figure 5). In addition, from Figure 5 we can also see that under identical pressure, the CO2 adsorption capacity is different in the two ZIFs. At low pressures, more CO2 molecules were adsorbed in ZIF-69 than those in ZIF-68 (parts a and b of Figure 5), especially in the small pores formed by nIM linkers. The reason is that the isosteric heat of CO2 adsorption at infinite dilution qst0 in ZIF-69 is larger than that in ZIF-68 (Table 1). At higher pressures, this is also the case that the amount of CO2 molecules adsorbed in the large pores in ZIF-69 are bigger than that in the large pores in ZIF-68 (parts e and f of Figure 5); this can be attributed to the fact that cbIM linker in ZIF-69 has a chlorine atom in phenyl rings as shown in Scheme 1, which enhances the binding energy and thus increases the CO2 adsorption. So, it can be anticipated that if the hydrogen atom in site a of cbIM (Scheme 1) is replaced by another chlorine atom, the CO2 uptake should be improved. On the other hand, if the hydrogen atoms in site b or c are replaced, the effect should be negative because the chlorine atoms in site b or c may block the preferential adsorption sites (the corners formed by phenyl rings) in the large pores as shown in Figure 5. The modification of ZIFs using computational method based on the information mentioned above is ongoing in our group, and the post-synthesis covalent functionalization of ZIFs has been carried out experimentally without altering the original structural integrity very recently.33 3.3. Influences of Electrostatic Interactions. Previous works have shown that electrostatic interactions presented in MOFs are important in CO2 adsorption, which can not only enhance the CO2 storage capacity11 but also change the shape of the adsorption isotherms.34 To investigate the contribution of the electrostatic interactions produced by the frameworks to the CO2 adsorption in ZIFs, additional GCMC simulations were performed by switching off the electrostatic interactions between the framework atoms and CO2 molecules. The percentage contributions of the electrostatic interactions from the frameworks are shown in Figure 6 as a function of pressure (loading). Obviously, the lower the pressure, the larger the contribution of electrostatic interaction produced by the frameworks; this

Zeolitic Imidazolate Frameworks

J. Phys. Chem. C, Vol. 113, No. 12, 2009 5007

Figure 5. Probability distribution plots of COM of CO2 in XY planes in ZIF-68 (a, c, e) and ZIF-69 (b, d, f) at different pressures (Zn, pink; O, red; C, gray; N, blue; Cl, green; H, white).

SCHEME 1: Structures of Linker bIM in ZIF-68 and cbIM in ZIF-69

contribution can be as large as more than 50%, which should not be ignored at low pressures although it decreases with increasing pressure. These results are similar to the situations in other MOFs.11,16 Furthermore, it can be seen from Figure 6

that the contribution is more evident in ZIF-69 at low pressures, whereas it becomes more obvious in ZIF-68 with increasing pressure. The reason may be that the presence of chlorine atoms in cbIM linkers in ZIF-69 leads to larger binding energy and smaller pore free volume, resulting in larger CO2 density under identical bulk pressure. Thus, a quicker decrease of electrostatic contribution from framework occurs in ZIF-69. 3.4. Influences of Temperature on CO2 Adsorption. Step phenomena have been observed in light gas adsorption in MOFs at temperatures lower than room temperature.34 Therefore, we performed additional GCMC simulations on CO2 adsorption isotherms in ZIF-68 and ZIF-69 at low temperatures down to

5008 J. Phys. Chem. C, Vol. 113, No. 12, 2009

Figure 6. Percentage contributions of the electrostatic interactions between ZIFs and CO2 to CO2 adsorption at 273K in ZIF-68 and ZIF69.

Liu et al. examined with MD simulations at 273 K. The results are shown in Figure 8. Obviously, the diffusion of CO2 in ZIF-69 is slower than that in ZIF-68, which is mainly due to the relatively smaller pore size and the steric hindrance effects of chlorine atoms in cbIM linkers in the former. Interestingly, it seems that the diffusivity of CO2 in ZIF-68 and ZIF-69 is much smaller than that in other MOFs, nearly an order of magnitude slower than that in IRMOF-10, IRMOF14, and MOF-177 studied in our previous work.11 Although the pore sizes of the ZIFs are somewhat smaller, another important reason may be due to the structural characteristics of these ZIFs that cause larger steric hindrance. Because no experimental data is available, and this work used a different force field from our previous work,11 a definite conclusion that the diffusion of CO2 in ZIFs is much slower than that in other MOFs cannot be derived at the moment. 4. Conclusions

Figure 7. Adsorption isotherms for CO2 in ZIF-68 at various temperatures.

The computational study in this work shows that CO2 molecules first adsorb in the small pores formed by the nIM linkers in ZIF-68 and ZIF-69. With increasing pressure, they begin to adsorb in the corners formed by phenyl rings in the large pores, followed by the filling of the pores. The presence of chlorine atoms in cbIM linkers in ZIF-69 leads to enhanced binding energy but reduced diffusivity for CO2. This work further shows that the electrostatic interactions produced by the frameworks are important that cannot be ignored, particularly at low pressures; at temperatures down to 180 K, no steps were found in isotherms. In addition, this work demonstrates that the diffusion of CO2 in ZIFs is likely to be much slower than that in other MOFs, caused mainly by the structural characteristics of this kind of materials. However, more detailed study on this topic should be performed to get a definite conclusion. Acknowledgment. The financial support of the National Natural Science Foundation of China (Nos. 20725622, 20876006 and 20821004) is greatly appreciated. References and Notes

Figure 8. Self-diffusivities of CO2 in ZIF-68 and ZIF-69 at 273 K as a function of pressure.

180 K to see if steps will occur in CO2 adsorption in ZIFs. Some typical results in ZIF-68 are taken as examples as shown in Figure 7. Obviously, no steps or inflections in isotherms have been observed, which is also the case for ZIF-69. The reason may be that the bIM or cbIM linkers in this kind of materials prevent the CO2 molecules to form monolayer in the pores and the pore size is also not very large (Table 1), which have been considered as the main reasons that result in the steps and inflections in isotherms.34 3.5. Diffusivity of CO2 in ZIFs. The diffusivity of gas molecules in porous materials is an important property in practical applications. Therefore, the self-diffusivities (Ds) of CO2 in ZIF-68 and ZIF-69 as a function of pressure were further

(1) Park, K. S.; Ni, Z.; Coˆte´, A. P.; Choi, J. Y.; Huang, R.; UribeRomo, F. J.; Chae, H. K.; O’Keeffe, M.; Yaghi, O. M. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 10186–10191. (2) Huang, X.-C.; Lin, Y.-Y.; Zhang, J.-P.; Chen, X.-M. Angew. Chem., Int. Ed. 2006, 45, 1557–1559. (3) Liu, Y.; Kravtsov, V. C.; Larsen, R.; Eddaoudi, M. Chem. Commun. 2006, 14, 1488–1490. (4) Hayashi, H.; Coˆte´, A. P.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Nat. Mater. 2007, 6, 501–506. (5) Tian, Y.-Q.; Zhao, Y.-M.; Chen, Z.-X.; Zhang, G.-N.; Weng, L.H.; Zhao, D.-Y. Chem.sEur. J. 2007, 13, 4146–4154. (6) Wu, H.; Zhou, W.; Yildirim, T. J. Am. Chem. Soc. 2007, 129, 5314– 5315. (7) Banerjee, R.; Phan, A.; Wang, B.; Knobler, C.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Science 2008, 319, 939–943. (8) Wang, B.; Coˆte´, A. P.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Nature 2008, 453, 207–211. (9) Keskin, S.; Liu, J.; Rankin, R. B.; Johnson, J. K.; Sholl, D. S. Ind. Eng. Chem. Res. 2008, (DOI: 10.1021/ie800666s). (10) Ramsahye, N. A.; Maurin, G.; Bourrelly, S.; Llewellyn, P. L.; Serre, C.; Loiseau, T.; Devic, T.; Fe´rey, G. J. Phys. Chem. C 2008, 112, 514– 520. (11) Yang, Q.; Zhong, C.; Chen, J.-F. J. Phys. Chem. C 2008, 112, 1562– 1569. (12) Yang, Q.; Zhong, C. J. Phys. Chem. B 2005, 109, 11862–11864. (13) Yang, Q.; Zhong, C. J. Phys. Chem. B 2006, 110, 655–658. (14) Liu, B.; Yang, Q.; Xue, C.; Zhong, C.; Chen, B.; Smit, B. J. Phys. Chem. C 2008, 112, 9854–9860. (15) Liu, B.; Yang, Q.; Xue, C.; Zhong, C.; Chen, B.; Smit, B. Phys. Chem. Chem. Phys. 2008, 10, 3244–3249. (16) Liu, D.; Yang, Q.; Zhong, C. Mol. Simulat. 2009, 35, 213–219.

Zeolitic Imidazolate Frameworks (17) Acceleys, Inc., Materials Studio, 3.0V; Accelrys Inc.: San Diego, CA, 2003. (18) Breck, D. W. Zeolite Molecular SieVes; John Wiley & Sons: New York, 1974. (19) Rappe´, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard III, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024–10035. (20) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244–13249. (21) Tsuzuki, S.; Lu¨thi, H. P. J. Chem. Phys. 2001, 114, 3949–3957. (22) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833–1840. (23) Ramsahye, N. A.; Maurin, G.; Bourrelly, S.; Llewellyn, P.; Loiseau, T.; Fe´rey, G. Phys. Chem. Chem. Phys. 2007, 9, 1059–1063. (24) Ramsahye, N. A.; Maurin, G.; Bourrelly, S.; Llewellyn, P. L.; Devic, T.; Serre, C.; Loiseau, T.; Fe´rey, G. Adsorption 2007, 13, 461–467. (25) Garberoglio, G.; Skoulidas, A. I.; Johnson, J. K. J. Phys. Chem. B 2005, 109, 13094–13103. (26) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2005, 109, 15760– 15768.

J. Phys. Chem. C, Vol. 113, No. 12, 2009 5009 (27) Frost, H.; Du¨ren, T.; Snurr, R. Q. J. Phys. Chem. B 2006, 110, 9565–9570. (28) Keskin, S.; Sholl, D. S. J. Phys. Chem. C 2007, 111, 14055–14059. (29) Frost, H.; Snurr, R. Q. J. Phys. Chem. C 2007, 111, 18794–18803. (30) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications. San Diego: Academic Press, 2002. (31) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117–1157. (32) Span, R.; Wagner, W. J. Phys. Chem. Ref. Data 1996, 25, 1509– 1596. (33) Morris, W.; Doonan, C. J.; Furukawa, H.; Banerjee, R.; Yaghi, O. M. J. Am. Chem. Soc. 2008, 130, 12626–12627. (34) Walton, K. S.; Millward, A. R.; Dubbeldam, D.; Frost, H.; Low, J. J.; Yaghi, O. M.; Snurr, R. Q. J. Am. Chem. Soc. 2008, 130, 406–407.

JP809373R