Article pubs.acs.org/JPCC
Influence of Zeolite Topology on CO2/N2 Separation Behavior: ForceField Simulations Using a DFT-Derived Charge Model Michael Fischer and Robert G. Bell* Department of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom S Supporting Information *
ABSTRACT: Force-field-based grand-canonical Monte Carlo (GCMC) simulations of carbon dioxide and nitrogen adsorption have been carried out for 18 different zeolite topologies, assuming a purely siliceous composition. The electrostatic potential in the pores of the zeolite framework is accurately modeled using point charges derived from density-functional theory calculations. These parameters are complemented by newly derived, empirical Lennard-Jones parameters, permitting an unambiguous decomposition of the total interaction energy into an electrostatic and a dispersive contribution. For each of the systems, such a decomposition is carried out, based on GCMC simulations of single-component adsorption. In the next step, mixture adsorption isotherms are computed, permitting the prediction of working capacities and adsorption selectivities, important quantities that determine the suitability of an adsorbent. In the final part, preferred carbon dioxide adsorption sites are calculated and analyzed for a subset of seven zeolite topologies. It is found that particular structural features, such as double-crankshaft chains and double eight-ring units, are common to those systems that exhibit the highest selectivities. While the majority of the systems have not been synthesized in their purely siliceous form, the conclusions can easily be generalized to other materials with zeolite topologies, such as (silico)aluminophosphates.
1. INTRODUCTION As rising atmospheric carbon dioxide levels have led to increasing environmental concerns, considerable research efforts are directed toward the development of technologies that prevent the release of carbon dioxide into the atmosphere. A promising route is the removal of carbon dioxide from exhaust gases (or other anthropogenic sources), with a subsequent sequestration of CO2 in suitable geological formations.1,2 Classically, carbon dioxide removal is achieved by amine scrubbing, that is, absorption of CO2 in solvents such as methanolamine.3 However, as absorption in amines is often considered to be too energy intensive, alternative approaches using adsorption-based or membrane-based technologies have been proposed.4 Because of the low partial pressure of carbon dioxide, membranes do not appear to be very suitable for the task of removing CO2 from exhaust gases, although they hold more promise for separation processes that operate at higher pressures.5,6 Therefore, adsorption-based separation appears to be the most promising route. The development of new adsorbent materials, such as metal−organic frameworks (MOFs) and porous polymers, has stimulated extensive research to assess the potential of these materials for CO2 removal. A variety of reviews dealing with carbon dioxide capture using solid adsorbents have been published.5,7−10 A typical exhaust gas from a coal-fired plant contains 10− 15% carbon dioxide (in molar units).11,12 The main constituent is nitrogen. Therefore, the model system for the removal of carbon dioxide from exhaust gases is the separation of binary © 2012 American Chemical Society
CO2/N2 mixtures at low CO2 partial pressures. Other components, especially water, may pose a problem for the separation process, as they tend to be strongly adsorbed, but they are normally neglected in laboratory-scale studies. The exhaust gas typically has a total pressure of approximately 1 bar. The considerable compression needed to implement a separation by conventional pressure swing adsorption (PSA) leads to an additional energy penalty. Therefore, vacuum swing adsorption, which uses an adsorption pressure of 1 bar, and a desorption pressure considerably below 1 bar, appears to be the preferred technology.8,11−13 Although the initial temperature of the exhaust gas is typically around 100 °C, it is usually assumed that the separation is performed at around room temperature. Siliceous zeolites, cation-exchanged zeolites, and related materials with zeolite-like structures have been characterized with regard to their potential for carbon dioxide removal, usually with more emphasis on CO2/CH4 separation (upgrade of natural gas/landfill gas) than CO2/N2 separation. Studies of all-silica zeolites include Silicalite (MFI framework type),14,15 clathrasil DD3R (DDR),16,17 and ITQ-29 (LTA).18 A much larger number of investigations have been reported for cationexchanged zeolites, including LTA zeolites with different Si/Al ratios or different cations,18,19 zeolite X and Y with various cation types (FAU),20,21 chabazites (CHA),22,23 mordenites Received: October 9, 2012 Revised: November 19, 2012 Published: November 20, 2012 26449
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
(MOR),24 and zeolite RHO.25 While these systems are all synthetic materials, other studies have assessed the potential of naturally occurring zeolites.26,27 Related materials exhibiting zeolite topologies, but incorporating other framework atoms, have also been considered, among them aluminophosphates28 and silicoaluminophosphates.29 Interestingly, a recent survey of experimental data implied that, while MOFs have beneficial properties for CO2/CH4 separation, cation-exchanged zeolites seem to be more promising for VSA CO2/N2 separation.8 Molecular simulation methods, principally grand-canonical Monte Carlo (GCMC) simulations and molecular dynamics (MD) calculations, have made an important contribution to the understanding of gas adsorption, diffusion, and separation in porous materials.30−32 A crucial point in these simulations is the use of a sufficiently accurate set of force-field parameters to describe solid−fluid and fluid−fluid interactions. While reference potentials for fluid−fluid interactions can be obtained relatively straightforwardly, the derivation of solid−fluid interaction parameters typically invokes a much higher degree of empiricism. The majority of existing force fields are based on an empirical fit to experimental adsorption isotherms, such as the widely used parameter set of Calero and co-workers.33,34 However, some efforts have been made to derive the parameters from electronic structure calculations.35,36 The D2FF force field, recently developed by Sholl and co-workers, based on a large set of dispersion-corrected density-functional theory (DFT) calculations, appears to be particularly promising in this respect.36 Initial studies of carbon dioxide adsorption in zeolites typically focused on one system, or a small group of structurally different materials.15,33,35,37,38 More recently, an increasing number of researchers has concentrated on the “screening” of very large sets of zeolites and zeolite-like structures. For example, Krishna and van Baten evaluated the potential of 27 all-silica zeolites (and several MOFs) for selective CO2 removal using GCMC and MD simulations for a wide range of conditions.39 Focusing on MOFs, Sholl and co-workers assessed the potential of ∼500 different systems, starting with a computation of the Henry constants for CO2 and N2 for all systems, and subsequently performing calculations of the full adsorption isotherms for the most promising systems.40 A much larger number of real and hypothetical zeolite structures, as well as zeolitic imidazolate frameworks (ZIFs), was screened by Smit and co-workers.41 Despite these screening studies covering a variety of structures, relatively little is known concerning the structural features of zeolites that are beneficial to achieve high adsorption selectivities. Krishna and van Baten provided an extensive analysis of simulation snapshots; however, their main focus was placed on the impact of structural properties (such as window size) on the diffusion properties.39,42,43 On the other hand, there have been only a few structural characterization studies of the preferential adsorption sites, although some neutron diffraction studies23 and DFT investigations44−46 have provided interesting insights. In this study, the CO2/N2 separation properties of 18 allsilica zeolites are investigated using GCMC simulations. The choice of all-silica systems is not based on their particular relevance for applications (cation-exchanged zeolites and, in some cases, aluminophosphates usually exhibit more favorable separation properties), but was made to compare a wide range of materials with different topologies, but identical pore wall composition. While only a few of these systems are actually
available in their purely siliceous form, the choice of topology was motivated by (a) the “relevance” of the structures, preferring common topologies over more exotic ones, and (b) the availability of the topologies in an aluminosilicate form. Given the number of structures typically included in recent screening studies,40,41 the current work cannot be considered in this category. Instead, the focus lies on a detailed case-by-case evaluation of the factors governing the CO2/N2 separation behavior for all 18 systems. This includes an investigation of the preferential adsorption sites to understand the microscopic origin of the observed properties. Another distinction of the present study in comparison to earlier works is the use of a semiempirical force field, which combines DFT-based point charges that have been specifically derived for each zeolite with new empirical Lennard-Jones parameters to model dispersive interactions. The use of point charges that are designed to accurately reproduce the DFT electrostatic potential permits an unambiguous decomposition of the interaction energy into an electrostatic and a dispersive part. While previous empirical force fields also include both interaction types, the relative contributions of the two are often represented in an ambiguous way, both charges and LennardJones parameters being derived by reference to the same experimental data. The remainder of this Article is organized as follows: After the technical details of the computations involved are summarized, the construction of the semiempirical force field is described. The agreement between simulations using these parameters and experimental data is discussed critically. The next section reports the results of GCMC simulations of singlecomponent adsorption in the 18 systems. Special attention is paid to the analysis of the different contributions to the total interaction energy. Results for mixture adsorption are discussed subsequently, focusing on both the adsorption selectivity and the working capacity. In the final part, the preferred carbon dioxide adsorption sites are analyzed for a subset of zeolite topologies of special interest, with the aim of identifying structural features that are beneficial for a high affinity toward CO2.
2. METHODS AND MODELS 2.1. Grand-Canonical Monte Carlo Simulations. Grandcanonical Monte Carlo simulations of carbon dioxide and nitrogen single-component and mixture adsorption were carried out using the Multipurpose Simulation Code MuSiC.47 Technical details of the simulations are provided in the Supporting Information. Long-range dispersive interactions and short-range repulsion were modeled using Lennard-Jones potentials, whereas an atom-centered point-charge representation was employed for electrostatic interactions. To derive a set of parameters that reproduces the experimental results most accurately, initial simulations were performed for a subset of four all-silica zeolites for which experimental data are available. In these calculations, the temperature and pressure ranges were matched to the experimental conditions, and the details are reported in section 3.2. Having established a suitable parameter set, further simulations were carried out for 18 zeolite structures: BEA, CHA, DDR, EDI, ERI, FAU, FER, GIS, GME, ISV, LTA, LTL, MER, MFI, MOR, MTT, MWW, and RHO. Single-component CO2 and N2 adsorption isotherms were computed at room temperature (T = 298 K). The results from these calculations 26450
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
optimization of the snapshot was carried out, using the Forcite module of Accelrys “Materials Studio”,49 and force-field parameters identical to those in the GCMC simulations. The nonequivalent minima were identified from the snapshot, and a second geometry optimization with only one CO2 molecule per simulation box was performed. For each minimum, the crystallographic coordinates, as well as the interaction energy Utot, were determined. 2.2. Models of Fluid Molecules and Adsorbents. The models of the fluid molecules used in the GCMC simulations were taken from the TraPPE force field, which is designed to reproduce fluid properties, especially vapor−liquid equilibria.50 As these models are widely used in the literature, their detailed description is provided in the Supporting Information. The structural data of the siliceous zeolites were taken from the IZA database.51 This database provides optimized geometries of all recognized zeolite frameworks, assuming an allsilica composition, the optimizations being performed using a distance least-squares refinement. These idealized structures were used to ensure that all frameworks are treated on an equal footing, using data that are easily and freely available. The structures are very similar to those that would be obtained from a force field optimization, for instance, using the program GULP,52 although they may differ slightly from experimental structures of all-silica polymorphs, where known. We note that a recent study has shown that the adsorption properties are relatively insensitive toward small changes in the structure, for example, when comparing IZA database structures to experimental structures, whereas the diffusion properties are affected to a much larger extent.53 Supercells of sufficient size (typical box length of 35−40 Å) were used in all calculations. Where necessary, inaccessible voids in the structures were blocked using noninteracting spheres. The zeolite structures were assumed to be rigid. While zeolite flexibility has been found to play a significant role on the CO2 adsorption properties in some systems, such as cation-exchanged zeolite RHO,25 a detailed investigation of this effect is beyond the scope of this study. This approximation is also in line with the assumptions made in earlier modeling studies, where framework flexibility was found to have little impact on the adsorption properties.54,55 In contrast to this, it is well established that framework deformations upon adsorption are a very interesting feature of many MOFs.56−58 In the GCMC simulations, electrostatic interactions were modeled using point charges placed on the atomic sites of the zeolites. The point charges were obtained from DFT calculations, described below, using the REPEAT method.59 Before calculating charges for all 18 systems, the dependence of the REPEAT charges on different computational parameters was assessed for a simple model system, as will be described in more detail in section 3.1. Dispersive interactions between the zeolite framework and the fluid molecules were represented by Lennard-Jones potentials. It was assumed that only the zeolite oxygen atoms contribute to this interaction. The attribution of all dispersive interactions to the oxygen atoms is a very common approach, used in many empirically derived force fields for zeolites.30 Some recently developed parameter sets, such as the DFTderived D2FF force field, take into account dispersive interactions with both the oxygen and the silicon atoms.36 However, such a distinction requires some means to unambiguously decompose the dispersive interactions into contributions from the two different atom types. While this is
were primarily used for the analysis of the adsorption energetics. In particular, the isosteric heat of adsorption qst was calculated as: qst = −⟨Utot /N ⟩ + RT
(1)
Here, ⟨Utot/N⟩ is the average interaction energy per adsorbed molecule (in kJ mol−1), and R is the gas constant. The total interaction energy Utot was decomposed into the different contributions: Utot = UvdW,sf + Ues,sf + Uff
(2)
UvdW,sf and Ues,sf correspond to the dispersive and electrostatic contribution to solid−fluid interactions, whereas Uff includes all fluid−fluid interactions. As these interactions are much weaker than the solid−fluid interactions, a further decomposition into a dispersive and an electrostatic part was not performed. GCMC simulations of CO2/N2 mixture adsorption were carried out for a pressure range from 0.1 to 1 bar, T = 298 K, and two different compositions of the bulk phase, an equimolar CO2/N2 mixture, and a mixture with a CO2/N2 ratio of 1:9. The latter case corresponds to the idealized scenario considered by Bae and Snurr for flue gas separation using VSA.8,48 The calculations were not extended to higher pressures for two reasons: First, one focus of this study is the analysis of solid− fluid interactions, which are dominant at low pressures, whereas fluid−fluid interactions increase in importance with increasing pressure. Second, vacuum swing adsorption (VSA), which uses pressures up to 1 bar only, currently appears to be the economically most viable option for CO2 removal from flue gases.8,11−13 From the binary mixture isotherms, the theoretical working capacity was calculated from the amount of CO2 adsorbed from the binary mixture, assuming an adsorption pressure of 1 bar, and a desorption pressure of 0.1 bar: Δn(CO2 ) = n(CO2 )p = 1bar − n(CO2 )p = 0.1bar
(3)
The ideal adsorption selectivity SCO2/N2 was calculated by dividing the molar ratio in the adsorbed phase by the molar ratio in the gas phase, which corresponds to the ratio of partial pressures. SCO2/N2 =
n(CO2 )/n(N2) p(CO2 )/p(N2)
(4)
In addition to these two quantities, the sorbent selection parameter Ssp was also calculated from the results for the 1:9 mixture. This parameter, which was defined by Bae and Snurr based on earlier work,8 combines the adsorption selectivity at both the adsorption and the desorption pressure with the working capacities for both mixture components. Ssp is a particularly useful criterion to compare the separation performance among different groups of materials. Further information on the calculation of Ssp, as well as the values calculated for all 18 zeolite topologies, are reported in the Supporting Information (Table S20). The trends observed for this quantity, in terms of the most promising topologies, agree with the conclusions derived from the working capacity and the adsorption selectivity. Therefore, Ssp is not discussed separately in the main body of this Article. Finally, an analysis of the carbon dioxide adsorption sites was performed. A simulation snapshot obtained at T = 298 K, p = 1 bar was used as a starting point of this analysis. To relocate the carbon dioxide molecules to local minima, a geometry 26451
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
moderate cutoff of 400 Ry was used in all future CP2K calculations, and no supercells were employed. However, as shown in Table 1, the charges were found to change
possible when exhaustive DFT data are available (as is the case for the D2FF force field), it is not a sensible approach in an empirical fitting procedure, as employed here. The details of the fitting procedure are described in section 3.2. Best agreement between simulation and experiment was obtained with the following zeolite oxygen parameters: σ(Ozeo) = 3.08 Å, ε(Ozeo)/kB = 107 K. The solid−fluid interaction potentials were calculated from these parameters, together with the TraPPE parameters for the fluid molecules, using Lorentz− Berthelot combination rules. 2.3. Density-Functional Theory Calculations. The majority of the density-functional theory calculations were carried out with the “Quickstep” part of the CP2K program, a DFT code that uses a mixed Gaussian and plane-wave approach.60 For the core electrons, relativistic Goedecker− Teter−Hutter pseudopotentials were employed.61 Some computational parameters, described fully in section 3.1, were varied in an initial study, to test their influence on the resulting REPEAT charges. The following setup was found to give charges of reasonable accuracy at moderate computational expense and was used in all following calculations: PBE exchange-correlation functional,62 plane-wave cutoff of 400 Ry, DZVP-SR (double-ζ plus polarization) basis set from the series of highly contracted MOLOPT basis sets developed by VandeVondele and Hutter.63 In addition to the CP2K calculations, a few reference calculations were carried out using the DMol3 all-electron code.64,65 These calculations used the PBE exchangecorrelation functional, numerical basis sets of different size (see section 3.1), and an orbital cutoff of 6.0 Å. To perform a calculation of the REPEAT charges, the DFT electrostatic potential was exported on a sufficiently fine grid (typically 0.1−0.2 Å). This was then used as input for the REPEAT code as provided by Woo and co-workers.59 No additional restraints were invoked in these calculations.
Table 1. REPEAT Charges of the Silicon Atom in All-Silica Sodalite Calculated with Different DFT Codes and Different Basis Setsa CP2K
DMol3
basis set
q(Si)/e
SZV SVZ-SR DZVP DZVP-SR TZVP TZV2P TZV2PX SN DN DNP TNP
2.098 2.091 1.544 1.487 1.477 1.436 1.429 2.087 2.376 1.609 1.397
a
Because there is only one symmetry-unique oxygen position, the charge of the oxygen atom is −1/2q(Si).
significantly with increasing size of the basis set, ranging from q(Si) = +2.1 e for a single-ζ (SZV) basis set to q(Si) = +1.43 e for the largest triple-ζ (TZV2PX) basis set. Interestingly, a “short-range” double-ζ (DZVP-SR) basis set designed for use with periodic solids resulted in a value of q(Si) = +1.49 e, thus being in very good agreement with the results for the much larger TZV2PX basis. Because calculations using the diffuse triple-ζ basis sets become computationally very demanding for larger systems, the DZVP-SR basis set was employed in all further calculations. With these computational settings, the use of different exchange-correlation functionals (LDA, BP, BLYP, PBE) caused a relatively insignificant variation of the charges. Therefore, the PBE functional was used in all CP2K calculations. In the additional all-electron calculations performed with the DMol3 codes, the basis set size was varied from a single-ζ (SN) to a triple-ζ (TNP) basis set, using the PBE functional and a 4 × 4 × 4 k-mesh. A trend similar to that in the CP2K calculations is observed, with decreasing charge on increasing basis set size (Table 1). The quantitative agreement between the results from the two different codes is very good. The results obtained using the largest basis sets also agree well with the previously reported value for q(Si) in sodalite from the original work introducing the REPEAT method.59 Taken together, the systematic variation of different input parameters in the DFT calculations reveals a high sensitivity of the REPEAT charges to the basis set size, whereas the sensitivity to most other computational parameters seems to be much less pronounced. For the specific case of the CP2K code, the DZVP-SR basis set delivers sufficient accuracy, permitting computations for relatively large systems at a very modest computational expense. Analogous CP2K calculations were carried out for the 18 zeolite structures considered in this study, and the REPEAT charges were derived for each structure. These charges are given in the Supporting Information. The charges for the silicon atoms typically range from +1.2 e to +1.7 e, and the charges for the oxygen atoms fall mostly in a range between −0.5 e and −0.9 e, although a few outliers are observed. It is noteworthy that the charges are roughly 20−40% smaller than
3. RESULTS AND DISCUSSION 3.1. Calculation of REPEAT Charges. The REPEAT method is a charge-derivation scheme designed to reproduce the electrostatic potential from an electronic structure calculation, in a fashion similar to the more traditional ESP method.66 However, while the traditional ESP method performs very poorly for periodic structures, the REPEAT method allows the treatment of these systems by introduction of a modified error functional.59 The electrostatic potential obtained from a DFT calculation is used as an input to the REPEAT code, which delivers a set of atom-centered point charges. The REPEAT method has been shown to deliver physically reasonable charges for a variety of crystalline systems. For example, GCMC simulations of adsorption in MOFs that employ REPEAT charges agree very well with calculations that directly incorporate the DFT electrostatic potential.67 To test the sensitivity of the resulting REPEAT charges to the DFT setup, calculations for all-silica sodalite, a nonporous zeolite with a small unit cell, were carried out, using two different DFT codes. DFT calculations were performed with the pseudopotential code CP2K, and with the all-electron code DMol3, varying several input parameters. In the CP2K calculations, the plane-wave cutoff, the size of the unit cell, the basis set size, and the exchange-correlation functional were varied. The variation of the plane-wave cutoff (in a range from 200 to 800 Ry) and the use of a 2 × 2 × 2 supercell had virtually no influence on the resulting charges. Therefore, a 26452
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
Figure 1. (a) CO2 adsorption isotherms at room temperature up to pressures of 1 bar. (b) CO2 adsorption isotherms for higher pressures, and different temperatures. LTA: blue ⬡ = experiment,18 blue ◆ = simulation. MFI: orange lines = Langmuir fit to experiment (solid, T = 303 K; dashed, T = 408 K),14 orange symbols = simulation (◆, T = 303 K; ▼, T = 408 K). DDR: green lines = Langmuir fit to experiment (solid, T = 298 K; dashed, T = 348 K),16 green symbols = simulation (◆, T = 298 K; ▼, T = 348 K).
the charges used in an empirical force field that has been frequently used for the modeling of CO2 adsorption (q(Si) = +2.05 e).33 On the other hand, they lie in a range similar to that of the charges that were used in a modeling study of liquid water adsorption in Silicalite (q(Si) = +1.40 e), where the outcome of the simulations was found to strongly depend on the magnitude of the point charges.68 The relatively broad distribution of the REPEAT charge values indicates that a replacement of these DFT-derived charges by generic values of q(Si) and q(O) would lead to a considerable loss of accuracy in the resulting electrostatic potential. Therefore, the use of REPEAT charges, rather than generic charge values, should deliver a considerably more accurate representation of electrostatic interactions in simulations of adsorption. 3.2. Derivation of Lennard-Jones Parameters. CO2 Adsorption. As described earlier, Lennard-Jones parameters were assigned only to the oxygen atoms of the framework. Because the solid−fluid parameters are calculated from the atomic parameters using Lorentz−Berthelot mixing rules, only two parameters were derived, the van der Waals radius of the zeolite oxygen atoms, σ(Ozeo), and the potential well depth ε(Ozeo). The parameter adjustment was performed, carrying out GCMC simulations for three all-silica zeolites for which experimental CO2 adsorption data are available (LTA, MFI, and DDR). For Si-LTA, only a room temperature measurement has been reported.18 For Silicalite (Si-MFI), experimental data for temperatures ranging from T = 303 to 408 K were considered.14 Finally, data for temperatures from T = 273 to 348 K were taken into account for all-silica DDR.16 The calculations were carried out up to pressures of 6 bar, as experimental data for higher pressures are available for DDR only. In the cases of MFI and DDR, the experimental data were recalculated from Langmuir fits provided in the original articles, which agreed very well with the actual measurements. This procedure permits a direct calculation of the relative error between simulation and experiment for each pressure point. The initial set of framework parameters was taken from the work of Dubbeldam et al.69 The parameters were then refined in several steps until a good match between experimental and simulated isotherms was attained. For MFI and DDR, the parameter refinement aimed at a minimization of the relative error between the simulated isotherms and the Langmuir fits
obtained from experiment, whereas a visual inspection of the match between calculation and experiment was employed for LTA. The refinement took advantage of the well-established relationships between σ(Ozeo) and ε(Ozeo) and the gas uptakes in different pressure regimes: While the radius of the framework atoms σ(Ozeo) determines the uptake at high pressures, the shape of the isotherm at low pressures depends on the potential well depth ε(Ozeo) once σ(Ozeo) has been fixed.69,70 Figure 1a visualizes the CO2 adorption isotherms for zeolites LTA, MFI, and DDR at temperatures around room temperature for pressures up to 1 bar, together with experimental data. Figure 1b includes pressures up to 6 bar, and selected isotherms for higher temperatures. For DDR and MFI, the isotherms for all temperatures, as well as the relative errors, are displayed in the Supporting Information (Figures S1−S4). For LTA, the simulated adsorption isotherm for T = 303 K exhibits excellent agreement with the experimental data throughout the whole pressure range, with only marginal deviations at higher pressures. The other two zeolites exhibit opposite trends at room temperature, with a slight overestimation of the carbon dioxide uptake of MFI throughout the whole pressure range, and a more noticeable but still minor underestimation of the uptake of DDR at pressures above 2 bar. These two results suggest that little further improvement of the parameters is possible: Any change that improves the description of one system is likely to make matters worse for the other. The overall agreement is very good, with the relative error between simulation and experiment exceeding 10% (Figures S2/S4) only at one single pressure point (absolute root-mean-square error MER > GME ≈ CHA. The amount of CO2 adsorbed in these systems is considerably lower than the carbon dioxide uptakes of many MOFs, which have a lower density, and in many cases possess localized sites of strong interaction (such as unsaturated metal sites). Here, very high CO2 uptakes have been reported for room temperature and a pressure of 1 bar. For example, MgMOF-74 takes up more than 8.5 mmol g−1 under these conditions,72 and uptakes around 5 mmol g−1 are attained in the commercially available Cu3(btc)2 and other MOFs with unsaturated copper sites, such as PCN-11 and PCN-16.73 Values near 5 mmol g−1 are also reached in Na-exchanged zeolite LTA.18 While these uptakes are not quite matched by the all-silica systems studied here, the predicted CO2 uptake of GIS is quite close to the capacities reported for two adsorbents employed in VSA pilot-scale studies, zeolite 13X (4.4 mmol g−1
a
For carbon dioxide (nitrogen), the simulated values correspond to a pressure of 0.05 bar (0.5 bar) and a temperature of 298 K. Where coverage-dependent data are available, the closest matching experimental data point is reported. Experimental data are from refs 14, 16−18, 36, 71.
adsorption is underestimated by 2−3 kJ mol−1, an excellent agreement between experiment and simulation is observed. Moreover, there is no systematic deviation toward too high or too low values. Taking together all of these results, we confidently expect that the Lennard-Jones parameters developed in this work, in conjunction with DFT-based REPEAT charges, will allow for an accurate prediction of the carbon dioxide uptake of all-silica zeolites for a variety of pressure and temperature conditions. Because the Lennard-Jones parameters are associated with the oxygen atoms only, they can also be employed for other systems with zeolite topology and homogeneous pore wall composition, such as aluminophosphates. To investigate this possible extension, ALPO-17 (ERI topology) was studied as an example, using REPEAT charges that were specifically calculated for the aluminophosphate structure. The simulated carbon dioxide adsorption isotherm is shown in the Supporting Information, together with experimental data (Figure S6).28 Very good agreement is observed, emphasizing the potential transferability of the present approach to other zeolite-like materials. 3.3. N2 Adsorption. Of the four all-silica zeolites used for the parameter derivation in section 3.2, experimental N2 adsorption data for temperatures around room temperature are available for MFI and DDR only.17,71 GCMC simulations of nitrogen adsorption were carried out for MFI for T = 305 K and T = 342 K, and for DDR for three different temperatures. The zeolite oxygen parameters derived in the previous section were used to calculate the solid−fluid interaction parameters from Lorentz−Berthelot combination rules. The simulation results are displayed in Figure 2, together with experimental data. Excellent agreement is observed for MFI for both temperatures. For DDR, there is a tendency to slightly 26454
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
Figure 3. CO2 and N2 adsorption capacity of the 18 systems at a pressure of 1 bar, T = 298 K.
sites, qst usually ranges between 25 and 40 kJ mol−1.72,73 Even higher values (at low coverage) are reached in a MOF with amino-decorated pores.75 While a high isosteric heat of CO2 adsorption is beneficial to reach a high selectivity, there is a trade-off with regard to applications, because the regeneration of the adsorbent becomes more energy-consuming with increasing qst. In all zeolite frameworks studied, the dispersive solid−fluid interactions constitute the major part of the total interaction energy. Electrostatic interactions usually account for about 10% of the total solid−fluid interaction energy, indicating that these interactions play a significant role in the adsorption of carbon dioxide in siliceous zeolites, which are relatively nonpolar adsorbents. Looking at the different structures, there is a considerable variation of the relative contribution of electrostatics, ranging from less than 5% in EDI, FER, and ISV to more than 15% in GIS, MER, MOR, and RHO. These quite pronounced differences highlight the dependence of electrostatic interactions on the pore topology, which further points to the necessity of using system-specific point charges to accurately model these contributions. Finally, the contribution of fluid−fluid interactions varies between 0.6 and 2.1 kJ mol−1. Unsurprisingly, this contribution is highest for those systems where the carbon dioxide adsorption isotherms approach saturation (GIS, MER), because the CO2 molecules are already quite closely packed in these structures at a pressure of 1 bar. The evolution of the isosteric heat of carbon dioxide adsorption on increasing loading is displayed in the Supporting Information. In the majority of the cases, the isosteric heat of adsorption increases with increasing amount of CO2 adsorbed, typically by approximately 0.5−2.0 kJ mol−1.This behavior can be straightforwardly explained with the increasing contribution of fluid−fluid interactions. There are, however, a few cases where the opposite trend is observed: In GME, MOR, and RHO, the isosteric heat of adsorption decreases by up to 2.8 kJ mol−1 with increasing CO2 loading. In the final part of this Article, an explanation for this decrease will be offered. The isosteric heats of nitrogen adsorption at 1 bar are summarized in Figure 5 (bottom). With the exception of FAU, all values of qst fall in a range between 12 and 18 kJ mol−1. Electrostatic interactions are practically negligible in the majority of the systems and are roughly one order of magnitude smaller than the corresponding values of Ues,sf observed for carbon dioxide. This difference in magnitude directly reflects the difference in quadrupole moment. Finally, because only
at 293 K),11 and a molecular sieve carbon material (4.2 mmol g−1 at 294 K).74 A closer inspection of the carbon dioxide adsorption isotherms, displayed in Figure 4, shows that GIS has nearly
Figure 4. CO2 adsorption isotherms of the five systems exhibiting carbon dioxide uptakes above 2 mmol g−1 at 1 bar and 298 K.
reached saturation at 1 bar, and that GME and MER also exhibit a visible decrease in the slope of the isotherm. In contrast, the isotherms of RHO and CHA show a more linear behavior, indicating that saturation will be attained at considerably higher pressures. The nitrogen uptakes of the 18 systems are roughly a factor of 10 lower, and there is no significant correlation between the CO2 and the N2 uptake. All N2 isotherms exhibit a linear shape. The isosteric heats of carbon dioxide adsorption at p = 1 bar are shown in Figure 5 (top), together with the decomposition of the average interaction energy into the different contributions. The highest isosteric heats are predicted for GIS and MER, with 34.4 and 31.3 kJ mol−1, respectively. Of the remaining systems, eight fall in a range between 26 and 30 kJ mol−1, five lie between 20 and 26 kJ mol−1, and three systems exhibit an isosteric heat of CO2 adsorption of less than 20 kJ mol−1. While these values lie in the same range as the isosteric heat of CO2 adsorption of a molecular sieve carbon (26.5 kJ mol−1),74 cation-exchanged zeolites and MOFs with accessible metal sites typically exhibit higher heats of adsorption. For NaLTA, qst = 25−48 kJ mol−1, depending on the amount of cations and the CO2 coverage.18 In MOFs with accessible metal 26455
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
Figure 5. Isosteric heats of CO2 and N2 adsorption, calculated at a pressure of 1 bar and 298 K. The isosteric heat of adsorption is decomposed into the (constant) RT term, fluid−fluid contributions to the total interaction energy Utot, and dispersive and electrostatic solid−fluid contributions to Utot.
Figure 6. Carbon dioxide working capacity of the 18 systems, assuming CO2/N2 mixtures of two different compositions, an adsorption pressure of 1 bar, and a desorption pressure of 0.1 bar. For clarity, the bar representing the working capacity of GIS is not shown to its full extent.
small amounts of N2 are adsorbed under these conditions, fluid−fluid interactions are virtually absent. 3.5. CO2/N2 Mixture Adsorption. The primary data of the CO2/N2 mixture adsorption isotherms at T = 298 K for all 18 systems are shown in the Supporting Information. The working capacities derived from these calculations for the equimolar mixture and the (more realistic) 1:9 mixture are compiled in Figure 6. For those systems that exhibit a nearly linear evolution of the carbon dioxide isotherm at the pressures considered, the working capacity drops by a factor of 4−5 when
going from the equimolar to the 1:9 composition. However, the drop is much less pronounced in the zeolites that show a sharp rise of the isotherm at low pressures, especially GIS, which exhibits the highest working capacities with 2.7 mmol g−1 (equimolar mixture) and 1.4 mmol g−1 (1:9 mixture). On the other hand, the working capacity of this system will be much lower for normal PSA conditions (i.e., an adsorption pressure of 5 bar or more, and a desorption pressure of 1 bar), because the CO2 adsorption isotherm flattens at higher pressures. This highlights that a given value of the working capacity is only 26456
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
Figure 7. CO2/N2 adsorption selectivity of the 18 systems, obtained for CO2/N2 mixtures of two different compositions at a total pressure of 1 bar. For clarity, the bars representing the selectivity of GIS are not shown to their full extent.
than all other zeolites, it is not sensible to establish a quantitative correlation. Nevertheless, it is obvious that a large difference in affinity is a necessary prerequisite to obtain a high CO2/N2 selectivity. As discussed above, the isosteric heats of nitrogen adsorption fall in a much narrower range than the values of qst(CO2). Therefore, systems with a high affinity for carbon dioxide hold particular promise with regard to reaching a high selectivity. To complement the macroscopic insights obtained at this stage with information at an atomistic level, the final section of this study will explore the structural origins of the affinity for carbon dioxide in selected zeolite topologies. 3.6. Comparison to Literature GCMC Studies. Because of the fact that no stable all-silica polymorphs have been reported for the best-performing systems in this study, a quantitative comparison with experimental adsorbents is not particularly meaningful. Nevertheless, a few benchmark values from a recent survey of experimental data8 should be mentioned to put the results into a broader context: In this compilation, which focused on MOFs, but also included some zeolites, the same conditions as in the present work were assumed, that is, CO2/N2 separation under VSA conditions. Zeolite 5A, zeolite 13X, and Ni-MOF-74 emerged as the materials with the highest working capacities, which amount to 2.4, 1.4, and 3.2 mmol g−1, respectively. CO2/N2 selectivities of 62, 86, and 41 were obtained for these systems. All three materials contain extra-framework cations or unsaturated metal sites that interact strongly with the carbon dioxide molecules. A comparison of the adsorption selectivities reported here with previous simulation studies reveals good agreement in many instances, but also some interesting discrepancies. An early investigation by Goj et al. reported Henry’s law selectivities for MFI and ISV.37 While the value for MFI is moderately higher than the value calculated here (26.9 as opposed to 17.7 at p = 0.1 bar), the selectivities obtained for ISV are in excellent agreement (7.6 as compared to 8.5). A more recent study by Liu and Smit delivered a selectivity of 22 at 1 bar for MFI, very close to the result of the present work of 20.9.77 For DDR, two previous modeling studies obtained S(CO2/N2) = 2615 and 3078 at 1 bar, values that are higher than the selectivity calculated here (20.8). The latter investigation also included a CO2/N2 selectivity value for LTA, which was considerably lower than the value of the present work (5.5 as compared to 8.7). The most exhaustive GCMC studies of CO2/N2 mixture separation in siliceous zeolites have been performed by Krishna and van Baten.39,43 A number of the
valid for the specific conditions considered, in particular the chosen adsorption and desorption pressures. While reporting a single value of the working capacity, as is done here and in previous studies,39 can serve as a useful “benchmark”, it is important to acknowledge that a change of conditions may qualitatively affect the behavior of an adsorbent. Therefore, as pointed out earlier by Bae and Snurr,8 conclusions concerning the suitability of a system should not be generalized prematurely. The adsorption selectivities S(CO2/N2) of all systems as a function of total pressure are shown in the Supporting Information. The values calculated for a pressure of 1 bar, and the two different mixtures, are summarized in Figure 7. Clearly, the exceptionally high selectivity of GIS, in the range of 100, is the most striking feature, and is more than twice the value calculated for the second-best system (MER). In total, there are four systems (GIS, MER, GME, MTT) that exhibit a selectivity S(CO2/N2) > 20 toward a 1:9 mixture. The selectivities toward an equimolar mixture are typically about 5−25% higher, an evolution that can be explained simply by the increased contribution of attractive fluid−fluid interactions, which are stronger for CO2 than for N2. Interestingly, GME is the only exception to this rule, the selectivity decreasing by about 20%. The evolution of the selectivity with total pressure can be inferred from the plots compiled in the Supporting Information. In the large majority of the cases, the selectivity increases on increasing pressure, by up to 45% in the case of GIS. Again, GME is the only case where the selectivity drops with pressure, although there are a few systems where the selectivity remains virtually constant (ERI, MOR, RHO). When comparing the data on the CO2/N2 selectivity with the isosteric heats of adsorption reported earlier, it is apparent that those systems that exhibit the highest values of qst(CO2) also show the highest selectivities (GIS, MER). Indeed, a correlation between the isosteric heat of CO2 adsorption and the selectivity was found in a survey of experimental results.8 A refined approach was proposed by Zhong and co-workers, based on extensive GCMC simulations, who found a correlation between the CO2/N2 selectivity and difference in isosteric heats of adsorption, Δqst = qst(CO2) − qst(N2).76 This observation is not altogether surprising, given that Δqst constitutes a measure of the difference in affinity toward the two species. Because of the limited number of systems treated in the present work, and the fact that one single system exhibits a much higher selectivity 26457
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
Table 3. Solid−Fluid Interaction Energies Usf for the Carbon Dioxide Adsorption Sites Determined from Geometry Optimizations (One CO2 Molecule per Supercell)a
topologies considered in the current study were included in these works. While the very compact presentation of the results in these articles makes a quantitative comparison difficult, some interesting qualitative differences are observable: • DDR was predicted to have a considerably higher CO2/ N2 higher selectivity than MFI (ratio of about 3:1); very similar selectivities were obtained in the present work. • The selectivity of ERI was found to be higher than the selectivity of MFI, whereas the opposite behavior is observed here. • At 1 bar, the following order of selectivities was computed: CHA > DDR > GME. Precisely the opposite trend is predicted in the present study. As the selectivities were calculated for comparable conditions, the remarkable differences highlight the dependence of the simulation results on the force-field parameters. Interestingly, the deviations between results of this study, and previous works (which used empirically derived force fields, in most cases the parameters of Calero and co-workers33) are not systematic: One parameter set may give a higher selectivity for one system, but a lower selectivity for another system. This leads to qualitative differences, for example, in the sequence of selectivities observed for a group of systems. Clearly, further careful studies are necessary, preferably in close conjunction with experimental measurements of single-component and binary mixture adsorption data, to unequivocally determine an optimized set of force-field parameters. In addition to the fully empirical parameter sets used in the literature, and the semiempirical force field proposed here, further improvements could be expected from parameters entirely based on electronic structure calculations, such as the D2FF force field.36 However, to date, this force field has not been extended to incorporate N2 as a guest molecule. 3.7. Carbon Dioxide Adsorption Sites − Results. To explain the high affinity of some zeolite structures toward carbon dioxide, the adsorption sites were analyzed for all 18 framework types considered. To highlight the most significant findings, a subset of seven systems is discussed in this section. In this subset, GIS and MER are included as the two systems exhibiting the highest affinity toward CO2, and consequently the highest CO2/N2 selectivity. Furthermore, GME, MOR, and RHO are considered as those zeolites that have high working capacities and/or adsorption selectivities, and moreover show a decrease of the isosteric heat of carbon dioxide adsorption with increasing loading, which is different from all other systems studied. Finally, the adsorption sites in FAU and BEA are analyzed. In these two systems, the CO2/N2 selectivities are among the lowest, and they are included to highlight the distinctly different nature of the carbon dioxide adsorption sites when compared to the other five systems. The interaction energies calculated for each geometry-optimized adsorption position are shown in Table 3. The assumption of local energy minima as static adsorption “sites” is of course a simplification of adsorption at finite temperature, where the whole accessible pore space is sampled. Nevertheless, it can be expected that most of the carbon dioxide molecules are adsorbed at the local minima and their direct surroundings at the temperature conditions considered. Therefore, these areas contribute most to the adsorption energetics, and an analysis of their structural environment can provide an explanation of the material’s adsorption behavior.
GIS MER
GME MOR
RHO FAU
site site site site site site site site site site site site site site site site site
A B A B C D E A B A B C A B A B C
Usf/kJ mol−1
Ues,sf/Usf
−34.5 −33.4 −34.1 −32.8 −32.5 −30.5 −30.0 −35.2 −27.8 −39.0 −36.1 −32.4 −35.0 −33.1 −19.7 −18.6 −17.9
0.22 0.26 0.23 0.29 0.22 0.21 0.00 0.22 0.16 0.29 0.18 0.21 0.29 0.23 0.11 0.19 0.06
a
The relative contribution of the electrostatic solid−fluid energy Ues,sf is also given. In the case of GME, MOR, and RHO, only the most favorable sites are included.
GIS. The accessible pore space of GIS consists of intersecting zigzag channels formed by eight-membered rings. The channel intersection, surrounded by four eight-ring windows, is designated as gis cavity.51 An analysis of the adsorption sites reveals two different sites that are very close in energy. The first site, shown in Figure 8a, is located inside the gis cavity,
Figure 8. (a) CO2 adsorption sites in the GIS structure. Three adjacent gis cavities are shown. For clarity, only the sites lying in one plane are displayed. The two overlapping, symmetry-equivalent molecules at site A are both shown. (b) CO2 adsorption sites in the MER structure. One mer cavity and two adjacent d8R rings are shown. Four symmetry-equivalent images of site B overlap with each other. Site E is outside the region of the structure displayed. Figures showing the full unit cell of both zeolites with all CO2 sites are given in the Supporting Information.
displaced from the cavity center toward either of the bowl-like units that constitute the double-crankshaft (dcc) chains (the site overlaps with its own image, meaning that only 50 percent of the positions can be occupied). Site B is located at the center of the eight-ring windows. While both positions are very close in energy, adsorption snapshots reveal that site A is practically exclusively occupied at a pressure of 1 bar. This can be explained by favorable molecule−molecule interactions: In site A, neighboring molecules are oriented parallel to each other. 26458
dx.doi.org/10.1021/jp3099768 | J. Phys. Chem. C 2012, 116, 26449−26463
The Journal of Physical Chemistry C
Article
The zigzag shape of the channels allows the molecules to orient in a way that attractive electrostatic CO2−CO2 interactions are maximized. On the other hand, neighboring molecules occupying site B would point toward each other with the oxygen atoms at a very close distance (below 3 Å), leading to strong electrostatic repulsion. It is therefore easy to understand why site A will dominate with increasing CO2 loading. An occupation of site A by 1/2 (which is the maximum due to the overlap with its own image) corresponds to 4 molecules per unit cell, a loading that is almost reached at 1 bar and 298 K (3.7 molecules). MER. The tetragonal MER framework contains ellipsoidal pores (mer cavities) that are connected by intersecting eightring channels, as well as double eight-rings (d8R building unit). Because of the increased structural complexity as compared to GIS, there is a larger number of local energy minima, visualized in Figure 8b. As all five minima are relatively close in energy, it can be expected that they all will play a relevant role in adsorption. Site A is located near the eight-ring windows of the d8R units, the carbon atom of the CO2 molecule being slightly displaced from the plane defined by the framework oxygen atoms in the eight-ring. Site B is located inside the d8R units, parallel to the windows. For geometric reasons, only one molecule can occupy either site A or B per d8R unit (2 molecules per unit cell). Site C is located centrally between two bowl-shaped units formed by three four-ring units of the dcc chains, much alike the first adsorption site observed in GIS. Only one-half of these sites can be occupied due to close contacts with their own images (4 molecules per unit cell). Sites D and E are located in the eight-ring windows connecting the mer cavities. In site D, the CO2 molecule assumes an orientation perpendicular to the window plane, leading to favorable electrostatic interactions. Conversely, attractive and repulsive electrostatic interactions cancel out for the case of site E, which is oriented parallel to the window. Again, there are close contacts with the other sites. Therefore, the loading of 5.6 molecules per unit cell at 1 bar and T = 298 K, which approaches saturation, is much lower than the sum over the adsorption sites. GME. The GME framework consists of straight 12-ring channels that are connected by smaller side cavities (gme cavities), with eight-ring windows forming the entrances to these cavities. Four local energy minima were obtained from the geometry optimizations. However, only the first two sites are included in Table 3, as the other two sites are energetically much less favorable. Site A is located inside the gme cavities, with the oxygen atoms of the CO2 molecule pointing toward the central four-ring of the dcc units that build up the cavity wall (Figure 9a). Site B is located at the eight-ring windows that connect these cavities to the main channels. Within one cavity, the images of both site A and site B point toward each other via their oxygen atoms, indicating that fluid−fluid interactions between adjacent molecules will tend to be repulsive. Moreover, no more than three molecules per gme cavity can be adsorbed at these sites for geometric reasons. Two other adsorption sites are located inside the 12-ring channels, the CO2 molecules orienting along the running direction of the channels. At less than −20 kJ mol−1, the interaction energies at these sites are considerably lower than the average solid−fluid interaction energy (Figure 5). Therefore, it can be concluded that the majority of the CO2 molecules under the conditions studied are adsorbed in sites A and B. Indeed, the observed loading at 1 bar and T = 298 K amounts to 6.2 molecules per
Figure 9. (a) CO2 adsorption sites in the GME structure. One gme cavity is shown. The two overlapping, symmetry-equivalent molecules at site B are both shown. (b) CO2 adsorption sites in the MOR structure. The side cavity and the adjacent part of the narrow eightring channel are shown. The two overlapping, symmetry-equivalent molecules at site B are both shown. Figures showing the full unit cell of both zeolites with all CO2 sites are given in the Supporting Information.
unit cell, which is in good correspondence with the six molecules per unit cell that sites A and B can contribute. MOR. The pore topology of the MOR framework includes 10-ring channels with side cavities formed by two eight-rings. The side cavities further connect to very narrow eight-ring channels. However, an analysis of the snapshots showed that these channels are not fully accessible to carbon dioxide molecules. Of the total of five local energy minima found, only three are included in Table 3. Similarly to GME, the other two sites, both located in the 10-ring channels, exhibit much lower interaction energies (
