ARTICLE pubs.acs.org/JPCC
Ab Initio, Physically Motivated Force Fields for CO2 Adsorption in Zeolitic Imidazolate Frameworks Jesse G. McDaniel,† Kuang Yu,† and J. R. Schmidt* Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin—Madison, Madison, Wisconsin 53706, United States
bS Supporting Information ABSTRACT: We present an entirely ab initio methodology, based on symmetry adapted perturbation theory (SAPT), for constructing force-fields to study CO2 adsorption in nanoporous zeolitic imidazolate frameworks (ZIFs). Our approach utilizes the SAPT energy decomposition to generate physically motivated force fields for the CO2-ZIF interaction, with explicit terms representing exchange, electrostatic, induction and dispersion interactions. Each of these terms is fit to the corresponding term in the SAPT energy decomposition, yielding a force field entirely free of empirical parameters. This approach was utilized to construct force fields describing the CO2 interaction with both ZIF-8 and ZIF-71. In conjunction with our existing CO2CO2 force field, parametrized in a consistent manner, we validate our force fields using grand canonical Monte Carlo simulations and obtain good agreement with the corresponding experimental CO2 adsorption isotherms. Furthermore, the explicit correspondence between force field terms and fundamental interaction types (dispersion, electrostatics, and induction) allows for an analysis of the underlying physics controlling ZIF gas adsorption that is far more direct and well-defined than with the generic force fields that had been previously utilized to study these systems. As our force fields are free from empirical parameters, these results demonstrate the potential for computationally screening novel ZIFs for flue gas separation applications with near quantitative accuracy.
’ INTRODUCTION Concern about the increasing production of carbon dioxide (CO2) and its role in global warming has led to significant interest in postcombustion CO2 capture and separation. As coal fired power plants are one of the major point-source emitters of CO2, they are the obvious targets for the implementation of carbon capture and sequestration (CCS) technologies. One inherent challenge is that the flue gas emitted from such a power plant is relatively dilute in CO2, with the major component being nitrogen gas. Nanoporous metalorganic frameworks (MOFs) have been shown to exhibit high CO2 uptake with excellent adsorption selectively over nitrogen and other gases, and have thus received a large amount of research interest as potential CO2 separation materials. Here zeolitic imidazolate frameworks (ZIFs),135 consisting of metal cations bridged by organic imidazolate anions, have received particular attention due to their selective CO2 adsorption,1,2,7,8 excellent chemical and thermal stability,6,36 and potential for tunability via functionalization of the organic linker groups.2,4,7,8 Uptake and selective adsorption of CO2, N2, and CH4 in ZIFs is relevant to many applications, including natural gas recovery and CCS,13,5,7,8,3740 and a wide variety of ZIFs with different topologies, pore/aperture sizes, and chemical functionality have previously been developed.111,15,16,18,20,24,25,35 This tunability has led to the description of ZIF synthesis as structural “design” as opposed to materials “discovery”.7 Yet in order to develop and synthesize ZIFs optimized for a particular application, it is essential to have a fundamental understanding of the detailed r 2011 American Chemical Society
interactions between different gas molecules and ZIF frameworks. Molecular dynamics (MD) and Monte Carlo (MC) computer simulations are important tools for studying gas adsorption isotherms and diffusion in porous materials such as ZIFs. The use of accurate and physically appropriate force fields in these simulations can provide valuable insight into the important interactions between gas and framework that dictate adsorption capacity and/or selectivity. There has been a significant amount of previous computational work modeling the uptake and diffusion of various gases in different ZIFs. Much of this computational work has utilized standard “generic” force fields such as UFF41 or Drieding42 to describe the interactions of gas molecules with the ZIF framework, yielding mixed results and raising concerns regarding the applicability of these standard force fields to these systems. Nonetheless, many previous grand canonical Monte Carlo (GCMC) and MD studies have examined the uptake and separation of various gases including CO2, N2, CH4, and H2 with ZIFs, characterizing adsorption isotherms of these gases,5,38,4357 separation and selective adsorption of gas mixtures,4446,49 diffusion,4345,48,49 and important interaction sites in the ZIFs.38,43,58 In particular, ZIF-68 and ZIF-69 have been the subject of numerous computational studies targeting CO2 adsorption and Received: September 27, 2011 Revised: December 9, 2011 Published: December 15, 2011 1892
dx.doi.org/10.1021/jp209335y | J. Phys. Chem. C 2012, 116, 1892–1903
The Journal of Physical Chemistry C separation due to their excellent capacity and selectively.2 Various force fields have been used in these studies with mixed conclusions. For example, describing both CO2CO2 interactions and CO2ZIF interactions using the UFF force field yields reasonable agreement with the experimentally measured CO2 adsorption,50 although it is well-known that the UFF force field does a poor job of describing bulk CO2CO2 interactions. It is therefore perhaps not surprising that utilizing either the UFF/ Drieding parameters for the CO2ZIF interactions in conjunction with the (accurate) TraPPE59 force field for CO2CO2 interactions led to dramatic over prediction of CO2 adsorption.48 Other authors resorted to scaling the UFF parameters describing CO2ZIF interactions so that interactions were less attractive, then yielding good agreement with experimental isotherms.46 Recently it was also noted that by blocking small pores that were deemed “inaccessible” to gas molecules, either a UFF or Drieding description of gas-framework interaction could reproduce experimental isotherms to reasonable accuracy.47 Other notable applications of standard force fields to ZIFs include the simulations by Perez-Pellitero et al., in which adsorption isotherms for CO2, CH4, and N2 in ZIF-8, ZIF-76, and ZIF69 were found to agree with experiment only after introducing empirical scaling factors of 0.69 for ε values and 0.95 for UFF σ values,38 indicating a very significant empirical correction. Battisti et al. utilized this same scaling of the UFF gas-framework parameters to systematically study the adsorption and dynamics of CO2, CH4, N2, and H2 and mixtures thereof in ZIF-2 through ZIF-10.49 Morris et al. examined ZIFs 25, 71, 93, 96, and 97,5 using the EPM2 force field60 for CO2 and the UFF force field to describe the framework and were for the most part able to reproduce the experimental trends in CO2 adsorption with functionalization of the imidazolate linker groups, with the exception of ZIF-96. In the latter case they found that by using OPLS-ua parameters61 to describe the nitrile and amine groups they were able to obtain better agreement with the experimental isotherm for this system. In contrast, the literature contains few examples of ab initiobased force fields to describe gas adsorption/diffusion processes in MOFs. Schmid and co-workers have developed ab initio force fields to describe framework flexibility in a variety of MOFs,6264 this being of potential importance in understanding gas diffusion in these systems. Goddard and co-workers have parametrized functional group-H2 interactions based on ab initio calculations,53,54,65,66 enabling first principles calculation of H2 uptake in MOFs doped with metal cations. They have also used this methodology to investigate H2 adsorption in ZIFs.53,54 Thus very little work has focused on developing force fields tailored specifically to ZIFs, and in nearly all cases force field construction has been limited to deriving atomic charges from ab initio calculations and coupling these with Lennard-Jones (LJ) parameters of a generic force field such as UFF or Dreiding.5,38,43,47,48,50 To our knowledge, there have not been any attempts to develop rigorous ab intio force fields for the interaction of ZIFs with species relevant to flue gas separation. Thus overall, previous simulations of gas uptake and selective adsorption have often been forced to resort to arbitrary scaling factors parametrized to prior experimental data and a patchwork of standard force fields to achieve semi-quantitative agreement with experiment. This is somewhat unsatisfying as this scaling procedure may not capture the correct underlying physics. In other words, there are likely many possible scaling procedures (and thus many force fields) that reproduce the same isotherm, leading to an underspecified fitting problem and the possibility of
ARTICLE
“unphysical” parameters. We illustrate this principle via a simple example, CO2 adsorption in ZIF-8. Using the EPM2 model60 for CO2 in conjunction with UFF ZIF parameters and density functional theory (DFT) derived charges for the ZIF-8 framework (from ref 38), we locate three dramatically different parameter sets that all reproduce the same experimental adsorption isotherm. For the first parameter set, we used the scaling developed by Perez-Pellitero et al.,38 a uniform scaling of σ1scale = 0.95σUFF and ε1scale = 0.69εUFF. For the second parameter set, we use a different uniform scaling of the LJ parameters of σ2scale = 0.9σUFF and ε2scale = 0.9εUFF, and artificially modify the role of = 0.25qoriginal , where electrostatics by scaling all charges by q2scale i i original are the original charges developed by Perez-Pellitero et al. qi Note that this extreme scaling is certainly unphysical and is merely meant to illustrate the diverse parameter sets that are consistent with the experimental data. The third parameter set utilizes different scalings for the imidazolate ring atoms and the methyl functional group attached to the ring, thus modifying the importance of the functional group in comparison to the ring 3scale atoms: σ3scale ring = 0.95σUFF and εring = 0.85εUFF for the ring 3scale atoms and σmeth = 0.75σUFF and ε3scale meth = 0.5εUFF for the methyl atoms. The complete example parameter sets used in these simulations are listed in Table 1. As shown in Figure 1, these parameter sets yield very similar isotherms as predicted by GCMC. However, they suggest dramatically different interpretations as to the importance of electrostatic and dispersive interactions, as well as the role of the methyl functional group in determining the CO2 uptake, depending on the scaling method utilized. Furthermore, these algorithms would likely yield very different and inconsistent results when applied to other functionalized ZIFs. This inherent ambiguity when empirically optimizing existing force fields motivates our present efforts to develop parameter-free, ab initio force fields for nanoporous systems; nonetheless, approximate generic force fields are still extremely valuable for use in computational screenings, where the availability of parameters for arbitrary systems is a principle concern. These considerations motivated us to develop entirely firstprinciples based force fields for CO2ZIF interactions. In this work, we present a methodology for the development of such force fields and, utilizing our recently developed CO2 model for adsorbateadsorbate interactions,67 have carried out GCMC simulations to compute CO2 adsorption isotherms. We validate our approach via calculations on ZIF-8 and ZIF-71, as these ZIFs contain relatively “simple” linker groups. However the methodologies developed here are equally applicable to a vast array of existing and novel ZIFs. As in our previous work,67 the force fields developed here rely on the use of symmetry adapted perturbation theory (SAPT)68 to compute the interaction energies between two species. SAPT has the advantage of not only high accuracy, but it also naturally decomposes the interaction energy into the fundamental components of electrostatic, exchange, induction and dispersion interactions. This allows for the utilization of separate, physically appropriate terms to represent each of these interaction types. As they include essentially all of the relevant “physics”, these physically motivated force fields yield not only extremely accurate results but also a natural energy decomposition allowing for greater interpretation of thermodynamic predictions (adsorption isotherms) in terms of the underlying fundamental interactions. In order to ensure robust, physical, and asymptotically valid parameters in our force field, we obtain parameters from monomer calculations (polarizabilities, charge 1893
dx.doi.org/10.1021/jp209335y |J. Phys. Chem. C 2012, 116, 1892–1903
The Journal of Physical Chemistry C
ARTICLE
Table 1. Scaled UFF Parameter Sets Used to Compute the CO2 Isotherms Shown in Figure 1a parameter set 1 q
σ (Å)
q
ε (kJ/mol)
parameter set 3 σ (Å)
q
ε (kJ/mol)
σ (Å)
C1
0.640
0.3033
3.259
0.160
0.3954
3.088
0.640
0.3734
3.259
C2
0.080
0.3033
3.259
0.020
0.3954
3.088
0.080
0.3734
3.259
C3
0.670
0.3033
3.259
0.168
0.3954
3.088
0.670
0.2197
2.573
H2
0.144
0.1271
2.440
0.036
0.1657
2.314
0.144
0.1565
2.440
H3
0.144
0.1271
2.440
0.036
0.1657
2.314
0.144
0.0920
1.928
0.540
0.1993
3.097
0.135
0.2598
2.935
0.540
0.2454
3.097
1.100
0.3582
2.338
0.275
0.4669
2.215
1.100
0.4410
2.338
N Zn a
ε (kJ/mol)
parameter set 2
The atom types are shown in the corresponding methyl imidazolate ring pictured in Figure 2.
Figure 2. Imidazolate atom types.
Figure 1. The experimental CO2/ZIF-8 isotherm at 303 K38 is shown as the solid black dots. The GCMC predicted isotherm resulting from the first set of parameters is depicted by the red squares and solid red line; second set of parameters depicted by the green diamonds and green dashed line; third set of parameters depicted by the blue triangles and blue dotted-dashed line.
distribution, etc.) whenever possible, and fitting remaining parameters that rely on charge-density overlap to the SAPT interaction energies. We should note that many previous interaction potentials have been derived based on SAPT interaction energy calculations,6975 and we briefly comment on the primary differences between our methodology and techniques used in previous works. Almost all previous works have fit the total interaction energy between molecules,69,70,72,74,75 while we choose to independently fit each individual interaction energy component. This not only allows for greater physical interpretation when such a force field is used in a computer simulation, but also reduces the possibility of error cancelation and helps to ensure physically meaningful parameters, especially in larger molecules. Our procedure is probably most similar to that of Misquitta et al.,73 with the major difference being that we abandon the use of explicit anisotropic force field terms so that our final force field has a form amenable to standard simulation packages. This paper is organized as follows: We first consider our choice of ZIF model system, examining the transformation from an infinite periodic ZIF crystal, to a medium sized model system fragment (Zn2Im7Li63+), and finally to a single Li-capped imidazolate type ring that we utilize to obtain force-field parameters.
Next we discuss the various monomer calculations and fitting procedures that we use to obtain the force-field parameters. Finally, we present the GCMC simulations and the computed CO2 adsorption isotherms and discuss possible applications of this approach to ZIF screening.
’ MODEL SYSTEM CONSTRUCTION SAPT was used to compute the interaction energies necessary to fit parameters in our force field. This technique is not amenable to periodic systems, and thus our first step was to develop a finite model system that was large enough to be representative of the bulk crystal, but small enough to enable tractable SAPT calculations. Our goal in constructing this model system was to reproduce the full three-dimensional electron density on the imidazolate rings and surrounding Zn2+ cations as compared to the corresponding bulk. In creating this model system, we considered a ZIF system with the simplest possible linker group, an unfunctionalized imidazolate ring. The periodic ZIF crystal used as a precursor was created by obtaining the XRD crystal structure of ZIF-8 from the Cambridge Structural Database (CSD), replacing the methyl groups with hydrogen. Bader charge analysis76 was employed to characterize the electron density of the imidazolate rings in both the periodic system and the model fragment computed utilizing the PBE density functional with a plane wave basis (500 eV energy cutoff). The electron density for both the model fragment and the periodic system was computed using the VASP77,78 package. From the periodic system, we constructed a Zn2Im7Li63+ fragment, where two Zn ions were each attached to three outer imidazolate rings, and connected through a common imidazolate ring (see the Supporting Information). The outer six imidazolate rings were capped with Li ions. This system was constructed by cutting out the appropriate fragment from the periodic system and capping the 1894
dx.doi.org/10.1021/jp209335y |J. Phys. Chem. C 2012, 116, 1892–1903
The Journal of Physical Chemistry C
ARTICLE
Table 2. Comparison of Bader Charges on Imidazolate Rings of Periodic System and Zn2Im7Li63+ Fragmenta fragment atom C1
a
periodic
central ring
outer rings (average)
1.062
1.075
1.061
N
1.292
1.305
1.274
N C2
1.292 0.588
1.300 0.596
1.353 0.596
C2
0.676
0.660
0.621
H1
0.024
0.027
0.025
H2
0.180
0.155
0.119
H2
0.180
0.168
0.149
Total
0.594
0.570
0.593
The ring atom types are shown in Figure 2.
dangling nitrogen lone pairs with Li ions to mimic dative ligand bonding. The LiN bond distances were adjusted until the Bader charges on the fragment rings matched those on the periodic system. Shown in Table 2 is a comparison of the Bader charges of the imidazolate rings in the Zn2Im7Li63+ fragment to the charges of the imidazolate rings in the corresponding periodic crystal. Comparing the electron density difference and Bader charges for ring atoms (see the Supporting Information for electron density difference isosurface), it is evident that the Zn2Im7Li63+ fragment does an excellent job of reproducing the electron density of the periodic system in the central region, although differences obviously arise near the terminating Li groups. Finally, a constrained geometry optimization was carried out on the fragment, where all metalorganic bonds were kept fixed. The main purpose for this optimization was to generate physical CH, CC, and CN bond distances and compensate for the finite resolution of the XRD structure. We subsequently carried out approximately 200 SAPT calculations for different configurations of CO2 interacting with the frozen Zn2Im7Li63+ fragment. The configurations used were taken from snapshots of an NVT MD simulation using a box with a CO2 molecule and the Zn2Im7Li63+ fragment. The GROMACS package79,80 was used to generate these conformations. A temperature of 3000 K was utilized to generate sampling in repulsive regions, and the intermolecular interactions were described using the OPLS force field61 for the fragment, and the TraPPE59 force field for CO2; the results do not depend sensitively on the choice of force field due to the high temperature, and our goal is to sample a wide variety of representative intermolecular conformations. The SAPT calculations were extremely time-consuming (45 days per point on 8 processors, without computing dispersion interactions), and are infeasible for more than one type of ZIF. However, the results provide benchmark interaction energies against which we evaluate smaller model systems. The details for the calculations were almost identical to those done in our previous work.67 The density fitting DFT-SAPT (DF-DFT-SAPT) method8191 was used as implemented in the MolPro 2009 package.92 The PBE exchange-correlation functional93,94 was used with an asymptotic correction, which is essential for accurate interaction energies.86,91 Ionization potentials required for the asymptotic correction were computed at the aug-ccpVTZ (AVTZ)/PBE0 level. We used a dimer-centered basis set consisting of Dunning style aug-cc-pVDZ (AVDZ) basis95 for organic atoms and def2-SVP basis96 for metal atoms (AVDZ/def2-SVP). For twelve of these points, dispersion interactions were calculated, and as
Figure 3. Structures of Li2mIm ring (left) and Li2dcIm ring (right).
done previously,67 mid bond basis functions were used, as these are essential to converge dispersion energies.97 Three sites were used for the mid bond functions, with the sites placed between the CO2 molecule and the three closest imidazolate rings, and the locations determined using the method of Podeszwa et al.74 We seek a smaller model system to reproduce the SAPT interaction energies between the Zn2Im7Li63+ fragment and CO2 molecule (and thus the bulk-CO2 interaction). We considered three different rings: an imidazole ring, imidazolate ring (anion), and a Li-capped imidazolate ring (net cation), where the LiN bond distance for this last ring was taken to be the same as that used in the Zn2Im7Li63+ fragment. As before, GROMACS was used to run NVT simulations at 3000 K to generate configurations for each case. Approximately 700 single-point DF-DFTSAPT calculations were carried out using the same methodology as before. In contrast to the large model system, these smaller SAPT calculations took only several processor-hours, thus allowing for many of these calculations to be done for various systems. Using the methodology described in the next section, force fields were developed for each of these three ring systems. We find that the force field generated from the small Liimidazolate (Li2Im) SAPT calculations, in combination with charges for the Zn ions determined using the methodology described below, was able to accurately reproduce the SAPT interaction energies between CO2 and the larger Zn2Im7Li63+ fragment (see the Supporting Information). We thus utilize Licapped imidazolate-type rings in all subsequent calculations to efficiently generate force fields for various bulk ZIFs.
’ FORCE FIELD DEVELOPMENT Based on the above results, we developed corresponding small model systems for ZIF-8 and ZIF-71. For ZIF-8, a Li-capped methyl-imidazolate ring (Li2mIm) was used and for ZIF-71 a Licapped dichloro-imidazolate ring (Li2dcIm) was employed (see Figure 3). As described above, ∼700 single-point DF-DFTSAPT calculations with different CO2 configurations were carried out for each imidazolate-type ring system. We utilized a dimer centered AVTZ/def2-TZVPP basis plus midbond functions, yielding SAPT interaction energies closely approximating the basis set limit. The midbond basis consisted of the same basis functions as before, located at the center of mass of the dimer. Unless explicitly stated, only parameters corresponding to imidazolate ring/substituent atoms were fit, while the CO2 parameters were taken from our previous work.67 Using a similar methodology as before,67 we separately fit each physically distinct term in our force field to the corresponding SAPT term in our calculations. We believe that this approach not only results in a more robust force field with physically accurate parameters, but also allows for greater insight and analysis when such a force field is used in simulations. All fitting in this paper, unless otherwise explicitly stated, utilized a least mean square error approach, in which parameters 1895
dx.doi.org/10.1021/jp209335y |J. Phys. Chem. C 2012, 116, 1892–1903
The Journal of Physical Chemistry C
ARTICLE
were explicitly constrained to exhibit the appropriate symmetry )(ESAPT,comp EMM,comp )2 of the ring, with χ2 = ∑ni W(ESAPT,tot i i i SAPT,comp where Ei is the desired SAPT interaction energy comis the corresponding energy ponent for data point i, EMM,comp i ) component predicted by the force field, and W(ESAPT,tot i is a weighting function for each data point depending on the total interaction energy for that configuration. The weighting function was used to make the fit insensitive to high-energy repulsive points (which are unlikely to be sampled at relevant temperatures) and was of the Fermi-Dirac form WðEÞ ¼
1 expððE μÞ=kB TÞ þ 1
where μ was typically chosen to be 5 mHartree and kBT was chosen to be 1 mHartree (∼300 K). We have taken many steps to reduce the correlation of parameters during force field fitting, which is necessary to ensure the resulting parameters are physically meaningful and transferable. Using physical considerations and obtaining parameters from monomer calculations whenever possible, we rely on SAPT interaction energies only for fitting terms explicitly involving electron overlap (e.g., exponential, Buckingham-like terms, vide infra). Nonetheless, we have found cases where there is coupling between Buckingham coefficients on neighboring atoms when fitting interaction energies. Therefore, while fitting all Buckingham coefficients on functionalized imidazolate rings, we have introduced weak harmonic constraints for carbon, nitrogen, hydrogen, and lithium atoms. The values around which parameters were constrained were determined from corresponding fits to SAPT calculations for a Licapped unfunctionalized imidazolate ring interacting with a CO2 molecule. For each energy component, only four atomic parameters were used, namely those for the four atom types mentioned above. In order to determine the uniqueness of the unfunctionalized imidazolate ring fits, we analyzed the fit of Buckingham terms for the exchange energy. We constructed a covariance matrix, describing the sensitivity of the fit of the data points to percentage changes in the parameters. Small eigenvalues of this covariance matrix indicate that the parameters are under-determined due to coupling, and the corresponding eigenvectors describe the coupling of these parameters. We therefore characterized the change in rms error of the fit with a change in parameter space along the eigenvector of the covariance matrix corresponding to the smallest eigenvalue. We found a 5% increase in rms error with an average change of 2025% of each parameter along this eigenvector. Therefore, we expect these parameters to be qualitatively accurate and physically meaningful for use as atomic base values, and thus all subsequent fits utilized these resulting parameters to center the harmonic constraints for similar atom types. The magnitude of the constraint is very weak, leading to rms errors in the energy fits that were negligibly different for constrained and unconstrained fits. The SAPT interaction energy decomposition is as follows: ð1Þ
ð1Þ
ð2Þ
ð2Þ
ð2Þ
Eint ¼ Epol þ Eexch þ Eind þ Eind exch þ Edisp ð2Þ
þ Edisp exch þ Eδhf To fit the first order exchange repulsion energy, E(1) exch, we utilized a pairwise additive Buckingham functional form ð1Þ
Eexch =
where the index i runs over all atoms on the imidazolate-type ring and j runs over the three atoms of the CO 2 molecule. The exponents,Bij, were determined using the procedure used previously.67 Briefly, initial pairwise exponents, B0ij, were determined using atomic ionization potentials which describe the asymptotic exponential decay of atomic electron densities. Then, along with a scaling factor, λ, such that the prefactors A exch ij B ij = λB 0ij, were fit to the SAPT first-order exchange interacprefactors (and all other tion energies. Note that for the Aexch ij similar Buckingham prefactors used in the other interaction energy were fit, and a standard terms), only atomic ring parameters Aexch ii exch 1/2 = (Aexch combination rule, Aexch ij ii Ajj ) , was used to generate the pairwise parameters. For CO2, the Ajj parameters were taken from our previous work.67 As before,67 a scale factor of λ ≈ 1.2 was found, and the final exponents Bij were used accordingly in all other Buckingham terms for the remaining energy components. The electrostatic term, E(1) pol , was fit using a combination of screened atom-centered point charges and Buckingham-type terms to account for charge penetration. The atom-centered point charges were obtained as follows. First, the electron density of the imidazolate-type ring was computed using the same level of theory as used for the SAPT calculations, and a distributed multipole analysis (DMA)98,99 was carried out using the Molpro package.92 Using the method proposed by Ferenczy,100 atom centered point charges were fit to represent these distributed multipole moments. Fitting atomic charges to DMA generated multipole moments rather than the electrostatic potential (ESP) leads to negligible differences in reproducing the ESP, but has the advantage of resulting in more transferable charges and prevents unphysical charges on “buried” atoms.101,102 The final functional form used to fit the electrostatic SAPT term was qi qj ð1Þ Epol = f1 ðBij , rij Þ þ Aelec ij expð Bij rij Þ rij i, j i, j
∑
∑
For the screening function, the TangToennies functional form103 was used fn ðλ, rÞ ¼ 1 eλr
n
∑ m¼0
ðλrÞm m!
where the parameters required for the screening function were taken to be the exponents used in all Buckingham terms, λ = Bij, as both screening and Pauli repulsion arise from electron cloud overlap. The atomic ring parameters Aelec ii , required to generated parameters, were then fit to reproduce the the pairwise Aelec ij SAPT first-order electrostatic energies. As the penetration terms for electrostatics were attractive in nature, the combination rule elec 1/2 = (Aelec was used. Aelec ij ii Ajj ) (2) The second order polarization/induction energy, E(2) pol = Eind + (2) Eindexch, was fit using an explicitly polarizable treatment for CO2. This explicit polarization was accomplished with the use of a shell model, and is described in our previous work.67 All intermolecular chargecharge interactions, whether static or shell, were screened using the Tang-Toennies function as described above. The screening parameters for the shells were taken to be the same as those used for the corresponding atoms. All intramolecular chargecharge interactions (due to polarization) were screened with a Thole-type function104 as described previously.67 No explicit polarization was used for the ZIF framework. The total functional form used to model the induction energy is
expð Bij rij Þ ∑i, j Aexch ij
ð2Þ
Epol = Ushell þ 1896
∑i, j Aindij expð BijrijÞ
dx.doi.org/10.1021/jp209335y |J. Phys. Chem. C 2012, 116, 1892–1903
The Journal of Physical Chemistry C Here, Ushell is defined as the difference between the energy of the system with optimized shell positions, and that of the corresponding nonpolarizable system. The atomic ring parameters, Aind jj , required for the pairwise Aind ij parameters were fit to reproduce the SAPT second-order polarization/induction energy, E(2) pol . As with electrostatics, it was found that the penetration terms for induction were ind ind 1/2 attractive in nature, and the combination rule Aind ij = (Aii Ajj ) was used to generate the pairwise parameters. Since the shell model can exhibit overpolarization in strong electric fields,105 and the capping Li atoms are formally cations, configurations of CO2 close to Li exhibit induction energies that are incapable of being fit with the shell model; we omit these data points (using a methodology described below) in our fits to yield more physical Aind ij parameters. It is important to note that this problem with the shell model only arises due to the ionic nature of the Licapped imidazolate-type rings that we use to generate parameters for the ZIF systems, and as the ZIF crystals contain only buried cationic centers we do not anticipate these problems in simulations of such systems. (2) (2) The second-order dispersion energy, E(2) VdW = Edisp + Edispexch, was fit using a combination of steps to arrive at an accurate and physically meaningful parameter set; here a na€ive fit directly to the SAPT energies yields ambiguous results due to the longrange nature of dispersion and thus strong coupling between atomic parameters. We thus employ the approach developed by Williams and Stone,106 as implemented in the CAMCASP program,107110 to compute distributed imaginary frequencydependent polarizabilities and obtain atomic C6, C8, and C10 dispersion parameters. For these monomer calculations we used an AVTZ basis with an asymptotically corrected PBE0 exchangecorrelation functional. The linear response employs an ALDAX + CHF kernel, meaning 25% exact HartreeFock exchange was used, consistent with the amount present in the PBE0 functional (see CAMCASP manual for further discussion). This methodology was used to compute atomic C6, C8, and C10 dispersion parameters for CO2 and the Li2dcIm and Li2mIm imidazolatetype rings (our previous bulk CO2 force field67 did not utilize C8 or C10 terms). The pairwise dispersion parameters, Cij0 n were ii0 jj0 1/2 , then generated using the combination rule Cij0 n = (Cn Cn ) where the atomic parameters Cii0 n were those computed using CAMCASP. A notable exception to this rule was for hydrogenheavy atom pairwise parameters. Here, consistent with literature recommendation,109 only atomic C6 parameters were computed for hydrogen atoms, and therefore the nonzero C8 parameter for = (CHyd0 CHeavy0 )1/2, which such pairs was given by CHydHeavy0 8 6 10 is a consistent definition when noting the rank of the polarizability tensors involved in computing each Cn term. As with the induction energy fit, we found it best to remove points from the SAPT data set in which the CO2 molecule was “too close” to either Li atom of the imidazolate type ring. The data points were removed by estimating the dispersion energy of each configuration using the parameters from the CAMCASP calculations; data points were removed if either Li atom of the imidazolate type ring was the main contributor to the dispersion interaction energy. For all other energy components, only Buckingham coefficients were fit, accounting for charge penetration effects. The dispersion energy, in addition to fitting Buckingham coefficients, was also used to fit two scale factors, one of which adjusts the magnitude of the asymptotic C6C10 terms. As the C6C10 treatment of dispersion becomes unphysical in charge penetration regions, we found it best to remove said data points for this fit.
ARTICLE
The final functional form used to reproduce the SAPT dispersion energies is ð2Þ
EVdW =
Cij
expð Bij rij Þ ∑ ∑ fn ðβij , rij Þ nn ∑i, j Adisp ij rij n ¼ 6, 8, 10 i, j
Here, we have included two scale factors such that βij = λdampBij and Cijn = λasympCij0 n . Both of these scale factors, as well as atomic ring parameters, were fit to the SAPT imidazolate-type Adisp ii ring/CO2 dispersion interaction energies and thus are not empirically determined. Similar to the treatments of electro= statics and induction, the pairwise combination rule Adisp ij disp 1/2 was used. Because of the noted sensitivity of (Adisp ii Ajj ) the CO2 isotherms to dispersion parameters, we added a penalty function to the typical least mean square error fitting approach. This penalty function was proportional to the absolute mean deviation of the fitted energies to the SAPT values, and the purpose of this penalty function was to prevent any systematic error in the fit. Also, to ensure that the correct asymptotic behavior was preserved, rather than constructing the χ2 function using the absolute energy differences, we used the percentage energy deviations. This ensured that small dispersion energies were fit accurately; although small in magnitude, these weak interactions contribute significantly to bulk adsorption due to the corresponding increase in the volume element at large distances. Because our previous CO2 model67 used a different treatment for dispersion, we obtained new atomic parameters for CO2 consistent with this new functional form. This was done by refitting the previously computed SAPT CO2CO2 dispersion interaction data67 using the functional form and methodology discussed above. The final dispersion energy fits utilized a scale factor λasymp that was typically ∼1.05, indicating that the Cijn parameters used in the force field were typically 5% greater than those predicted by the CAMCASP calculations (possibly due to basis set size as well as order of polarizability tensors utilized109). The scale factor λdamp used to scale the “exponents” for the TangToennies damping function was typically ∼0.9, which in combination with the original ∼1.2 scaling of the exponents from their asymptotic values, led to damping exponents very close to those predicted by the exact asymptotic atomic densities. The delta HartreeFock terms, Eδhf, have been described previously67,86 and were fit using exponential functions Eδhf =
∑i, j Aδhf ij expð Bij rij Þ
were fit to reproduce the The atomic ring parameters Aδhf jj calculated Eδhf interaction energies, using a slightly different combination rule to generate the pairwise parameters Aδhf ij . The δhf δhf 1/2 , where combination rule used here was Aδhf ij = (|Aii ||Ajj |) δhf δhf the sign of Aδhf ij was determined to be positive if Aii Ajj was positive, and negative otherwise. All fits for the Li2mIm and Li2dcIm ring systems are given in the Supporting Information (sections 3 and 4) as are the complete set of parameters generated by these fits. Finally, we note one possible potential improvement in our force field fitting methodology. For intermolecular configurations that are closer than the potential minimum (i.e., on the repulsive “wall”), the total interaction energy results from cancellation of substantially larger individual interaction energy components. Thus small errors (
