Adsorption and Diffusion of Light Gases in ZIF-68 and ZIF-70: A

Sep 8, 2009 - William Morris , Ning He , Keith G. Ray , Peter Klonowski , Hiroyasu ... Yao A. Houndonougbo , Mark Asta , Omar M. Yaghi , and Brian B. ...
0 downloads 0 Views 5MB Size
Download PDF
16906

J. Phys. Chem. C 2009, 113, 16906–16914

Adsorption and Diffusion of Light Gases in ZIF-68 and ZIF-70: A Simulation Study Rees B. Rankin,†,‡ Jinchen Liu,†,‡ Anant D. Kulkarni,‡ and J. Karl Johnson*,†,‡ National Energy Technology Laboratory, Pittsburgh, PennsylVania, 15236, Department of Chemical and Petroleum Engineering, UniVersity of Pittsburgh, Pittsburgh, PennsylVania, 15261 ReceiVed: April 23, 2009; ReVised Manuscript ReceiVed: August 10, 2009

We have computed the adsorption and diffusion of CO2, N2, CH4, and H2 in zeolitic imidazolate framework (ZIF) materials ZIF-68 and ZIF-70 from atomistic simulations. These simulations were performed using geometries obtained from density functional theory (DFT) optimization of the experimental crystal structure of ZIF-68 and on four structures of ZIF-70, one based on the experimental crystal structure and three having different imidazole/nitroimidazole substitution ratios. The framework charges for charge-quadrupole interaction (CQI) terms in our simulations were parametrized with charges obtained from Bader charge decomposition (periodic DFT calculations) and ChelpG charges (cluster DFT calculations). The adsorption and diffusion of the quadrupolar fluids can be dramatically different when using the Bader and ChelpG charges. Agreement between simulations and experiments for the N2 adsorption isotherms in ZIF-68 and 70 is very good when CQI terms are included. In contrast, simulations overpredict the amount of CO2 adsorbed at 298 K compared with experiments. 1. Introduction Rational design of functionalized nanoporous materials for various applications in gas adsorption and separation is of tremendous importance for developing highly selective, lowcost, energy-efficient membranes. The realization of such membranes could make possible the capture of CO2 from coal-fired power plants, for example, at a reasonably small increase in the cost of electricity. The broad group of materials known as metal-organic frameworks (MOFs) is especially promising for rational design, because they are highly tailorable. Studies regarding the performance of this class of materials have focused on the storage of gases, such as methane and hydrogen, and separation of gas mixtures.1-14 In recent years, a new class of materials has been developed that in some ways bridges the gap between MOFs and zeolites, namely, zeolitic imidazolate frameworks (ZIFs).15-19 A simple way to compare zeolites and ZIFs is to note that the traditional zeolite Al(Si)O2 formula is replaced with a T(Im)2 formula unit in ZIF materials. Here, T is a transition metal (tetrahedrally bound), typically zinc, and Im is an imidazole (or imidazole derivative). The structural topologies of ZIFs resemble the traditional zeolite cage topologies, while the chemical functionalization potential of other MOFs is retained in the ability to form heterolinkers using different Im building units. An enormous number of different possible MOFs and ZIFs can be constructed by varying the linker groups and metallic building blocks. The ability to quantitatively predict adsorption and transport properties of fluid mixtures in these materials would be a tremendous benefit to the design of MOFs or ZIFs for specific applications. There are very few molecular simulation studies of gas adsorption and diffusion in ZIFs. Liu et al.20 have reported adsorption isotherms and self-diffusivities for CO2 in ZIF-68 and ZIF-69. They used the UFF force field for the framework ZIF atoms. They also used the UFF potential parameters for * To whom correspondence should be addressed. E-mail: [email protected]. † National Energy Technology Laboratory. ‡ University of Pittsburgh.

the C and O atoms in the CO2 molecule, instead of using a well-established potential model for CO2, such as the one developed and tested by Potoff and Siepmann.21 Liu et al. employed the Mulliken charge partitioning method to calculate partial charges for the solid atoms and for the C and O atoms in the CO2 molecule. The calculated charges on the C and O atoms in CO2 molecules are about 20% smaller than those used in the model of Potoff and Siepmann,21 indicating that the CO2 quadrupole moment of their model is significantly smaller than that of Potoff and Siepmann’s model,21 and likewise smaller than the experimental CO2 quadrupole moment. The bulk properties of their CO2 model were not presented or compared with experiments. Using their CO2 model, the calculated adsorption isotherms at 273 K agree well with experiments for pressure up to 0.8 bar. We here report results of atomic-level modeling of adsorption and diffusion of CO2, H2, N2, and CH4 in ZIF-68 and ZIF-70. These materials were selected for study due to their reported high CO2/CO selectivity.16 We begin with application of periodic density functional theory (DFT) calculations to determine the minimum energy ideal crystal structures of these ZIFs starting from their reported X-ray diffraction (XRD) structures. From these optimized crystal structures, we further utilize DFT calculations to characterize the atomic point charges on the framework atoms of these two ZIFs. Previous work has shown that inclusion of atomic point charges (via charge-quadrupole interactions, for example) can significantly impact predicted adsorption and diffusion of light gases in materials such as CuBTC.22 Using the crystal structures and atomic charges from these DFT calculations, adsorption isotherms are computed from grand canonical Monte Carlo simulations. Additionally, selfdiffusivities and transport diffusivities are computed using equilibrium molecular dynamics simulations. 2. Simulation Techniques and Computational Details 2.1. Simulation Models. For ZIF-68 we employed both experimental and DFT-optimized structures in atomistic simulations. Because the experimental crystal structure from ZIF-70

10.1021/jp903735m CCC: $40.75  2009 American Chemical Society Published on Web 09/08/2009

Light Gases in ZIF-68 and ZIF-70

J. Phys. Chem. C, Vol. 113, No. 39, 2009 16907

TABLE 1: Interaction Potential Parameters for Adsorbate Molecules Used in This Work adsorbate

ref

ε/k (K)

σ (Å)

L (Å)

q (e)

H2-H2 CH4-CH4 N-N (in N2) center of N-N C-C (in CO2) O-O (in CO2)

25 26 27 27 21 21

34.2 148.2 36.4

2.960 3.812 3.318

1.098

27.0 79.0

2.800 3.050

-0.4048 +0.8096 +0.7 -0.35

had too many ambiguous multiple occupancies, we have only used the DFT-optimized structure of this material. The atomic positions of ZIF-68 using the experimental structure were obtained from XRD data,16 from which we removed the atoms belonging to the residual solvent molecules. Our simulations used the universal force field (UFF)23 for the framework atoms in the simulations. We also investigated the DREIDING24 potential for selected systems. The charges on framework atoms are required for simulations of N2 and CO2 for inclusion of the charge-quadrupole interactions (CQIs). We used the framework partial charges calculated as described in section 3.1. A spherical Lennard-Jones (LJ) 12-6 potential was used to model adsorbed H2 and CH4 molecules.25,26 We have ignored H2-framework CQIs because these have been shown to have a very small effect for adsorbed H2 at room temperature.22 A two-center LJ model with three point charges to represent the quadrupole moment was used for N2.27 CO2 was modeled as a three-site LJ with partial point charges located at the center of each LJ site to represent the quadrupole moment of CO2.21 The C-O bond length is 1.16 Å. The adsorbate-adsorbate interaction potential parameters used in our simulations are reported in Table 1. Lorentz-Berthelot mixing rules were employed to calculate the fluid-solid LJ cross-interaction parameters. The fluid-fluid and fluid-solid intermolecular LJ potentials were truncated at 17 Å for adsorption simulations and no long-range corrections were applied. Fluid-fluid and fluid-solid intermolecular LJ potentials were truncated at 13 Å for diffusion simulations, with longrange corrections applied in this case. We have verified that diffusivities calculated using a cutoff radius of 13 Å with longrange corrections gives results that are indistinguishable from calculations using a truncation of 17 Å without long-range corrections. Fluid-fluid electrostatic interactions were truncated at 25 Å with no long-range corrections. Charge-quadrupole interactions were calculated by summation of the charge-charge interactions between each framework atom and each charge site of the adsorbate molecule. These charge-charge interactions were pretabulated by direct calculation of the atomic charge-charge interactions, carried out to long-range.28 2.2. Classical Simulation Techniques. We have used the conventional GCMC technique in this work to compute adsorption isotherms.29,30 We specified the temperature and fugacity of the adsorbing gases and calculated the number of adsorbed molecules at equilibrium. Simulations at the lowest fugacity for each system were started from an empty ZIF matrix. Each subsequent simulation at higher fugacity was started from the final configuration of the previous run. Simulations consisted of a total of 107 trial configurations, with only the last half of the configurations used for data analysis. A configuration is defined as an attempted translation, reorientation (for N2 and CO2), creation, or deletion of an adsorbate molecule. The Lennard-Jones 12-6 equation of state31 was used to obtain the relationship between the bulk pressure and the fugacity for CH4 and H2; this is appropriate because we are using LJ models for these fluids in our simulations. Equations of state based on experimental data were used for CO232 and N233 because no

equations of state are available for the potential models used in our simulation. We have converted the total adsorption obtained from simulation to excess adsorption in order to compare with experimental data34 (total adsorption results for each figure are presented in the Supporting Information). We used equilibrium molecular dynamics (EMD) simulations to compute the self-diffusivities and corrected diffusivities of the pure fluids. EMD simulation details have been described in previous studies of zeolites and carbon nanotubes.28,35-38 A Nose´-Hoover thermostat was employed in NVT-MD simulations. We performed 40 independent MD simulations, each having a simulation length of 8 ns (4 × 106 steps), at each loading considered. Sampling of a large number of independent trajectories was vital in order to compute the corrected diffusivities. After creating initial states with the appropriate loading using GCMC, each system was first equilibrated with EMD for ∼20 ps (1 × 104 steps), prior to taking data points. We have calculated the transport diffusivity from

∂[ln(f)] ( ∂[ln(c)] )

Dt(c) ) D0(c)

T

(1)

where Dt(c) is the transport diffusivity (cm2/s), D0(c) is the corrected diffusivity (cm2/s), f is the fugacity of the adsorbed fluid (in equilibrium with the bulk), c is the concentration of the adsorbate in the ZIF, and [∂[ln(f)]/∂[ln(c)]]T is the thermodynamic correction factor; the thermodynamic correction factor can be calculated from adsorption isotherms obtained from GCMC simulations. 2.3. Density Functional Theory Computational Details. We have performed DFT calculations to quantitatively characterize and refine the atomic crystal structure of ZIF-68 and ZIF70. Existing experimental crystal structures for these materials unfortunately included residual solvent atoms, site double occupancy, and ring orientation uncertainties. Both of these ZIF materials belong to hexagonal crystal systems, specifically possessing space group symmetry P63/mmc.16 Their cagelike topologies are akin to that of the zeolite gmelinite (GME). The reported lattice constants (a/c) for these materials are 26.64 Å/18.49 Å and 27.01 Å/18.02 Å, for ZIF-68 and -70, respectively.16 The crystal structure files of these two materials possess 600 and 480 atoms per unit cell, respectively. The chemical formula of ZIF-68 is C7.06H4.94N3.53O1.59Zn0.71, and the chemical formula of ZIF-70 is C6H5.13N4.87O1.74Zn; the primary difference in the formula between these two ZIF materials is attributed to the benzyl groups that protrude into the large c-axis pores of the ZIF-68 structure.16 We have utilized DFT calculations with the Vienna ab initio simulation package (VASP)39-42 to account for the fully periodic crystal structures. In our calculations we employed the PW9143,44 generalized gradient approximation functional supplied in VASP. We have used the projector-augmented wave (PAW)45,46 method to account for inner core electrons. Plane-wave cutoff energies up to 430 eV were used. Minimization of forces on the atomic coordinates of the ZIF structures was carried out using the standard conjugate gradient algorithm in the VASP package until forces on all atoms were less than 0.03 eV/Å. A Monkhorst-Pack47 scheme for k-point sampling was used with grids corresponding to 2 × 2 × 2 k-points. Atomic charges on the framework atoms were computed from our periodic DFT calculations using the Bader charge decomposition method.48 However, we note that Bader charges computed for the MOF known as CuBTC (or HKUST-1) were shown to be substantially in error when compared to cluster-

16908

J. Phys. Chem. C, Vol. 113, No. 39, 2009

based charge calculations.22 Therefore, we have computed cluster-based DFT charges on representative fragment moieties of these ZIFs using the geometries of the fully periodic crystal structure calculations. ZIF-like clusters extracted from the periodic structures possess artificial dangling bonds. We have saturated these bonds with CH3 groups or H atoms, as appropriate. We have performed calculations on ZIF clusters using the Gaussian 03 software suite.49 The ModRedundant option was employed to perform selective optimization of the moieties associated with saturating the dangling bonds on the ZIF clusters at the HF/6-31G level of theory. We expect this level to be adequate to relax the clusters. Fitting of the electrostatic charges on the atoms was performed at the B3LYP/ 6-311++G(d,p) level of theory using the ChelpG50 algorithm integrated into the Gaussian 03 suite. A benchmarking study on electrostatic potential fitted charge calculation schemes via Chelp, ChelpG, and Merz and Kollman (MK) was carried out by Maciel and Garcia51 for 89 compounds. Their investigations demonstrated the reliability of ChelpG and MK charge computation schemes over the Chelp scheme. They also reported that the MP2 and DFT levels of theory behave in a similar fashion for these schemes. The reliability of the ChelpG scheme has been suggested by Martin and Zipse52 for their work on comparison of methods for charge distribution calculation of the water molecule. It is to be noted that the ChelpG method is generally accepted in similar studies to be the most reliable charge fitting scheme for these types of materials.10,22,53 3. Results and Discussion 3.1. DFT CalculationssZIF-68, -70 Bulk Structure and Charges. In our fully periodic VASP calculations we have constructed the 600 and 480 atom unit cells of ZIF-68 and ZIF70 from the available experimental crystal structures. In each material, residual solvent atoms were removed, site double occupancies were removed, and ring orientation uncertainties were averaged out. We began by geometrically optimizing the full atomic structure of the unit cell by minimizing the forces on the atoms to a force criteria of less than 0.12 eV/Å. We subsequently performed calculations allowing a further optimization of the cell shape/volume with accompanying optimization of the atomic coordinates to the previously mentioned force convergence criterion of 0.03 eV/Å. The optimized lattice constants (a/c) of ZIF-68 and ZIF70 from our fully periodic DFT calculations are 26.73 Å/18.83 Å and 27.21 Å/18.12 Å, respectively. These are close to, but slightly larger than, the experimentally reported values of 26.64 Å/18.49 Å and 27.01 Å/18.02 Å. The optimized structures of ZIF-68 and ZIF-70, from the available crystal structure files, are shown in Figures 1 and 2, respectively. In allowing relaxation of cell shape, deviations away from the reported hexagonal crystal system vectors were not pronounced; the largest changes observed were on the order of 0.09°. The lattice constants calculated in the optimization of the crystal structure are in good agreement with those reported experimentally. The most interesting feature in the optimization of the lattice constant is that, in ZIF-68, the expansion is less symmetric between the a- and c-axes (∼0.3% vs ∼1.8%), whereas in ZIF-70 it is more symmetric about each axis (∼0.7% and ∼0.6%). We postulate that the nature of this relaxation is attributed to ZIF-68 being much more rigid in expansion about the a- and b-axes due to symmetry confinement effects in the benzyl rings protruding into the large pores. We have optimized the structures of three additional ZIF-70 materials. The available experimental crystal structure file has

Rankin et al.

Figure 1. The optimized periodic crystal structure of ZIF-68 is shown looking down the c-axis. Solid black lines denote one the boundaries of one unit cell. Carbon, hydrogen, nitrogen, oxygen, and zinc are colored gray, white, blue, red, and green, respectively.

Figure 2. The optimized periodic crystal structure of ZIF-70 is shown looking down the c-axis. Solid black lines denote one the boundaries of one unit cell. The atom types are the same as in Figure 1.

a IM:nIM (IM ) imidazole, nIM ) nitroimidazole) ratio of 1:3, which is in disagreement with the published stoichiometric ratio of ∼2:2.16 We have therefore manually altered (via the appropriate highest symmetry deletions) this structure to be in the nominal IM:nIM ratio of 2:2. We have also generated two different 3:1 structures (see Figure S1, Supporting Information). We have optimized these structures and generated N2 adsorption isotherms for each of the four different structures in order to assess how the IM:nIM ratio affects adsorption. A discussion of the effect of the differing ZIF-70 structures on adsorption is presented in section 3.2. Atomic point charges corresponding to the framework atoms of ZIF-68 and ZIF-70 were computed using the fully periodic structures mentioned above. As we have done in our previous work with CuBTC,22 we have computed these charges with two methods. The first method used was the Bader charge decomposition method48 within VASP. The second method used was to construct cluster-based fragment moieties of ZIF-68 and ZIF70 and to fit their charges using the ChelpG50 electronic partitioning scheme in Gaussian 03. The fragment moieties of ZIF-68 and ZIF-70 used for the cluster DFT calculations are

Light Gases in ZIF-68 and ZIF-70

J. Phys. Chem. C, Vol. 113, No. 39, 2009 16909

TABLE 2: Normalized Atomic Partial Charges in This Work for the Framework Atoms of ZIF-68 and ZIF-70a ZIF-68

ZIF-70

atomic species

small cluster

large cluster

Bader charges (periodic)

small cluster

large cluster

Zn N1 N2 N3 N4 N5 O C1 C2 C3 C4 C5 C6 C7-C10(av) H(av)

0.550 -0.250 -0.315 -0.250 -0.315 0.570 -0.500 0.145 0.145 0.145 0.145 0.030 0.230 -0.150 0.110

0.668 -0.302 -0.257 -0.302 -0.257 0.658 -0.572 0.148 0.148 0.148 0.148 0.128 0.138 -0.152 0.098

1.200 -1.160 -1.250 -1.160 -1.250 0.310 -0.480 0.350 0.390 0.350 0.390 0.960 1.310 0.040 0.070

0.677 -0.283 -0.353 -0.283 -0.353 0.577 -0.603 0.107 0.107 0.107 0.107 0.177 0.177 N/A 0.087

0.689 -0.231 -0.341 -0.231 -0.341 0.509 -0.531 0.009 0.009 0.009 0.009 0.264 0.264 N/A 0.089

Figure 3. The adsorption isotherms for CH4 in ZIF-68 at 298 K using experimental and DFT-optimized structures are shown.

a All results are in atomic charge units (e). The as-computed atomic partial charges are listed in Table S1 in the Supporting Information. The atomic species labels are identified in Figure S2 in the Supporting Information.

presented in the Supporting Information. From these results, charges computed from nonperiodic clusters must be renormalized in order to produce a charge neutral unit cell for the periodic structure; this is because the clusters do not possess the same stoichiometry as the bulk material. We renormalized the charges using

(2)

Figure 4. Comparison of simulated and experimental adsorption isotherms for N2 in ZIF-68 and ZIF-70 at 77 K. The different IM:nIM ratios correspond to the structures shown in Figure S1 (Supporting Information). The charge-quadrupole interactions between N2 and the framework charges are included.

where q′i is the renormalized charge for species i, qi is the charge computed from the cluster calculation, xi is the fraction of species i having charge qi in the fully periodic material, and N is the number of distinct species. This is the same method we have previously applied in our study of charges in the framework of CuBTC.22 The normalized atomic charges are summarized in Table 2. We find reasonably consistent agreement between the charges computed for ZIF-68 and ZIF-70; this result is not surprising given the many similarities in their structures. However, results from the Bader charge calculations are problematic in ZIF-68 and ZIF-70, much as they were observed to be for CuBTC.22 In particular, the charges on Zn, most of the N atoms, and a few of the C atoms are much larger in magnitude as computed from Bader analysis than from ChelpG. Labels identifying the atomic species in Table 2 are presented in Figure S2 (Supporting Information). Our results show that in order to approach convergence in the charges on Zn and equivalent N atoms in the cluster charge calculations, the clusters must be very large in size, approaching ∼200 atoms or more. Hence, there is need for a more accurate and reliable charge partitioning method that can be used within periodic DFT calculations for these materials. 3.2. Simulation Results: Adsorption in ZIF-68 and ZIF70. Classical GCMC simulations for the adsorption of the light gases CO2, N2, CH4, and H2 in the ZIF-68 and ZIF-70 materials have been performed. We begin our discussion by examining the effects of using the as-is (albeit “cleaned up”) experimental crystal structure and the DFT-optimized crystal structure of ZIF68 with respect to changes in the adsorption of CH4 at 298 K.

We do this to assess the effects due to the uncertainty in the experimental crystal structure. The results of these simulations up to 55 bar are shown in Figure 3. The difference in the CH4 isotherms for the two structures is small, although the DFToptimized structure consistently gives slightly less adsorption. Similar results were observed for CO2 and H2 (not shown). This result was qualitatively the same as observed in our previous work on adsorption and diffusion of light gases in CuBTC.22 On the basis of this result, it may seem apparent that one can simply use the available crystal structures of these ZIF materials to obtain simulation results without requiring the time-consuming DFT calculations to optimize the structure. However, our analysis of the adsorption of these gases in ZIF-70 provides an important counterexample. The results of the N2 adsorption isotherms in ZIF-68 and ZIF70 at 77 K are plotted in Figure 4. For ZIF-70, there are various simulation results depending on the framework structure used in the simulation (see section 3.1 for discussion). As the stoichiometric ratio of IM to nIM units in the ZIF-70 framework decreases, the adsorption of N2 decreases by over 20%. The closest agreement between experimental and simulated adsorption of N2 is at the composition of IM:nIM equal to 2:2, in agreement with the published stoichiometry16 but in contrast to the published XRD data.16 The differences in the experimental and simulated isotherms in Figure 4 for ZIF-70 are not due to the utilization of a specific force field, since UFF and DREIDING force fields present similar adsorption isotherms. Comparison of the simulated and experimental N2 isotherms in Figure 4 indicates that, if the experimentally determined

N

q′i ) qi -

∑ xiqi i)1

16910

J. Phys. Chem. C, Vol. 113, No. 39, 2009

Rankin et al.

TABLE 3: Unit Cell Masses and Pore Volumes for ZIF-68 and ZIF-70a material ZIF-68-exp ZIF-68-opt ZIF-70-opt (IM:nIM ) 1:3) C144H108N132O72Zn24 ZIF-70 (IM:nIM ) 2:2) C144H120N120O48Zn24 ZIF-70-1 (IM:nIM ) 3:1) C144H132N108O24Zn24 ZIF-70-2 (IM:nIM ) 3:1) C144H132N108O24Zn24

mass, g/mol of UC

helium pore volume, cm3/g

7069.862 7069.862 6408.398

0.495 0.484 0.629

5868.427

0.729

5328.456

0.828

5328.456

0.854

a Pore volumes computed from experimental crystal structures are compared with optimized/restructured crystal structures. Two ZIF-70 structures both having an IM:nIM ratio of 3:1 but different site occupancies are given. The helium pore volumes were computed from simulations according to the method of Myers and Monson.54

Figure 5. Adsorption isotherms for CO2 in ZIF-68 at 273 K using framework charges from Bader analysis (diamonds) and large (circles) and small (open squares) cluster calculations. See Table 2 for the charge values.

crystal structure is free from ambiguities and structural inconsistencies, then using those crystal structures in adsorption simulations appears to be justified. However, for materials where the crystal structure is ambiguous, manual manipulation of the available experimental crystal structure or careful application of DFT calculations to optimize the crystal structure may still be required as an input into the simulations. Table 3 gives the mass and pore volumes of ZIF-68 and the different ZIF-70 structures examined. Adsorption simulations have been performed both with and without charges on the framework atoms for computing CQI between N2 or CO2 and the framework charges. Our previous work on CuBTC has shown that framework charges can significantly affect the adsorption and diffusion of light gases in pores of such materials.22 Note that we have not included CQI for H2 because we only examine adsorption near room temperature for H2, where CQI have been shown to have negligible effect on adsorption.22 The adsorption isotherms for CO2 in ZIF-68 at 273 K using point charges derived from the periodic Bader charge analysis and the large and small representative clusters (Table 2) are plotted in Figure 5. These simulations reveal that use of the charges derived from the smaller cluster does not significantly impact the adsorption properties. This is important because of the very high computational cost required to obtain these charges. This result is surprising because the charges on most of the structurally

Figure 6. Isotherms from experiments and simulations for CO2 in ZIF68 at 273 K. Experiments16 are shown as points. Simulations were performed using the UFF potential with charges (solid), DREIDING with charges (dashed), UFF without charges (dot dashed), and DREIDING without charges (dotted).

equivalent atoms in the smaller cluster differ from those on the larger cluster by about 10-30%. Adsorption simulations using charges from the Bader charge analysis method give results that are unphysical, showing essentially pore-filling at extremely low pressures of CO2. The unphysical behavior is a direct result of the large magnitude charges on many of the atoms from the Bader analysis. Adsorption isotherms for CO2 in ZIF-68 at 273 K computed from simulations using different potentials and from experiments are plotted in Figure 6. We have used both the UFF and DREIDING potentials, each with and without framework charges. We note that simulations using the UFF and DREIDING potentials give significantly different results, with the UFF potential resulting in higher adsorption than the DREIDING potential. Addition of charges substantially increases the adsorption for both potentials, especially at very low pressures (below about 0.1 bar). Note also that the increase due to charges is about the same for both potentials. The isotherm computed using the DREIDING potential without charges is in fairly good agreement with experiments; addition of charges overestimates the amount adsorbed rather dramatically, as seen in Figure 6. In contrast, Liu et al.20 obtained very good agreement between simulations and experiments for CO2 adsorption in ZIF-68 at 273 K for pressures up to 0.8 bar. As noted above, they developed their own CO2 potential rather than using an existing model from the literature. It is instructive to compare the pure fluid properties predicted by their model with experimental data because a physically reasonable model must give both good fluid-fluid and fluid-solid interactions. We have therefore computed the bulk phase isotherm at 273 K for models of Liu et al.20 and of Potoff and Siepmann21 and have compared both sets of calculations with experimental data. The results are plotted in Figure 7. The model developed by Potoff and Siepmann agrees very well with experiments for pressures up to the saturation pressure. However, the model used by Liu et al.20 significantly underestimates the bulk density compared with experiments (∼30% lower at the saturation pressure). We have also carried out Gibbs ensemble55,56 calculations at 273 K for both potential models using the Towhee program.57,58 Our calculations with the Potoff-Siepmann model agree with their published results,21 and are in good agreement with experiments for the saturated vapor and liquid densities. In contrast, we were unable to locate the vapor-liquid equilibrium state for the CO2 model of Liu et al. This indicates that 273 K lies above the

Light Gases in ZIF-68 and ZIF-70

Figure 7. Bulk gas phase isotherms from experiments and simulations for CO2 at 273 K. Experiments59 are shown as points. Simulations using the models of Liu et al.20 and Potoff and Siepmann21 are shown as dashed and solid, respectively.

Figure 8. Isotherms from experiments and simulations for CO2 in ZIF70 at 273 K. Experiments16 are shown as points. Simulations were performed using the UFF potential with charges (solid), UFF without charges (dot dashed), and DREIDING without charges (dotted).

critical temperature for this model, due to the fluid-fluid interactions in the Liu et al. model being too weakly attractive. Note that the weak fluid-fluid interactions lead to significant errors in the isotherm at fluid densities higher than about 20 kg/m3 (see Figure 7). The experimental density of adsorbed CO2 in the pore of ZIF-68 is ∼260 kg/m3 at 1 bar,16 which is significantly higher than the density of bulk CO2 gas at saturation. Thus, the fluid-fluid CO2 interactions in the model of Liu et al. must be substantially weaker under these conditions than they should be. Hence, the good agreement between experiments and simulations of Liu et al. for CO2 adsorption in ZIF-6820 may be due to a fortuitous cancellation of errors. Adsorption isotherms of CO2 in ZIF-70 are plotted in Figure 8. All simulation results for ZIF-70 from this point on will be for the 2:2 IM:nIM ratio structure, unless otherwise stated. UFF and DREIDING potentials give isotherms that are very similar in this case and both are below the experimental data. The larger difference between UFF and DREIDING for ZIF-68 compared with ZIF-70 is a result of the former having smaller pores than the latter. Differences in the solid-fluid interaction potentials will be magnified for a material with small pores. Addition of charges again results in predicted values that are much larger than those from experiments. Note that we have not included the DREIDING potential with charges in Figure 8 because the results are almost identical to the data from UFF with charges. It is tempting to conclude that inclusion of charges results in a solid-fluid potential that is too attractive to be physically realistic, since the experimental data are overpredicted by a

J. Phys. Chem. C, Vol. 113, No. 39, 2009 16911

Figure 9. Simulated adsorption isotherms for CO2 in ZIF-68 and ZIF70 using the UFF framework potential with CQI.

substantial amount. However, there is no physical basis for ignoring framework charges. The only reliable basis for ignoring charges is to perform calculations with and without framework charges and verify that the results are essentially the same in either case. Likewise, empirically adjusting potential parameters to match a single experimental isotherm is unjustified. Moreover, our simulations are in good agreement with the experimental N2 adsorption isotherms, at least for P/P0 > 0.2 (see Figure 4). However, the simulations overpredict the amount of N2 adsorbed at lower values of P/P0 (see Figure S4, Supporting Information), similar to that shown in Figures 6 and 8 for CO2. We note that the vapor pressure of CO2 at 273 K is 34.7 bar, so that the entire pressure range covered in Figures 6 and 8 corresponds to 0 < P/P0 < 0.03 (see Figures S5 and S6 of the Supporting Information for a plot for the excess adsorption vs reduced pressure). We therefore predict that experimental isotherms for CO2 at high pressures will be in better agreement with simulations than the low pressure data. We have performed adsorption isotherm simulations up to pressures of 30 and 50 bar at 273 and 298 K, respectively, for CO2, in ZIF-68 and ZIF-70 (2:2) to examine the adsorption behavior at high pressures. The 273 K isotherms for both ZIF68 and ZIF-70 (Figure 9) exhibit a maximum loading in the excess adsorption around 25 bar; the maximum loading at 298 K is shifted to higher pressures, and the maximum amount is decreased compared with the 273 K isotherm, as expected. The total adsorption isotherms (see Supporting Information, Figure S7) monotonically increase in loading; isotherms at both temperatures asymptotically approach the same loading. The ratio of the total adsorption capacities for ZIF-70:ZIF-68 (see Figure S7, Supporting Information) is approximately 1.4:1 at both temperatures, in accordance with the ratio of their respective pore volumes (∼1.5:1 from Table 3). Adsorption isotherms have also been computed for CH4 and H2 at high pressures and are plotted in Figures 10 and 11. Only the 298 K results are shown for clarity. The CH4 adsorption curves in ZIF70 and ZIF-68 cross at approximately 10 bar. The crossing is a result of a competition between adsorption potential and pore volume. At lower pressures the solid-fluid interaction potential dominates, and ZIF-68, having smaller diameter pores, exhibits a larger uptake than ZIF-70. At higher pressures pore volume is the dominant factor in the adsorption. The expected ∼1.4:1 saturation adsorption capacity behavior is observed for CH4 at about 55 bar (see Figure S8, Supporting Information), as was seen for CO2 at lower pressures. H2 is found to have a higher adsorption capacity in ZIF-68 than in ZIF-70 up to 100 bar. The result for this is attributable to the same effects seen at

16912

J. Phys. Chem. C, Vol. 113, No. 39, 2009

Rankin et al.

Figure 10. Simulated adsorption isotherms for CH4 in ZIF-68 and ZIF-70 at 298 K.

Figure 13. The average density of CO2 in ZIF-68 and ZIF-70 normalized to the saturated liquid CO2 density of the bulk as a function of the adsorption loading at 298 K.

Figure 11. Simulated adsorption isotherms for H2 in ZIF-68 and ZIF70 at 298 K.

Figure 14. The simulated diffusivities of CO2 in ZIF-68 at 298 K.

Figure 12. The simulated isosteric heats of adsorption for CO2 in ZIF68 and ZIF-70 at 298 K.

low pressure with CH4, where again the dominant factor determining the adsorption is solid-fluid interaction potential. We have calculated the isosteric heats of adsorption for CO2 in ZIF-68 and ZIF-70 as a function of coverage as an additional way to characterize the adsorption behavior. The calculated values for the isosteric heats are plotted as a function of the excess adsorption in Figure 12. These results are from simulations at 298 K. The isosteric heat, qst, initially declines steeply as a function of coverage due to the highest energy sites being occupied. The small cages in the ZIFs are the energetically most favorable adsorption sites. The isosteric heats for both materials then exhibits a plateau region followed by an increase. The increase in qst at high loading is due to the increase in fluid density, giving rise to enhanced fluid-fluid interactions. This can be seen by comparing Figures 12 and 13, the latter showing the average density of CO2 inside the pores divided by the bulk

density of saturated liquid CO2 at 298 K, as a function of excess adsorption. We note that the increase in qst occurs about where the average density of CO2 in the pore exceeds the liquid density of CO2. Similar results have been observed in other nanoporous materials.60 In the ZIF-70 pores, the isosteric heats decline again after reaching a density significantly higher than in the bulk fluid, indicating that fluid-fluid repulsion is responsible for the decrease in qst at the highest loadings. 3.3. Simulation Results: Diffusion in ZIF-68 and ZIF-70. We have performed equilibrium molecular dynamics simulations to compute the self-diffusivities (Ds) and transport diffusivities (Dt) of the light gases H2, CH4, and CO2 in ZIF-68 and ZIF-70. We have examined the effect of CQI on diffusivity of CO2 in both ZIFs. All results presented utilize the UFF potential for the framework atoms. The self-diffusivity and transport diffusivity of CO2 in ZIF-68 at 298 K, both with and without CQI, are plotted in Figure 14. The results from these simulations clearly reveal the importance of including CQI terms in the diffusion calculations. At low CO2 loading the addition of CQI results in approximately 1 order of magnitude drop in the diffusivities, a result that reflects the significantly enhanced adsorption on high energy sites due to CQI. As loading increases, Ds computed with CQI first increases and then decreases. This is due to higher loading of CO2 in the small pores decreasing the binding energy, as seen by the decrease in qst in Figure 12, and hence decreasing the diffusion barrier height. At high loadings Ds computed with and without CQI converge toward the same value. This indicates that CQIs have a moderate impact on diffusion barriers at high loadings, where Ds is apparently dominated by fluid-fluid collisions, as indicated by the drop in Ds at high loadings. The Dt values with and without CQI converge toward the same values as well, but not as closely as Ds. Note also that Dt monotonically increase,

Light Gases in ZIF-68 and ZIF-70

J. Phys. Chem. C, Vol. 113, No. 39, 2009 16913

Figure 15. Calculated diffusivities of CO2 in ZIF-70 at 298 K using the UFF potential with framework charges.

Figure 17. Calculated diffusivities of H2 in ZIF-68 and ZIF-70 at 298 K.

decrease linearly and Dt increases in a roughly linear fashion. The low loading diffusivities of H2 in these materials are seen to be about an order of magnitude larger than for CH4, or about 2 orders of magnitude larger than CO2; the weaker interactions with the framework and the lighter mass both act to enhance diffusivity of H2. 4. Conclusions

Figure 16. Diffusivities of CH4 in ZIF-68 and ZIF-70 at 298 K calculated using the UFF framework potential.

irrespective of inclusion of CQI, but do so more dramatically when CQI are included. The diffusivities presented, both Ds and Dt are averages over the three principal crystallographic axes, as given by

D ) 1/3(Dx + Dy + Dz)

(3)

In ZIF-68, the diffusion is significantly closer to 1-D due to the pore structure. The diffusivity along the c-axis is typically 1-3 orders of magnitude larger than in the other two directions. The self-diffusivity and transport diffusivity of CO2 in ZIF70 are plotted in Figure 15. The zero-coverage diffusivity for CO2 in ZIF-70 is about 1 order of magnitude larger than in ZIF-68. The lower diffusivities in ZIF-68 are due to the narrower pores, the rougher potential energy surface that results from the benzyl groups pointing into the pores, and the nearly 1-D nature of the ZIF-68 pores. Only the simulations including CQI are shown in Figure 15. Both Ds and Dt are qualitatively similar in ZIF-68 and ZIF-70. The self-diffusivities and transport diffusivities of CH4 and H2 in ZIF-68 and ZIF-70 at 298 K are plotted in Figures 16 and 17, respectively. CH4 diffusivities in ZIF-68 and ZIF-70 are about an order of magnitude higher than those for CO2. Moreover, Ds for CH4 does not increase substantially with increased loading, but rather remains relatively flat over a range of loadings before decreasing. Similar to CO2, Dt for CH4 monotonically increases, with the greatest increase occurring in ZIF-68. At high loadings the ratio of Dt in ZIF-70 and ZIF68 is about 3 and is due to the factor of 3 appearing in eq 3, since ZIF-68 is essentially 1-D. The self-diffusivity and transport diffusivity of H2 in ZIF-68 and ZIF-70, plotted in Figure 17, are comparatively constant with loading, although Ds does

Our calculations have shown that using the structural information for ZIF-70 reported from XRD data gives results that are in poor agreement with experimental N2 adsorption data. Importantly, the error is due to the X-ray data having the wrong stoichiometry, a result of partial occupancy of sites containing nIM groups. We have found that, when the XRD-derived structure is manually adjusted to have the reported 2:2 IM:nIM ratio, then the simulation results are in good agreement with the experimental N2 adsorption isotherm. In contrast, the adsorption isotherms for ZIF-68 using the raw XRD structure and the DFT optimized structure were nearly identical. Hence, if a reliable, unambiguous crystal structure from experiment is available, computationally expensive DFT calculations to refine the structure may be unnecessary. Inclusion of framework charges to account for N2-ZIF CQI gives isotherms in good agreement with experiments, for both ZIF-68 and -70. However, simulations are in poor agreement with experiments for CO2 adsorption in these materials at low pressures when charges are incorporated. The origin of this discrepancy is not currently understood. We suggest that CO2 adsorption at high pressures be measured on carefully synthesized and activated ZIF samples and compared with our predictions for high pressure adsorption. Note that activation of some MOFs is known to have a large impact on adsorption.34 The effects of the inclusion of charges on diffusivities are seen to be far more pronounced at low loadings than at high loadings; at low loadings the diffusivities may differ by more than 1 order of magnitude when charges are included versus when they are ignored. Acknowledgment. This work was performed in support of the National Energy Technology Laboratory’s ongoing research in the area of carbon management under the RDS contract DEAC26-04NT41817 Supporting Information Available: Structures of various IM:nIM ratios in ZIF-70, fragment moieties used in charge calculations, experimental and simulation isotherms for N2 and CO2, total adsorption isotherms, and diffusivities as a function of total loading. This material is available free of charge via the Internet at http://pubs.acs.org.

16914

J. Phys. Chem. C, Vol. 113, No. 39, 2009

References and Notes (1) Eddaoudi, M.; Li, H.; Yaghi, O. M. J. Am. Chem. Soc. 2000, 122, 1391. (2) Wang, Q. M.; Shen, D. M.; Bulow, M.; Lau, M. L.; Deng, S. G.; Fitch, F. R.; Lemcoff, N. O.; Semanscin, J. Microporous Mesoporous Mater. 2002, 55, 217. (3) Du¨ren, T.; Sarkisov, L.; Yaghi, O. M.; Snurr, R. Q. Langmuir 2004, 20, 2683. (4) Du¨ren, T.; Snurr, R. Q. J. Phys. Chem. B 2004, 108, 15703. (5) Rowsell, J. L. C.; Spencer, E. C.; Eckert, J.; Howard, J. A. K.; Yaghi, O. M. Science 2005, 309, 1350. (6) Millward, A. R.; Yaghi, O. M. J. Am. Chem. Soc. 2005, 127, 17998. (7) Rowsell, J. L. C.; Yaghi, O. M. Angew. Chem., Int. Ed. 2005, 44, 4670. (8) Wong-Foy, A. G.; Matzger, A. J.; Yaghi, O. M. J. Am. Chem. Soc. 2006, 128, 3494. (9) Pan, L.; Olson, D. H.; Ciemnolonski, L. R.; Heddy, R.; Li, J. Angew. Chem., Int. Ed. 2006, 45, 616. (10) Yang, Q.; Zhong, C. ChemPhysChem 2006, 7, 1417. (11) Yang, Q.; Xue, C.; Zhong, C.; Chen, J.-F. AIChE J. 2007, 53, 2832. (12) Babarao, R.; Hu, Z.; Jiang, J.; Chempath, S.; Sandler, S. I. Langmuir 2007, 23, 659. (13) Keskin, S.; Sholl, D. S. J. Phys. Chem. C 2007, 111, 14055. (14) Wang, S.; Yang, Q.; Zhong, C. Sep. Purif. Technol. 2008, 60, 30. (15) Hayashi, H.; Coˆte´, A. P.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Nat. Mater. 2007, 6, 501. (16) Banerjee, R.; Phan, A.; Wang, B.; Knobler, C.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Science 2008, 319, 939. (17) 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. (18) Baburin, I. A.; Leoni, S.; Seifert, G. J. Phys. Chem. B 2008, 112, 9437. (19) Garberoglio, G. Chem. Phys. Lett. 2009, 467, 270. (20) Liu, D.-H.; Zheng, C.-C.; Yang, Q.-Y.; Zhong, C.-L. J. Phys. Chem. C 2009, 113, 5004. (21) Potoff, J. J.; Siepmann, J. I. AIChE J. 2001, 47, 1676. (22) Liu, J.; Rankin, R. B.; Johnson, J. K. Mol. Simul. 2009, 35, 60. (23) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (24) Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III J. Phys. Chem. 1990, 94, 8897. (25) Buch, V. J. Chem. Phys. 1994, 100, 7610. (26) Jiang, S. Y.; Gubbins, K. E.; Zollweg, J. A. Mol. Phys. 1993, 80, 103. (27) Murthy, C. S.; Singer, K.; Klien, M. L.; McDonald, I. R. Mol. Phys. 1980, 41, 1387. (28) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2005, 109, 15760.

Rankin et al. (29) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (30) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: London, 1996. (31) Johnson, J. K.; Zollweg, J. A.; Gubbins, K. E. Mol. Phys. 1993, 78, 591. (32) Span, R.; Wagner, W. J. Phys. Chem. Ref. Data 1996, 25, 1509. (33) Span, R.; Lemmon, E. W.; Jacobsen, R. T.; Wagner, W.; Yokozeki, A. J. Phys. Chem. Ref. Data 2000, 29, 1361. (34) Liu, J.; Culp, J. T.; Natesakhawat, S.; Bockrath, B. C.; Zande, B.; Sankar, S. G.; Garberoglio, G.; Johnson, J. K. J. Phys. Chem. C 2007, 111, 9305. (35) Skoulidas, A. I.; Bowen, T. C.; Doelling, C. M.; Falconer, J. L.; Noble, R. D.; Sholl, D. S. J. Membr. Sci. 2003, 227, 123. (36) Chen, H. B.; Sholl, D. S. J. Membr. Sci. 2006, 169, 152. (37) Sanborn, M. J.; Snurr, R. Q. Sep. Purif. Technol. 2000, 20, 1. (38) Sholl, D. S. Acc. Chem. Res. 2006, 39, 403. (39) Kresse, G.; Hafner, J. Phys. ReV. B 1993, 47, 558. (40) Kresse, G.; Hafner, J. Phys. ReV. B 1994, 49, 14251. (41) Kresse, G.; Furthmu¨ller, J. Comput. Mater. Sci. 1996, 6, 15. (42) Kresse, G.; Furthmu¨ller, J. Phys. ReV. B 1996, 54, 11169. (43) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. ReV. B 1992, 46, 6671. (44) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. ReV. B 1993, 48, 4978. (45) Blo¨chl, P. E. Phys. ReV. B. 1994, 50, 17953. (46) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758. (47) Monkhorst, H. J.; Pack, J. D. Phys. ReV. B 1976, 13, 5188. (48) Henkelman, G.; Arnaldsson, A.; Jonsson, H. Comput. Mater. Sci. 2006, 36, 354. (49) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; et al. In Gaussian 03; Gaussian, Inc.: Pittsburgh, 2003. (50) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (51) Maciel, G. S.; Garcia, E. Chem. Phys. Lett. 2005, 409, 29. (52) Martin, F.; Zipse, H. J. Comput. Chem. 2005, 26, 97. (53) Heinz, H.; Suter, U. W. J. Phys. Chem. B 2004, 108, 18341. (54) Myers, A. L.; Monson, P. A. Langmuir 2002, 18, 10261. (55) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813. (56) Panagiotopoulos, A. Z.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527. (57) http://towhee.sourceforge.net. (58) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508. (59) NIST. NIST Standand Reference Database Number 69; June 2005 release ed., 2005. (60) Zhu, W.; Kapteijn, F.; Moulijn, J. A.; Jansen, J. C. Phys. Chem. Chem. Phys. 2000, 2, 1773.

JP903735M