Article pubs.acs.org/JPCC
Accurate Computation of Gas Uptake in Microporous Organic Molecular Crystals Wenliang Li,†,‡ Stefan Grimme,*,‡,§ Helge Krieg,‡ Jonas Möllmann,‡,§ and Jingping Zhang*,† †
Faculty of Chemistry, Northeast Normal University, Changchun, China Theoretische Organische Chemie, Organisch-Chemisches Institut der Universität Münster, Germany § Mulliken Center for Theoretical Chemistry, Institut für Physikalische und Theoretische Chemie der Universität Bonn, Germany ‡
S Supporting Information *
ABSTRACT: Microporous organic molecular crystals (MOMCs) are materials composed of discrete organic molecules interacting noncovalently. The gas uptake of CO2, CH4, N2, and Xe in one of the most widely investigated MOMCs, trispiro(benzodioxole[2′.2:2″.4:2‴.6]1,3,5,2,4,6-triazatriphosphinine) (1), was computed by grand canonical Monte Carlo (GCMC) simulations based on our lately developed force field method vdW3, which is fitted to accurate ab initio potentials (double hybrid functional B2PLYP with a D3 dispersion correction using def2-TZVPP basis sets). The B2PLYP-D3 results compare very well with CCSD(T)/CBS interaction energies for benzene−gas model complexes. Our multiscale modeling approach to accurate interaction potentials is found to be essential in order to obtain reasonable adsorption isotherms and heats of adsorption. The good agreement between the simulation results and experimental data in all cases gives us the opportunity to design novel complex materials for gas adsorption based solely on first-principles calculations.
1. INTRODUCTION Microporous materials are solids containing interconnected pores of less than 2 nm in size. They are widely used for heterogeneous catalysis, adsorption, separation, gas storage, and a number of emerging technologies.1−3 Inspired by the intriguing properties of metal organic frameworks (MOFs), the preparation of purely organic porous materials (including covalent organic frameworks (COFs),4,5 polymers of intrinsic microporosity (PIM),6,7 porous aromatic frameworks (PAF),8,9 hypercrosslinked polymers (HCPs),10 crystalline triazine-based organic frameworks (CTFs),11 and conjugated microporous polymers (CMPs)12,13) is currently of intense interest. In addition to extended networks and frameworks, an alternative way to form porosity is to use discrete molecules between which there are only noncovalent interactions. Because microporous organic molecular crystals (MOMCs, sometimes termed “organic zeolites”14,15) do not possess an extended framework with covalent or coordination bonds, they can be dissolved and then reassembled in suitable solvents. However, because molecules pack efficiently in the solid state and therefore are nonporous, MOMCs are comparatively rare. MOMCs that are stable to desolvation include trispiro(benzodioxole[2′.2:2″.4:2‴.6] 1,3,5,2,4,6-triazatriphosphinine) (1; with the short name of TPP in previous works),16−21 calixarenes,22 cucurbit[n]uril (n = 6, 7),23,24 some dipeptides,25,26 3,3′,4,4′-tetra(trimethylsilylethynyl)biphenyl,27 and, recently, imine-based organic cages.28−32 Computational modeling is a powerful tool to complement experiment and design new porous materials.33−37 Grand © 2012 American Chemical Society
canonical Monte Carlo (GCMC) simulations are able to predict gas uptake and the heats of adsorption (Qst), which can be compared directly with experimental data, and therefore offer an alternative way to laborious experimental work. Empirical force fields, such as the universal force field (UFF),38 DREIDING,39 and the optimized potential for liquid simulations-all atom (OPLSAA) models,40 have often been used to describe the adsorbate−adsorbate interactions as well as adsorbate−adsorbent interactions in many previous GCMC simulations of gas adsorption in porous materials. However, because the empirical force fields use simple mathematical models (e.g., of Lennard-Jones (LJ) type) to describe the potentials, GCMC simulations performed with these force fields can potentially produce inaccurate gas uptake in porous materials.41 The lack of a sophisticated force field and accurate force field parameters for specific systems significantly limits the investigation of new porous materials.42 Recently, a significant number of theoretical works which combine first-principles calculations with GCMC simulations have been explored in order to compute the uptake of gases in porous materials. Düren and her co-workers directly used a DFT/CC-based potential energy surface implemented in GCMC simulations to predict the adsorption of methane in CuBTC.43 Froudakis and co-workers investigated hydrogen adsorbed in COF44−46 by fitting the LJ potential to an ab initio Received: November 22, 2011 Revised: April 5, 2012 Published: April 5, 2012 8865
dx.doi.org/10.1021/jp2112632 | J. Phys. Chem. C 2012, 116, 8865−8871
The Journal of Physical Chemistry C
Article
reduce the computation time of the MP260 part in the functional. All CCSD(T)/CBS calculations were performed with the Gaussian 03 program.61 The MP2 correlation energies were extrapolated to the CBS limit by applying the two-point extrapolation scheme of Helgaker et al.57
calculated potential energy surface (calculated with DFT and/ or MP2). Similarly, Sun and co-workers evaluated hydrogen uptake in MOF47 and PAF.48 Besides the LJ potential, the Morse potential has been used extensively to describe the relevant noncovalent interactions. A multiscale simulation study of hydrogen adsorption in MOFs was performed by Han et al. in 2007 using the Morse potential derived from RI-MP2/ TZVPP and RI-MP2/QZVPP.49 In a similar way, the adsorption of hydrogen50 and methane51 in COF was predicted. More work about multiscale simulation of adsorptive processes can be found in recent reviews.42,52 Many multiscale simulation works gave pretty good results in comparison to experimental data, but there still remain some challenges for the development of this method. A sufficient amount of accurate reference data is needed for fitting of the force fields. Additionally, developing a more sophisticated force field which can produce accurate interaction energies for a wide range of systems and intermolecular distances still remains a main challenge. In this work, we compute the gas uptake (including carbon dioxide, methane, nitrogen, and xenon) and Qst in crystal of 1, which represents a realistic material to investigate the main features of MOMCs and their potential applications. The molecular structure and packing topology of 1 are illustrated in Figure 1. Several sets of potential energy curves are calculated
CBS EMP2 =
X3·EXMP2 − Y 3·EYMP2 X3 − Y 3
where X and Y are the cardinal numbers of the two basis sets. In the present work, we used aug-cc-pVDZ and aug-cc-pVTZ62 for the extrapolation of MP2. The CCSD(T)/aug-cc-pVDZ correlation energies were then estimated for the complete basis set limit (i.e., CCSD(T)/CBS) by applying the following additivity scheme:63 CBS CBS ECCSD(T) = EMP2 + (ECCSD(T) − EMP2)aug ‐ cc ‐ pVDZ
Corrections for basis set superposition error (BSSE) are not applied for the B2PLYP-D3/def2-TZVPP calculation considering the following aspects. Though it is typically 10−20% (which is not larger than basis set effects for many other chemical properties for which typically also no corrections are applied) of the interaction energy at the B2PLYP/def2-TZVPP level, the good agreement between estimated CCSD(T)/CBS results (which is more or less BSSE free) and the B2PLYP-D3/def2TZVPP data in selected cases indicates that BSSE effects are even smaller. Note that the D3 dispersion correction has been developed without any BSSE correction so that remaining effects are partially absorbed by the D3 short-range parametrization. Because counter-poise corrections require a huge amount of computational resources and are also not free of error, they are not considered for mass production of reference data for the force field parametrization. This procedure is reasonable because it can be judged from a comparison of theoretical and experimental thermochemical observables.55 2.2. Force Field. The vdW3 force field developed by one of the authors is a purely intermolecular potential that is designed to accurately describe the noncovalent interactions between distant molecules. We model the total interaction energy between a pair of molecules as a sum of four contributions that are computed from atom-pairwise and charge-pairwise interactions. The terms included are Pauli-exchange repulsion, Coulomb interaction, induction, and dispersion. Repulsion is modeled in the spirit of the Buckingham potential,64 whereas the latter three follow the known physical interatomic distance dependencies of R−1, R−3, and R−6. We also included higher order terms for the dispersion energy. The key ingredients of vdW3 are the use of off-site charges computed from semiempirical monomer wave functions for electrostatic and induction terms and the application of ab initio computed accurate system-dependent dispersion coefficients.55 Details can be found in the Supporting Information. The proposed force field model is applicable to all 94 elements up to Pu though containing only a small number of 10 empirical parameters in total. In order to increase the accuracy, we optimized the parameters specifically for the systems of interest in this work. 2.3. Fitting Procedure. In total, 30 potential curves for different conformations of model complexes of one molecule of 1 in combination with one gas molecule (CO2, CH4, N2, or Xe) as well as model systems of gas−gas dimers have been used for the fitting of the force field. These are 312 different points of reference in total. The cuts of the potential energy surface are
Figure 1. Molecular structure and packing topology of 1: (a and b) molecule of 1; (c) 2*2*2 supercell packing of the pseudohexagonal framework.
for each adsorption system using density functional theory (B2PLYP53/def2-TZVPP54 with the D355 dispersion correction). These are then taken as reference data for fitting of our recently developed force field vdW3, which is specifically designed for the description of intermolecular noncovalent interactions. Finally, the GCMC simulations are performed based on the energy evaluation of the fitted vdW3.
2. METHODOLOGY 2.1. Reference Data. The molecular structure of 1 is extracted directly from crystal data.16 Gas−gas and gas−1 interaction energies are calculated at the B2PLYP-D353,55/def2TZVPP54 level, which is verified by the CCSD(T)/CBS method.56,57 The double hybrid functional calculations were carried out with TURBOMOLE6.0.58 An SCF convergence criterion of 10−7 EH and a grid size of “m4” for the numerical integration of the exchange-correlation energy was used. The resolution-of-the-identity (RI) approximation59 was applied to 8866
dx.doi.org/10.1021/jp2112632 | J. Phys. Chem. C 2012, 116, 8865−8871
The Journal of Physical Chemistry C
Article
Figure 2. Comparison of calculated interaction energies of benzene−gas complexes at the B2PLYP-D3/def2-TZVPP level (black triangles) with the ones estimated at the CCSD(T)/CBS level (red circles). The lines are drawn to guide the eye.
selected on the basis of how gas molecules locate in the channel of the porous material. The channel is mainly composed of benzene groups of TPP (see Figure 1C), and for this reason, the potential curves of gases surrounding benzene rings are taken into consideration. Figures 3 and S2−S5 (in the Supporting Information) present all of the generated potential curves including information about the coordinate that is varied. The vdW3 force field was optimized by minimizing the root-mean-square deviation (rms) between the interaction energies computed with the reference method B2PLYP-D3/ def2-TZVPP and with the applied force field. We minimized the rms using numerically computed gradients following the BFGS algorithm.65−68 2.4. GCMC Procedure. The gas loadings and heats of adsorption in a fixed 2*2*2 supercell of 116 were computed from the GCMC simulations. All simulations were performed with the Monte Carlo simulation suite of the MUSIC code69 where we replaced the standard LJ potential with the fitted vdW3 potential. For CO2, CH4, and N2, we used rotation, translation, insertion, and deletion as possible Monte Carlo steps, whereas only the latter three were possible steps for Xe. Each simulation started with a one-million step equilibration period followed by a one-million step production run. Noncovalent interactions of two molecules with the distance of their centers of mass less than 9.9 Å were taken into consideration. GCMC simulations give the absolute number of molecules loaded, while experiments typically produce the excess amount adsorbed. Thus, all absolute loadings were converted to excess loadings to allow direct comparison of the results with experimental data.70
The isosteric heat of adsorption (Qst) was calculated for each gas as the difference of the partial molar enthalpy of the sorbate in the bulk phase and the partial molar internal energy in the adsorbed phase.71 This is given by the following equation Q st = RT −
⎛ ∂U ⎞ ⎜ ⎟ ⎝ ∂N ⎠T , V
#tab;where T is the temperature, R is the universal gas constant, N is the adsorption loading, and U is the internal energy of the sorbate in the adsorbed phase that includes contributions from both adsorbate−adsorbent and adsorbate−adsorbate interactions. To compare quantitatively the obtained theoretical trends with experimental results, the Langmuir−Freundlich (LF) adsorption model is adopted.72 N=
Mkf β 1 + kf β
where f is the fugacity of the bulk gas at equilibrium with the adsorbed phase and M, k, and β are model parameters of the maximum adsorption amount, the affinity constant, and the deviation from the simple Langmuir equation, respectively.
3. RESULTS AND DISCUSSION 3.1. Force Field Fitting. Producing the reference data is a fundamental step in multiscale simulations. For the porous coordination frameworks, the system has to be divided into molecular units to perform quantum mechanical calculations,42 which means some covalent bonds have to be cut and dangling bonds must be saturated. However, for the current molecular crystal, taking the molecule of 1 as one monomer in the gas− 8867
dx.doi.org/10.1021/jp2112632 | J. Phys. Chem. C 2012, 116, 8865−8871
The Journal of Physical Chemistry C
Article
molecule complex is straightforward. The B2PLYP-D3/def2TZVPP method was employed in this work because it gives very accurate interaction energies for a vast range of systems73 and, especially, for the gas−benzene system. In order to evaluate the accuracy of the density functional for the systems of interest in this work, we compare B2PLYP-D3/def2-TZVPP results with highly accurate CCSD(T)/CBS reference energies. We restrict our comparison to relatively small gas−benzene complexes because the computational effort for generating CCSD(T)/CBS reference energies for gas−1 systems is currently too demanding and the gas−benzene system is a suitable model for the major types of interactions of gases in the channels of 1. Although lots of potential curves are needed to describe the whole potential surface of the gas−benzene complexes, only three special symmetric orientations are used to represent the potential surface for each gas because the computational effort can be reduced by symmetry. This is known to have very little influence on the comparison between the two methods.74,75 We generated potential energy curves for the nine model systems by changing the distance between the fragments’ centers of mass while keeping the geometry of the fragments frozen. As shown in Figure 2, the interaction energies obtained by the dispersion corrected double hybrid density functional with a triple-ζ basis set are in good agreement with the interaction energies computed with the reference method for all of the gas−benzene model systems. We found a very low rms value of
