Robust, Transferable, and Physically Motivated Force Fields for Gas

May 31, 2012 - Muthuramalingam Prakash , Navid Sakhavand , and Rouzbeh Shahsavari. The Journal of Physical Chemistry C 2013 117 (46), 24407-24416...
0 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCC

Robust, Transferable, and Physically Motivated Force Fields for Gas Adsorption in Functionalized Zeolitic Imidazolate Frameworks Jesse G. McDaniel and J. R. Schmidt* Theoretical Chemistry Institute and Department of Chemistry, University of WisconsinMadison, Madison, Wisconsin 53706, United States S Supporting Information *

ABSTRACT: We extend our existing methodology for generating physically motivated, tailored ab initio force fields via symmetry-adapted perturbation theory (SAPT). The revised approach naturally yields transferable atomic exchange, charge penetration, and dispersion parameters, facilitating the creation of versatile, optimized force fields; this approach is general, applicable to a wide array of potential applications. We then employ this approach to develop a force field, “ZIF FF”, which is tailored to accurately model CO2/N2 adsorption in zeolitic imidazolate frameworks (ZIFs). In conjunction with our previous “SYM” force field used to model adsorbate−adsorbate interactions, we compute adsorption isotherms for both CO2 and N2 in nine different ZIFs, yielding results that are in excellent accord with the corresponding experimental results. We find that ZIF FF accurately predicts isotherms for three different topologies of ZIFs (RHO, SOD, GME) and reproduces gas adsorption trends for varying functionalization across an isoreticular series of ZIFs of the GME topology. Because ZIF FF is free of empirical parameters, it presents the opportunity for computationally screening novel ZIFs that have not yet been synthesized and/or characterized.

■

INTRODUCTION Recent interest in carbon capture and sequestration has spurred an associated push toward the development of CO2 separation technologies. Here the challenge lies in the separation of dilute CO2 from a gas mixture, typically the flue gas originating from combustion at a coal-fired power plant. While monoethanolamine-based sorbents are currently the best developed technology,1 alternatives such as ionic liquids,2−7 membranes,8−10 and nanoporous sorbents11−17 are currently under investigation. In the latter case, metal−organic frameworks (MOFs) have been shown to exhibit excellent CO2 adsorption capacities and selectivity. Specifically, zeolitic imidazolate frameworks (ZIFs), composed of zinc or cobalt cations and imidazolate-based linkers, have shown particular promise in such applications due to their enhanced chemical and thermal stability compared to other MOFs.18,19 Like all MOFs, the functionality of the associated organic linker groups can be varied, yielding a fascinating array of ZIF structures with associated variations in CO2 adsorption, selectivity, and transport properties. Recent experimental developments in the high-throughput synthesis of ZIFs17,20 have led to the characterization of many new materials, some of which show excellent CO2 adsorption capacity and selectivity over N2 gas.16 However, such experimental synthesis and characterization of new materials is expensive and timeconsuming, and thus it is essential to have a physical understanding of the properties of these systems to direct future experimental efforts. Here, computer simulations such as © 2012 American Chemical Society

molecular dynamics (MD) and Monte Carlo (MC) are invaluable tools that can aid both in the understanding of gas adsorption properties of known materials21−34 and even in the screening of new materials not yet characterized by experiment.35−38 Yet all such simulations rely on the employed empirical force fields to accurately describe the solute−solute and adsorbate− framework interactions in the system. For ZIFs in particular, the ability of generic force fields such as the ubiquitous universal force field (UFF)39 to describe ZIF−gas interactions across different crystal topologies is questionable. For example, UFF parameters were found to dramatically overpredict gas adsorption in SOD (ZIF-8)23,31 and GME (ZIF-68, ZIF69)23,27 topologies, and these parameters had to be significantly scaled to obtain agreement with experimental isotherms. However, for RHO (ZIF-71, ZIF-25) topologies, unadjusted UFF parameters were found to satisfactorily predict CO2 adsorption isotherms.22 Such inconstancies should not be surprising in that such “generic” force fields (UFF, Drieding, etc.) have not been tailored for any specific application. Furthermore, even in cases where agreement with measured isotherms is seemingly good, we suggest that in many cases this agreement may arise from spurious error cancelation (vide infra). Received: April 19, 2012 Revised: May 29, 2012 Published: May 31, 2012 14031

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

Bii = (Bi + Bj ) × ((Bi Bj )/(Bi2 + Bj2 ))

In response, we previously reported a methodology which utilizes ab initio, symmetry-adapted perturbation theory (SAPT)40 calculations to construct physically motivated force fields for ZIF systems.21 We also demonstrated that these force fields lead to accurate predictions of CO2 adsorption isotherms. However, a significant limitation of that approach was parameter transferability, requiring the construction of a unique force field for each individual functionalized ZIF. In this paper, we extend our previous methodology to construct a set of transferable force-field parameters that can be utilized to study ZIFs of differing topologies and functionalities, and even different adsorbates, yielding a tailored ZIF force field, “ZIF FF”. In conjunction with our previously developed CO2 and N2 force fields,41 we use this force field to compute adsorption isotherms of these gases in nine different ZIFs with various topologies (SOD, RHO, GME) and functionalities. We find excellent agreement with experimental isotherms, thus validating the transferability and accuracy of the force-field parameters.

for pairwise exponents. Such a procedure leads to a well-defined set of pairwise exponents for every pair of elements (or ions). The physical basis for the use of these exponents and our combination rule is detailed in Supporting Information. As previously, the prefactors, Aexch ii , are then fit for specific atom types to SAPT exchange energies between CO 2 and appropriately chosen systems, with cross-terms given by a simple combining rule, as in our previous work. We also make modifications to our methodology related to the atomic dispersion parameters. It is extremely challenging to obtain physical dispersion parameters directly from a SAPT dispersion energy fit due to the long-range nature of these interactions and thus intrinsic coupling between atomic parameters. Thus, in the present work, we develop and employ a novel extension of the Williams−Stone method57 to obtain physical and transferable C6, C8, and C10 atom-centered dispersion parameters. The Williams−Stone method involves calculating the frequency-dependent, point-to-point response of a molecule (e.g., the electrostatic potential at grid point Q produced by the molecular response to a point charge at a different grid point P) and then fitting all of these responses to a set of distributed polarizability tensors,

■

FORCE-FIELD DEVELOPMENT The general methodology for our physically motivated forcefield development, based on fitting SAPT intermolecular interaction energies, has been described in our previous works.21,41 Briefly, we use density functional theory-based SAPT (DFT-SAPT)42−53 to calculate the intermolecular interactions for a variety of adsorbate/adsorbent geometries. We fit the inherent SAPT energy decomposition (exchange, electrostatic, induction, dispersion, and a small δ-Hartree−Fock term) on a term-by-term basis to physically appropriate functional forms to mitigate error cancelation, yielding a robust force field. We obtain atomic charges by fitting a distributed multipole analysis (DMA)54,55 for a molecule, and C6, C8, and C10 dispersion coefficients from monomer linear-response calculations (vide infra). Remaining exchange and charge penetration parameters are fit to appropriate SAPT dimer interaction energies (see ref 21 for complete details). Induction energies are fit using Drude oscillator models previously developed for CO2 and N2 in our SYM force field.41 Our previous ZIF force fields are extremely accurate and robust but may have limited transferability (because parametrization was performed for a specific adsorbate/adsorbent pair). Here we present several modifications to our methodology to enable the generation of transferable atomic force-field parameters. As in our previous work, the first-order exchange repulsion (1) , was fit to a pairwise additive exponential energy, Eexch functional form, (1) Eexch ≅

αPQ ̃ = − ∑ T0Pat αtuabTubQ 0 abtu

In the above equation (eq 2 of ref 57), α̃ PQ is the modelpredicted point-to-point response of the molecule at grid points 58 P and Q. TPa describing the 0t is an interaction function electrostatic potential or its derivatives (depending on index t) at site a, due to a point charge at point P; similarly, TbQ u0 is an interaction function describing the electrostatic potential at point Q due to a moment of rank u at site b. The polarizability tensors, αab tu , of rank t and u at sites a and b, respectively, are the quantities to be fit to best reproduce the ab initio molecular response.57 Here, we are only concerned with a local, isotropic polarizability model, as we want atom-centered, isotropic C6, C8, and C10 coefficients, and thus our frequency-dependent “tensors” are described by just two indices αai , where a refers to the site and i refers to the rank. The dispersion parameters are related to a quadrature of the imaginary frequency polariability tensor components by the Casimir−Polder formula.58 For large molecules, however, coupling between the various αai can lead to unphysical atomic parameters, and therefore constraints must be introduced.59 We propose a novel way of constraining the fits to the point-to-point molecular responses, and we show that this method gives physical and transferrable atomic C6, C8, and C10 coefficients for large molecules. In our approach, physical, zeroeth-order values for the αia(ω) (explicitly writing the frequency dependence) can be obtained through unconstrained implementation of the Williams−Stone method on small molecules (i.e., N2, CH4). Next, specific atom types can be defined, and point-to-point responses can be calculated for a group of molecules containing these desired atom types (e.g., linear and branched alkanes, with different atom types of carbon and hydrogen atoms depending on their hybridization). Then, an iterative fitting procedure is carried out, where all molecules are fit independently, but the αai (ω) are only explicitly fit for one atom type at a time, with the remaining αai (ω) for all other atom types explicitly constrained to the values obtained from previous iterations. After the αai (ω)

∑ Aijexchexp(−Bijrij) i,j

where the indices i and j cover all atoms of the respective monomers; the same exponents Bij are also utilized in all charge-penetration terms for other interaction energy components. In our previous work, both the prefactor, Aij, and a scale factor for the exponents were obtained by directly fitting the SAPT exchange energy, leading to a strong coupling between the magnitude of the prefactors and exponent scale factor. Utilizing the fact that the asymptotic decay of atomic electron densities is rigorously exponential, with an exponent equal to 2 × (2 × IP)1/2, where IP is the atomic ionization potential,56 we set Bii = 2 × (2 × IPi)1/2, and we use the combination rule 14032

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

parameter fitting) basis. These test SAPT calculations were performed using the CAMCASP program.62 At very large separations (small dispersion energies), the monomer-centered dAVTZ basis gives dispersion energies slightly lower (more attractive) than that of the dimer-centered basis set, which is to be expected as midbond and partner monomer basis functions are useless at large separations, and the extra diffuse functions in the monomer dAVTZ basis makes this basis more complete. However, at separations of 1.1−1.4 times the van der Waals surface, we find that the dimer-centered + midbond basis set almost always yields systematically more attractive dispersion energies, sometimes by as much as 10% for certain configurations at 1.1−1.2 van der Waals separations. We fit a single scale factor to correct for this basis-set deficiency in our monomer point-to-point response (and thus C6, C8, C10 coefficients) calculations. As our Li2−imidazolate and CO2 interaction SAPT calculations form the basis for our ZIF parameter fitting, we fit the dispersion energy for those single points having CO2-ring separations of 1.25 times the van der Waals surface or greater. We use the C6, C8, C10 coefficients from our self-consistent fitting procedure and fit a single scale factor for these parameters to reproduce the SAPT dispersion energies. We find a scale factor of ∼1.07 from this fit, consistent with the observed basis set discrepancies found in the above test calculations. This scale factor, along with all C6, C8, C10 coefficients from the self-consistent point-to-point response calculations, and the exponents, Bij, for the damping function, form the complete set of parameters needed to model the dispersion interactions for all systems studied: E(2) VdW ≅ −∑n=6,8,10∑i,j f n(Bij,rij)(Cijn/rnij). We find that these parameters do an excellent job of reproducing the dispersion interactions from SAPT calculations for all systems (see Supporting Information). The remaining imidazolate electrostatic, polarization, and δHF terms were fit to ∼1500 CO2/Li2−imidazolate SAPT calculations, as in our previous work. These parameters are then used, without modification, for all functionalized ZIFs. However, functionalized ZIFs require additional functional group parameters to describe the unique chemical functionalities. These transferable functional group parameters were obtained from additional SAPT calculations on dimer systems of the specified model system and a CO2 molecule, as detailed in Table 1. For these additional calculations, the density fitting DFT-SAPT (DF-DFT-SAPT) method42−53 was used as implemented in the MolPro 2009 package.63 Note that for parameter fitting in some cases constraints were employed, using previous calculated atomic parameters to decouple parameters, in the same spirit as discussed above. Additional details are given in Supporting Information.

are fit for a given atom type, these values are averaged over all like atom types in all of the molecules considered, and these average values are then used in the subsequent iteration. Upon convergence, this procedure yields self-consistent αia(ω) parameters for all desired atom types in a specific class of molecules. Using this procedure, a “library” of αia(ω) parameters for different atom types can be constructed, and when fitting the response of a large molecule, one can use previously determined αai (ω) values for specific atom types as explicit constraints in the point-to-point response fits so as to greatly reduce the number of parameters fit and ensure physical resulting parameters. We use the CAMCASP program59−62 to calculate point-topoint responses for our various molecules of interest. The monomer response is computed at the coupled-perturbed Kohn−Sham (CPKS) level of theory (consistent with the DFTSAPT calculations), using an asymptotically corrected PBE0 exchange-correlation functional for the SCF calculation, and an ALDAX + CHF kernel (25% exact Hartree−Fock exchange) for the linear response kernel. We use a dAVTZ basis for all calculations, as this was previously shown to give almost identical dispersion parameters as an AVQZ basis, but with much lower computational cost.59 The CAMCASP defaults of 2000 grid points at distances between 2× and 4× the van der Waals surface were used, and responses were calculated at 10 different imaginary frequencies for every molecule. Polarizability tensors (in this case one value for isotropic treatment) were considered up to rank three for heavy atoms, and rank one for hydrogen atoms as recommended.59 We use the approach outlined above to fit αai (ω) parameters for various atom types of C, H, and N using both alkane and aromatic groups of molecules. Subsequently, additional αai (ω) parameters were fit for functional groups, using explicit constraints of appropriate atom types from the aromatic and alkane “library”. Atomic C6, C8, and C10 parameters, Ciin, can then be constructed using the known relations.57,58 The atomic dispersion parameters that are generated in this manner all have physically sensible values, exhibiting trends across different elements and varying molecular environments. Most importantly, the αai (ω) parameters (and thus dispersion parameters) are shown to be transferable between molecules, as we show that the αai (ω) parameters fit to small molecules are directly able to reproduce the molecular response of larger molecules, with negligible error compared to explicitly fitting parameters for these larger systems (see Supporting Information). We apply a single, system-independent basis set correction to the Ciin coefficients obtained from the monomer point-to-point response calculations. This scale factor accounts for the significant basis set dependence of the dispersion parameters for both the monomer and SAPT (dimer interaction energy) calculations. The equivalence (or lack thereof) of basis sets used in monomer response calculations and those used in a dimer interaction energy calculation has been previously discussed,59 and in general for a given cardinal number of basis set, a dimer-centered plus midbond basis used in the dimer interaction calculation will be more complete than the monomer-centered basis used in the monomer response calculation. We have tested the basis set effects on dispersion energy by comparing DFT-SAPT calculations for Li2− imidazolate and CO2 interactions, using both a monomercentered dAVTZ basis (consistent with the monomer point-topoint response calculations) and a dimer-centered AVTZ + midbond (consistent with the DFT-SAPT calculations used in

Table 1. Model Systems and the Corresponding Atom-Type Force-Field Parameters Fit to Molecule/CO2 SAPT Dimer Calculations model system Li2−imidazolate methyl−Li2− imidazolate dichloro−Li2− imidazolate benzene nitromethane methyl bromide 14033

atom types fit C, N, H of imidazolate rings H of methyl group (methyl carbon explicitly constrained) Cl C, H of benzene N, O of nitro group Br dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

Table 2. Atom-Type Parameters Used for All ZIF Simulationsa atom type

Aexch

Aelec

Aind

Adhf

Bb

C6

C8

C10

C_CO2 O_CO2 N_N2 C1_imidazolate C2_imidazolate N_imidazolate H_imidazolate C_benzene H_benzene C_methyl H_methyl N_nitro O_nitro Cl Br Zn

7.4377 35.4373 69.2851 52.6213 43.6947 113.3368 0.8751 42.3935 1.0600 43.6947 0.9218 42.8024 37.7135 302.9496 456.4617 0.0000

2.4130 10.8209 18.7757 12.7137 12.7137 44.4084 0.1770 20.4172 0.2749 12.7137 0.1780 18.2844 14.0482 109.6660 225.9380 0.0000

11.2513 0.2545 0.9393 0.0000 0.0000 6.4396 0.0973 1.8741 0.0892 0.0000 0.0972 0.9404 1.2079 8.3313 29.9575 0.0000

0.0795 −1.8179 −3.9712 4.7903 4.7903 34.7902 0.1390 4.4968 0.0971 4.7903 0.1307 3.9620 4.8749 24.6430 35.9849 0.0000

3.4384 3.7795 3.9061 3.4384 3.4384 3.9061 3.7776 3.4384 3.7776 3.4384 3.7776 3.9061 3.7795 3.6899 3.5219 6.4579

11.4742 8.6728 10.7006 15.7559 15.7559 10.5033 1.6428 15.8465 1.6428 10.6710 1.5657 10.7337 10.4874 55.2918 99.1300 0.0000

6.3294 4.2665 7.0484 7.1796 7.1796 7.8413 0.0000 6.8553 0.0000 10.0158 0.0000 7.1966 2.9248 29.1647 72.0298 0.0000

2.9660 2.8761 4.8132 6.2594 6.2594 12.2944 0.0000 17.1885 0.0000 17.8584 0.0000 5.3127 2.4801 29.6625 66.4232 0.0000

Units are as follows: Aexch, Aelec, Aind, Adhf are in kJ/mol/104; B are in Å−1; C6 are in [kJ/(mol·Å6)]/102; C8 are in [kJ/(mol·Å8)]/103; C10 are in [kJ/(mol·Å10)]/104 . Here C1_imidazolate is always the carbon atom at the top of the imidazolate ring, bonded to two nitrogen atoms. C2_imidazolate is therefore the two carbon atoms at the bottom of the imidazolate ring, bonded to each other and each to one other nitrogen atom. b Note that the “B” exponents reported here are given explicitly by Bii = 2 × (2 × IPi)1/2 and are the values used in all CO2−ZIF simulations. For N2− ZIF simulations, all other parameters are the same, but these exponents are scaled by 0.98 as described below. a

expense of overpredicting repulsive interactions (shown in Supporting Information), which in general should not present a problem, as these configurations do not contribute significantly to the partition function. We utilized these scaled exponents in all N2−ZIF simulations.

Note that the imidazolate parameters obtained directly from the SAPT fits yield total interaction energies which are systematically too attractive at very short distances. This is most likely due to a breakdown of the “isotropic, singleexponential” treatment of exchange and charge penetration terms at these short distances. Because such configurations are locally quite repulsive, they typically do not contribute to the partition function at relevant temperatures. Nonetheless, in RHO topologies (e.g., ZIF-25, -71) there exist pore apertures framed by four functionalized imidazolate rings. At these sites, charge penetration becomes increasingly important, as gas molecules can experience close interactions with conjugated-πelectron density from four different aromatic rings, leading to locally repulsive interactions (2−10 mH, based on SAPT estimates) that are strongly stabilized by long-range dispersive interactions. As such, the C1_imidazolate exchange parameter was scaled by 20%, to compensate for the systematic error in total energy for these configurations; the modified parameter is reflected in Table 2. This scaling negligibly affects the rms errors of the original SAPT fits for the imidazolate systems, indicating that almost all of the CO2−imidazolate single points utilized for parameter fitting are insensitive to this 20% scaling of the C1_imidazolate exchange parameter. By exploiting the transferability of the derived atomic exchange, electrostatic, polarization, and dispersion parameters, we are able to trivially generate a corresponding force field for the N2−ZIF interaction. We fit molecular N2 parameters to our previous N2−N2 SAPT calculations41 and use combination rules with the previously fit ZIF parameters to generate crossterms. We validate the results against ∼1000 SAPT N2/Li2− imidazolate calculations. The force-field-predicted interaction energies are in relatively good agreement with SAPT energies, with an rms error in the total energy of ∼0.5 mH, which is very satisfactory considering that no parameters were explicitly fit to this data set (see Supporting Information). However, the forcefield-predicted total energy is slightly systematically too attractive compared to the SAPT energies. Simply scaling the exponents, Bij, by 0.98 corrected this systematic error, at the

■

RESULTS AND DISCUSSION We have focused our attention to ZIFs that (we believe) have well-determined structures, composed of either single linker types or mixtures of linkers that, due to steric constraints, crystallize in a well-defined pattern. As such, we are able to guarantee that our simulated systems are commensurate with the corresponding experimental data. Out of the nine ZIFs that we investigate, ZIF-8, ZIF-71, ZIF-25, ZIF-68, and ZIF-70 all have symmetrically functionalized imidazolate rings and therefore well-defined crystal structures (in ref 16, the imidazolate to nitroimidazolate ratio is reported to be 1.13/ 0.87, but the CIF file from the Cambridge Structural Database (CSD) contains a 1:1 ratio; we use the latter without modification). The remaining ZIFs of the GME topology (ZIF-69, -78, 79, and -81) all have nonsymmetrical functionalization of the benzene−imidazolate ring. However, steric considerations essentially determine the choice of occupation for these functional groups. We utilize staggered configurations for functional groups on adjacent benzene rings to avoid unfavorable steric interactions, consistent with the unit cell of ZIF-69 used in previous simulations and depicted in refs 27 and 33. In contrast, for other ZIFs (not investigated in this work) ,such as ZIF-82 (GME topology) and ZIF-93, ZIF-96, and ZIF-97 (RHO), there seems no definitive way to choose unique occupancies for the functional groups. Grand canonical Monte Carlo (GCMC) simulations were used to compute CO2 and N2 adsorption isotherms in the nine ZIFs mentioned above, using the force fields developed in this work. For CO2−CO2 and N2−N2 interactions, we utilize our previous “SYM” force fields which have been developed using a similar methodology and are extensively benchmarked to bulkphase thermodynamic data.41,64 All ZIF−adsorbate interactions 14034

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

exceptions, as discussed below), the predicted isotherms are in good agreement with experiment, and it is evident that our transferable force field is robust across different topologies (RHO, SOD, GME) and functionalizations (methyl, benzene, nitro, halides). We have demonstrated previously that accurately reproducing an experimental gas adsorption isotherm is not sufficient to guarantee “physical” force-field parameters.21 However, our ZIF FF was explicitly constructed to contain physical parameters that reproduce individual interaction energy components both asymptotically and in charge penetration regions, and therefore the accuracy of the predicted isotherms (in the absence of any empirical, adjustable parameters) testifies to the quality of these parameters. The accuracy and transferability of ZIF FF in predicting CO2 and N2 adsorption isotherms opens the door to a truly predictive computationally screening of novel, previously unsynthesized ZIFs. In Figure 3 and Figure 4 we show the average adsorption energies of CO2 and N2 in the different ZIFs (at ∼0.8−1.0 bar,

are treated with ZIF FF, described above. A rigid treatment of ZIF atom positions is employed in all simulations. We use a hybrid MC/MD approach, with simulation methodology analogous to our previous work.21 For the long-range part of the Buckingham potential, cutoffs are used as dictated by the size of the unit cell (or super cell), with values of 16 Å for ZIF-8 (2 × 2 × 2 super cell), 14 Å for ZIF-71 and ZIF-25 (unit cell), and 16 Å for all ZIFs of the GME topology (2 × 2 × 2 super cells). The Peng−Robinson equation of state was used to correct for nonidealities in the chemical potential for both CO2 and N2. The GCMC-computed CO2 isotherms at 298 K for six ZIFs in the isoreticular series (ZIF-68, ZIF-69, ZIF-70, ZIF-78, ZIF79, ZIF-81) of the GME topology as well as ZIF-71 and ZIF-25 of the RHO topology are shown in Figure 1 (the N2 isotherms

Figure 1. CO2 isotherms at 298 K, solid lines are experimental data, symbols are our simulations. Experimental data for ZIF-25 and ZIF-71 are taken from ref 22, and experimental data for GME ZIFs are taken from ref 16.

Figure 3. Energy decomposition reflecting ratio of average dispersion interactions to electrostatic interactions in different ZIFs from simulations at ∼0.8−1.0 bar, and ∼300 K. Blue columns are for CO2, and red columns are for N2.

for these ZIFs are given in Supporting Information). Both CO2 and N2 isotherms computed at 303 K over a different pressure range (so as to compare with experimental data) in ZIF-8 (SOD topology) are shown in Figure 2. In general (with several

Figure 4. Average magnitude of adsorption energies in different ZIFs from simulations at ∼0.8−1.0 bar, and ∼300 K. Blue columns are for CO2, and red columns are for N2.

∼300 K) as well as the relative contributions of dispersion and electrostatic interactions to these total adsorption energies. The last six ZIFs in the figures are the isoreticular series of the GME topology, and they are listed in order of increasing CO2 adsorption at 298 K and 0.8 bar (ZIF-70 to ZIF-78). We can see that for this isoreticular series, CO2 uptake generally

Figure 2. CO2 (red) and N2 (blue) isotherms in ZIF-8 at 303 K. Experimental data are shown as solid lines, and our simulation results are shown as symbols. The experimental isotherms are taken from ref 23. 14035

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

errors in our force field occur in the treatment of short-range charge penetration, and such errors should be revealed in these calculations. As a control, we also conduct this analysis for ZIF8, ZIF-71, and ZIF-68, which all exhibit excellent isotherms. In those cases the CO2 interaction energies are accurately reproduced by our force field (see Supporting Information). In contrast, for ZIF-78 and -25, our force field overpredicts CO2 adsorption, and we expect some error in our force-fieldpredicted energies for configurations sampled from this GCMC simulation. Indeed we see errors in predicting large exchange energies, as well as significant scatter in predicting the electrostatic energy, resulting in larger errors in total energy compared to the three previously mentioned systems (see Supporting Information). As these individual interaction components are usually significantly larger than the total interaction energy, the prediction of total interaction energies is more susceptible to error in charge penetration regions where there is in general greater cancelation between individual components. It seems that these challenging cases may require a more sophisticated treatment of exchange and charge penetration, rather than just pairwise, atom−atom isotropic exponential terms. However, it is important to realize that our treatment of exchange and charge penetration is successful in most cases, with the significant advantage that it provides a methodology to generate transferable parameters for potentially any pairwise atom types. As the induction energy is a relatively small contribution to the total adsorption energy of CO2 and N2 in ZIFs (∼5−10% of the total interaction energy), it is useful to examine whether it is necessary to use the expensive Drude oscillator models in the SYM force field for these gas solutes. One might imagine that for computational screening of ZIFs for gas adsorption, it would be sufficient to use ZIF FF in conjunction with SYM without Drude oscillator polarizability. To test this, we computed CO2 isotherms for three ZIFs exhibiting different contributions of induction energy, each isotherm computed with and without Drude oscillators for the CO2 molecules. In Figure 6, we see that contributions to the isotherms from the Drude oscillators range from a small amount for ZIF-71, to slightly more for ZIF-79, to a relatively significant amount

increases with increasing average adsorption energy (Figure 5). This is not, in general, true across different topologies, as ZIF-

Figure 5. Correlation between total CO2 uptake (simulation results) and the average adsorption energy (magnitude) per CO2 molecule for the six members of the isoreticular GME series at 298 K and 0.8 bar.

25 has an average CO2 adsorption energy higher than that of any GME ZIF, and yet it has a CO2 uptake much lower than that of most of the GME series. This reflects the fact that the number of adsorption sites, surface area, and free volume are also very important factors determining gas uptake, and these properties can dramatically vary over different topologies. Using the intrinsic energy decomposition of our force fields, we can analyze the contribution of different interaction components to the total adsorption energies. We find that gas adsorption in these ZIFs is dominated by dispersion interactions, as was found previously,21 with the ratio of average dispersion energies to electrostatic energies being 2:1 to 3:1 for CO2 and 4:1 to 5:1 for N2. It might be expected that as electrostatic interactions become more important, CO2/N2 adsorption selectivity should increase due to the larger quadrupole moment of CO2, and this does hold for ZIF-78 which has the highest selectivity (∼50 experimentally)16 and the smallest dispersion/electrostatic interaction ratio (∼2.3 for CO2) of the GME ZIFs. However, the trend is far from perfect, as ZIF-70 also exhibits a relatively low dispersion/electrostatic ratio (∼2.4 for CO2) and yet has a low selectivity (∼17 experimentally)16 for GME ZIFs. This implies that the magnitude of the adsorbate−ZIF interactions is also important for selectivity, and ZIF-70 probably exhibits low selectivity due to its relatively low average adsorption energies of CO2 and N2. There are a few cases where agreement with the experimental isotherms is imperfect, and it is instructive to analyze the nature of this disagreement for future force-field improvements. For instance, although the simulated CO2 isotherms in ZIF-81 and ZIF-78 agree very well with experiment at lower pressures, the predicted CO2 uptakes at higher pressures in both of these systems are systematically too high. Similarly, we see that with ZIF-25 our force field predicts systematically too high adsorption in all pressure regimes. To understand these errors, we sampled representative CO2 configurations from the GCMC trajectories for various ZIFs, and conduct DF-DFTSAPT calculations between these CO2 molecules and the closest functionalized imidazolate ring (Li-capped as previously) for each configuration. We anticipate that the dominant

Figure 6. Comparison of predicted CO2 adsorption isotherms in ZIF71, ZIF-78, and ZIF-79 using ZIF FF combined with the SYM CO2 model with (triangles, solid lines) and without (circles, dashed lines) Drude oscillator polarizability. The lines are only meant to guide the eye. 14036

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

tions. We find that the GCMC-predicted CO2 and N2 isotherms using this force field are in good agreement with experimental results. Although force-field parameters were only generated for a limited number of functional groups, our methodology is general and can easily be extended to develop parameters for other functional groups as needed for future ZIF studies. As such, this new force field has potential applications in computational screening of novel ZIFs. These results also help us explain the otherwise scattered and inconclusive results of utilizing the UFF force field to predict gas adsorption in ZIFs. The significant overprediction of gas adsorption in ZIF-8 (SOD topology) and ZIF-68 and ZIF-69 (GME topology) is due to the systematic overestimation of asymptotic dispersion, requiring significant parameter scaling for agreement with experimental isotherms.23,27,31 Furthermore, we conclude that the seemingly good performance of UFF in predicting gas adsorption in ZIFs of the RHO topology (ZIF-25, ZIF-71)22 is only achieved through error cancelation between this too-attractive treatment of asymptotic dispersion and overrepulsion in charge penetration regions. These errors are likely not (to a large extent) due to a deficiency in the UFF parametrization but rather are “intrinsic” in the unphysical LJ functional form and the use of “universal” (environmentindependent) atomic dispersion parameters. We suggest that a tailored, physically motivated force field such as our ZIF FF offers a more accurate and robust alternative to modeling gas sorption in ZIFs.

(∼30%) for ZIF-78. However, for simulations with the goal of fast prescreening of many systems for gas adsorption, this error in isotherm prediction may be acceptable, yielding approximately 1 order of magnitude in simulation performance. To understand the inconsistent performance of previous UFF ZIF/CO2 simulations, we compare our accurate SAPT ab initio energies to UFF predictions for a variety of representative CO2/imidazolate dimers. Overall, we find that UFF, combined with our accurate monomer point charges, generally predicts energies that are too attractive asymptotically and too repulsive at close charge penetration regions compared to SAPT interaction energies (see Supporting Information). This is unsurprising, as the Lennard−Jones (LJ) functional form does not have the ability to describe either charge penetration or asymptotic dispersion accurately, while still predicting van der Waals minima energies correctly. It is well-known that a Buckingham-like exponential term, as we use in our force fields, is much better suited for physically describing charge penetration than is the R12 repulsive term in a LJ potential.65 Also, because the LJ force field truncates the multipoleexpanded dispersion energy at the R−6 term, the C6 terms in UFF are likely artificially large to compensate for the neglect of attractive R−8 and R−10 force-field terms. As such, the asymptotic dispersion energy is predicted to be too large by the UFF parameters. For systems dominated by long-range dispersion interactions, UFF should overpredict gas adsorption, and this explains the previously mentioned simulation results in SOD (ZIF-8)23,31 and GME (ZIF-68,ZIF-69)23,27 topologies. However, because UFF is too repulsive at close charge penetration regions, there is the opportunity for fortuitous error cancelation in ZIFs with a delicate balance of attractive and steric interactions. As we have found in our simulations of ZIF-25 and ZIF-71, it indeed seems as though this balance occurs in ZIFs of the RHO topology, and this error cancelation could explain the fact that the UFF force field is able to qualitatively predict gas adsorption in these systems.22 This analysis of the UFF force field is consistent with the results of Chen et al., who found that for ligand sites in the MOF CuBTC, the UFF force field predicted stronger framework− methane interactions compared to high-level ab inito calculated energies.66

■

ASSOCIATED CONTENT

S Supporting Information *

Figure showing the correlation between exponents fit to atomic density overlaps and the exponent-predicted asymptotically by the ionization potential. Figure showing the correlation between exponents fit to the overlaps of different atomic densities and the exponents generated using our combination rule formula for exponents. Tables showing the C6, C8, C10 coefficients generated by the self-consistent fits to point-topoint responses of alkanes and aromatics. Table showing the rms error to low-frequency point-to-point responses of larger molecules modeled using polarizability tensor components from our previously fit library. Figures showing UFF-predicted dispersion energies and total energies compared to SAPT calculations for CO2 and Li2−imidazolate configurations. Figures showing the SAPT energy fits used to generate forcefield parameters for CO2 and Li2−imidazolate interactions. Figures showing the SAPT energy fits used to generate forcefield parameters for CO 2 and methyl−Li2 −imidazolate interactions. Figures showing the SAPT energy fits used to generate force-field parameters for CO2 and dichloro−Li2− imidazolate interactions. Figures showing the SAPT interaction energy fits for dimethyl−Li2−imidazolate and CO2 for configurations pulled from a ZIF-25 simulation used to modify the C1_imidazolate exchange parameter. Figures showing the SAPT interaction energy fits for N2 and Li2−imidazolate configurations used to justify the 0.98 scale factor for exponents used in N2 GCMC simulations. Figures showing charges used for all functionalized imidazolate rings in ZIF simulations. Figures showing N2 computed isotherms for ZIF-25, ZIF-71, ZIF-68, ZIF-69, ZIF-70, ZIF-78, ZIF-79, and ZIF-81. Figures showing force-field predictions compared to SAPT energies for closest-ring CO2 configurations pulled from GCMC simulations for ZIF-8, ZIF-71, ZIF-68, ZIF-25, and ZIF-78. This

■

CONCLUSION We have developed a transferable force field that can be used to study CO2 and N2 adsorption in functionalized ZIFs, yielding “ZIF FF” force field. In the process, we generalized our existing SAPT-based approach to yield transferable atomic exchange, charge penetration, and dispersion parameters. In the latter case, we utilize a self-consistent, iterative fitting approach based on the Williams−Stone method57 of fitting frequency-dependent point-to-point molecular responses to generate a library of frequency-dependent polarizability tensor components for particular atom types. We have demonstrated the transferability of these atom-type frequency-dependent polarizability tensors by using them to model the point-to-point response of larger molecules for which they were not parametrized and have demonstrated their ability to very accurate reproduce SAPT dispersion energies. We believe that this general methodology for creating accurate, robust, physically motivated, and transferable force fields is widely applicable to a host of systems. We apply the resulting force fields to nine ZIFs spanning different topologies (RHO, SOD, GME) as well as an isoreticular series (GME topology) with varying functionaliza14037

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

(23) Perez-Pellitero, J.; Amrouche, H.; Siperstein, F. R.; Pirngruber, G.; Nieto-Draghi, C.; Chaplais, G.; Simon-Masseron, A.; Bazer-Bachi, D.; Peralta, D.; Bats, N. Chem.Eur. J. 2010, 16, 1560−1571. (24) Keskin, S. J. Phys. Chem. C 2011, 115, 800−807. (25) Sirjoosingh, A.; Alavi, S.; Woo, T. K. J. Phys. Chem. C 2010, 114, 2171−2178. (26) Krishna, R.; van Baten, J. M. J. Membr. Sci. 2010, 360, 323−333. (27) Liu, B.; Smit, B. J. Phys. Chem. C 2010, 114, 8515−8522. (28) Guo, H. C.; Shi, F.; Ma, Z. F.; Liu, X. Q. J. Phys. Chem. C 2010, 114, 12158−12165. (29) Zhou, M.; Wang, Q.; Zhang, L.; Liu, Y. C.; Kang, Y. J. Phys. Chem. B 2009, 113, 11049−11053. (30) Pantatosaki, E.; Pazzona, F. G.; Megariotis, G.; Papadopoulos, G. K. J. Phys. Chem. B 2010, 114, 2493−2503. (31) Battisti, A.; Taioli, S.; Garberoglio, G. Microporous Mesoporous Mater. 2011, 143, 46−53. (32) Rankin, R. B.; Liu, J. C.; Kulkarni, A. D.; Johnson, J. K. J. Phys. Chem. C 2009, 113, 16906−16914. (33) Liu, D. H.; Zheng, C. C.; Yang, Q. Y.; Zhong, C. L. J. Phys. Chem. C 2009, 113, 5004−5009. (34) Assfour, B.; Leoni, S.; Seifert, G. J. Phys. Chem. C 2010, 114, 13381−13384. (35) Han, S. S.; Choi, S. H.; Goddard, W. A. J. Phys. Chem. C 2010, 114, 12039−12047. (36) Han, S. S.; Choi, S. H.; Goddard, W. A. J. Phys. Chem. C 2011, 115, 3507−3512. (37) Wilmer, C. E.; Leaf, M.; Lee, C. Y.; Farha, O. K.; Hauser, B. G.; Hupp, J. T.; Snurr, R. Q. Nat. Chem. 2012, 4, 83−89. (38) Haldoupis, E.; Nair, S.; Sholl, D. S. J. Am. Chem. Soc. 2012, 134, 4313−4323. (39) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024−10035. (40) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Chem. Rev. 1994, 94, 1887−1930. (41) Yu, K.; McDaniel, J. G.; Schmidt, J. R. J. Phys. Chem. B 2011, 115, 10054−10063. (42) Bukowski, R.; Podeszwa, R.; Szalewicz, K. Chem. Phys. Lett. 2005, 414, 111−116. (43) Hesselmann, A.; Jansen, G. Chem. Phys. Lett. 2002, 362, 319− 325. (44) Hesselmann, A.; Jansen, G. Chem. Phys. Lett. 2002, 357, 464− 470. (45) Hesselmann, A.; Jansen, G. Chem. Phys. Lett. 2003, 367, 778− 784. (46) Hesselmann, A.; Jansen, G. Phys. Chem. Chem. Phys. 2003, 5, 5010−5014. (47) Hesselmann, A.; Jansen, G.; Schutz, M. J. Chem. Phys. 2005, 122, 14103−14119. (48) Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Phys. Rev. Lett. 2003, 91, 33201−33204. (49) Misquitta, A. J.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. J. Chem. Phys. 2005, 123, 214103−214116. (50) Misquitta, A. J.; Szalewicz, K. Chem. Phys. Lett. 2002, 357, 301− 306. (51) Misquitta, A. J.; Szalewicz, K. J. Chem. Phys. 2005, 122, 214109− 214127. (52) Podeszwa, R.; Bukowski, R.; Szalewicz, K. J. Chem. Theory Comput. 2006, 2, 400−412. (53) Williams, H. L.; Chabalowski, C. F. J. Phys. Chem. A 2001, 105, 646−659. (54) Stone, A. J. Chem. Phys. Lett. 1981, 83, 233−239. (55) Stone, A. J.; Alderton, M. Mol. Phys. 1985, 56, 1047−1064. (56) Levy, M.; Perdew, J. P.; Sahni, V. Phys. Rev. A 1984, 30, 2745− 2748. (57) Williams, G. J.; Stone, A. J. J. Chem. Phys. 2003, 119, 4620− 4628. (58) Stone, A. J. The Theory of Intermolecular Forces; Clarendon Press: Oxford, 1996. (59) Misquitta, A. J.; Stone, A. J. Mol. Phys. 2008, 106, 1631−1643.

material is available free of charge via the Internet at http:// pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank DOE for their support through grant DE-FG0209ER16059. Additional computational resources were provided via the Center for High Throughput Computing at the University of Wisconsin67 and by National Science Foundation Grant CHE-0840494.

■

REFERENCES

(1) Metz, B. IPCC Special Report on Carbon Dioxide Capture and Storage; Cambridge University Press: Cambridge, 2005. (2) Dang, L. X.; Wick, C. D. J. Phys. Chem. B 2011, 115, 6964−6970. (3) Gurkan, B.; Goodrich, B. F.; Mindrup, E. M.; Ficke, L. E.; Massel, M.; Seo, S.; Senftle, T. P.; Wu, H.; Glaser, M. F.; Shah, J. K.; Maginn, E. J.; Brennecke, J. F.; Schneider, W. F. J. Phys. Chem. Lett. 2010, 1, 3494−3499. (4) Gurkan, B. E.; de la Fuente, J. C.; Mindrup, E. M.; Ficke, L. E.; Goodrich, B. F.; Price, E. A.; Schneider, W. F.; Brennecke, J. F. J. Am. Chem. Soc. 2010, 132, 2116−2117. (5) Shiflett, M. B.; Drew, D. W.; Cantini, R. A.; Yokozeki, A. Energy Fuels 2010, 24, 5781−5789. (6) Wick, C. D.; Chang, T.-M.; Dang, L. X. J. Phys. Chem. B 2010, 114, 14965−14971. (7) Zhao, X.; Xing, H. B.; Li, R. L.; Yang, Q. W.; Su, B. G.; Ren, Q. L. Prog. Chem. 2011, 23, 2258−2268. (8) Choi, S.; Gray, M. L.; Jones, C. W. ChemSusChem 2011, 4, 628− 635. (9) Li, J.-L.; Chen, B.-H. Sep. Purif. Technol. 2005, 41, 109−122. (10) Mavroudi, M.; Kaldis, S. P.; Sakellaropoulos, G. P. Fuel 2003, 82, 2153−2159. (11) Britt, D.; Furukawa, H.; Wang, B.; Glover, T. G.; Yaghi, O. M. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 20637−20640. (12) Caskey, S. R.; Wong-Foy, A. G.; Matzger, A. J. J. Am. Chem. Soc. 2008, 130, 10870−10871. (13) McDonald, T. M.; D’Alessandro, D. M.; Krishna, R.; Long, J. R. Chem. Sci. 2011, 2, 2022−2028. (14) Sumida, K.; Rogow, D. L.; Mason, J. A.; McDonald, T. M.; Bloch, E. D.; Herm, Z. R.; Bae, T.-H.; Long, J. R. Chem. Rev. 2011, 112, 724−781. (15) Babarao, R.; Jiang, J. W. J. Am. Chem. Soc. 2009, 131, 11417− 11425. (16) Banerjee, R.; Furukawa, H.; Britt, D.; Knobler, C.; O’Keeffe, M.; Yaghi, O. M. J. Am. Chem. Soc. 2009, 131, 3875−3877. (17) Phan, A.; Doonan, C. J.; Uribe-Romo, F. J.; Knobler, C. B.; O’Keeffe, M.; Yaghi, O. M. Acc. Chem. Res. 2010, 43, 58−67. (18) Park, K. S.; Ni, Z.; Cote, A. P.; Choi, J. Y.; Huang, R. D.; UribeRomo, F. J.; Chae, H. K.; O’Keeffe, M.; Yaghi, O. M. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 10186−10191. (19) Kusgens, P.; Rose, M.; Senkovska, I.; Frode, H.; Henschel, A.; Siegle, S.; Kaskel, S. Microporous Mesoporous Mater. 2009, 120, 325− 330. (20) Banerjee, R.; Phan, A.; Wang, B.; Knobler, C.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Science 2008, 319, 939−943. (21) McDaniel, J. G.; Yu, K.; Schmidt, J. R. J. Phys. Chem. C 2011, 116, 1892−1903. (22) Morris, W.; Leung, B.; Furukawa, H.; Yaghi, O. K.; He, N.; Hayashi, H.; Houndonougbo, Y.; Asta, M.; Laird, B. B.; Yaghi, O. M. J. Am. Chem. Soc. 2010, 132, 11006−11008. 14038

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039

The Journal of Physical Chemistry C

Article

(60) Misquitta, A. J.; Stone, A. J. J. Chem. Theory Comput. 2008, 4, 7− 18. (61) Misquitta, A. J.; Stone, A. J.; Price, S. L. J. Chem. Theory Comput. 2008, 4, 19−32. (62) See also http://www-stone.ch.cam.ac.uk/programs. html#CamCASP. (63) H.-J.Werner; Knowles, P. J.; Lindh, R.; Manby, F. R.; Schutz, M.; Celani, P.; Korona, T.; Mitrushenkov, A.; Rauhut, G.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hetzer, G.; Hrenar, T.; Knizia, G.; K̈ oppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; Palmieri, P.; Pfluger, K.; Pitzer, R.; Reiher, M.; Schumann, U.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Wolf, A. MOLPRO, version 2009.1, a package of ab initio programs, see http:// www.molpro.net. (64) Yu, K.; Schmidt, J. R. J. Chem. Phys. 2012, 136. (65) Stone, A. J.; Misquitta, A. J. Int. Rev. Phys. Chem. 2007, 26, 193− 222. (66) Chen, L. J.; Grajciar, L.; Nachtigall, P.; Duren, T. J. Phys. Chem. C 2011, 115, 23074−23080. (67) Litzkow, M.; Livney, M.; Mutka, M. Condor - A Hunter of Idle Workstations; IEEE: New York, 1988.

14039

dx.doi.org/10.1021/jp303790r | J. Phys. Chem. C 2012, 116, 14031−14039