Accurate Treatment of Electrostatics during Molecular Adsorption in

ARTICLE pubs.acs.org/JPCC

Accurate Treatment of Electrostatics during Molecular Adsorption in Nanoporous Crystals without Assigning Point Charges to Framework Atoms Taku Watanabe, Thomas A. Manz, and David S. Sholl* School of Chemical & Biomolecular Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0100, United States

bS Supporting Information ABSTRACT: Molecular simulations have become an important complement to experiments for studying gas adsorption and separation in crystalline nanoporous materials. Conventionally, these simulations use force fields that model adsorbate-pore interactions by assigning point charges to the atoms of the adsorbent. The assignment of framework charges always introduces ambiguity because there are many different choices for defining point charges, even when the true electron density of a material is known. We show how to completely avoid such ambiguity by using the electrostatic potential energy surface (EPES) calculated from plane wave density functional theory (DFT). We illustrate this approach by simulating CO2 adsorption in four metal-organic frameworks (MOFs): IRMOF-1, ZIF-8, ZIF-90, and Zn(nicotinate)2. The resulting CO2 adsorption isotherms are insensitive to the exchange-correlation functional used in the DFT calculation of the EPES but are sensitive to changes in the crystal structure and lattice parameters. Isotherms computed from the DFT EPES are compared to those computed from several point charge models. This comparison makes possible, for the first time, an unbiased assessment of the accuracy of these point charge models for describing adsorption in MOFs. We find an unusually high Henry’s constant (109 mmol/g 3 bar) and intermediate isosteric heat of adsorption (34.9 kJ/mol) for Zn(nicotinate)2, which makes it a potentially attractive material for CO2 adsorption applications.

1. INTRODUCTION Molecular simulations of crystalline nanoporous materials such as zeolites and metal-organic frameworks (MOFs) have come to play an important role in understanding the adsorption and diffusion of molecules in these materials.1-11 These materials are potentially important for many different types of gas adsorption and separation applications, including CO2 capture,12,13 hydrogen storage and purification,5,14 and natural gas sweetening.11,15 Accurate simulation of adsorption isotherms in these materials requires a force-field with multiple components including, at a minimum, (i) electrostatic interactions between the adsorbing molecule and the framework, (ii) dispersion interactions between the adsorbing molecule and the framework, (iii) electrostatic interactions between two adsorbing molecules, and (iv) dispersion interactions between adsorbed molecules. Polarization effects and framework flexibility can also play a role for some materials and adsorbates. The present paper concentrates on the electrostatic interactions between the adsorbing molecule and the framework within the rigid framework approximation. Depending on the material, these electrostatic interactions may be either larger or smaller than the dispersion interactions. To reproduce experimental isotherms, it is necessary to model the largest of these effects accurately. Comparison to experiment may involve further difficulties if the experimental material contains a significant number of defects or has some remaining solvent molecules or ions trapped inside it. r 2011 American Chemical Society

The potential energy of an adsorbed molecule due to the nanoporous framework in simulations of this kind is most commonly represented as a sum of van der Waals (vdW) interactions, typically pairwise Lennard-Jones potentials, and electrostatic interactions. As in many other kinds of molecular modeling, these electrostatic interactions are typically treated by first assigning point charges to atoms in the material and then summing up the pairwise interactions among these point charges. Although they are conceptually convenient, point charge models have important disadvantages. Most importantly, they suffer from the long-standing problem that the continuous electron density of a real material cannot be exactly represented by a finite collection of point charges. For example, in atoms-in-molecules (AIM) approaches, the total electron density, F( B r ), is partitioned between the various atoms of a material, and this allows the electrostatic potential, V( B r ), in the material’s pores to be represented as a distributed multipole expansion. Truncating this expansion at monopole order gives a point charge model in which each atom’s net atomic charge is assigned to its nuclear position. Different methods for partitioning the material’s electron density give different point charges.16-18 This leads to a more subtle problem. Because the vdW interactions that are required for these simulations must be fixed using mixing rules or Received: February 1, 2011 Published: February 28, 2011 4824

dx.doi.org/10.1021/jp201075u | J. Phys. Chem. C 2011, 115, 4824–4836

The Journal of Physical Chemistry C by fitting to experimental data, it is not possible to use calculations for adsorption in nanopores based on point charge models alone to determine which set of point charges is the most (or least) appropriate. If results from one set of point charges are in better agreement with experimental data than another set of point charges, it is not currently possible to exclude the possibility that the latter point charges are more realistic but the vdW interactions used in the calculations are imperfect. A key objective of this paper is to demonstrate that it is not necessary to assign point charges to framework atoms to simulate adsorption isotherms in nanoporous materials like MOFs. Instead, we use the electrostatic potential energy surface (EPES) computed from plane wave DFT calculations, which rigorously includes information from the full electron density of the periodic adsorbing material. A second objective of our work is to examine the accuracy of existing point charge models. Although we anticipate that our DFT-based approach will find widespread use, there will certainly be instances in which using point charges is more convenient. An obvious and important example is the simulation of flexible framework materials,19,20 where it is not practical to perform a plane wave DFT calculation to assess the EPES during every time step of a classical MD simulation. Because our DFT-based method gives an unambiguous description of electrostatic interactions for a given crystal structure, results from this method can be used to judge the validity of approximate point charge methods. Multiple methods are available for assigning point charges from a known electron density, and we examine some of the most widely used below. The Hirshfeld method has been widely used in the simulation of adsorption isotherms for MOFs;21,22 however, it generally underestimates the magnitudes of atomic point charges.23,24 Electrostatic potential (ESP) charge methods minimize the rms electrostatic potential difference between the point charge model and V( B r ) for a grid of points outside the material’s vdW surface. Until recently, ESP charges in periodic materials were obtained from nonperiodic cluster models,17,25 but in calculations of this kind the choice of optimal cluster model is unclear.26 The repeating electrostatic potential extracted atomic (REPEAT) charge method is a recently developed ESP method for assigning point charges in porous periodic materials without building nonperiodic cluster models.27 More chemically meaningful charges are given by the density derived electrostatic and chemical (DDEC) method.18 DDEC simultaneously optimizes the net atomic charges (NACs) to reproduce realistic chemical states and to reproduce V( B r ) from quantum chemistry calculations.18 Like the REPEAT method, the DDEC method does not require nonperiodic cluster models for computing NACs in periodic materials. Efforts have also been made to develop transferable point charge models for families of MOFs without explicitly calculating the electron density for new materials. We examine one approach based on connectivitybased atom contribution (CBAC) charges derived from point charges assigned to clusters representing ∼30 MOFs.28 To illustrate our methods, we examined four MOFs of interest for adsorption and separation applications: IRMOF-1, ZIF-8, ZIF-90, and Zn(nicotinate)2.29-31 These materials have orthogonal lattice vectors and are nonmagnetic. These two features simplify the calculations but are not required for our approach to be applicable. Figure 1 shows the crystal structures of these materials.29-32 Among the four MOF’s examined, we found that Zn(nicotinate)2 has a Henry’s constant, isoteric heat of

ARTICLE

Figure 1. Crystal structures of (a) IRMOF-1 viewed along [100], (b) ZIF-8 viewed along [111], (c) ZIF-90 viewed along [111], and (d) Zn(nicotinate)2 viewed along [100]. Atoms are shown in light blue for C, red for O, gray for Zn, white for H, or dark blue for N.

adsorption, and saturation uptake that make it a potentially attractive material for adsorption-based CO2 separations.

2. METHODS Plane wave DFT calculations were performed using the Vienna ab initio simulation package (VASP) with the projector augmented wave (PAW) formalism.33 An energy cutoff of 500 eV was used in all calculations. MOF geometries were optimized using a conjugate gradient method until the force on every atom reached less than 0.03 eV/Å. Calculations were performed at the Γ point for IRMOF-1, ZIF-8, and ZIF-90 and using a Monkhorst-Pack grid of 1 1 2 k-points for Zn(nicotinate)2. Test calculations indicated that this k-space sampling gave total energies converged to better than 0.01 eV. The exchangecorrelation functionals studied were the generalized gradient approximation (GGA) by Perdew and Wang,34 denoted PW91, and the local density approximation (LDA) by Ceperly and Alder.35 We first compared the DFT-optimized geometry of each MOF using two different exchange-correlation functionals with the reported X-ray diffraction (XRD) crystal structures to check the fidelity of the DFT calculations. The optimized and experimental lattice parameters and unit cell volumes are shown in Tables 1-4. We used the XRD structures from Li et al. for IRMOF-1,29 Wu et al. for ZIF-8,36 Morris et al. for ZIF-90,37 and Rather et al. for Zn(nicotinate)2.32 The DFT results are generally consistent with the XRD structures within a few percent in volume. For IRMOF-1, ZIF-8, and ZIF-90, PW91 (LDA) overestimated (underestimated) the lattice parameters by 1.6-2.9% (2.4-8.8%). These results are consistent with the observation that GGA (LDA) tends to overestimate (underestimate) bond lengths.38-40 For Zn(nicotinate)2, the lattice parameter change 4825

dx.doi.org/10.1021/jp201075u |J. Phys. Chem. C 2011, 115, 4824–4836

The Journal of Physical Chemistry C

ARTICLE

Table 1. Lattice Parameters, Zn-O Bond Lengths, and Unit Cell Volume of IRMOF-1 from Experimental Data29 and DFT-Optimized Structuresa XRD

PW91

LDA

a, b, c (Å)

25.83

26.06 (þ0.9%)

25.63 (þ0.8%)

Zn-O1 (Å)

1.923

1.965 (þ2.2%)

1.913 (-0.5%)

Zn-O2 (Å)

1.941

1.966 (þ1.3%)

1.923 (-0.9%)

volume (Å3)

1724

1770 (þ2.6%)

1683 (-2.4%)

a

Oxygen atoms of different coordination types are represented by O1 and O2. Values in parentheses are the differences from the XRD structure.

Table 2. Similar to Table 1, but for ZIF-8, with Experimental Data from XRD31 XRD

PW91

LDA

a, b, c (Å)

16.99

17.16 (þ1.0%)

16.69 (-1.8%)

Zn-N (Å)

1.982

2.019 (þ1.8%)

1.948 (-1.7%)

volume (Å3)

4904

5050 (þ2.9%)

4650 (-5.5%)

Table 3. Similar to Table 1, but for ZIF-90, with Experimental Data from XRD37 XRD

PW91

LDA

a, b, c (Å)

17.27

17.36 (þ0.5%)

16.75 (-3.0%)

Zn-N (Å)

1.995

2.019 (þ1.3%)

1.958 (-1.8%)

volume (Å3)

5152

5232 (þ1.6%)

4701 (-8.8%)

Table 4. Similar to Table 1, but for Zn(nicotinate)2, with Experimental Data from XRD32 XRD

PW91

LDA 21.04 (-1.5%)

a, b (Å)

21.35

21.14 (-1.0%)

c (Å)

6.918

6.907 (-0.2%)

6.743 (-2.6%)

Zn-N (Å)

2.057

2.040 (-0.8%)

2.012 (-2.2%)

Zn-O1 (Å)

2.056

2.079 (þ1.1%)

2.031 (-1.2%)

Zn-O2 (Å) volume (Å3)

2.344 3154

2.327 (-0.7%) 3086 (-2.2%)

2.335 (-0.4%) 2985 (-5.6%)

was anistropic and the volume change was small. Overall, PW91 reproduced the experimental lattice parameters within ∼3% accuracy. The potential energy of an adsorbed CO2 molecule defined by (C, O, O) atoms at positions B r i (i = 1, 2, 3) was given by 3

VCO2 ð B r iÞ ¼

3

∑ qjVe ð Br jÞ þ j∑¼ 1 ∑k VCO - MOF ð Br j, Br k Þ j¼1 2

þ

∑

1 VCO2 - CO2 2 CO2

where B r k are the position of atoms in the MOF framework. In all r j,rBk) are the dispersion interacof our calculations, VCO2-MOF( B tion defined by the universal force field (UFF),41 and VCO2-CO2 is the EPM2 force field42 for CO2-CO2 interactions. The EPM2 model for CO2 is a point charge model derived to accurately describe CO2 vapor-liquid phase equilibria. In principle, one could use more sophisticated models for CO2 which include its electron density distribution more directly, but this would be

more computationally expensive and is beyond the scope of the present work. The main aim of this paper is to examine different ways of describing the electrostatic potential created by the r j). MOF, Ve( B We tabulated the DFT electrostatic potential energy surface (EPES) on a finely spaced grid in the unit cell of each material. The electrostatic energy of a charged adsorbate atom at an arbitrary position can then be efficiently computed via interpolation from the pretabulated grid. To perform CO2 adsorption simulations, we chose grid sizes of 256 256 256 for IRMOF1, 252 252 252 for ZIF-8 and ZIF-90, and 320 320 100 for Zn(nicotinate)2. These grid sizes are sufficiently fine to provide accurate calculation of the charge density during the plane wave DFT calculation and to give results for the overall EPES via interpolation that are well-converged. The grids of the DFT EPES output from VASP were directly imported into our in-house molecular simulation program used to compute the adsorption isotherms. This DFT EPES is an approximation to the real EPES in the MOFs, because (i) the exact exchangecorrelation potential is unknown and (ii) we used the rigid framework approximation. These computations are based on the Born-Oppenheimer approximation, which assumes that the electronic and nuclear motions are separable. A further approximation is that the PAW method incorporates relativistic effects for frozen-core electrons but not for valence electrons. As more accurate ways to solve for the true electron density within the Born-Oppenheimer approximation become available, our method can be applied to obtain increasingly more accurate EPES for benchmarking point charge models within the rigid framework approximation. Because exactly the same electron density was used to generate the DDEC, REPEAT, and Hirshfeld charges as was used to generate the DFT EPES, our method provides a rigorously correct comparison between these methods. The local electrostatic potential for all point charge models was calculated using Ewald summation.43,44 The decay constant, R, and reciprocal space cutoff, kcut, for the Ewald summation were chosen to simultaneously optimize computational efficiency and accuracy. The values of R and kcut were 0.2, 0.5, 0.5, and 1.0 Å-1 and 2.0, 3.0, 3.0, and 10.0 Å-1, respectively, for IRMOF-1, ZIF-8, ZIF-90, and Zn(nicotinate)2. Point charges for the framework atoms were calculated using the PW91-optimized crystal structures. The REPEAT, DDEC, and Hirshfeld methods were applied to all four MOFs considered. DDEC and Hirshfeld charges were computed using the software available at ddec.sourceforge.net as described by Manz and Sholl.18 Since the DDEC method is based on the charge density distribution, it can provide multipole components. We refer to the atom centered monopole model as DDEC(m) and the dipole corrected model as DDEC(d). Unrestrained REPEAT charges were calculated from the code developed by the original authors of the REPEAT method,27 using the vdW multiplier of 1.3 found to give good results by Chen et al.45 CBAC charges for IRMOF-1 and ZIF-8 were taken from the literature.27,28,46 By definition, the DDEC, Hirshfeld, and REPEAT charges sum over the unit cell to the correct total charge. The CBAC charges, which do not necessarily sum to zero for a neutral unit cell, were corrected to give a neutral unit cell by dividing the sum of uncorrected atomic charges over the unit cell by the number of atoms and subtracting this value from each atom. For example, corrections of -0.012 (0.090) e/atom for IRMOF-1 (ZIF-8) were required for the CBAC charges. 4826



ARTICLE

Table 5. Values of the Atomic Charges for IRMOF-1a REPEAT

a

DDEC(m)

Hirshfeld

CBAC28 0.8095

C1

0.5187

0.7949

0.1759

C2

0.1387

-0.1105

-0.0051

0.0535

C3

-0.1813

-0.0766

-0.0191

-0.1265

Zn

1.2787

1.2149

0.4229

1.5955

H

0.1687

0.1227

0.0379

0.1125

O1

-0.6113

-0.6764

-0.2171

-0.7895

O2

-1.5713

-1.4023

-0.3381

-1.9175

All values are in units of electron charge.

Table 6. Values of the Atomic Charges for ZIF-8a

a

REPEAT

DDEC(m)

Hirshfeld

CBAC28

C1 C2

-0.1013 -0.3623

0.4291 -0.0839

0.0683 -0.0504

0.6460 -0.0620

C3

-0.7617

-0.4526

-0.1023

-0.0410

Zn

-0.9807

0.6920

0.3591

0.6970

H1

0.2367

0.1128

0.0421

0.0100

H2

0.2463

0.1325

0.0391

0.0100

H3

0.2792

0.1306

0.0402

0.0100

N

0.4164

-0.3879

-0.1237

-0.4400


Table 7. Values of the Atomic Charges for ZIF-90a REPEAT

Hirshfeld

C1

0.5232

0.1730

0.0411

C2

-0.1031

-0.0115

-0.0123

C3

0.0692

0.2801

0.0862

Zn

0.5895

0.6695

0.3867

H1

0.1573

0.1123

0.0510

H2

0.0494

0.0380

0.0341

-0.3578 -0.3434

-0.4258 -0.3009

-0.2291 -0.1015

O N a

DDEC(m)


Table 8. Values of the Atomic Charges for Zn(nicotinate)2a

a

REPEAT

DDEC(m)

Hirshfeld

C1

-0.5395

-0.2023

-0.0283

C2

-0.1881

-0.1734

-0.0090

C3

0.1538

0.0432

-0.0005

C4 C5

0.2165 0.5682

0.1363 0.1448

0.0357 0.0360

C6

0.7713

0.7342

0.1496

Zn

1.3097

0.8256

0.3878

H1

-0.3635

0.0852

0.0417

H2

-0.0088

0.0971

0.0425

H3

0.0961

0.1181

0.0466

H4

0.2906

0.1394

0.0471

O1 O2

-0.6375 -0.4394

-0.6429 -0.6115

-0.2485 -0.2372

N

-0.5747

-0.2810

-0.0697


Atomic charges for the four MOFs using the methods listed above are shown in Tables 5-8, with the location of each charge species shown in Figure S1 of the Supporting Information. Here we make two key observations. First, the REPEAT charges for ZIF-8 are drastically different from those of ZIF-90, IRMOF-1, or Zn(nicotinate)2. Specifically, for ZIF-8 the Zn atoms have negative charges, which seems unphysical. In the unrestrained REPEAT method, atomic charges are mathematical constructs optimized to fit the EPES without any specific chemical meaning. We return to this point later in section 3. Second, the Hirshfeld charges are quite small compared to the rest of the charge models, which is consistent with the well-known fact that Hirshfeld charges typically underestimate the magnitudes of net atomic charges.23,24 The accuracy of each point charge model for reproducing the DFT EPES was quantified by the root-mean-square error (RMSE), RMSE ¼ f

∑i ½Vi q þ Vof fset - Vi DFT 2g1=2

ð1Þ

where Viq (ViDFT) is the potential energy at grid point i due to the point charges (DFT EPES), and Voffset equals the average of ViDFT - Viq over the entire set of valid grid points. Voffset was recomputed for each point charge model and for the case when all charges are set to zero. Valid grid points lie between the surfaces defined by γinner and γouter times the vdW radii. For nanoporous materials like MOFs, it is desirable to exclude the volume within 1.3 the vdW radius of each atom (i.e., γinner = 1.3, γouter = limited only by the pore size),45 while for molecular systems it is desirable to set γinner = 1.4 and γouter = 2.0.47 We used the vdW radii of the REPEAT code, which are listed in the Supporting Information (Table S1).48 The relative root-mean-square error (RRMSE) was defined as the RMSE of the point charge model divided by the RMSE when all charges are set to zero. Since this RRMSE is defined to give a value of 1 when all the point charges are set to zero, RRMSE values for different materials can be directly compared. Code for computing RMSE and RRMSE is available at ddec. sourceforge.net. CO2 adsorption isotherms were calculated with grand canonical Monte Carlo (GCMC) simulations.2,6,44 All framework atoms were fixed at their crystallographic positions during these calculations. All GCMC simulations were done at 303 K with 5 106 MC moves for equilibration and 5 106 MC moves for data collection. Since CO2 has a nonzero quadrupole moment,49 both vdW and electrostatic interactions are necessary to correctly describe its adsorption in MOFs. We modeled the vdW interactions using the Lennard-Jones parameters of the UFF for framework atoms,41 the EPM2 force field42 for atoms in CO2, and the Lorentz-Berthelot mixing rules to define CO2-MOF interactions. Although these force fields are not specifically tailored to model MOF-CO2 interactions, they have been found to reproduce isotherms of several gas mixtures in MOFs.2,50,51 The Lennard-Jones interaction cutoff was 17 Å for all elements. The electrostatic interaction was computed by multiplying the charge of each CO2 atom (0.6512 e for C and -0.3256 e for O) by the local electrostatic potential, V, interpolated from the grid as described in the next paragraph. Henry’s constants were obtained from the slope of the linear regression to the data points shown in Figures 2-4 without any constraints. Isosteric heats of CO2 adsorption in MOFs were obtained during GCMC 4827



ARTICLE

Figure 2. Computed CO2 adsorption isotherms at 303 K using the DFT EPES from PW91 and LDA functionals for IRMOF-1, ZIF-8, ZIF-90, and Zn(nicotinate)2. XRD indicates calculations using the XRD structures, while results without this label used the DFT-optimized structure.

simulations using52 qst ¼ RT -

ÆNV æ - ÆNæÆV æ ÆN 2 æ - ÆNæ2

ð2Þ

where R is the gas constant, T is the temperature, N is the number of molecules present, V is the potential energy of the adsorbed molecules, and the angled brackets indicate ensemble averages. During GCMC calculations, the energy of the inserted CO2 molecules was calculated by interpolating the potential energy grid. Because accurate interpolation is essential, we examined the effect of the interpolation scheme. Trilinear interpolation has the desirable property that it cannot introduce any extra (i.e., spurious) maxima or minima in the EPES. Other interpolation schemes like cubic splines53 provide the advantage of a continuous differentiable potential at the cost of potentially introducing spurious maximum or minimum54 in the potential. We tested both trilinear and cubic spline interpolation for calculating the CO2 isotherms in each material. We found several instances where cubic splines introduced spurious maxima when interpolating the DFT EPES, even with very fine meshes. Artifacts of this kind were not observed when using cubic splines to interpolate the potentials associated with point charge models. The DFT EPES is computed in VASP via a fast Fourier transform (FFT) such that the smallest wavelength components are comparable to the grid spacing. As a result, cubic spline interpolation does not work well unless the FFT-generated EPES potential has been numerically smoothed to damp such high-frequency oscillations.

Since these adsorption isotherm computations did not require a differentiable interpolation, trilinear interpolation was clearly the more convenient option and was used for all the results that follow.

3. RESULTS AND DISCUSSION 3.1. Effect of the Exchange-Correlation Functionals. In section 2, we examined the reliability of the PW91 and LDA functionals for optimizing the geometries of four MOF materials. A more critical issue for our calculations is the effect of the exchange-correlation functional on the CO2 adsorption isotherms computed from the DFT EPES. We focused primarily on the low fugacity (i.e., low loading) regime, since the net energy of adsorbed molecules is dominated by adsorbateframework interactions in this regime. A comparison of the resulting CO2 adsorption isotherms from LDA and PW91 for IRMOF-1, ZIF-8, ZIF-90, and Zn(nicotinate)2 at low loadings is shown in Figure 2. Isotherms computed using the reported XRD structures are labeled XRD. Data from the PW91 (LDA) optimized structure with the PW91 (LDA) EPES are simply labeled PW91 (LDA). For every material, the PW91- and LDAbased isotherms using the XRD structures were almost indistinguishable. This important observation indicates that for a fixed framework geometry the DFT EPES is relatively insensitive to the exchange-correlation functional. That is, the DFT calculated EPES can be defined without ambiguity once the geometry of the adsorbent is determined. 4828



ARTICLE

Figure 3. CO2 adsorption isotherms (left) and corresponding Henry’s constant for ZIF-8 at -4%, -2%, 0%, and þ2% lattice strains using the PW91 EPES. 0% strain corresponds to the XRD structure.

Figure 4. CO2 adsorption isotherms at 303 K using the point charge models for IRMOF-1, ZIF-8, ZIF-90, and Zn(nicotinate)2. PW91-optimized structures are used in all of the calculations. Results using the PW91 EPES are included for comparison, denoted by the solid curve.

Tables 9 and 11 list our computed CO2 Henry’s constant values for the four MOFs. These Henry’s constants varied among the four MOFs by 2 orders of magnitude. The high Henry’s constant of Zn(nicotinate)2 for CO2 may make this material interesting for applications involving low CO2 pressures. To the best of our knowledge, these are the first computations of CO2 adsorption isotherms for Zn(nicotinate)2 and no gas adsorption or separation experiments have been reported for this material. Keskin et al. recently reviewed the experimental heats of adsorption and loadings of CO2 in MOFs and found that many of the MOFs with high heats of adsorption did not display high loadings near STP conditions.2 Zn(nicotinate)2 appears to be highly

unusual, because our calculations show that it should display both a high loading near STP (∼4 mmol/g), a high Henry’s constant (∼109 mmol/g 3 bar), and a moderate isoteric heat of adsorption (∼34.9 kcal/mol). This unusual combination of properties may make Zn(nicotinate)2 a potentially attractive material for adsorption-based CO2 separations. Figure 2 also shows the adsorption isotherms computed from DFT-optimized geometries, and it is clear that in some cases these isotherms differ noticeably from the results using XRD structures. Examination of Figure 2 and Tables 1-4 suggests that these isotherm variations are closely related to lattice strain. To test the effect of lattice strain, we performed isotherm 4829



ARTICLE

calculations with strained ZIF-8 using the PW91 functional. In these calculations, the crystal was isotropically contracted by -4% or -2% or expanded by þ2% in terms of its lattice parameters. In each case, the DFT EPES was calculated and used in simulating the CO2 adsorption isotherm with GCMC. Results of these calculations are shown in Figure 3. It is clear that the Henry’s constant for this example is a strong function of the lattice parameter. This effect accounts for essentially all of the variation in the isotherms shown in Figure 2. In particular, the strong lattice contraction seen upon optimizing Zn(nicotinate)2 with LDA accounts for the significantly higher Henry’s constant in this case than in the other calculations for this same material. An open question that cannot be resolved by our calculations is whether it is better to use an experimental structure or a DFToptimized structure when attempting to make quantitative predictions about a specific material. The results above highlight the key contribution of this paper, namely, that it is possible to compute the adsorption isotherm of molecules such as CO2 in MOFs without having to assign point charges to the atoms in the MOF framework. We suspect that this concept will find widespread applications in future simulations of MOFs and other crystalline materials. 3.2. Assessment of Point Charge Models. As described in the Introduction, even though the DFT-derived EPES is a rigorous model of electrostatic interactions in nanoporous crystalline materials, point charges are still necessary to perform molecular simulations that include dynamic framework motions. Thus, it is important to assess the relative accuracy of different point charge models. Morris et al. applied the REPEAT method to ZIF-25, -71, -93, -96, and -97 and compared the results with the experimental CO2 sorption data at 298 K.55 Agreement between the theoretical and experimental isotherms varied significantly over the different materials in their study. Since they did not compare any other charge methods for the same set of materials, the relative quality of the REPEAT method was not clear. Vaidyanathan et al. recently published an experimental and computational study of CO2 adsorption geometry and energetics in the MOF Zn2(3-amino-1,2,4,-triazole)2(oxalate). GCMC simulations they performed using REPEAT charges gave good

agreement with the experimental data, but comparisons to other charge methods were not made.56 We calculated CO2 adsorption isotherms in the four MOFs using the REPEAT, DDEC(m), Hirshfeld, and CBAC charge models. In all calculations, we fixed the MOFs in their DFTPW91-optimized geometry to avoid the influence of structural differences. For this geometry, we take the results from the DFTderived EPES to be the correct representation of the electrostatic potential. An isotherm computed with no electrostatic interactions (i.e., no charges on framework atoms) is also included for each material. Figure 4 shows the CO2 isotherms for all five materials at 303 K in the Henry’s regime. The corresponding Henry’s constants from these calculations are listed in Table 10. The most complete comparison of point charge methods is possible for IRMOF-1, since this MOF has been examined in multiple previous calculations. Neglecting charges entirely underestimates the Henry’s constant by ∼14%. Applying the Hirshfeld charges did not improve the isotherm significantly. The CBAC charges of Zhong and co-workers overestimate the Henry’s constant by 12%. The DDEC(m) and REPEAT charges give the best results among the point charges we tested. This trend agrees with the isosteric heat shown in Table 11 calculated for the Henry’s regime. It is not surprising that the REPEAT charges do the best at reproducing the electrostatic potential, since they are specifically optimized to fit this potential with no other constraints. Previously reported computational and experimental CO2 adsorption isotherms for IRMOF-1 and ZIF-8 are summarized in Table S2 of the Supporting Information.57-64 The experimental Henry’s constants for IRMOF-1 ranged from 0.91 to 1.67 mmol/g 3 bar, while the previously reported theoretical Henry’s constants for this material ranged from 1.36 to 1.68 mmol/g 3 bar. Our PW91 EPES Henry’s constant was 0.83 mmol/g 3 bar, which is slightly below the experimental values. Differences between our Henry’s constant and the previously reported GCMC simulations can be explained by the different simulation conditions. Previous calculations were performed at lower temperatures (293 and 295 K) using the XRD structures. More importantly, the electrostatic charges and vdW interaction models also differed between these calculations: Skoulidas and Sholl used the Hirshfeld charges with EPM2 and UFF for vdW interactions, Walton et al. used no framework charges with the TraPPE and DREIDING force fields, and Yang and Zhong used ChelpG charges with the TraPPE and the OPLS-AA force fields. Differences between the point charge models were significant in ZIF-8. The largest Henry’s constant was found for the PW91 EPES. Hirshfeld and no charges gave the smallest Henry’s constants. The best agreement with the PW91 EPES was obtained with the REPEAT charges. This was followed by the CBAC and DDEC(m) methods. The experimental Henry’s constants for ZIF-8, listed in Table S2 of the Supporting Information, range from 0.60 to 2.98 mmol/g 3 bar. Our PW91 EPES Henry’s constant lies within this range, although with the

Table 9. Henry’s Constants from the Adsorption Isotherms Calculated Using the DFT EPESa DFT Optimized

a

XRD

PW91

LDA

PW91

LDA

IRMOF-1

0.83

0.84

0.85

0.86

ZIF-8

1.75

2.46

1.93

1.91

ZIF-90

3.04

5.13

1.88

1.92

Zn(nicotinate)2

109

157

102

102

All values are in mmol/g 3 bar.

Table 10. Similar to Table 9 but for Adsorption Isotherms Calculated Using Point Charge Modelsa PW91

a

REPEAT

DDEC(m)

Hirshfeld

No Charge

CBAC

IRMOF-1 ZIF-8

0.83 1.75

0.84 (þ2.0%) 1.65 (-6.1%)

0.85 (þ3.0%) 1.44 (-17.9%)

0.73 (-11.8%) 1.34 (-23.7%)

0.71 (-13.6%) 1.35 (-22.9%)

0.92 (þ11.8%) 1.47 (-16.4%)

ZIF-90

3.04

2.73 (-10.3%)

2.54 (-16.3%)

1.30 (-57.1%)

0.91 (-70.0%)

—

Zn(nicotinate)2

109

107 (-1.8%)

114 (þ3.8%)

112 (þ2.8%)

108 (-0.9%)

—

Percentage differences relative to the values from the PW91 EPES are shown in parentheses. 4830



ARTICLE

large variation in the experimental data this is not a strong statement. The only available computed CO2 isotherm by PerezPellitaro et al. is 0.60 mmol/g 3 bar,59 from calculations that used modified UFF parameters for vdW interactions and ESP charges. Among all the MOFs considered, the effects of electrostatic interactions were most pronounced in ZIF-90. Similar to ZIF-8, the PW91 EPES gave the largest Henry’s constant in ZIF-90 and the REPEAT charges showed the best agreement with the DFT Table 11. Isosteric Heats of Adsorption from the Isotherm Calculations in the Henry’s Regime for Each MOFa PW91 REPEAT DDEC(m) Hirshfeld No Charge CBAC

a

IRMOF-1

13.8

13.7

14.0

12.8

12.5

14.6

ZIF-8 ZIF-90

19.2 22.2

19.1 21.7

18.4 21.4

18.1 18.6

18.0 17.6

18.6 —

Zn(nicotinate)2 34.9

34.8

35.1

34.9

35.2

—

All values are given in kJ/mol.

Table 12. RRMSE (dimensionless) and RMSE (kJ/mol, in brackets) between V from Different Point Charge Models and the PW91 EPES for a Uniform Grid of Points Outside the Surface Defined by 1.3 vdW Radii REPEAT DDEC(d) DDEC(m) Hirshfeld No Charge CBAC IRMOF-1 ZIF-8 ZIF-90 Zn(nicotinate)2

0.1 [1.3] 0.2 [1.7] 0.1 [3.2] 0.1 [1.1]

0.2 [2.3] 0.4 [3.4] 0.1 [4.0] 0.2 [1.9]

0.4 [4.1] 0.5 [3.9] 0.1 [4.0] 0.5 [4.7]

0.6 [6.3] 0.9 [7.2] 0.5 [15] 0.7 [6.2]

1.0 [11] 1.0 [7.7] 1.0 [32] 1.0 [9.3]

1.0 [11] 2.1 [17] — —

EPES. This was followed by DDEC(m). The Hirshfeld and no charge models significantly underestimated the Henry’s constant. The structures of ZIF-8 and ZIF-90 are similar, except ZIF90 contains an aldehyde side group instead of the methyl side group in ZIF-8. The O atom of the aldehyde group in ZIF-90 is more charged than the H atom of the methyl group in ZIF-8, which gives rise to stronger electrostatic interactions in the pores of ZIF-90 than in ZIF-8. This resulted in a Henry’s constant approximately 75% higher for ZIF-90 than for ZIF-8. Both the REPEAT and DDEC(m) methods accurately captured this increase in Henry’s constant, while the Hirshfeld and no charge methods incorrectly predicted a slightly lower Henry’s constant for ZIF-90 than for ZIF-8. For the no charge case, the Henry’s constant in ZIF-90 was only 67% of that for ZIF-8 due to the slightly larger lattice parameter and cell volume. To the best of our knowledge, no experimental or computational CO2 adsorption isotherms for ZIF-90 have been previously reported. The highest Henry’s constant among the four MOFs was observed for Zn(nicotinate)2. The isotherms and Henry’s constant values for different point charge methods were virtually identical for this material. The differences in the Henry’s constants for different charge methods and the DFT EPES were

Accurate Treatment of Electrostatics during Molecular Adsorption in

Recommend Documents