Ab Initio Evaluation of Henry Coefficients Using Importance Sampling.

We present a new algorithm that allows for an efficient evaluation of the Henry coefficient of a guest molecule inside a porous material, which permits to use ab initio energy calculations. The Widom insertion method, which is currently used to compute these Henry coefficients, typically requires millions of energy evaluations. Our new methodology reduces this number by more than 1 order of magnitude, enabling the use of an ab initio potential energy surface. The methodology we propose is reminiscent of the well-known importance sampling technique which is frequently used in Monte Carlo integrations. First, a conventional Widom insertion simulation is performed using a force field. In the second step, the Widom results are used to select a limited number of configurations and only for these configurations the ab initio evaluation of the energy is required. Finally, by appropriately reweighting the latter energies, an accurate estimation of the ab initio Henry coefficient is possible at a moderate computational cost. We apply our methodology to the adsorption of CO2 in Mg-MOF-74, a prototypical case where interactions of a polar guest molecule with unsaturated metal sites dominate the adsorption mechanism. In this case generic force fields such as UFF or Dreiding are inappropriate and the use of ab initio methods is indispensable. In a second case study, we compute Henry coefficients of methane in UiO-66 using different levels of theory. We pay particular attention to the influence of the dispersion corrections and the role of many-body effects. For the final example, we qualitatively investigate adsorption features for a series of functionalized UiO-66 frameworks. Overall the cases we present show that accurate computations of Henry coefficients is extremely challenging, as different levels of theory provide strongly varying results. At the same time ab initio calculations have added value compared to force fields, as they provide a physically more sound description of the adsorption mechanism and in some cases clearly improve correspondence with experiment.

Introduction

The separation and storage of gases are important industrial processes in modern-day society. Example applications include capture of the greenhouse gas CO2, storage of methane or hydrogen gas for cleaner energy sources, and separation of toxic gases such as CO or ammonia to reduce environmental pollution. Porous materials (activated carbons, silicas, or zeolites) are often used in this respect, but recently much attention has been devoted to a new class of materials built from inorganic metal nodes and organic linkers: metal–organic frameworks (MOFs).[1] Because of their tunable structure and modifiable functionality, MOFs are very promising materials for the gas storage and separation applications mentioned earlier.[2] This paper deals with the computational analysis and prediction of the gas adsortion properties of MOFs, although the presented methodology is applicable to any porous medium.

There are many studies in which the adsorption of small guest molecules in MOFs has been simulated, often in a high-throughput fashion where calculations are performed for all structures in a database. Chung et al. constructed the CoRE MOF database and performed Grand-Canonical Monte Carlo (GCMC) simulations of methane adsorption for all of the more than 4700 structures.[3] Li et al. considered the same database and calculated the Henry coefficient (which completely characterizes adsorption at low pressures) for CO2, H2O, and N2 in order to perform an initial screening for MOFs with a high selectivity.[4] In a study by Banerjee et al., the Henry coefficient of noble gases Xe and Kr in 125 000 MOF structures has been computed.[5] All of these high-throughput studies have in common that a generic force field (such as UFF[6] or Dreiding[7]) is used in order to keep the computational effort under control.

It has been demonstrated that generic force fields provide results in agreement with experimental adsorption isotherms for many MOFs.[8−10] At the same time there are some specific cases in which these classical force fields fail to properly describe the adsorption mechanism.[11] Grajciar et al. performed a systematic investigation of different guest molecules adsorbing on Cu2+ and Fe3+ clusters, which are models of coordinatively unsaturated metal sites (CUS).[12] The generic UFF potential performed surprisingly well in the description of interactions with organic linkers. The adsorption at CUS was however underestimated by more than 10 kJ mol–1, and in such cases the use of ab initio methods is required. In the same study it was noted that the accuracy of often-used levels of theory, such as dispersion-corrected density functional theory (DFT), depends on the specific combination of adsorbate and adsorbent. This has been linked to the inaccuracy of DFT for the electronic structure of (open-shell) metal cations.[11]

In the context of CUS, the MOF-74 structure is an often studied example. For instance, Valenzano et al. investigated the adsorption of CO, N2, and CO2 in Mg-MOF-74 both computationally and experimentally.[13] Adsorption enthalpies are calculated by adding thermal corrections and zero-point energies to the electronic energy of the adsorbed complex extracted from periodic, dispersion-corrected B3LYP calculations. The experimental adsorption enthalpy of CO2 is −47 kJ mol–1, in reasonable agreement with the calculated value of −38 kJ mol–1. Similarly, Kundu et al. studied CO and N2 in the same MOF making use of ab initio calculations.[14] Even in this case study, in which only one structure is considered, performing ab initio GCMC simulations is computationally not feasible due to current hardware limitations. Those authors therefore propose a model in which the Gibbs free energy of a single adsorption site is computed using ab initio methods and combine this with the empirical Bragg-Williams/Langmuir model in order to compute an adsorption isotherm, which was in close agreement with experiment. Similarly Alonso et al. employed a combination of periodic DFT calculations with the dual-site Langmuir model in order to investigate the effect of SO2 poisoning on the adsorption of CO2 in Mg-MOF-74.[15]

Such an approach however is only appropriate on condition that there are a few clearly identifiable adsorption sites. Additionally, Lee et al. showed that, even if this is the case, the correlation between heats of adsorption and the Henry coefficient is not perfect.[16] This demonstrates the necessity for a complete ab initio evaluation of the Henry coefficient, and hereafter we give a short discussion on models aiming to do so. Peirs et al. reported ab initio computed Henry constants, separation constants and heats of adsorption for N2 and O2 in faujasite.[17] By using a cluster model for the cage and exploiting symmetry, the necessary integrals were computed on a cubic grid at the HF/6-31G* level of theory. Chen et al. also computed adsorption energies on a cubic grid, using methane in CuBTC as a case study.[18] To keep the computational cost under control, no rotational averaging was performed, or in other words, methane was approximated as a spherical molecule. By combining the grid (containing guest-framework interactions) with an empirical model for guest–guest interactions, an isotherm was constructed in better agreement with experiment than corresponding results obtained from classical force fields. Lee et al. recently proposed a biased Widom insertion method and applied it to compute the Henry coefficient for CO2, N2, CH4, and C2H2 in Zn-MOF-74 and Mg-MOF-74.[16] This methodology relies on partitioning the framework into three regions, representing strongest, intermediate, and weakest adsorption domains based on the CO2/MOF binding energies. Insertions are still random, as in the conventional Widom insertion scheme, but biased into the strongly adsorbing regions to quickly achieve good statistics on the required integrals.

Clearly, many methodologies have been developed and applied successfully in the context of gas adsorption in porous materials, specifically MOFs. Most often classical force fields are used because the large number of GCMC or Widom insertion simulations are computationally too expensive to rely on an ab initio potential energy surface (PES). Yet in many cases it would be very interesting, and sometimes indispensable, to study adsorption properties using an approach based on ab initio principles to enhance the accuracy of energy predictions.[19] Next to increased accuracy, an ab initio approach offers the additional advantage that it is more generally applicable and is able describe a wide range of chemical environments. Contrary to many force fields, there is no need for a time-consuming calibration specific to the specific combination of adsorbate and adsorbent. It should also be noted that empirical force fields can sometimes provide results in agreement with experiment, but without correctly capturing the underlying physical interactions.[20] This is another point where an appropriate ab initio level of theory provides added value.

In this work, we present a novel methodology to compute Henry coefficients, based on the concept of importance sampling that is generally applicable. We show that the number of required energy evaluations to compute the Henry coefficient is more than 1 order of magnitude smaller than using conventional Widom insertion. This drastic reduction makes the use of an ab initio PES realistic. First we present the theoretical development of our methodology. Its efficiency will then be demonstrated on some selected test cases.

Methodology

Theoretical Development

The adsorption of guest molecules in a porous material is described by an adsorption isotherm, which relates the number of adsorbed guest molecules to the external pressure (or the chemical potential, which is tightly coupled to the external pressure) at a certain temperature. Such an adsorption isotherm can be computed using for example GCMC simulations. At sufficiently low pressures, however, the number of adsorbed guest molecules is proportional to the external pressure P and the adsorption isotherm is described bywith ρa the volumetric density of adsorbed guest molecules. This is essentially Henry’s law, and therefore KH is called the Henry coefficient. The Henry coefficient KH at a certain temperature completely characterizes the adsorption in the low-pressure regime and is related to the excess chemical potential μex:[21]where . The excess chemical potential of the guest molecule, which is considered to be rigid, in the porous material can be computed using the Widom insertion method.[22]

In this formula ⟨···⟩host denotes a canonical ensemble average over configurations of the porous material, the integration ∫ds runs over all configurations of the guest molecule and ΔU = Uhost+guest – Uhost – Uguest is the adsorption energy. In this work we will employ the often-used rigid framework approximation, which means that the average over host configurations is reduced to a single configuration. The Henry coefficient can then be calculated using the Widom insertion method as follows:In practice, the integral is approximated by generating a random unbiased set of N configurations {s}, calculating the corresponding set of adsorption energies {ΔU} and averaging the integrand:The hat symbol on K̂H indicates that this represents a Monte Carlo estimation of the true value. The presented Widom insertion method requires a high number of configurations that need to be sampled in order to obtain an accurate estimation of K̂H: typically N should be of the order of millions. For pairwise-additive force fields this incurs an acceptable computational cost. For ab initio calculations on the other hand, this is prohibitive (except perhaps for some very small systems). Below we present a novel approach that renders the ab initio estimation of Henry coefficients computationally feasible.

In order to gain more insight into why a large number of configurations is required in the Widom insertion method, we study the sampling variance of the Monte Carlo estimator, which is related to the error bar on the computed integral:The term between brackets will fluctuate over many orders of magnitude because of the presence of the Boltzmann factor, which results in a large variance. Another way to look at this problem is as follows: most of the randomly generated configurations will not be favorable (i.e., ΔU is a positive number that is much larger than ) and will make a negligible contribution to K̂H. On the other hand, the most favorable adsorption sites (i.e., ΔU is a negative number) that contribute the most to the integral, will only be sampled ever so often. This is a problem that is often encountered in random sampling and a well-known solution is to use importance sampling. To this end we introduce a biasing potential energy surface Ũ, of which the adsorption energies are represented by ΔŨ. The expression for the Henry coefficient can now be rewritten aswhere K̃H is the Henry coefficient for the PES Ũ. This approach is only useful if the biasing potential Ũ is computationally cheap (in practice a force field) and the original potential U is computationally expensive (in practice an ab initio method). In this case, the Henry coefficient K̃H can be accurately computed using the Widom insertion method by using a large value for the number of samples N. Ab initio calculations are only required to evaluate the expectation value λ = ⟨exp (− β[ΔU – ΔŨ])⟩. This expectation value is estimated by generating a set of n configurations {s} where the probability to include a configuration is given by p ≈ exp (− βΔŨ):A crucial difference with eq is that, on condition that the biasing potential Ũ and the original potential U are sufficiently similar, the integrand of λ will always be close to unity and therefore show a much smaller variance. The advantage is that here ab initio computation are performed only for a limited number of configurations which are truly important. In other words, the number of evaluations of the ab initio PES U is much smaller than in the Widom insertion method.

The method we present can be interpreted and summarized as follows. We start by computing the Henry coefficient K̃H with a force-field potential Ũ using the Widom insertion method. As explained, this requires a large number of energy evaluations which is however feasible for the computationally cheap force field. Next, we select a fraction of configurations from this simulation where the probability to include configuration α is proportional to exp(− βŨ). In other words, the biasing potential selects favorable adsorption sites (which contribute most to the Henry coefficient) and ensures that the computationally expensive ab initio calculations are only performed for this limited set of configurations that are most important. By computing the average of exp(− β[ΔU – ΔŨ]) for these biased configurations, we can finally estimate the ab initio value of the Henry coefficient using eq .

To conclude, we note that the method outlined above can also be used to compute other quantities of interest. For example the average adsorption energy (related to the adsorption enthalpy as ) can be calculated as follows:Again the expectation values will converge quickly with increasing number of samples on condition that U and Ũ are sufficiently similar. Note that our procedure to calculate the ab initio value of the Henry coefficient also provides all necessary ingredients to compute the average adsorption energy as expressed in eq .

Validation of the Methodology on CH4 Adsorption in UiO-66

We validated our approach by comparing with the numerical evaluation of the integral in eq on a regular grid, considering CH4 adsorption in UiO-66 as a test case. Details about this regular grid, that samples translations of the center-of-mass as well as rotations of the guest molecule, are provided in the Supporting Information. This numerical integration requires a lot of computational power and is therefore not a recommended approach: we use it here for one particular case in order to obtain a reference value for the Henry coefficient, which allows to validate the derivation and implementation of our importance sampling method.

Next to ab initio energies, force-field energies for all configurations considered in the regular grid were computed, thus constructing a database of ab initio and corresponding force-field adsorption energies to which we applied our methodology as outlined above. The ab initio level of theory is PBE+D3(BJ)[23,24] and the necessary periodic single-point DFT calculations are performed using VASP[25−28] version 5.4 using the projector augmented wave (PAW) method,[29,30] with the supplied PAW–PBE potentials. Different force fields are considered for the biasing potential Ũ: MEDFF,[31,32] MM3,[33] Dreiding[7]/TraPPE,[34,35] and UFF6/TraPPE. Note that only host–guest force-field terms are required, as the frameworks and guest molecules are considered rigid and furthermore only one guest molecule is considered so guest–guest interactions are absent. All force-field calculations are performed using the in-house code YAFF.[36] The results are shown in Figure , where the estimated Henry coefficient at room temperature as a function of the number of samples (n in eq ) is plotted in color. We use the notation PBE+D3(BJ)@MEDFF to indicate for instance that the Henry coefficient is calculated for the PBE+D3(BJ) level of theory, using our importance sampling approach with MEDFF as the biasing potential. For each sample size n, the simulation was repeated 1000 times: the dots indicate the average of these 1000 simulations and the error bar indicates the standard deviation. Irrespective of which force field is used as the biasing potential, the average values are always very close to the reference value from the grid-based calculation (shown as a full black line). Concerning the error bars, however, there are striking differences between the force fields. If MEDFF is used as the biasing potential, the error bar is only a few percent even if only about 1000 ab initio calculations are considered. For MM3, slightly larger sample sizes are required to reach similar accuracy, while for UFF and Dreiding the error is still about 10% even if 10 000 ab initio calculations are included. The reason for this has been mentioned in the theoretical derivation: it is required that the original potential U and the biasing potential Ũ are sufficiently similar. Indeed, in our previous work we have shown that MEDFF adsorption energies for this system are in good agreement with ab initio values, while this is not the case for generic force fields such as Dreiding or UFF.[20]

Figure 1

Comparison of our method with a grid-based approach for the Henry coefficient KH of CH4 in UiO-66 at room temperature.

The number of ab initio calculations that need to be considered to reach convergence, will obviously depend on the system of interest. It is not possible to estimate this number a priori, as it depends for instance strongly on the biasing force field. It could however be suspected that the symmetry of the adsorbent will play an important role. For the UiO-66 structure considered here, there are 24 symmetry operators. For a more complicated structure where this symmetry is completely absent, a significantly larger computational effort may be required.

The previous findings would suggest that pure force-field predictions of the Henry coefficient using MEDFF would also give good agreement with the PBE+D3(BJ) value, as both PES agree rather well. However, this is contradicted by the numerical values for the Henry coefficients taken up in Table , which show that on the contrary the UFF Henry coefficient is very close to the ab initio PBE+D3(BJ) result. This apparent paradox is resolved by comparing ab initio and force-field adsorption energies. The RMSD between PBE+D3(BJ) and MEDFF is 2.7 kJ mol–1, considering all configurations of the regular grid for which the PBE+D3(BJ) adsorption energy is negative. The RMSD between PBE+D3(BJ) and UFF on the other hand is 14.5 kJ mol–1, indicating a much worse correlation for this force field with the ab initio reference. Scatter plots (shown in the Supporting Information) also reveal that UFF sometimes provides adsorption energies much more attractive than the reference value, while in many other cases it is severely more repulsive. The reason for the apparent good agreement (concerning the Henry coefficient) is due to a fortuitous cancellation of errors (sometimes overbinding, sometimes underbinding) when averaging is done as in eqs or 12. As indicated in our earlier study on the prediction of adsorption isotherms, this is another case where apparent good results of generic force fields must be treated with caution.[20]

Table 1

Henry Coefficient and Average Adsorption Energy of CH4 in UiO-66 at Room Temperature: Comparison of Ab Initio and Force-Field Results

level of theory	KH [mol kg–1]	ΔHads [kJ mol–1]
PBE+D3(BJ) (grid)	2.23	–19.9
PBE+D3(BJ) (importance sampling)	2.22	–19.9
MEDFF (Widom)	1.17	–18.0
MM3 (Widom)	1.15	–18.4
Dreiding (Widom)	1.72	–19.4
UFF (Widom)	2.24	–19.9

An interesting question is whether a limited number of ab initio calculations, which are needed anyway to apply the importance sampling method, could be used to reparametrize the force field that was used as the biasing potential. This reparametrized force field could then be used to calculate the Henry coefficient and should provide a result close to the ab initio value. We have investigated this approach and report the results in the Supporting Information. In short, we find that it can give acceptable approximations to the ab initio result (although less precise than the importance sampling method), provided the original force field is fairly accurate and features a lot of force-field parameters that are optimized. If this is not the case or one is interested in a very precise ab initio value, our importance sampling methodology is usually preferred because it can be fully automated and systematically converges to the correct value.

From the results presented in this section we conclude, by comparing with a grid-based approach, that our method provides correct results for the ab initio prediction of Henry coefficients. In principle any force field can be used as the biasing potential, but we showed that MEDFF leads to a faster convergence because it predicts adsorption energies that correlate well with ab initio values. Therefore, in the remainder of this work we will always employ MEDFF as the biasing potential. These conclusions were drawn for the PBE+D3(BJ) level of theory, but in the Supporting Information we show that they are valid for other ab initio methods as well.

Results and Discussion

We now apply the developed procedure to two relevant test cases. First we discuss CO2 adsorption in Mg-MOF-74, where we show that in general DFT dispersion corrected calculations improve the agreement with experiment compared to generic force fields. The second case is the adsorption of CO2, CH4, and N2 in UiO-66 and its functionalized variants. There we mainly focus on the comparison of different DFT levels of theory, discussing among others the influence of many-body dispersion correction schemes.

Adsorption of CO2 on Mg-MOF-74

The adsorption of CO2 in MOF-74 frameworks is considered very challenging from a computational point of view. It is a prototypical case where the main adsorption mechanism is due to interactions between coordinatively unsaturated sites of the framework and a polar guest molecule. This typically results in generic force fields underestimating adsorption at low pressure by a few orders of magnitude, because such force fields cannot describe the interactions with open metal sites.[14]

Here we consider the MOF-74 variant with magnesium as metal element, which is referred to as Mg-MOF-74, Mg-CPO-27, or Mg2(dobdc). The Mg2+ metal nodes are linked by 2,5-dioxido-1,4-benzenedicarboxylate (dobdc4–) ligands (as shown in Figure ), with each magnesium atom 5-fold coordinated by oxygen atoms. Variants containing transition metals (Mn, Fe, Co, Ni, and Cu) exhibit multiple magnetic configurations.[37] Additionally, DFT provides a poor description of interactions concerning localized d electrons, which means corrections such as the Hubbard U model are required.[38] To avoid these additional computational complications, we do not consider frameworks containing transition metals.

Figure 2

Graphical representation of the primitive unit cell of Mg-MOF-74, viewed along the ⟨1,1,1⟩ direction. Insets show the dobdc4– ligand and the unsaturated MgO5 polyhedron. Color codes: Mg (orange), O (red), C (brown), and H (white).

We constructed a database of 50 × 106 randomly generated configurations of CO2 in Mg-MOF-74 (geometry from a full optimization of the empty framework, see the Supporting Information) and computed the corresponding MEDFF adsorption energies. Next, in order to test the convergence of the conventional Widom method, 100 different subsets were chosen randomly, each subset containing N samples. For each subset the Henry coefficient and adsorption enthalpy were estimated using eqs and 12. By repeating this procedure for different values of N, we can plot the convergence of the Widom insertion method as a function of the number of considered samples as shown in Figure on the left (in these violin plots the thickness indicates the probability to find the corresponding value, while bars indicate minimal, mean, and maximal value encountered). To obtain a Henry coefficient that is within a range of 10% of the converged result, more than 106 Widom insertions are required. Concerning the adsorption enthalpy, the same number of Widom insertions is sufficient to obtain a result with less than 1 kJ mol–1 uncertainty (with respect to the “correct” result for the given PES). We also applied our importance sampling methodology to this example with MEDFF as the biasing potential. As mentioned earlier, a generic force field is probably a poor option to choose as the biasing potential for this application. We assume that MEDFF will be a better choice, as it will provide a better correspondence with the ab initio PES because it has been shown to be very robust and particularly suitable for the prediction of interaction energies.[20,31] This assumption is confirmed by the obervation that MEDFF predicts Henry coefficients and adsorption enthalpies closer to ab initio results for CO2 in Mg-MOF-74, which will be discussed further on in this section. Then, 20 000 configurations were selected, with probability to include configuration α proportional to exp(− βŨ) where Ũ is the MEDFF adsorption energy. For each of the selected 20 000 configurations the ab initio adsorption energy was calculated using various levels of theory. The electronic energy is calculated using VASP, while Grimme’s dispersion corrections are calculated using the dftd3(39) program (see the Supporting Information for more details). Again, 100 different subsets of n samples were chosen and for each subset the Henry coefficient and adsorption enthalpy were estimated, now using eqs and 15. The results obtained using this importance sampling scheme are shown in Figure on the right for the PBE+D3(BJ) ab initio method. Clearly the results converge much faster than for the conventional Widom simulation: about 10 000 samples are sufficient to obtain an accuracy of 10% on the Henry coefficient and 1 kJ mol–1 on the adsorption enthalpy (again with respect to the fully converged value for the given PES). This demonstrates again that the importance sampling method reduces the required number of calculations by about a factor 100, making it computationally feasible for ab initio methods. Note that we compared the convergence behavior with two different methods (Widom insertions vs importance sampling) but also with two different PES (MEDFF vs PBE+D3(BJ)). Although the shape of the PES can also influence the convergence, we suspect that this has a minor effect considering that MEDFF and PBE+D3(BJ) adsorption energies correlate well. This means that the main difference in convergence is indeed due to the choice of the sampling methodology.

Figure 3

Convergence of the Henry coefficient and adsorption enthalpy of CO2 in Mg-MOF-74 as a function of the number of samples. Left: MEDFF using conventional Widom insertion, right: PBE+D3(BJ) using our importance sampling method. Dotted lines indicate deviations with respect to the converged value.

We now turn our attention to a comparison of the adsorption properties from force-field and ab initio simulations as well as from experiments for CO2 in Mg-MOF-74. An overview of results is provided in Table , where experimental enthalpies of adsorption are obtained as ΔH = −Q st with Q st the isosteric heat of adsorption. Ab initio results are obtained using the importance sampling method with MEDFF as biasing potential while force-field results are obtained using Widom insertion. Each simulation was performed on two different framework geometries: the DFT optimized structure of the empty framework and the DFT optimized structure of the framework loaded with 6 CO2 molecules per unit cell, with each CO2 located at a primary adsorption site. The presence of the CO2 molecules induces subtle changes in the framework geometry as visualized in Figure , where a cutout of the open-metal site in the empty framework is shown, overlaid with the same atoms from the loaded framework shown in gray. A notable difference is the increase of the distance between the Mg atom and the apical oxygen (this is the oxygen opposite of the open side of the metal atom), which is 2.035  Å for the empty framework and 2.062  Å in the loaded framework. The cell parameters change ever so slightly, with the volume going from 1351.5 Å3 in the empty framework to 1350.6 Å3 for the loaded framework. In previous work it was found that the difference in binding energy changes by about 2.5 kJ mol–1 when allowing the framework to relax in the presence of CO2, which we deem significant.[40]

Table 2

Henry Coefficients in mol kg–1 bar–1 and Enthalpies of Adsorption in kJ mol–1 of CO2 in Mg-MOF-74 at 298 K

	geometry empty		geometry loaded
level of theory	KH	ΔHads	KH	ΔHads
PBE+D2	34	–34.2	124	–38.1
PBE+D3(BJ)	48	–34.1	168	–38.2
PBE+D3M(BJ)	112	–37.5	434	–41.4
vdW-DF2	130	–37.8	522	–42.1
MEDFF	62	–37.8	488	–44.8
MM3	28	–31.9	71	–35.1
Dreiding	27	–30.6	45	–32.5
UFF	31	–31.2	47	–32.8

Figure 4

Cutout of the open-metal site in the empty framework (Mg, orange; O, red; C, brown). The same atoms in the framework loaded with CO2 are shown in gray. The arrow indicates the distance between the metal atom and the oxygen atom opposite of the open side of the metal atom.

The computed Henry coefficients and enthalpies of adsorption for both geometries are compared in Table (upper). Despite the seemingly small differences in atomic positions, the impact on adsorption properties is considerable: the Henry coefficients are a factor 3 to 4 larger in the geometry of the loaded framework and adsorption enthalpies are about 4 kJ mol–1 lower, so the effect is slightly larger than in previous work.[40] The values obtained for the loaded framework tend to be in better agreement with experiment as discussed below. The previous observations imply that the rigid framework approach is not valid for this material, despite its application in earlier studies.[16,19,41] An extension of the approach followed here using a hybrid MC/MD approach (where the MD part allows changes in framework geometry due to the presence of guest molecules) might be interesting for future computational studies but is out of the scope of this work.

Finally we discuss the comparison with experiment: several experimental results are compiled in Table (lower). As noted before in the literature, generic force fields (UFF, Dreiding, MM3) underestimate the experimental Henry coefficient by more than 1 order of magnitude and the absolute value of the adsorption enthalpy by about 10 kJ mol–1, irrespective of which framework geometry is considered. MEDFF fares slightly better, and for the loaded framework geometry it is very close to experimental values. All ab initio methods predict higher Henry coefficients than the generic force fields and are thus “closer” to experiment: again, for the loaded framework geometry the agreement with experiment is acceptable. The spread on the results from different levels of theory is also remarkable. For instance, PBE+D3M(BJ)[42] only differs from PBE+D3(BJ)[24] in the data set that was used to fit the damping parameters of the dispersion correction scheme. Yet, the former Henry coefficient is more than twice the value of the latter while the adsorption enthalpies differ by about 3 kJ mol–1. On the other hand, the results obtained with another functional vdW-DF2[43−46] are relatively close to those from PBE+D3M(BJ). The PBE+D2[47] and PBE+D3(BJ) predictions of the adsorption enthalpy of CO2 in Mg-MOF-74 are systematically underestimated with respect to vdW-DF2. This finding has been confirmed by a recent work[40] where the performance of several van der Waals corrected functionals was compared.

It has been suggested in the literature that computed isotherms (and Henry coefficients) should be rescaled to account for the limited availability of adsorption sites in experimental samples, compared to ideal crystal structures used in computations.[14,48] Based on experimental adsorption isotherms, a scaling factor of 0.765 is suggested for the case of Mg-MOF-74.[49] By rescaling the Henry coefficient predicted by vdW-DF2 and MEDFF for the loaded framework, the agreement with experiment improves. In all other cases howevers, the agreement is worse. Note that results in Table are not rescaled.

We conclude that the Henry coefficient and enthalpy of adsorption are rather sensitive quantities, both to the PES and to the framework geometry. This makes a consistent quantitative agreement with experiment very challenging.

Influence of (Many-Body) Dispersion Corrections on CH4 Adsorption in UiO-66

We used our methodology to study the adsorption of CH4 in the pristine UiO-66[53] MOF using various ab initio methods, always with MEDFF as the biasing potential. The force-field value is determined using N = 10 000 000 Widom insertions and the ab initio value using n = 5000 importance sampled configurations. All results reported in this section are calculated using the PBE functional. The electronic part of the adsorption energy is calculated using VASP, Grimme’s dispersion corrections using the dftd3 and many-body dispersion corrections using the Pymbd[54] program. Using pure PBE energies, nearly no CH4 adsorption is predicted because of the lack of (long-range) dispersion interactions in GGA functionals. Many dispersion correction schemes (which add an energy contribution on top of the PBE energy) are available in the literature, and we now present a comparison of such schemes paying special attention to the impact of many-body effects.

The most widely used dispersion schemes are variants of DFT-D as originally proposed by Grimme,[47] which add an r–6 interaction term between all pairs of atoms, appropriately damped at short-range. The method has been refined by adjusting the damping (for instance Becke-Johnson damping D3(BJ)[24]) or by simply reparametrizing on a larger data set (D3M(BJ)[42]). Additionally, three-body interactions have been added using the expression derived by Axilrod–Teller–Muto (ATM).[55] The Henry coefficient and adsorption enthalpy of CH4 in UiO-66 is reported in Table for three variants of the D3 dispersion scheme. For the first column only two-body dispersion corrections are included, while for the second column the three-body ATM term is taken into account. The difference in adsorption enthalpy between PBE+D3 and PBE+D3M(BJ) is 2.5 kJ mol–1, which can be considered significant considering that this means a deviation of more than 10%. The inclusion of the three-body term consistently leads to more repulsive adsorption energies, and the absolute value of the adsorption enthalpy decreases by about 2 kJ mol–1. As a result, the Henry coefficient decreases by about 50% which can be rationalized by considering that exp(− β2.0 kJ mol–1) ≈ 0.45 at room temperature. The experimental value for KH is 0.6  mol  kg–1 bar–1,[56] so the inclusion of three-body dispersion interactions improves the agreement. We conclude that, even for a small guest molecule such as CH4, three-body dispersion interactions have a non-negligible impact on adsorption.

Table 3

Henry Coefficients [ mol  kg–1 bar–1] and Enthalpies of Adsorption [ kJ mol–1] of CH4 in UiO-66 at 298 K Using Various Levels of Theory

	without ATM		with ATM
level of theory	KH	ΔHads	KH	ΔHads
PBE+D3	3.7	–21.6	1.8	–19.5
PBE+D3(BJ)	2.2	–19.9	1.1	–17.8
PBE+D3M(BJ)	1.7	–19.1	0.8	–17.1

The many-body dispersion (MBD) method was developed by Tkatchenko et al.[57] and also adds a van der Waals contribution to the underlying electronic structure calculation. Here, many-body effects up to infinite order are treated using the coupled fluctuating dipole model, which can be much better motivated from a physical point of view than the three-body corrections employed in DFT+DATM schemes. The adsorption energies resulting from the MBD calculations were partitioned into contributions from second up to sixth order, which allowed to calculate the Henry coefficient and adsorption enthalpy including up to Nth order many-body effects for N = 2, ..., 6. The resulting graph shown in Figure reveals that, by adding three-body interactions, the Henry coefficient as well as the absolute value of the adsorption enthalpy decreases, indicating that three-body interactions are generally repulsive. By additionally including four-body interactions, the Henry coefficient increases again which means that this has a generally more attractive effect, in line with derivations from perturbation theory. Finally we note that omitting higher than two-body dispersion effects leads in this case to an overestimation of the Henry coefficient by a factor 2 with respect to the value obtained from the full inclusion of many-body effects up to infinite order. On the other hand, the value quickly converges as soon as three-body and certainly four-body interactions are treated. Again we note that the Henry coefficient computed with many-body dispersion interactions is closer to the experimental value of KH = 0.6  mol  kg–1 bar–1 than the Henry coefficient computed with only two-body dispersion interactions.

Figure 5

Henry coefficient (left) and adsorption enthalpy (right) of CH4 in UiO-66 computed with PBE with MBD corrections, as a function of the order of included many-body effects (dotted lines are only present to guide the eye).

Linker Functionalizations in UiO-66

In the previous section we demonstrated that adsorption properties are rather sensitive to the underlying PES that is used in the simulation. In many cases, especially high-throughput studies, one is more concerned with trends in properties rather than the absolute value of these properties. We therefore now turn our attention to a comparison of the Henry coefficients of CO2, CH4, and N2 in a series of functionalized UiO-66 frameworks. The pristine UiO-66 MOF features benzene dicarboxylate (bdc) organic linkers. By functionalizing the benzene rings with amino (−NH2), nitro (–NO2), methoxy (−2,5-OMe2), and naphthyl (−1,4-Naph) groups (as shown in Figure ), we obtain four additional frameworks.

Figure 6

From left to right: bdc linker of pristine UiO-66, –NH2, –NO2, −2,5-OMe2, and −1,4-Naph functionalized.

Location of Functional Groups

We considered functionalized versions of the UiO-66 by adding a functional group on each organic linker. Because these functional groups can be added on multiple locations of the linker, it is possible to construct many different functionalized frameworks with the same chemical formula. For each functional group, we constructed five different frameworks by placing the functionalization on a random position for each linker. This allows us to investigate whether the precise placement of functional groups has an important effect on the adsorption properties of the framework. We calculated the Henry coefficient of CO2 in all frameworks with randomly placed functional groups using the MM3 force field as shown in Figure (similar plots for other force fields and for CH4 as the guest molecule are included in the Supporting Information). For each of the functional groups, there are five dots in this figure with each dot representing the Henry coefficient in one randomly constructed version of that specific functionalized framework. These results show that the precise placement of the functional groups certainly influences the adsorption properties, as there is a certain spread on the Henry coefficients for the different versions of each functionalized framework. At the same time, it should be realized that the spread is in general smaller than differences between different functional groups. In other words, if one is mostly interested in a comparison of different functionalized groups it suffices to study one specific placement for each group, as we have done here. Finally note that, when we compare results from various levels of theory, the same framework is considered, which enables a consistent comparison of different PES. A proper comparison with experiment is however more difficult, as it is often uncertain whether the simulated functionalized frameworks faithfully represent the ones used in experiment. Although Cmarik et al.[56] performed adsorption experiments on functionalized UiO-66 frameworks, these results will not be discussed here for this reason.

Figure 7

Henry coefficients for CO2 at 298 K in different versions of the functionalized UiO-66 materials, computed with the MM3 force field.

Results

As mentioned before, we focus on qualitative features in this section, and a comparison of force-field and ab initio results already reveals insights in this respect. We calculated the Henry coefficient of CO2, CH4, and N2 at room temperature in all frameworks at the PBE+D3(BJ) and PBE+MBD level of theory using our importance sampling method with MEDFF as the reference potential (N = 10 000 000 Widom insertions and n = 5000 importance samples for CH4 and N2, N = 50 000 000 Widom insertions and n = 20 000 importance samples for CO2). Additionally, the UFF values (a popular choice of PES for high-throughput studies) are calculated using Widom insertion and shown together with the ab initio results in Figure . The PBE+D3(BJ) and PBE+MBD results share a similar pattern, with the PBE+MBD Henry coefficients being generally lower. Following the discussion in the previous section, this can be attributed to the lack of (generally repulsive) many-body dispersion effects in PBE+D3(BJ). For CO2 and CH4, the UFF results are in good agreement with the PBE+D3(BJ) results for UiO-66, UiO-66-NH2, and UiO-66-NO2. For UiO-66–2,5-(OMe)2 and UiO-66–1,4-Naph, however, UFF predicts a steep increase in the adsorption which is not observed in the ab initio results. This indicates that some caution is warranted when using such a generic force field, even if only a qualitative comparison of some similar frameworks is made. MEDFF on the other hand is generally quite close to the ab initio values, especially if the PBE+MBD level of theory is considered as reference. Also trends predicted by MEDFF are in better agreement with the ab initio results.

Figure 8

Henry coefficients at room temperature of CO2, CH4, and N2 in functionalized versions of the UiO-66 MOF.

We finally turn our attention to the selectivities of equimolar mixtures of CO2/CH4, CO2/N2, and CH4/N2, which in the low-pressure regime are simply provided as the ratios of the respective Henry coefficients. If a certain PES is systematically more attractive than another for all guest molecules, such deviations would compensate when looking at selectivities, which makes this an interesting property to look for qualitative differences between PES. Figure reveals for instance that UiO-66-NH2 and UiO-66-NO2 show a higher CO2/CH4 selectivity than UiO-66–2,5-(OMe)2 and UiO-66–1,4-Naph, according to both PBE+D3(BJ) and PBE+MBD. The generic UFF predicts exactly the opposite. Again this a clear example where one should be very careful when dealing with force-field results, even when simply extracting qualitative adsorption properties. For these cases, MEDFF provides patterns that are more similar to the ab initio curves, although some small inconsistencies are still present. This suggests that it might be worth considering ab initio based force fields (such as MEDFF) in favor of generic force fields (such as UFF) in high-throughput screening applications, despite the additional computational cost of the former PES.

Figure 9

Selectivities of equimolar mixtures at room temperature in functionalized versions of the UiO-66 MOF.

Conclusions

In this work, a novel methodology to compute the Henry coefficient of small guest molecules in porous materials was developed, which allows to use ab initio evaluations of the PES. It comprises two stages: the first stage is a conventional Widom insertion simulation (using a computationally cheap force field PES), in the second stage ab initio adsorption energies are calculated only for configurations selected from the first stage using importance sampling. This approach enables an efficient evaluation of ab initio Henry coefficients, provided that force-field and ab initio adsorption energies correlate fairly well. It was demonstrated that MEDFF fulfills this condition in all investigated cases. MEDFF is a recently derived force field which is very robust because only a limited number of parameters are fitted to interaction energies from highly accurate CCSD(T)/CBS calculations. This robustness renders the current method applicable to a wide range of systems.

The adsorption of CO2 in Mg-MOF-74 was investigated as a prototypical case where generic force fields fail to properly describe interactions between a polar guest molecule and open metal sites. Indeed, ab initio Henry coefficients were generally closer to experimental values than force-field results, but the detailed geometric features of the framework were shown to be important, thus hindering a fair comparison between experiment and simulation. Also the influence of the choice of level of theory was remarkably high, confirming the sensitivity of the Henry coefficient. More specifically it turned out that many-body dispersion effects have a significant impact, as shown for adsorption of CH4 in the UiO-66 framework. Finally, qualitative features of adsorption in functionalizationed UiO-66 materials are discussed and here it is observed that generic force fields may predict dissimilar rankings compared to more advanced methods.

In conclusion, we estimate that our methodology will help to investigate adsorption in porous materials whenever force fields are either not available or are not sufficiently accurate. Furthermore, it can also be a useful tool for benchmarking force fields to higher-level methods enabling a more judicious selection of a certain potential for the system at hand.