Gas Membrane Selectivity Enabled by Zeolitic ... - ACS Publications

Jun 2, 2014 - ... Jessica M. R. Burton†, Yao Houndonougbo⊥, Brian B. Laird#, and ... Hiroyasu Furukawa , Margarita Russina , Alberto Albinati , Pe...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/cm

Gas Membrane Selectivity Enabled by Zeolitic Imidazolate Framework Electrostatics Keith G. Ray,*,†,‡ David L. Olmsted,†,‡ Jessica M. R. Burton,†,§ Yao Houndonougbo,⊥ Brian B. Laird,# and Mark Asta† †

Department of Materials Science and Engineering, University of California, Berkeley, California, United States, Department of Chemistry and Biochemistry, Eastern Washington University, Cheney, Washington, United States # Department of Chemistry, University of Kansas, Lawrence, Kansas, United States ⊥

ABSTRACT: CO2 and CH4 diffusivities are calculated using molecular-dynamics simulations for a set of five RHO topology zeolitic imidazolate frameworks (ZIF-25, -71, -93, -96, and -97), which differ only in linker functionalization. It is found that CO2 diffuses more readily than CH4, by a factor of 3.25, through the channels of ZIF-93, despite being much more strongly bound to the framework. Using adsorption selectivity values from previous studies, we calculate the membrane selectivities at 1 bar for all five materials and find this quantity to be highest in ZIF-93, with a value of 26.7. First-principles calculations based on van der Waals density functional methods are used to investigate the transport barriers for CO2 in ZIF-93 and ZIF-97. These two materials have similar channel diameters, but have very different CO2/CH4 diffusion selectivities. The excellent molecular-sieving performance of ZIF-93 is determined to be enhanced by electrostatic interactions, that result in smaller energy barriers for CO2 transport than in ZIF-97. Thus, for channels small enough to act as molecular sieves, the electrostatic field in the channel is found to be an important factor in the design of ZIFs for the separation of gas molecules possessing non-negligible electrostatic interactions with the framework. CH4 CO2 2 where SD = DCO /CCH4 are the ratios of trans/Dtrans and SA = C CO2 and CH4 transport diffusivities and concentrations, respectively.8−10 The diffusivities are both given at the concentrations, CCO2 and CCH4, and the concentrations are both given at partial pressures, pCO2 and pCH4 of the feed gas. Individually, we will call SD the diffusion selectivity and SA the adsorption selectivity. In eq 1, we have assumed that pCO2 = pCH4 and that the downstream loadings are vanishingly small. In addition to selectivity, a good membrane material for naturalgas upgrading should display reasonably high permeabilities for CO2, defined as

1. INTRODUCTION The burning of natural gas for power generation produces less of the greenhouse gas (GHG), CO2, per unit energy than the burning of many other fossil fuels, including coal and gasoline.1 Furthermore, the capture and use of biogas and landfill gas for energy prevents the release of methane, a more potent GHG than CO 2 , into the atmosphere. 2,3 Consequently, the processing of natural gas, biogas, and landfill gas into a form useful for energy generation is important for the reduction of GHG emissions in any foreseeable energy portfolio. The separation of CO2 from CH4, a crucial step in that processing, increases gas energy density and reduces acidity.4 Emerging membrane technologies are particularly attractive for gas separations because these systems potentially consume less energy and are more compact than alternative methods.3,5,6 Metal−organic framework (MOF) materials, including zeolitic imidazolate frameworks (ZIFs),7 have been actively investigated for membrane-based gas separation applications. These materials display a rich variety of chemistries and structural topologies depending on the type of metal atom and organic linker. There is thus a large parameter space available to optimize the performance of these materials for membrane separations. Typically, this performance is measured by the membrane selectivity, which in its simplest form is given for CO2/CH4 separations as S = SDSA

CO2 CO2 CO2 P CO2 = Dtrans C /f

where f CO2 is the fugacity, for which we will use the partial 2 pressure, pCO2, for our purposes, DCO trans is the CO2 transport CO2 diffusivity, and C is the concentration of CO2 per unit volume of the membrane.11 MOFs possess a regular array of well-defined pores and channels, which enable larger permeabilities and may enable larger selectivities than many polymer matrices.3,12,13 Specifically, MOFs may act as molecular sieves, in which the diameters of the channels are small enough to restrict the passage of larger gas molecules while allowing the transport of smaller gas molecules (steric Received: April 29, 2014 Revised: May 30, 2014

(1) © XXXX American Chemical Society

(2)

A

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

separation).13 In the following we show that this picture is incomplete and that the detailed electrostatic field in the channel, which is influenced by both channel morphology and chemical composition, can significantly affect the membrane selectivity and thus should be a design consideration for future MOF membranes. Previous experimental and simulation studies have been performed on ZIFs and MOFs for the diffusion and permeability of CO2 and CH4, as well as other gases. Several diffusion measurements have been made on ZIF-8 alone,8,14 showing only a small diffusion selectivity, around 1.5−2, for CO2 over CH4 at atmospheric pressure. MOF-177 has a diffusion selectivity of only 1.6 for CO2 over methane and CH4 has been found to have a higher diffusivity than CO2 in MOF5.15 In these materials, the promise of high membrane selectivity, eq 1, must come from the adsorption selectivity contribution, SA. Simulations on CO2/CH4 diffusion in MOFs and ZIFs largely show a similar pattern of low diffusion selectivity. Molecular-dynamics simulations on MOF-5,16 74 (both Zn and Mg versions),11 and 177,11 as well as ZIF-2,17,18 -3,18,19 -6,18 -10,18,19 -60,18 -65,18 -68,20 -69,18 -70,20 -79,18 and -8118 all show CO2/CH4 diffusion selectivities between 0.07 and 1 at pressures between 1 and 10 bar. In contrast, adsorption selectivities for these materials are between 1 and 20 at the same pressure. Higher CO2 diffusivities compared to CH4 have been calculated from simulation in ZIF-4,17 -8,17 -11,21 -12,21 -69,6 -77,6 and -90.6,18 However, calculations from refs 6 and 17 lack framework charges, which are important to the framework interaction with CO2, due to this molecule’s quadrupole moment.22 A CO2/CH4 diffusion selectivity of 2.34 is found in ZIF-90 in ref 18 along with an adsorption selectivity of 9.53 at 10 bar, making this potentially one of the better performing ZIFs for CO2/CH4 membrane separation applications. Also, an extraordinary diffusion selectivity of greater than 100 is predicted for ZIF-11 and -12 in ref 21; however, both of these materials have adsorption selectivities of less than 3. Thus, high adsorption selectivity is not found along with high diffusion selectivity in any of the materials reviewed above. Only ZIF-90 combines a diffusion selectivity greater than 2 with an adsorption selectivity greater than 5, within simulations taking electrostatic interactions into account. In studies to date, enhanced diffusion selectivity is often anticorrelated with adsorption selectivity. As one specific example, MOF-74 can be made with both Zn and Mg metal centers. The Mg version has greater CO2/CH4 adsorption selectivity11 through stronger electrostatic interactions between the metal center and the CO2 molecule,23 but the Zn version has higher diffusion selectivity.11 In the following, we examine CO2 and CH4 diffusion in a set of five RHO topology ZIFs. For membrane separation applications ZIFs are of interest relative to other classes of MOFs because of their chemical stability.6,12,24,25 The ZIFs considered in this work share the same RHO topology and differ only in the functionalization on the 4 and 5 sites of the imidazolate linker. Simulation and experimental studies have focused on gas adsorption on the same set of ZIFs,26−29 as well as other sets differing only in functionalization17,18,30−35 and topology.36 We focus on this set of ZIFs to examine the effect that the variation of the chemical functional group has on the diffusion selectivity. In particular, we find and discuss a material that combines good adsorption and diffusion selectivity, ZIF93, and discuss the chemical factors underlying these properties. The electrostatic field in the framework channels

of this material is shown to be critical to its CO2/CH4 diffusion and membrane selectivity, as well as CO2 permeability. The paper is organized as follows. In Section 2, we describe the materials studied and the computational methods employed. Section 3.1 contains results from moleculardynamics (MD) including the tracer and transport diffusivity for CO2 and CH4 in each of the five ZIFs and the hopping rates between the pores present in the RHO topology ZIFs. The effect of a flexible framework on diffusion barriers in these materials is analyzed in Section 3.2 utilizing the climbing-image nudged-elastic-band (cNEB) method37 in combination with van der Waals density functional (vdW-DF) methods.38−40 In Section 3.3, we apply the same methods to compare the CO2 transport barriers through the ⟨001⟩ channel in ZIF-93 and ZIF-97. In Section 3.4, we present membrane selectivities and permeabilities calculated utilizing the diffusivities from this work and adsorption isotherms from the literature.26,28 The results and insights provided from this work are summarized in Section 4.

2. MATERIALS AND METHODS 2.1. ZIF Structures. We examine a series of five RHO topology ZIFs, each with a different set of chemical functional groups bound to the 4 and 5 sites of the imidazolate ring. The chemical functionalities, depicted in Figure 1, are as follows: two −CH3 groups (ZIF-25), two −Cl atoms (ZIF-71), one −CHO and one −CH3 group (ZIF-93), one −CN and one −NH2 group (ZIF-96), and one −CH2OH and one −CH3 group (ZIF-97). The latter three ZIFs, -93, -96, and -97, each have two different functional groups on each imidazolate linker. Channel diameters for the two distinct channels in each of the five RHO topology ZIFs are presented in Table 1. The channel diameters were determined by the largest cylinder that will fit in the channel without coming within a Lennard-Jones σ/2 parameter of a framework atom. This measurement is also called the pore aperture size.41 These channels correspond to transport through the double ring of 8 zinc atoms and single ring of six zinc atoms along the ⟨100⟩ and ⟨111⟩ directions, respectively, that connect the Linde Type A cavities. Also listed are the minimum and maximum widths, with respect to orientation, of the CO2 and CH4 molecules, as determined by the Lennard-Jones σ parameters used in our MD simulations. 2.2. Molecular Dynamics Simulations. MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)42 package starting with configurations calculated with the grand canonical Monte Carlo (GCMC) code, Towhee.43 The configurations were chosen such that the number of gas molecules is equal to the average loading computed by the GCMC simulations at a pressure of 1 bar. The force fields employed in these simulations are parametrized by Lennard-Jones (L-J) and Coulomb potentials. The L-J parameters and charges for CO2 were taken from the Elementary Physical Model 2 (EPM2) of Harris and Young.44 CH4 was treated as a united atom at the carbon position using the Transferable Potentials for Phase Equilibria-United-Atom (TraPPEUA) force field.45 The L-J parameters for the framework were chosen from the universal force field (UFF)46 and optimized potentials for liquid simulations (OPLS)47 force field sets, as listed in ref 28 where they were validated by a comparison between simulated and experimentally measured adsorption isotherms. Lorentz−Berthelot mixing rules were utilized. The framework charges were derived using the REPEAT algorithm, involving a fitting of DFT-calculated electrostatic potentials in open regions of the framework, as described in an earlier publication.26 The tracer and corrected diffusivities were calculated by a linear least-squares fit, from 100 to 900 ps, to the mean square displacement and mean square of N times the center of mass displacement, respectively, in order to capture the diffusive regime tDtracer + B ≃ B

1 ⟨||ri(t + t 0) − ri(t 0)||2 ⟩ 6

(3)

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

Table 1. Diameters for the Channels in the ⟨001⟩ and ⟨111⟩ Directions for Each RHO ZIF ZIF

N

1 ⟨||∑ ri(t + t0) − ri(t0)||2 ⟩ 6N i = 1

(4)

The transport diffusivity was calculated using the following standard relation:

⎛ ∂ln f ⎞ ⎟ Dtrans = Dcorr⎜ ⎝ ∂ln c ⎠T

⟨111⟩ channel (Å)

4.55 5.40 3.58 5.69 3.65 min L-J width (Å)

2.91 4.04 3.05 4.62 3.08 max L-J width (Å)

CH4 CO2

3.73 3.03

3.73 5.35

obtained the thermodynamic factor, ((∂ln f)/(∂ln c))T, by fitting concentration as a function of chemical potential from the Monte Carlo simulations described in refs 26 and 28. In the MD simulations, we made use of the SHAKE algorithm48 to fix the bond lengths of the CO2 molecule, while keeping the bond angle as flexible. The ZIF coordinates were held rigid at their experimentally derived positions.26 All results reported below were obtained from simulations in an NVT ensemble, using a chain of 5 Nose−Hoover49,5051 thermostats. The simulation cells employed in these MD runs have dimensions of 2 × 2 × 1 ZIF unit cells. To test the sensitivity of the calculated diffusivities to the details of the simulation setup, we performed additional calculations with larger simulations cells (4 × 4 × 4) and different equations of motion (NVE), finding effects on the order of 10% or smaller. Statistical errors in the tracer diffusion coefficients are estimated to be less than 6%. 2.3. van der Waals Density Functional Calculations. To complement the classical simulations described above, and investigate the potential role of framework flexibility on gas transport in the ZIFs studied, we have undertaken additional quantum-mechanical calculations based on the formalism of van der Waals density functional (vdW-DF) theory. Specifically, transport barriers in ZIF-93 and ZIF-97 were calculated using two versions of the vdW-DF:38 vdW-DF239 and vdW-optB88.40 In both methods, the exchange correlation energy is divided into the local correlation energy, a semilocal exchange energy, and a nonlocal correlation energy.38 The nonlocal correlation energy gives rise to the attractive dispersion interaction.38 These two methods differ in a single parameter in the nonlocal correlation kernel39 and the choice of semilocal exchange functional.39,40 The Vienna ab initio simulation package (VASP)52−56 was used to perform the vdW-DF calculations. For these calculations we employed the projector-augmented wave (PAW) method57,58 with Perdew− Burke−Ernzerhof (PBE)59 based potentials from the VASP library. The brillouin zone sampling contained a single k-point (Γ) due to the large primitive cell and the plane-wave basis set cutoff is 400 eV. The vdW-DF nonlocal correlation energy is calculated using the computationally efficient algorithm of Roman-Perez and Soler.60 We used the van der Waals corrected density functionals to compute transport barriers within the climbing image Nudged Elastic Band (cNEB) scheme.37 In this approach a minimum energy path was found, between two CO2 positions for our purposes, using intermediate images effectively connected by springs. The climbing modification caused the highest energy image to converge toward the saddle-point. We used 4 intermediate images for these calculations and the activation energies are estimated to be converged to better than 5 meV. The cNEB method has been used in calculations of gas transport barriers for a few MOFs and ZIFs previously, including CO2, H2, and H2O in MOF-74-Mg.61 In order to quantify the effect of framework flexibility on the barrier through the ⟨001⟩ channel for CO2, we compared calculations allowing for full framework flexibility with those holding the framework atoms rigid. In such comparisons an important consideration is the fact that the framework coordinates relaxed with vdW-optB88 or vdW-DF2 differ from the experimental atomic positions. Therefore, we calculated the barriers for three cases with each density functional. These are the rigid framework with atomic

Figure 1. (Top) Representation of the ZIF-93 molecular structure viewed along the ⟨001⟩ direction. In the center of this top image is the ⟨001⟩ channel, which is surrounded by a ring of eight zinc atoms, shown as gray tetrahedra. At the corners of this top image are rings of six zinc atoms that encircle the ⟨111⟩ channels. (Middle) Same structure viewed along the ⟨111⟩ direction, illustrating the geometry of the ⟨111⟩ channel in its center. (Bottom) ZIF-25, 71, 96, and 97 imidazolate linkers with atoms labeled.

tDcorr + C ≃

⟨001⟩ channel (Å)

25 71 93 96 97 gas

(5)

In the above equations, t is time, B and C are constants, ri is the centerof-mass position of gas molecule i, N is the number of gas molecules, f is the fugacity, c is the concentration, and T is the temperature. We C

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

Table 2. Molecular-Dynamics Calculated Values of the Tracer Diffusion, Dtracer, Corrected Diffusion, Dcorr, and Transport Diffusion Coefficient, Dtrans, in Units of 10−9 m2 s−1; Concentration, C, from refs 26 and 28, in kmol m−3; Permeability, Where P = CDtrans/Pressure, in 104 barrers; Adsorption Selectivity, SA, for an Equimolar Mixture of CO2 and CH4 as Derived from the IAST in Previous Work,28 is Reproduced Here along with the Ratio of CO2/CH4 Transport diffusivities, SD; Membrane Selectivity Is Given by S = SASDa CO2 ZIF 25 71 93 96 97

Dtracer 2.3 3.7 1.2 3.0 0.18

± ± ± ± ±

0.06 0.11 0.06 0.09 0.01

2.5 3.6 1.4 3.5 0.17

CH4

Dcorr

Dtrans

C

P

± ± ± ± ±

4.5 4.5 1.8 4.4 0.22

0.94 0.72 1.59 2.08 1.04

12.4 9.6 8.3 27 0.68

0.13 0.25 0.08 0.19 0.016

Dtracer 3.7 6.1 0.40 7.6 0.79

± ± ± ± ±

0.11 0.17 0.02 0.15 0.04

CO2/CH4

Dcorr

Dtrans

C

P

SD

SA

S

± ± ± ± ±

4.1 8.0 0.54 9.5 0.88

0.46 0.32 0.42 0.37 0.29

5.7 7.5 0.67 10.4 0.76

1.08 0.56 3.25 0.46 0.25

2.5 2.7 8.2 10.2 6.1

2.7 1.5 26.7 4.7 1.5

3.7 6.5 0.42 7.5 0.76

0.21 0.27 0.03 0.43 0.05

a

All quantities are obtained from simulations for which the gas pressure is 1 bar and the temperature is 298 K. The statistical errors shown are the expected error in the mean.

Figure 2. CO2 minimum energy positions in the ⟨001⟩ channel of ZIF-97. The leftmost image shows the double 8-ring of Zn atoms, the middle image shows only the functional groups closest to the bound CO2 molecules, and the rightmost image shows a side view. The two horizontal CO2 molecules in the rightmost image are bound at equivalent binding sites, which we call the 8-ring site. The vertical CO2 in the middle position of the ⟨001⟩ channel (shown in all three images) is more weakly bound, but also at a local minimum in the potential energy, which we call the 8-ring mid site. positions set to experimental coordinates (refined to X-ray measurements), the rigid framework corresponding to vdW-DF2 or vdWoptB88 relaxed coordinates (forces relaxed to better than 0.01 eV/Å), and the flexible framework described within vdW-DF2 or vdWoptB88. In all cases, the experimentally measured value for the unit-cell volume was used. Finally, in order to investigate finite-temperature effects on gas transport through the ZIF structures considered, we performed an ab initio molecular dynamics (AIMD) simulation on ZIF-93 at 300 K, in order to analyze the distribution of the ⟨001⟩ channel diameters, as described in Section 3.2. The simulations were performed in a canonical ensemble using a Nosé−Hoover thermostat49,50 with a temperature of 300 K; the time-step was chosen as 1 fs, and hydrogen atoms were replaced by deuterium to lower vibrational frequencies. The total simulation time was 1.2 ps after a 0.4 ps equilibration.

these ZIFs are already known to have a significant CO2/CH4 adsorption selectivity,27,28 making them good candidate membrane separation materials. Next, we quantify the gas transport through each of two symmetry distinct channels in the RHO topology ZIFs. The channels which are found to be important for CO2/CH4 diffusion selectivity are examined more closely in later sections with first-principles cNEB calculations using the vdW-DF2 and vdW-optB88 methods. 3.1. CO2 and CH4 Diffusivities. Calculated tracer, corrected, and transport diffusivities for CO2 and CH4 obtained from MD simulations are presented in Table 2 along with the ratios of transport diffusivities. A correlation is observed between the channel diameters, listed in Table 1, and the overall diffusivities for both gases. For instance, the CH4 diffusion coefficient is highest in ZIF-25, -71, and -96, which also have the widest channels. CH4 diffusion is smallest through ZIF-93, which has the smallest channels. ZIF-93 stands out as being the only one of the frameworks considered for which the tracer diffusion coefficient for CO2 is larger than that for CH4. The origin of this effect is not clear from the channel widths shown in Table 1 alone. Specifically, the larger size of CH4

3. RESULTS AND DISCUSSION An ideal membrane separation material will both preferentially adsorb and enable rapid transport of the gas that is to be separated. However, these are often conflicting requirements, as described in the introduction.62 In the following, we first identify members of the set of RHO topology ZIFs considered here for which CO2 has a higher diffusivity than CH4. Some of D

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

the hopping rate for CO2 through the ⟨001⟩ channel in ZIF-97 is slower by at least a factor 6 and as much as a factor of 123, when compared to the other RHO topology ZIF materials considered in this work. It was found, both in classical simulations and in vdW-DF calculations, that there are two equivalent binding sites for each ⟨001⟩ channel of ZIF-97, as illustrated in Figure 2 by the two horizontally oriented CO2 molecules in the rightmost image. These are the 8-ring binding sites. Also shown in Figure 2 is a third minimum energy position for CO2 in the double 8-ring ⟨001⟩ channel of ZIF-97, that we call the 8-ring mid site. CO2 is more weakly bound here, than in the other two degenerate 8ring sites. However, the presence of this binding position allows us to calculate the barrier for transport through the ⟨001⟩ channel as that for transport from a single 8-ring site to the 8ring mid site. As described in the methods section, we applied the cNEB method to the rigid framework with experimentally measured positions (RFE), the rigid framework with calculated coordinates (RFC) relaxed at constant volume with vdWDF2 or vdW-optB88, and the flexible framework with calculated coordinates (FFC) under vdW-DF2 or vdWoptB88, also at constant volume. The atomic coordinates surrounding the ⟨001⟩ channel, for each of these systems, is depicted in Figure 3. The adsorption energies in the 8-ring site

relative to CO2 might be expected to lead to slower diffusivities for the former. To understand the diffusion results in more detail we have undertaken an analysis of the trajectories of the molecules as they hop from one pore in the framework to another. The analysis establishes that the primary channels for gas transport between Linde Type A pores lie in the ⟨111⟩ directions through the rings of six Zn atoms and in the ⟨001⟩ directions through the double rings of 8 Zn atoms. These channels and rings of Zn atoms are illustrated in Figure 1. The ⟨001⟩ channel of ZIF-97 is shown in more detail in Figure 2, from two different perspectives and with the most relevant linkers highlighted. By tracing the movement of gas molecules through two or more pores, the relative importance of transport through the two distinct channels has been quantified. These data are presented in Table 3, which gives the number of jumps through each type of channel per molecule per ns, for both CO2 and CH4. Table 3. CO2 and CH4 Jumps through Channels in the ⟨001⟩ and ⟨111⟩ Directions, Per Molecule, Per ns at a Temperature of 298 K, Obtained from Simulations Lasting a Total of at Least 120 nsa CO2

CH4

ZIF

⟨001⟩

⟨111⟩

⟨001⟩

⟨111⟩

25 71 93 96 97

1.4 2.2 0.69 1.8 0.036

0.055 0.40 0.029 0.61 0.081

2.2 3.6 0.24 4.4 0.42

0 0.61 0.016 1.6 0.013

a The ⟨001⟩ and ⟨111⟩ directions correspond to transport through the zinc double 8-ring and zinc 6-ring respectively, as shown in Figure 1.

For almost every combination of framework and gas, transport through the ⟨001⟩ channel dominates over transport through the ⟨111⟩ channel. In particular, this is true of ZIF-93, and we can see that the CO2/CH4 selective diffusion in ZIF-93 is mostly governed by transport through the ⟨001⟩ channels. Despite having similar ⟨001⟩ channel diameters, the CO2 jump frequency in this channel in ZIF-93 is higher than that in ZIF-97, by a factor of 19. However, the CH4 transport through this channel is higher in ZIF-97 than ZIF-93. Differences in CO2 transport through this channel are responsible for the large diffusion selectivity of ZIF-93 compared to ZIF-97 (3.25 vs 0.25), which results in a large effect on the overall membrane selectivity. 3.2. vdW-DF Calculations for Flexible Frameworks. All of the MD results described in the previous section were obtained from simulations where the framework was held rigid. We will first assess the effect of the rigid framework approximation using the cNEB method to calculate the energy barrier through a channel, with and without framework flexibility. Because these calculations are relatively expensive computationally, we limit the analysis to one gas molecule through one channel. The system is chosen to represent the largest possible effect of flexibility, i.e., transport in the material and gas combination with the smallest tracer diffusivity. The tracer diffusivity is smallest for CO2 through ZIF-97. This is the unique gas and RHO framework combination, in this study, for which the diffusion through the ⟨001⟩ channel is smaller than diffusion through the ⟨111⟩ channel. Specifically,

Figure 3. CO2 in double 8-ring ⟨001⟩ channel of ZIF-97 and 93. Top left: ZIF-97 experimental coordinates (RFE), bottom left: ZIF-97 coordinates relaxed with vdW-DF2, but held rigid in the presence of CO2 (RFC), bottom right: ZIF-97 coordinates relaxed with vdW-DF2 and kept flexible in the presence of CO2 (FFC). Note the positions of the H atoms in ZIF-97 nearest the CO2 molecule. Top right: ZIF-93 with experimentally measured atomic coordinates.

and the 8-ring mid site, and the barrier height from the former to the latter, calculated using the vdW-DF2 and vdW-optB88 functionals for these three different cases, are presented in Table 3. Adsorption energies are typically larger in magnitude for the vdW-optB88 method than for the vdW-DF2. Also, the RFC atom positions bind CO2 more strongly than the RFE atom positions for both functionals. The effect of flexibility on E

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

the adsorption of CO2 can be seen by comparing the adsorption energies for RFC and FFC structures. We see that the −415 to −567 meV adsorption energies are enhanced through flexibility by 8 to 20 meV for the two sites with the two methods. The energy barriers, calculated with cNEB, are presented on the rightmost column of Table 4. Also, the energy vs reaction Table 4. Adsorption Energies (AE) for CO2 in the 8-Ring Site and 8-Ring Mid Site in ZIF-97 As Well As the Energy Barrier into the 8-Ring Mid Site from the 8-Ring Site, in meV, Utilizing the vdW-DF2 and vdW-optB88 Methodsa ZIF

method

framework

flexible?

8-ring AE

mid AE

barrier to mid

97

vdW-DF2

97

vdWoptB88

experiment calculated calculated experiment

no no yes no

−463 −532 −552 −502

−399 −415 −426 −438

243 166 133 237

calculated calculated experiment

no yes no

−567 −582 −409

−470 −478 −332

137 118 96

93

vdW-DF2

Figure 4. Calculated vdW-DF2 cNEB results for CO2 transport through the ⟨001⟩ channel in ZIF-97 are plotted for three different framework systems, all referenced to the 8-ring binding energy. the framework systems include the Rigid Framework with Experimentally measured positions (RFE), the rigid framework with calculated coordinates (RFC) relaxed at constant volume with vdW-DF2, and the flexible framework with calculated coordinates (FFC) under vdWDF2, also at constant volume.

Different framework coordinates and their flexibility are considered. Also, the barrier between the 8-ring site and 8-ring mid site in ZIF-93 is calculated for the framework with experimentally measured coordinates. a

results suggest an increase by a factor of 5.4 in the hopping rates due to the effect of framework flexibility on the energetic and entropic contributions to k+TST. The calculated rates above are for CO2 transport through the ⟨001⟩ channel in ZIF-97. From the results in Table 3, if we increase the value of the CO2 jump frequency in the ⟨001⟩ direction by a factor of 5.4 and keep the jump frequency in the ⟨111⟩ direction constant, the tracer diffusivity would be expected to increase by a factor of ∼2.1. Considering the values in Table 2, such a change would not cause this tracer diffusivity to surpass values for other gas-ZIF combinations. So, if the effect of flexibility calculated here is an upper bound for these systems, this analysis based on harmonic transition-state theory suggests that framework flexibility would not change significantly the trends in the results obtained by the classical MD simulations, as reported in the previous section, i.e., the relative ordering of the jump rates for the different channels for the different gases and ZIFs considered in this work. The results from harmonic TST involved comparing the lowest energy transport trajectory for the entire frameworkCO2 system to the transport trajectory for CO2 through a rigid framework. The relatively small difference in calculated hopping rate for CO2 in the ZIF-97 ⟨001⟩ channel contrasts with the large effect of framework flexibility found in ZIF-8.64 To investigate the applicability of the harmonic limit of TST for analyzing flexibility effects on CO2 hopping rates, we examined the magnitude of anharmonic effects on channel widths, employing ab initio molecular dynamics (AIMD) simulations for ZIF-93 at 300 K. From this simulation we assess the ⟨001⟩ channel size distributions, based on data from the coordinates of 1200 frames, sampled at each time step. In the primitive cell there are 3 ⟨001⟩ channels, each consisting of two sets of four linkers with functional groups closest to the center of the channel. Focusing on the distance between pairs of linkers across the channel, we obtain 4 distances per channel corresponding to “x” and “y” directions in both sets of 4 linkers. A total of 4 distances per channel, 3 channels per primitive cell,

coordinate (representing CO2 center of mass position) results with the vdW-DF2 are presented in Figure 4. The paths plotted in this figure are spline fits to the energies of the two end points at the 8-ring site and 8-ring mid site, the energies of the four intermediate images, and the forces at each point. With both vdW-DF2 and vdW-optB88, the barrier is lower for the RFC structure than for the RFE structure. Furthermore, as illustrated in the plots, the adsorption energy in the RFC structure of the 8-ring mid site becomes higher relative to the barrier and the energy at the 8-ring site. Flexibility has the effect of lowering the barrier by 33 meV in the vdW-DF2 formalism and 19 meV in the vdW-optB88 formalism. These numbers amount to 20 and 14% reductions in the calculated activation energies, respectively. The reason for this relatively small value is that while the energy of the highest energy image does, in fact, decrease for both methods with flexibility, so does the adsorption energy in the 8-ring site. The effect of framework flexibility on the hopping rates can be estimated using the numbers given above within the framework of harmonic transition state theory (TST).63 Within this theory the hopping rate is given by the expression + k TST =

kBT ΔS / kB −E b / kBT e e h

(6)

where ΔS denotes the entropy difference between the CO2 molecule at the saddle and energy-minimum positions, and Eb is the activation energy. In harmonic TST, ΔS is estimated from a vibrational-frequency analysis under the harmonic approximation. In the present work we have estimated ΔS by computing vibrational frequencies from force−displacement calculations for the CO2 molecule and the four nearest linkers at the binding and saddle point configurations identified by the cNEB calculations. Using the resulting vibrational frequencies and the calculated activation energies we obtain hopping frequencies of 1.53 and 8.25 GHz at 300 K in RFC and FFC calculations using vdW-DF2 functionals, respectively. These F

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

and 1200 snapshots yield 14400 distances produced by the AIMD run. We find that the average width of the ⟨001⟩ channel at 300 K in ZIF-93 increases from the vdW-DF2 relaxed, zerotemperature value by 0.11 ± 0.09 Å. This is 1.6% of the atom to atom channel width of 6.54 Å measured experimentally. It is also close to the value of 0.10 Å that the channel expands in the presence of CO2 in the zero-temperature flexible framework considered above, which in that case produced a change in the hopping rate that is small compared to the difference in hopping rates for the materials in this work. Thus, we find the ZIF-93 framework to be relatively stiff, i.e., with a relatively small expansion of the channel width due to anharmonicity. 3.3. Transport Barriers and Electrostatic Potentials in ZIF-93 and ZIF-97. As mentioned in Section 3.1, the CO2 jump frequency through the ⟨001⟩ channel is 19 times higher in ZIF-93 than ZIF-97. This is unexpected since the ⟨001⟩ channel in ZIF-97 is slightly wider than in ZIF-93 with our MD potentials. To better understand this large difference in CO2 transport, we perform a cNEB calculation on carbon dioxide in ZIF-93 and compare the minimum-energy pathway for CO2 transport produced with the one for ZIF-97 considered in the previous section. In all of the calculations that follow, we hold framework atoms fixed at their experimentally determined positions and utilize only the vdW-DF2, to be consistent with previous work on gas binding in ZIFs.27,29 In ZIF-93 there are two symmetrically equivalent binding sites in the ⟨001⟩ channel and a mid channel energy minimum, as in ZIF-97, however the mid site is much more shallow. We find the barrier between binding and mid sites in ZIF-93 to be 96 meV, compared to 243 meV for ZIF-97. Consistent with the above calculations, transport through the ⟨001⟩ channel in ZIF-97 is observed in the classical MD simulations to be a two step process, into the mid site, and out the other side of the channel. In contrast, for ZIF-93 the mid channel binding site is weak enough that much less time is spent there during transit between the binding site on either side of the ⟨001⟩ channel. Both the 147 meV higher barrier in ZIF-97, and the 2-step transport process (which is expected to lead to a lower transmission coefficient through the channel) are consistent with fewer CO2 jumps through this channel compared to ZIF-93, as found in the classical MD simulations and reported in Table 3. Interestingly, in both ZIF-93 and ZIF97, the binding sites at the entrances to the channel are characterized by a CO2 that lines up the molecule parallel to the axis of the channel. To further understand the above calculated results, we investigate the role of electrostatic interactions in governing the relatively slower hopping rates for CO2 in ZIF-97 versus 93. Figure 5 shows a slice of the framework electrostatic potential calculated by vdW-DF2 for the ⟨001⟩ channels in ZIF-93 and -97. The darkest red color indicates regions where the potential is −6 V or lower and the darkest blue color indicates potentials that are +6 V and higher, whereas the green color corresponds to regions where the electrostatic potential is at 0 V. Positions of the CO2 molecule during transport through the channel are also shown. The orientation of the images in this figure is similar to that of the rightmost image in Figure 2. For ZIF-93 and -97 the CO2 molecules on the left and right are at the ⟨001⟩ channel binding site and mid site energy minima, respectively, and the center CO2 is close to a saddle point. The classical MD simulations also predict CO2 orientations parallel and perpendicular to the channel in the binding and mid site, respectively.

Figure 5. Slices of the electrostatic potential in the (020) plane showing the [001] channel in ZIF-93 and -97. Intermediate CO2 positions taken from cNEB vdW-DF2 calculations are shown, where red spheres represent oxygen and brown spheres represent carbon. Irregularities in the CO2 sphere and cylinder models come from the molecular coordinates passing through the plane of the plotted electrostatic potential.

The electrostatic potentials plotted for ZIF-93 and -97 in Figure 5 are observed to be qualitatively different. In ZIF-93 the region between the binding sites consists of negative potentials that are roughly homogeneous. In ZIF-97, the region between the left and right binding sites in Figure 5 consists of both positive and negative potentials, such that there exist positive potential regions away from the center of the channel, shown as four teal colored spots in the figure. The positive potentials attract the negatively charged oxygen atoms in CO2 which favors the perpendicular orientation and creates the moderate binding energy of the mid site. In ZIF-93, the binding of this site is much weaker and the barriers between the sites are relatively smaller. The difference in electrostatic potentials for ZIF-93 versus ZIF-97 is most likely due to the oxygen atoms on the −CHO functional groups in the former. These oxygen atoms protrude into the channel, as shown in Figure 3, and are responsible for the homogeneous negative electrostatic potential. By contrast, G

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

Figure 6. Calculated CO2 permeability and CO2/CH4 selectivity of ZIF-25, 71, 93, 96, and 97 at 1 bar and 298 K obtained in the present work are plotted along with values taken from the literature for MOFs, ZIFs, and MOF/polymer hybrid materials. The open symbols represent calculated values from molecular simulations and filled symbols represent experimental measurements. The values from the literature are taken from the following refs: (a) 21, (b) 6, (c) 65, (d) 17, (e) 19, (f) 16, (g) 66, (h) 67, (i) 68, (j) 15, (l) 69, (l) 70, (m) 71, (n) 72, (o) 73, and (p) 74. The solid black line represents the Robeson 2008 upper bound for polymer membrane CO2/CH4 separation performance.75

ZIF-97 has −CH2OH functional groups, for which the hydrogen atoms are closest to the center of the channel, producing the regions of positive potential, and for which the oxygen atoms are closer to the entrances of the channel and create regions of negative potential. Overall, the analysis presented in this section highlights that for the ⟨001⟩ channel in RHO topology ZIFs, the electrostatic potential greatly influences the energetics associated with the transport of CO2 and can contribute to the diffusion selectivity of CO2 over CH4. 3.4. Membrane Selectivity. The permeability for CO2 (CH4) is estimated by the product of the CO2 (CH4) transport diffusivity, calculated in this work, and the experimental CO2 (CH4) concentration at 1 bar and 298 K, taken from refs 26 and 28. To estimate the membrane selectivity, previously reported adsorption selectivities for equimolar CO2/CH4 mixtures,28 calculated with ideal adsorbed solution theory (IAST) from individual gas adsorption isotherms, are multiplied by the ratio of CO2/CH4 transport diffusivities, which were calculated in the current work with single component gases. The results are listed in Table 2. The CO2/CH4 membrane separation performance of ZIF-93 is predicted to be much higher than the other ZIFs studied in this work, with a selectivity of 26.7 and a CO2 permeability of 8.32 × 104 barrers. In addition, the performance of this material compares favorably with that of many other MOFs, ZIFs, and polymer/ MOF composites studied experimentally or computationally, as illustrated in Figure 6. The origin of the separation performance of ZIF-93 differs from that of two other ZIF materials predicted to combine high

selectivity and permeability: ZIF-11 and -12.21 ZIF-11 and -12 differ from each other in the transition metal atom (Zn and Co), but share the same RHO topology with the ZIFs considered in this work. The linkers in ZIF-11 and -12 are benzimidazoles, so that a benzene ring takes the place of the functional groups described in Section 2.1. These benzene rings protrude into the ⟨001⟩ and ⟨111⟩ channels. Yilmaz and Keskin calculated that the narrow channel diameters (more narrow than the ZIFs considered here) slow the diffusion for CH4, but not CO2, such that diffusion coefficients for these gases differ by about 2 orders of magnitude.21 However, Thornton et al. calculated a much lower CO2/CH4 membrane selectivity and CO2 permeability for ZIF-11.6 The calculations in the work of Thornton et al. differ from those performed by Yilmaz and Keskin in the potentials used, inclusion of framework partial charges, and simulation pressure. (We note that for other systems where independent simulations have been published, namely ZIF-2, -3, -69, -71, and -90, the different results plotted for each in Figure 6 show better agreement.) In contrast to ZIF-11 and -12, the wider channels of ZIF-93 do not hinder CH4 diffusion as greatly, however, electrostatic interactions that line up the CO2 molecule for transport at the channel entrance, along with a weak mid channel binding site, help retain diffusion selectivity for CO2. In addition, there is a larger CO2/CH4 adsorption selectivity in ZIF-93 than in ZIF11 or -12. Because the molecular sieving function is maintained for a larger channel diameter in ZIF-93, the results should be less sensitive to the simulation details that may affect the channel constriction. Also, for realistic applications, a membrane with wider channels should function better at H

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

(3) Basu, S.; Khan, A. L.; Cano-Odena, A.; Liu, C.; Vankelecom, I. F. Chem. Soc. Rev. 2010, 39, 750. (4) Kidnay, A. J.; Parrish, W. R. Fundamentals of Natural Gas Processing; CRC Press: Boca Raton, FL, 2006; Vol. 200. (5) Merkel, T.; Freeman, B.; Spontak, R.; He, Z.; Pinnau, I.; Meakin, P.; Hill, A. Science 2002, 296, 519. (6) Thornton, A.; Dubbeldam, D.; Liu, M.; Ladewig, B.; Hill, A.; Hill, M. Energy Environ. Sci. 2012, 5, 7637. (7) Phan, A.; Doonan, C.; Uribe-Romo, F.; Knobler, C.; O'Keeffe, M.; Yaghi, O. Acc. Chem. Res. 2010, 43, 58. (8) Chmelik, C.; Voß, H.; Bux, H.; Caro, J. Chem. Ing. Tech. 2011, 83, 104. (9) Krishna, R.; Van Baten, J. Chem. Eng. J. 2007, 133, 121. (10) Wijmans, J.; Baker, R. J. Membr. Sci. 1995, 107, 1. (11) Krishna, R.; van Baten, J. M. Phys. Chem. Chem. Phys. 2011, 13, 10593. (12) Caro, J. Curr. Opin. Chem. Eng. 2011, 1, 77. (13) Li, J.; Kuppler, R.; Zhou, H. Chem. Soc. Rev. 2009, 38, 1477. (14) Bux, H.; Chmelik, C.; van Baten, J. M.; Krishna, R.; Caro, J. Adv. Mater. 2010, 22, 4741. (15) Saha, D.; Bao, Z.; Jia, F.; Deng, S. Environ. Sci. Technol. 2010, 44, 1820. (16) Keskin, S.; Sholl, D. S. Ind. Eng. Chem. Res. 2008, 48, 914. (17) Battisti, A.; Taioli, S.; Garberoglio, G. Microporous Mesoporous Mater. 2011, 143, 46. (18) Atci, E.; Keskin, S. J. Phys. Chem. C 2012, 116, 15525. (19) Keskin, S. J. Phys. Chem. C 2010, 115, 800. (20) Liu, J.; Keskin, S.; Sholl, D. S.; Johnson, J. K. J. Phys. Chem. C 2011, 115, 12560. (21) Yilmaz, G.; Keskin, S. Ind. Eng. Chem. Res. 2012, 51, 14218. (22) Liu, J.; Rankin, R. B.; Karl Johnson, J. Mol. Simul. 2009, 35, 60. (23) Yu, D.; Yazaydin, A. O.; Lane, J. R.; Pascal, D.; Snurr, R. Q. Chem. Sci. 2013, 4, 3544. (24) Park, K. S.; Ni, Z.; Côté, A. P.; Choi, J. Y.; Huang, R.; UribeRomo, F. J.; Chae, H. K.; O’Keeffe, M.; Yaghi, O. M. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 10186. (25) Banerjee, R.; Phan, A.; Wang, B.; Knobler, C.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Science 2008, 319, 939. (26) 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. (27) Ray, K. G.; Olmsted, D. L.; Houndonougbo, Y. A.; Laird, B. B.; Asta, M. D. J. Phys. Chem. C 2013, 117, 14642. (28) Houndonougbo, Y.; Signer, C.; He, N.; Morris, W.; Furukawa, H.; Ray, K. G.; Olmsted, D. L.; Asta, M.; Laird, B. B.; Yaghi, O. M. J. Phys. Chem. C 2013, 117, 10326. (29) Ray, K. G.; Olmsted, D.; He, N.; Houndonougbo, Y.; Laird, B. B.; Asta, M. Phys. Rev. B 2012, 85, 085410. (30) Banerjee, R.; Furukawa, H.; Britt, D.; Knobler, C.; O'Keeffe, M.; Yaghi, O. M. J. Am. Chem. Soc. 2009, 131, 3875. (31) Liu, Z.; Grande, C. A.; Li, P.; Yu, J.; Rodrigues, A. E. Sep. Purif. Technol. 2011, 81, 307. (32) Liu, B.; Smit, B. J. Phys. Chem. C 2010, 114, 8515. (33) Li, B.; Wei, S.; Chen, L. Mol. Simul. 2011, 37, 1131. (34) Rankin, R. B.; Liu, J.; Kulkarni, A. D.; Johnson, J. K. J. Phys. Chem. C 2009, 113, 16906. (35) Amrouche, H.; Aguado, S.; Pérez-Pellitero, J.; Chizallet, C.; Siperstein, F.; Farrusseng, D.; Bats, N.; Nieto-Draghi, C. J. Phys. Chem. C 2011, 115, 16425. (36) Morris, W.; He, N.; Ray, K. G.; Klonowski, P.; Furukawa, H.; Daniels, I. N.; Houndonougbo, Y. A.; Asta, M.; Yaghi, O. M.; Laird, B. B. J. Phys. Chem. C 2012, 116, 24084. (37) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. J. Chem. Phys. 2000, 113, 9901. (38) Dion, M.; Dydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92, 246401. (39) Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Phys. Rev. B 2010, 82, 081101.

higher pressures and over a larger temperature range, as well as be more robust against membrane fouling.

4. CONCLUSION Classical MD simulations have been employed to compute the diffusivities of CO2 and CH4 in a series of RHO topology ZIFs: ZIF-25, -71, -93, -96, and -97. ZIF-93 was found to be unique in its diffusion selectivity, with a CO2 transport diffusivity that is 3.25 times greater than the value for CH4. Our MD simulations predicted a larger jump frequency for CO2 through the ⟨001⟩ channel in ZIF-93 than ZIF-97. However, this trend is reversed for CH4, which contributes to the large difference in diffusion selectivities in these two materials, despite their similar channel diameters. It is determined with cNEB calculations performed with the vdW-DF2 method that differences in the electrostatic potential in the channel contribute to the difference in CO2 transport barriers. Also, the electrostatic potentials in ZIF-93 and 97 are consistent with different binding strengths of the CO2 in a midchannel site, with a strong binding unfavorable for transport present in ZIF-97. In both ZIF-93 and -97, the binding site at the entrance to the channel line up the CO2 for parallel entrance. Thus, for future engineering of membrane materials, the electrostatic potential through the channels should be considered an important design consideration. With the MD calculated CO2 /CH4 diffusivities and previously reported adsorption data,26,28 permeabilities and membrane selectivities for ZIF-25, -71, -93, -96, and -97 were computed. In ZIF-93, diffusion selectivity is predicted to be combined with good adsorption selectivity. Accordingly, this material is predicted to have membrane selectivity of 26.7 for an equimolar mixture at 1 bar and 298 K, while having a permeability for CO2 of 8.32 × 104 barrers. This selectivity is higher than any ZIF the authors have found experimental data for in the literature and the permeability also compares favorably with other high-performing ZIF materials predicted from computational studies.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address §

Cornell University, Ithaca, New York, United States

Author Contributions ‡

K.G.R. and D.L.O. contributed equally to this work

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported as part of the Molecularly Engineered Energy Materials (MEEM), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award DESC0001342. This work made use of resources of the National Energy Research Scientific Computing Center, supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC02-05CH11231.



REFERENCES

(1) Burnham, A.; Han, J.; Clark, C. E.; Wang, M.; Dunn, J. B.; PalouRivera, I. Environ. Sci. Technol. 2011, 46, 619. (2) DeSutter, T. M.; Ham, J. M. J. Environ. Qual. 2005, 34, 198. I

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article

(40) Klimeš, J.; Bowler, D. R.; Michaelides, A. J. Phys.: Condens. Matter 2010, 22, 022201. (41) Deng, H.; Grunder, S.; Cordova, K. E.; Valente, C.; Furukawa, H.; Hmadeh, M.; Gándara, F.; Whalley, A. C.; Liu, Z.; Asahina, S.; Kazumori, H.; O'Keeffe, M.; Terasaki, O.; Stoddart, J. F.; Yaghi, O. M. Science 2012, 336, 1018. (42) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comp Phys. 1995, 117, 1−19, http://lammps.sandia.gov. (43) Martin, M. G.; Chen, B.; Wick, C. D.; Potoff, J. J.; Stubbs, J. M.; Siepmann, J. I. MCCCS Towhee http://towhee.sourceforge.net. (44) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (45) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2004, 108, 17596. (46) Rappe, A. K.; Casewit, C. J.; Colwell, K.; Goddard Iii, W.; Skiff, W. J. Am. Chem. Soc. 1992, 114, 10024. (47) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (48) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. J. Comput. Phys. 1977, 23, 327. (49) Nosé, S. J. Chem. Phys. 1984, 81, 511. (50) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (51) Martyna, G. J.; Klein, M. J.; Tuckerman, M. J. Chem. Phys. 1993, 97, 2635. (52) Kresse, G.; Hafner, J. Phys. Rev. B 1993, 47, 558R. (53) Kresse, G.; Hafner, J. Phys. Rev. B 1994, 49, 14251. (54) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15. (55) Kresse, G.; Furthmuller, J. Phys. Rev. B 1996, 54, 11169. (56) Klimeš, J.; Bowler, D. R.; Michaelides, A. Phys. Rev. B 2011, 83, 195131. (57) Blöchl, P. Phys. Rev. B 1994, 50, 17953. (58) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 50, 17953. (59) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (60) Román-Pérez, G.; Soler, J. M. Phys. Rev. Lett. 2009, 103, 096102. (61) Canepa, P.; Nijem, N.; Chabal, Y. J.; Thonhauser, T. Phys. Rev. Lett. 2013, 110, 026102. (62) Watanabe, T.; Keskin, S.; Nair, S.; Sholl, D. Phys. Chem. Chem. Phys. 2009, 11, 11389. (63) Hänggi, P.; Talkner, P.; Borkovec, M. Rev. Mod. Phys. 1990, 62, 251. (64) Haldoupis, E.; Watanabe, T.; Nair, S.; Sholl, D. S. ChemPhysChem 2012, 13, 3449. (65) Atci, E.; Erucar, I.; Keskin, S. J. Phys. Chem. C 2011, 115, 6833. (66) Song, Q.; Nataraj, S.; Roussenova, M. V.; Tan, J. C.; Hughes, D. J.; Li, W.; Bourgoin, P.; Alam, M. A.; Cheetham, A. K.; Al-Muhtaseb, S. A.; Sivaniah, E. Energy Environ. Sci. 2012, 5, 8359. (67) Zou, X.; Zhang, F.; Thomas, S.; Zhu, G.; Valtchev, V.; Mintova, S. Chem.Eur. J. 2011, 17, 12076. (68) Huang, A.; Dou, W.; Caro, J. J. Am. Chem. Soc. 2010, 132, 15562. (69) Li, Y.-S.; Liang, F.-Y.; Bux, H.; Feldhoff, A.; Yang, W.-S.; Caro, J. Angew. Chem. 2010, 122, 558. (70) McCarthy, M. C.; Varela-Guerrero, V.; Barnett, G. V.; Jeong, H.K. Langmuir 2010, 26, 14636. (71) Liu, Y.; Hu, E.; Khan, E. A.; Lai, Z. J. Membr. Sci. 2010, 353, 36. (72) Bae, T.-H.; Lee, J. S.; Qiu, W.; Koros, W. J.; Jones, C. W.; Nair, S. Angew. Chem., Int. Ed. 2010, 49, 9863. (73) Perez, E. V.; Balkus, K. J.; Ferraris, J. P.; Musselman, I. H. J. Membr. Sci. 2009, 328, 165. (74) Zhang, Y.; Musselman, I. H.; Ferraris, J. P.; Balkus, K. J., Jr. J. Membr. Sci. 2008, 313, 170. (75) Robeson, L. M. J. Membr. Sci. 2008, 320, 390.

J

dx.doi.org/10.1021/cm5015477 | Chem. Mater. XXXX, XXX, XXX−XXX