Methane Adsorption in Zr-Based MOFs: Comparison and Critical Evaluation of Force Fields.

The search for nanoporous materials that are highly performing for gas storage and separation is one of the contemporary challenges in material design. The computational tools to aid these experimental efforts are widely available, and adsorption isotherms are routinely computed for huge sets of (hypothetical) frameworks. Clearly the computational results depend on the interactions between the adsorbed species and the adsorbent, which are commonly described using force fields. In this paper, an extensive comparison and in-depth investigation of several force fields from literature is reported for the case of methane adsorption in the Zr-based Metal-Organic Frameworks UiO-66, UiO-67, DUT-52, NU-1000, and MOF-808. Significant quantitative differences in the computed uptake are observed when comparing different force fields, but most qualitative features are common which suggests some predictive power of the simulations when it comes to these properties. More insight into the host-guest interactions is obtained by benchmarking the force fields with an extensive number of ab initio computed single molecule interaction energies. This analysis at the molecular level reveals that especially ab initio derived force fields perform well in reproducing the ab initio interaction energies. Finally, the high sensitivity of uptake predictions on the underlying potential energy surface is explored.

Introduction

For the past two decades, metal–organic frameworks (MOFs) have received considerable attention in scientific literature. This is in part due to the seemingly unlimited number of possible combinations of metal nodes coordinated by organic ligands, with each combination resulting in the unique charateristics of that MOF. The mixing and matching of these building blocks lead to a wide range of chemical environments, topologies, and pore sizes, with practical applications including heterogeneous catalysis[1,2] and gas storage or separation.[3] From the onset, MOFs were mostly used and studied for their gas adsorption properties.[4−6] Also from a computational side, MOFs attracted extensive attention, as evidenced by many high-throughput studies that aim to identify the most promising candidate materials for a given application.[7,8]

Many studies on the computational prediction of adsorption isotherms are available; however, the correct prediction of isotherms is particularly challenging. To compute the adsorption isotherm of a guest molecule in a rigid framework, one can use Grand-Canonical Monte Carlo (GCMC) simulations which are able to unveil the influence of the underlying potential energy surface (PES) on the uptake of the guest molecules. Typically, millions of energy evaluations are required in GCMC simulations, making them intractable with ab initio methods. Given the current limitations on computing power, only force-field simulations are feasible. There is an abundance of force fields in literature that can be used for adsorption simulations (note that we only consider the noncovalent part of the force fields in this study, as covalent terms are not relevant in rigid-framework simulations). An important distinction is the difference between so-called generic force fields on the one hand and ab initio derived force fields on the other hand. In brief, generic force fields are typically fitted to reproduce certain experimental results such as crystal structures, sublimation energies, or fluid properties and are then applied to other systems and/or properties of interest. Parameters of generic force fields are available for many chemical environments, which makes them (in principle) widely applicable.[9] In contrast, ab initio derived force fields are typically more tailored toward one specific application. Previous classification gives only a brief distinction between the two major classes of force fields (generic versus ab initio derived). In reality many force fields were developed which have much more subtle nuances. A more complete description of force fields used in this work is given in section . For an in-depth review, we refer to the work by Coudert et al.[10]

Hereafter, we particularly highlight some computational studies on gas adsorption which are relevant for this paper, in which methane adsorption in Zr-based MOFs is taken as a case study. Vasanth Kumar et al. studied H2 and CH4 adsorption in UiO-66, UiO-67, and UiO-68, showing that linker–guest interactions are the main driving force for adsorption.[11] Their study was performed using GCMC simulations with a combination of different generic force fields (Dreiding, UFF, and TraPPE). In contrast to the current work, the obtained results were not validated using experimental adsorption data or ab initio computed adsorption energies. Snurr et al.[12] simulated the effect of missing linker defects on water and CO2 adsorption in UiO-66, using also a combination of different generic force fields (Dreiding, UFF, TraPPE, and TIP4P) and with additional account of framework atomic charges which were derived from ab initio calculations using the REPEAT method. It was shown that defect sites render the material more hydrophilic and that the location of the defects has an appreciable impact on uptake. Düren et al.[13] presented an adsorption study of methane in CuBTC in which direct information from the true potential energy surface (calculated with the DFT/CC[14] ab initio method) was used. For this MOF with coordinatively unsaturated sites, the presented approach provided quantitative agreement with experiment over a wide range of temperatures and pressures. However, it was observed that using other ab initio methods as a reference may provide deviating results, which gives an indication of the sensitiviy of the adsorption isotherm on the underlying potential energy surface. Finally, Schmidt et al.[15] compared a generic force field with an ab initio derived force field in a screening study involving 424 MOFs. Significant differences in CO2 and CH4 gas adsorption isotherms were observed; however, both types of force fields predict a similar ranking (with respect to uptake at a certain pressure and temperature) of the frameworks. The authors also found that the generic force field may benefit from a compensation of errors. Consequently, the physical interpretation of results obtained with generic force fields should be done with caution. Obviously, computational studies of gas adsorption in MOFs have covered a much wider array of frameworks and guest molecules than the few examples mentioned before.

In this work, we study methane adsorption in Zr-based MOFs, which proves an ideal case study for a comparison and evaluation of force fields to predict adsorption isotherms. At first sight it might be regarded as a rather easy problem from a computational point of view for two reasons. First, for the study at hand, the methane guest molecule does not feature a dipole or quadrupole moment, rendering electrostatic interactions with the framework unimportant (as shown in section S3), and polarization effects can be assumed negligible. Second, the frameworks that will be studied here do not feature coordinatively unsaturated sites (note that such sites can be present in defective UiO-66 structures, but those cases are not considered in this work). Such open metal sites can have a dramatic impact on adsorption properties and often require additional force-field terms not included in generic force fields.[16,17] Generic force fields perform better in defect-free materials with fully saturated metal centers. This is convincingly demonstrated in ref (18), where adsorption isotherms have been computed with generic force fields for CO2 and CH4 in dehydroxylated UiO-66. An excellent agreement with experiment is even obtained if Dreiding, UFF, and TraPPE parameters are selectively used for the CH4-MOF and CO2-MOF interactions. In another recent study, it was concluded that generic force fields are suitable for the qualitative prediction of methane adsorption in MOFs in the low-loading regime.[19] In the same paper, it is demonstrated that these force fields are capable of providing detailed molecular-level information, although in some systems such conclusions should be approached with caution.

In this paper, we study methane adsorption in Zr-based MOFs using five different force fields. Three of them are so-called generic force fields while we also included two force fields based on an ab initio description of the potential energy surface. Adsorption isotherms are calculated using GCMC simulations. Insight at the atomic level is gained from single molecule adsorption energies, obtained from density functional theory (DFT) calculations on the periodic lattice. The main goal of this paper is to critically assess which factors contribute to the overall reproduction of the adsorption isotherms, both quantitatively and qualitatively. The paper is organized as follows. section provides an elaborate description of the five used force fields. In section , one reads computational details. The main results of the GCMC simulations and the comparison with ab initio adsorption energies are presented in sections and 3.2, respectively. The conclusions are formulated in section .

Methods

Description of the Force Fields

Many different force fields are available in literature, which can be distinguished by their functional form, parametrization method, and design philosophy. This section discusses these features of the five force fields that will be used in this work.

The first two force fields are so-called generic force fields, based on generic, nonsystem-specific parameters. The first variant will be further designated UFF/TraPPE. The potential energy is a pairwise additive sum of Lennard-Jones potentials:where the sum runs over all pairs of atoms where atom i is part of the framework and atom j is part of the guest molecule. The parameters σ and ϵ are determined using Lorentz–Berthelot mixing rules respectively:The atomic parameters σ and ϵ of the MOF frameworks are taken from the Universal Force Field (UFF)[20] and are in this case uniquely determined by the atomic number of atom i. The MOFs we will study (UiO-66, UiO-67, DUT-52, MOF-808, and NU-1000) consist of Zr, O, C, and H atoms, leading to 8 parameters. The UFF parameters have been developed by making use of a combination of ab initio calculations, literature values, and empirical rules. Methane is described using the transferable potential for phase equilibria (TraPPE),[21] where a united-atom model with only one site is used for methane. Because methane is a neutral molecule, this site has zero charge which means that no electrostatic interactions are taken into account. The σ and ϵ of methane are parametrized to reproduce the experimental critical properties and coexistence densities. The second force field, designated as Dreiding-UFF/TraPPE, is very similar and shares the same functional form and number of parameters. The only difference with the UFF/TraPPE force field is that the framework parameters σ and ϵ are obtained from the Dreiding force field.[22] Again, parameters are only based on the atomic number of the atom and take modified values that were based on experimental crystal structures and sublimation enthalpies. Unfortunately, no Dreiding parameters are available for the Zr atom and therefore the ones available from UFF are taken. Both the Dreiding-UFF/TraPPE and the UFF/TraPPE force fields have been used before to model methane adsorption in UiO-66.[18]

The third force field, MM3-MBIS, uses generic parameters for the van der Waals contribution combined with system-specific charges to model electrostatic contributions:The van der Waals contribution is modeled using the Buckingham potential. The parameters σ and ϵ are again determined from atomic values using following mixing rules:Note that σ is defined as the sum of the atomic radii, rather than as the average of atomic radii (which is the case for the UFF and Dreiding force fields, see eq ). The definition of σ as the sum of the atomic radii is in line with the convention used by the authors of the MM3 force field.[23] The atomic parameters in turn are parametrized to reproduce experimental results such as molecular geometry, crystal unit cell parameters, and heat of sublimation.[23−25] For some elements of the periodic table, parameters are obtained by interpolation. Note that for MM3, the values of σ and ϵ are not uniquely determined by the atomic number of atom i, for instance a sp3 hybridized C atom has different parameters than a sp2 hybridized C atom. In the original MM3 work, electrostatic interactions are described using bond dipoles, but for simplicity, we will employ point charges instead. These atomic partial charges are obtained using the minimal basis iterative stockholder (MBIS) partitioning scheme.[26] It should be mentioned that for the nonpolar methane molecule electrostatic interactions have a modest impact, as shown in section S3.

The fourth force field is a purely ab initio based force field introduced by Schmidt et al. and is called SAPTFF.[27,28] The potential energy of this force field is decomposed into terms that correspond to the physically distinct intermolecular interaction energies from SAPT:where Eexch, Eelec, Epol, Edisp, and Eδhf correspond to the exchange-repulsion energy, electrostatic energy, polarization energy, dispersion energy, and higher-order contributions to polarization/exchange, respectively. For a methane molecule, the polarizability is taken to be zero and point-charge electrostatic interactions are not included.[15] Taking these simplifications into account, the above expression can be written asThe dispersion interactions are damped using the Tang-Toennies damping function and combination rules are as follows:where Atype is always positive for exchange, always negative for electrostatics and induction, and positive for δhf if and only if both Aδhf and Aδhf are positive. This functional expression is more complicated than the Lennard-Jones or Buckingham potentials used in the previously discussed generic force fields, but there is also an important difference with respect to the way parameters are determined. The SAPTFF parameters are determined based on monomer properties and dimer interactions, all computed using ab initio calculations. To generate the necessary reference data, cluster models of the frameworks under study are used. Different atom types are introduced of which 11 are necessary to describe the systems investigated in this paper. For each atom type, 9 parameters are introduced which leads to a total of 99 parameters as reported in section S6. SAPTFF has been used before in a screening study, which included methane adsorption in Zr-based MOFs.[15]

The fifth and final force field considered in this work is the monomer electron density based force field (MEDFF), which was proposed by the present authors.[29] It is also based on ab initio data and shares some common features with SAPTFF, although much less parameters need to be fitted to dimer interaction energies. The functional form is as follows:The first term, Eelst, represents electrostatic interactions including the penetration effect. The second term, (Uexch–rep – Uind)S, describes exchange-repulsion and short-range induction interactions based on the overlap of the monomer electron densities. The final terms, and , are dispersion terms that are damped at short-range using the Tang-Toennies damping function. Contrary to most other force fields, no atom types are defined in MEDFF. Instead the atomic parameters are obtained directly from an ab initio computed electron density through the MBIS partitioning scheme.[26] For instance, S is the overlap of the model electron densities, which can be accurately computed from the MBIS parameters of the atoms involved, as MBIS provides an analytically simple expression for electron density. Next to these atomic parameters, MEDFF also features three so-called interaction parameters: Uexch–rep, Uind, and Us8. The values for these parameters can not be obtained from the MBIS partitioning but are fitted to reproduce ab initio interaction energies (both SAPT and CCSD(T) reference data are used in the fitting procedure). In this work, we use the interaction parameters as determined in the original MEDFF paper[29] for a data set of dispersion-dominated complexes. This data set contains dimers of small molecules, such as benzene, ethene, and neopentane. As dispersion interactions are the driving force for methane adsorption, we expect these parameter values to result in a proper description of the systems at hand.

Computational Details

Monte Carlo Simulations

Monte Carlo simulations are an important tool for the computational study of adsorption. In this work, adsorption isotherms are calculated using GCMC with the RASPA software package.[30,31] Each run consists of 50000 equilibration cycles and at least 100000 production cycles, where a cycle consists of max(20, N) move attempts (with N the current number of adsorbed molecules). The error bars, computed using block averaging, indicate that increasing the number of cycles (to obtain better sampling of the visited parts of phase space) is unlikely to change the results by more than 1%. The applied chemical potential is computed from the fugacity, which is converted from the input pressure using the Peng–Robinson equation of state. Experimental values for the critical temperature, critical pressure, and acentric factor (which are input parameters of the Peng–Robinson equation) of methane are used.[32] Translation, rotation, reinsertion, deletion, and insertion moves are all attempted with a probability of 20%. For the all-atom force fields (MM3-MBIS, SAPTFF, and MEDFF) methane is considered as a rigid molecule with a geometry optimized at the B3LYP[33,34]/aug-cc-pVTZ[35] level-of-theory. For the frameworks UiO-66, UiO-67, and DUT-52, a 2 × 2 × 2 supercell of the conventional unit cell is employed. For NU-1000, a 1 × 1 × 2 supercell is used, while for MOF-808, the conventional unit cell is employed. During all RASPA simulations, the framework geometries are taken from experiment, except for the defective UiO-66 frameworks where the geometries are DFT optimized. In any case, the framework geometry is considered to be rigid during the simulations. The impact of this approximation is limited for the MOFs studied in this work, as we show in section S1. The vdW energy of each atom of the guest molecule and the electrostatic potential in the framework (and derivatives of these quantities) are tabulated on a cubic grid with a spacing of about 0.15 Å. During simulations, a tricubic interpolation scheme is used to obtain the host–guest energy using these grids.[36] Henry coefficients are also calculated with RASPA using 50000 Widom insertions.

The surface area (calculated with RASPA) is obtained geometrically based on the overlap of a helium probe with the framework atoms. Atoms are considered to overlap when they are closer than the minimum of the corresponding Lennard-Jones interaction, where the necessary radii are obtained from the Dreiding and UFF force fields. The available pore volume is calculated by RASPA as the amount of argon that is adsorbed at 87 K and a pressure of 0.5p0 using UFF parameters for the framework atoms and the parameters from Maitland et al. for argon.[37] The uptake (in moles) is converted to a volume by multiplying by the molar volume of bulk argon under the same conditions.

Force-Field Settings

All Monte Carlo simulations in this work are performed using force fields to evaluate the required energies, and therefore we now briefly discuss some settings related to the force-field energy evaluations. An important consideration for force fields in periodic systems is the truncation of the nonbonded interactions. Because of the slow decay of point-charge electrostatics (present only in MM3-MBIS and MEDFF), we employ the Ewald summation,[38,39] which allows one to exactly compute the electrostatic interaction of a periodic system under tinfoil boundary conditions. Similar techniques exist to exactly compute the contribution of dispersion terms decaying as r– (where typically n = 6,8,10, ...) but this is not implemented in RASPA. We therefore truncate the van der Waals interactions at a radius of 14 Å but employ tail-corrections to approximate the contributions beyond this cutoff. The influence of this setting is studied in more detail in the section S2. Note that this is the recommended treatment of long-range van der Waals interactions in the TraPPE model.

The MM3-MBIS and MEDFF force fields require MBIS parameters of both the frameworks and the methane guest molecule. For the MOF frameworks, the MBIS parameters are computed from the partitioning of the PBE[40] electron density computed with GPAW.[41,42] These calculations are performed using a real-space uniform grid with spacing of about 0.20 Å, with Γ-point-only sampling of the Brillouin zone, employing the supplied default PAW data set and standard convergence settings for the electronic self-consistent loop. The MBIS parameters for methane are obtained from the B3LYP[33,34]/aug-cc-pVTZ[35] electron density computed using Gaussian 09.[43]

Ab Initio Adsorption Energies

To gain more fundamental insight into the adsorption mechanisms that govern the GCMC simulations, we also compared force-field and ab initio computed single molecule adsorption energies. Adsorption energies of methane in UiO-66 were calculated with the PBE[40] functional with modified D3 dispersion corrections with Becke-Johnson damping (D3MBJ),[44] including the three-body contribution.[45] The three-body contribution consists of an Axilrod–Teller–Muto (ATM) term, and therefore, this level-of-theory will be referred to as PBE-D3MBJ-ATM. Periodic DFT calculations are performed using VASP.[46−48] The conventional unit cell of UiO-66 (a cubic cell with sides 20.75 Å containing 4 inorganic bricks) is used, which allows accurate sampling of the Brillouin-zone with the Γ point only. GW PAW potentials are used, the plane-wave basis set kinetic cutoff energy is set to 800 eV, and the electronic self-consistent loop is considered converged as soon as the change in energy drops below 10–6 eV. D3MBJ corrections together with the three-body term are calculated using the dftd3 program and are added to the PBE energies.

The Henry coefficient K of methane in UiO-66 is calculated using the same level-of-theory as follows. First, a cubic grid with spacing 0.262 Å (80 grid points in each direction, 512000 grid points in total) is constructed. For each grid point, the interaction energy is computed using Dreiding-UFF/TraPPE and if it is more repulsive than 1000 kJ mol–1 then this point is not considered for the ab initio calculations as it will not contribute to the Henry coefficient. This allows one to discard 327780 grid points. Of the remaining 184220 grid points, only 2401 points are unique because of the symmetry of the framework. For each of these grid points, several random rotations are considered and the number of rotations for grid point i is indicated as N. At least 5 random points are included in the final set, but additional samples are generated “on the fly” to ensure proper convergence. The Henry coefficient is finally simply calculated by numerical integration over the cubic grid:

Results and Discussion

Adsorption Isotherms of Methane in Zr-Based MOFs

Comparison with Experiment for Three Isoreticular MOFs

First of all, we study the adsorption of methane in three isoreticular, defect-free zirconium-based MOFs for which experimental data is available: UiO-66,[49] UiO-67,[49] and DUT-52,[50] which are illustrated in Figure .

Figure 1

Depiction of the fcu topology which is shared by UiO-66, DUT-52, and UiO-67. Also the inorganic brick is common between these materials, which differ only in the organic linker. Zirconium atoms are shown in cyan, oxygen atoms in red, carbon atoms in gray, and hydrogen atoms in white. Reprinted from ref (51). Copyright 2016 American Chemical Society.

All these materials are composed of Zr6(μ3-O)4(μ3–OH)4 bricks connected by ditopic organic ligands to form a network with fcu topology. The three MOFs differ only in the organic linker, which is benzene-1,4-dicarboxylate (BDC) in UiO-66, 2,6-naphthalenedicarboxylate (NDC) in the case of DUT-52, and biphenyl-4,4′-dicarboxylate (BPDC) in UiO-67. The influence of the linker on methane adsorption was studied experimentally by Cavka et al.[52] It was found that the maximum CH4 gravimetric loading (i.e., uptake per mass of the framework) is ordered as UiO-66 < DUT-52 < UiO-67, which is the same ordering as for the surface areas and pore volumes of these materials.

We simulated the methane uptake for pressures between 0.1 and 80 bar (the experimental range) using the force fields described in the previous section. All simulations are performed with a rigid framework, where the geometry is obtained from single-crystal X-ray diffraction by Lillerud et al.[53] for UiO-66 and UiO-67 and by Senkovska et al.[50] for DUT-52. The experimentally measured excess uptake is converted to an absolute amount of adsorbed gas as done in the original work reporting the experimental results. The resulting absolute adsorption isotherms are shown in Figure together with the experimental curves. Hereafter, we discuss these isotherms and their agreement with the experiment both at low and higher pressures.

Figure 2

Comparison of experimental and simulated gravimetric adsorption isotherms of methane in a series of isoreticular MOFs. Experimental data for (a), (c), (d), (e), and (f) from Cavka et al.,[52] experimental data for (b) from Walton et al.[54]

We start with a discussion of methane adsorption in UiO-66 at relatively high pressures (from 30 to 80 bar), which is reported in Figure a. As reported in earlier work, the Dreiding-UFF/TraPPE model offers a good agreement with experiment while the UFF/TraPPE model leads to an overestimation of adsorbed molecules.[18] The MM3-MBIS uptake is substantially lower, while both ab initio force fields (MEDFF) and (SAPTFF) predict uptakes substantially higher than the experiment. Perhaps surprisingly, the ab initio force fields are outperformed by a generic force field such as Dreiding-UFF/TraPPE, even though the latter is in no way fitted specifically for the system studied here. It has however been shown before that DFT-derived force fields are not necessarily superior to a generic force field for describing noble gas adsorption in MOFs.[55] A possible explanation for the apparent poor performance of MEDFF and SAPTFF lies in the description of the guest–guest interactions, which play an important role at higher pressures and thus higher loadings. Calculation of the vapor–liquid coexistence curve of bulk methane shows that both SAPTFF and MEDFF overestimate the liquid-phase density, whereas TraPPE offers a very good agreement with the experiment (see section S9). To investigate the influence of the guest–guest interactions on adsorption isotherms at higher pressures, we also performed simulations where the host–guest interactions are described using MEDFF, while the guest–guest interactions are described using TraPPE. The only difference between Dreiding-UFF/TraPPE and MEDFFF/TraPPE is thus in the description of the host–guest interactions. As seen in Figure a, the uptake predicted by MEDFF/TraPPE is indeed lower than the uptake predicted by MEDFF; however, changing the guest–guest interactions does not fully resolve the discrepancy with respect to the experiment.

The results presented at higher pressures already show the high sensitivity of the predicted isotherms of the used force field. We now focus on the low loading regime, where uptake is almost exclusively determined by host–guest interactions. In Figure b, the methane adsorption in UiO-66 is plotted for pressures lower than 1 bar. Experimental data are taken from Walton et al.[54] because Cavka et al.[52] only report one data point below 1 bar. Note that both experimental sources agree reasonably well for common pressures. Dreiding-UFF/TraPPE reproduces fairly well isotherms at higher pressures but, as all the other force fields, predicts a too large uptake at low pressures. A possible explanation is that the Dreiding-UFF/TraPPE model does not offer a fundamentally proper description of the host–guest interactions, as evidenced by the overestimation at low pressure. At higher pressures, however, this is compensated by deficiencies in the description of guest–guest interactions. In other words, the Dreiding-UFF/TraPPE model might, in this case, benefit from a compensation of errors. Also UFF/TraPPE and SAPTFF largely overestimate the experimental values at low pressure. The MM3-MBIS and MEDFF isotherms are in slightly better agreement with the experiment for pressures below 1 bar but still overestimate the experimental values. The systematic overestimation of all force fields suggests it could be worthwhile to investigate error bars on the experimental data points. Indeed, an ab initio calculation of the Henry constant of methane in UiO-66 is higher than the experimental value (as discussed later in section and Table ). Several pitfalls are present when comparing experimental and computional adsorption,[10] but this discussion will not be pursued in this paper.

Table 3

Henry Coefficient K (mol g–1 bar–1) of Methane in UiO-66 at T = 298 K

Dreiding-UFF/TraPPE	1.78
UFF/TraPPE	2.36
MM3-MBIS	1.21
SAPTFF	1.87
MEDFF	1.08
experiment	0.60
PBE-D3MBJ-ATM	0.91

To investigate whether similar features hold for the other isoreticular Zr-based MOFs, a similar analysis is performed for DUT-52 and UiO-67 for which also experimental data are available. The results are shown in Figure (panels c and e) for higher pressures. All force fields now predict too high uptakes in this pressure regime. Despite the similarity of these MOFs with UiO-66, the Dreiding-UFF/TraPPE model also considerably overestimates methane uptake at pressures above 30 bar compared to the experiment. A possible explanation might be the incomplete activation of the samples used to measure adsorption or the presence of defects in these samples. A procedure that has been suggested in literature to empirically remedy this problem is the rescaling of the entire adsorption isotherm. The rescaling factor can be determined as the ratio of the theoretical and experimental surface area[56] or as the ratio of the theoretical and experimental pore volume.[57] We have investigated both of these methods in section S11. The analysis reveals that the suggested rescaling procedure does not systematically improve the correspondence with simulated results. It has been noted before that a simple rescaling is rather artificial, as the rescaling factor can depend on the loading.[58]

Thus, far our assessment reveals that there is no systematic trend in the prediction of the adsorption trends among the three isoreticular MOFs. The only observed feature is that the absolute CH4 uptake at higher pressures is predicted to be the highest for the SAPTFF, followed by MEDFF, UFF/TraPPE, Dreiding-UFF/TraPPE, and MM3-MBIS for all materials. This is not true at low pressures, where in general UFF/TraPPE is the highest, followed by SAPTFF, Dreiding-UFF/TraPPE, MEDFF, and MM3-MBIS. There is however an exception to this rule, as for UiO-66 MEDFF predicts a lower uptake than MM3-MBIS.

Qualitative Comparison of Zr-Based MOFs

An alternative way to assess the various force fields is to compare absolute CH4 uptake for a series of materials on a qualitative basis. Therefore, we have extended our study with two other materials belonging to the Zr-based family (i.e., sharing the same inorganic brick) but with distinct other features, namely NU-1000 and MOF-808. For these latter two materials, no experimental values are available for CH4 adsorption. The topology and linkers of NU-1000 (csq topology and tetratopic linker) and MOF-808 (spn topology and tritopic linker) are shown in Figure and contrasted with the UiO-66 framework. In Table , we report geometric properties of these Zr-based MOFs: the density, surface area, pore volume, and pore diameters. The following order is obtained when the MOFs are ranked according to pore volume per framework mass:This order is reversed when considering the framework density, which is to be expected as all frameworks are chemically rather similar. When the surface area is considered, the order is basically the same as the order for the pore volume:We conclude that the surface area and pore volume show the same order for the five MOFs we consider. The adsorption isotherms of the five Zr-based MOFs are compared to each other and shown in Figure , for each force field separately. First we study the absolute methane uptake at a relatively high pressure of 30 bar. At these pressures, both host–guest and guest–guest interactions are important to compute the adsorption. We focus on qualitative features predicted by the force fields and therefore determine the ranking of the MOFs at the given pressure. If we order the MOFs according to their high-pressure uptake, Dreiding-UFF/TraPPE and MM3-MBIS predict:while for UFF/TraPPE, SAPTFF, and MEDFF, a slightly different order is obtained:If we exclude MOF-808 from the list, all force fields predict the same order, and this order also corresponds to the ranking with respect to the pore volume and surface area. Note that MOF-808 shows a considerable slope and is far from saturation even at 80 bar, which explains why it appears to be an outlier. At even higher pressures, the uptakes are strictly ordered according to the pore volume as shown in section S10. Such high pressures (above 100 bar) are however not relevant for many applications. The ranking at relatively high pressures is a qualitative feature on which all force fields agree, although it should be mentioned that even a very simple geometric calculation (surface area or pore volume) can suffice to predict this qualitative feature.

Figure 3

Node connectivity, linker and topology for UiO-66, NU-1000, and MOF-808. Reproduced with permission from ref (59). Copyright 2017 Royal Society of Chemistry.

Table 1

Geometric Properties of the Five Zr-Based MOFs Considered in the Qualitative Comparison

	UiO-66	DUT-52	MOF-808	UiO-67	NU-1000
density (g cm–3)	1.238	0.955	0.840	0.725	0.499
pore volume (cm3 g–1)	0.40	0.62	0.85	0.87	1.62
surface area (m2 g–1)	1113	2040	2042	2949	3217
small pore diameter (Å)	7.3	8.6	4.8	10.1	9.8
large pore diameter (Å)	8.8	9.3	18.4	13.0	29.1

Figure 4

Comparison of methane adsorption in different Zr-based MOFs. The dashed vertical line indicates P = 30 bar.

Influence of Defects in UiO-66

The importance of missing linker defects in UiO-66 has been studied both experimentally and computationally.[12,60,61] Recently, Lillerud et al. demonstrated that missing cluster defects are predominant and have a large impact on nitrogen adsorption.[62] Computations confirm that reo structures (with missing cluster defects) profoundly affect CO2 adsorption.[63]

To assess the influence of defects on methane adsorption computationally, we considered four defect structures of UiO-66. All adsorption simulations are performed using rigid frameworks, with geometries optimized using DFT.[64] Defect structures of UiO-66 can be created by removing BDC linkers from the pristine material, and this can be done in many ways depending on the number and position of removed BDC linkers.

The first defect structure (“1 missing linker”) has one missing linker in the conventional 4-brick unit cell and this leaves two bricks 11-fold coordinated while the two other bricks remain 12-fold coordinated. This structure is classified as (11,11,12,12), according to the nomenclature recently introduced in literature.[64] The second structure (“2 missing linkers”) has two missing linkers in the conventional 4-brick unit cell making all four bricks 11-fold coordinated and is classified as the (11,11,11,11)4 structure. The third structure (“3 missing linkers”) has three missing linkers in the conventional 4-brick unit cell, leaving two bricks 9-fold coordinated and two bricks 12-fold coordinated (9,9,12,12)224. Next to these linker defect structures, we also consider a “1 missing cluster” defect. Here, one of the four Zr-bricks of the unit cell is removed together with all connected linkers, which is also referred to as the reo phase. All defects are terminated using a formate group, which ensures that no coordinatively unsaturated metal sites are introduced. A schematic representation of the defects structures considered here is shown in Figure S12.

In Figure , we plot the gravimetric methane adsorption isotherms of these defect structures. This allows one to show that all force fields share the same qualitative features. More specifically, in all cases, the introduction of linker defects has a small influence on methane adsorption, with deviations with respect to the pristine material being lower than 10%. In general, all force fields predict a slightly lower methane uptake at lower pressures compared to the pristine material. The introduction of a missing cluster defect has a bigger impact and leads to a considerably lower gravimetric uptake at low pressures and considerably higher gravimetric uptake at higher pressures, again relative to the pristine UiO-66 structure. Although there are important differences between the force fields (in absolute amount of methane adsorbed and in the pressure at which the “1 missing cluster” isotherm crosses the “pristine” isotherm), in this case all computations lead to qualitatively similar conclusions.

Figure 5

Adsorption isotherms of methane in defective UiO-66 structures.

Single Molecule Adsorption Energies in UiO-66 at the ab Initio and Force-Field Level

In the previous section, we have shown that distinct force fields can predict very dissimilar isotherms. We now try to obtain more fundamental insight at the atomic level by studying the interaction of a single methane molecule with the framework making use of ab initio calculations. It is expected that such analysis yields insight into the adsorption features at low pressures, because in the low-pressure regime, the concentration of guest molecules is so low that host–guest interactions completely determine the uptake. In this section, we focus on methane in pristine UiO-66.

We computed the adsorption energy Eads of a methane molecule in UiO-66 using periodic DFT at the PBE level-of-theory with D3MBJ-ATM dispersion corrections. The adsorption energy Eads is defined as the energy of the host+guest system minus the energy of the separate host and guest systems:To validate the used level-of-theory, we performed a series of CCSD(T)/CBS calculations for a set of methane-terephthalic acid dimers, which is a good model for the interaction between the guest molecule and the framework BDC linkers. The set of 80 dimer configurations is extracted from a GCMC simulation employing MEDFF. Counterpoise-corrected PBE-D3MBJ-ATM/aug-cc-pVTZ interaction energies show an RMSD of 0.39 kJ mol–1 with respect to the CCSD(T)/CBS reference. This RMSD is much smaller than errors in the force-field interaction energies, which justifies the use of PBE-D3MBJ-ATM as the reference level-of-theory for the periodic DFT calculations. More details on the dimer calculations are provided in section S4.1.

An important point to discuss is the choice of configurations of the methane molecule in the periodic UiO-66 framework, as it is necessary to sample all energetically favorable sites. An efficient method to generate such configurations, is to extract snapshots from a GCMC simulation. It has been shown that generic force fields might not sample all relevant portions of the potential energy surface (PES).[65] We therefore extracted 100 configurations from snapshots of GCMC simulations (at 298 K and 1 bar) using all five force fields considered in this work, leading to a total set of 500 configurations. For the united-atom models, only the position of the carbon atom of methane can be extracted from the GCMC simulations: the orientation of the hydrogen atoms is in these cases determined randomly.

In Figure , we show scatter plots of the adsorption energies for the five force fields. It is important to stress that none of the considered force fields are fitted specifically to reproduce the ab initio data presented here, not even the ab initio derived force fields SAPTFF and MEDFF. In each case, the adsorption energies obtained with the target force field are compared with the PBE-D3MBJ-ATM data. Each plot shows the 100 configurations sampled from a Dreiding-UFF/TraPPE GCMC simulation as blue dots as well as the 100 configurations sampled from a SAPTFF GCMC simulation as red dots. The RMSD for all five force fields is also indicated in the plots and varies from 1.31 kJ mol–1 for MEDFF to 3.24 kJ mol–1 for UFF/TraPPE for the configurations sampled from the Dreiding-UFF/TraPPE GCMC simulation (blue dots). When considering the configurations sampled from the SAPTFF GCMC simulation (red dots), the deviations with respect to the ab initio data are notably larger, which indicates that SAPTFF indeed samples different regions of the PES. The RMSD values are the lowest for the two ab initio force fields. In the case of UFF/TraPPE, Dreiding-UFF/TraPPE, and SAPTFF, a large fraction of points is found below the diagonal, meaning that the configurations are too much stabilized compared to the ab initio adsorption data, and this will lead to an overestimation of the adsorption isotherm at low pressure. In Figure b, it is indeed indicated that these three force fields predict the highest methane uptake in UiO-66 at low pressures. Similar figures for all other force fields and a table summarizing all errors are provided in section S5.

Figure 6

Scatter plots of adsorption energies for 200 configurations of methane in UiO-66 (blue dots, sampled from Dreiding-UFF/TraPPE GCMC; red dots, sampled from SAPTFF GCMC). Force-field energies are compared with the PBE-D3MBJ-ATM results.

We visualize the potential energy surface of methane in UiO-66 by plotting the isosurface at Eads = −8 kJ mol–1 for Dreiding-UFF/TraPPE (left) and SAPTFF (right) in Figure . Regions enclosed by this isosurface show adsorption energies that are more favorable than −8 kJ mol–1. For SAPTFF, the adsorption energy at a point r is calculated as a rotational average:where N = 100 random rotations are considered and with T = 298 K. Because TraPPE describes methane as a single site, this rotational averaging is not necessary for Dreiding-UFF/TraPPE. Clearly, SAPTFF predicts a larger portion of the tetrahedral pores to be favorable adsorption sites. This can be quantified by considering the volume fraction of adsorption sites that are more stabilized than Eads = −8 kJ mol–1 as shown in Table . The difference in potential energy surface, as quantified in Table , explains for example why MEDFF/TraPPE predicts a higher uptake than Dreiding-UFF/TraPPE at pressures higher than 30 bar, as observed in Figure a. Indeed, because these two force fields share a common description of the guest–guest interactions, the difference in the number of attractive adsorption sites completely explains the different uptake at higher pressures.

Figure 7

Isosurface (Eads = −8 kJ mol–1) of methane in UiO-66 for Dreiding-UFF/TraPPE (left) and SAPTFF (right).

Table 2

Volume Fraction of Adsorption Sites More Stable than Eads = −8 kJ mol–1 in UiO-66

force field	volume fraction
Dreiding-UFF/TraPPE	10.9%
UFF/TraPPE	11.8%
MM3-MBIS	11.5%
SAPTFF	14.3%
MEDFF	12.9%

Henry Coefficients

At low pressures, only host–guest interactions determine the uptake as the guest–guest interactions are unimportant at low loadings. This means there is a close correlation between single molecule adsorption energies and uptake at low pressure, which will be investigated hereafter by means of the Henry coefficient. At sufficiently low pressures, the uptake is proportional to the pressure, and the proportionality factor is called the Henry coefficient K. The value of K for methane computed with the five force fields is given in Table as well as the value estimated from low-pressure experimental data[54] and the ab initio computed value.

The results can be linked to the comparison of ab initio and force-field single molecule adsorption energies, which are summarized in Figure by showing the RMSD (left) for each force field for the data set of 500 configurations of methane in UiO-66. The RMSD values for the ab initio derived force fields MEDFF and SAPTFF are the smallest followed by MM3-MBIS, while the generic force fields Dreiding-UFF/TraPPE and UFF/TraPPE perform significantly worse. When considering adsorption, it is also important to study the mean deviation (MD): a systematic overbinding will lead to a much higher predicted uptake and Henry coefficient than a systematic underbinding, while both scenarios can give rise to the same RMSD value. The MD is shown on the right-hand side of Figure , and from this we conclude that MEDFF and MM3-MBIS offer the smallest MD for the periodic data set, which is in line with the observations about the Henry coefficients.

Figure 8

Errors of force-field single molecule adsorption energies with respect to ab initio reference data for a data set of 500 configurations of methane in the UiO-66 framework.

Although the ab initio computed Henry coefficient is closer to the experimental value than any of the force fields we studied, there is still a significant discrepancy. Next to possible deficiencies of the ab initio method, there are also experimental uncertainties on the reported values. These originate from imperfections or shortcomings in the experimental set up, such as slow adsorption kinetics, incomplete sample activation, external surfaces of the crystal sample, and the presence of defects.[10] A more extensive discussion on experimental uncertainties is not within the scope of the current work.

The large variations in the Henry coefficients predicted by the five force fields merit special attention, in view of the relatively small errors noticed in the force-field single molecule adsorption energies (Figure ), which are below of 1 kcal mol–1 (smaller than the threshold for chemical accuracy). We therefore propose to perform a sensitivity analysis of the Henry coefficient (or uptake at low pressures). We investigate the influence on the Henry coefficient by slightly varying the parameters of the Dreiding-UFF/TraPPE force field. The applied procedure for the sensitivity analysis is outlined in section S8. We illustrate the outcome of this analysis in Figure , where we correlate the change in the adsorption energies with the relative change on the Henry constant (in both cases with the original Dreiding-UFF/TraPPE model as a reference). By changing the Dreiding-UFF/TraPPE parameters in such a way that adsorption energies change with an RMSD of kJ mol–1, the Henry coefficient (and thus the uptake predicted at low pressures) changes by as much as 40%. For slightly larger deviations on the RMSD (but still less than 2 kJ mol–1), this number can rise above 80%. This extreme sensitivity is not completely unexpected, as the uptake at low pressures is proportional to the Boltzmann factor e–. In other words, small changes in the potential energy surface are exponentially amplified in the corresponding predicted uptake.

Figure 9

Correlation plot to investigate sensitivity of the Dreiding-UFF/TraPPE force-field parameters.

Conclusions and Outlook

Within this paper, we studied the adsorption of methane in a series of Zr-based MOFs with the aim to critically assess the sensitivity of a diverse set of force fields to produce isotherms and single molecule adsorption energies. The selected materials include UiO-66, UiO-67, DUT-52, NU-1000, and MOF-808 and show distinctly different properties in pore volume and surface area, which is reflected in the uptake of methane. As generally known in literature, isotherms are very sensitive to the applied force field. However, to further unravel the physical origin of the observed correspondence between experimental and theoretical isotherms, we performed a systematic investigation of single methane adsorption energies using five different force fields and compared them with adsorption energies produced with periodic density functional theory data. To this end, 500 different configurations of methane adsorbed in the UiO-66 material were taken from the GCMC calculations. We find that some generic force fields such as UFF/TraPPE give an acceptable agreement with the experiment in the UiO-66 framework for pressure between 30 and 80 bar. However, these force fields fail to reproduce accurately single molecule adsorption energies (errors larger than 4 kJ mol–1 are found), and the good correspondence between theory and experiment for the isotherm at these higher pressures should be ascribed to a fortuitous cancellation of errors. The two ab initio derived force fields, SAPTFF and MEDFF, yield a remarkable accuracy of the individual adsorption energies with deviations of less than 2 kJ mol–1 on an overall adsorption energy (RMS value) of 15 kJ mol–1 for methane in UiO-66. Still, such accuracy does not guarantee a quantitative reproduction of adsorption isotherms, since the uptake at low pressures is proportional to the Boltzmann factor e–, and thus the errors are exponentially amplified in the corresponding predicted uptake. The required accuracy for single molecule adsorption energies in order to quantitatively reproduce the isotherms in the low pressure limit is very hard to achieve with current available force fields and even ab initio methods. Apart from these quantitative differences between various methods, we find that all five force fields yield overall similar trends for the reproduction of isotherms of methane in the five different materials, which is rewarding since this validates the common usage of force field based GCMC calculations for the study of adsorption isotherms in the field of nanoporous materials. Yet when constructing force fields for computational simulation of adsorption data, it is advisible to start from ab initio derived force fields as they succeed in reproducing single molecule adsorption energies with high accuracy, which underlines the proper inclusion of host–guest interactions in these models. For higher pressures, a specifically designed force field such as TraPPE may be useful to describe guest–guest interactions, as it was specifically designed for the description of phase equilibria. In order to generalize the conclusions found here, it might be useful to extend the current study to polar adsorbates and frameworks with coordinatively unsaturated sites, which introduce specific host–guest interactions.