Prediction of Carbon Dioxide and Methane Adsorption on UiO-66 Metal–Organic Framework via Molecular Simulation

Abstract

The adsorption equilibrium of methane (CH₄) and carbon dioxide (CO₂) on the metal–organic framework (MOF) UiO-66 is studied via molecular simulation. UiO-66 is a versatile MOF with vast potential for various adsorption processes, such as biogas upgrading, CO₂ capture, and natural gas storage. The molecular simulations employ the grand canonical Monte Carlo (GCMC) method, covering a temperature range of 298–343 K and pressures up to 70 bar for CH₄ and 30 bar for CO₂. The accuracy of different forcefields in describing the adsorption equilibria is evaluated. Two modelling approaches are explored: (i) lumping each hydrogen atom in the MOF framework to the heavy atom it is bonded to (united atom approximation) and (ii) considering explicit hydrogen atoms. Additionally, the influence of electrical charges on CO₂ adsorption is also evaluated. The findings indicate that the most effective forcefield to describe the adsorption equilibrium is a united atom forcefield based on the TraPPE parametrization. This approach also yields an accurate calculation of the isosteric heat of adsorption. In the case of CO₂, it is observed that the use of electrical charges enhances the prediction of the heat of adsorption, especially in the low-coverage region.

1. Introduction

Metal–organic frameworks (MOFs) are highly ordered crystalline materials composed of metal ion cluster nodes connected by organic linkers. These frameworks offer versatility through customization by altering the metal-center/organic-linker combination and incorporating specific functionalities. MOFs have been studied for several applications, including gas storage, catalysis, sensing, and drug delivery. Among these applications, the study of adsorption-based separations stands out as one the most interesting fields of study [1,2,3,4,5,6,7,8].

UiO-66, like ZIF-8, ZIF-67, or HKUST-1, stands out as an interesting material within the vast array of MOFs due to its suitability for separation processes. Notably, UiO-66 boasts excellent adsorption properties, alongside remarkable thermal and mechanical stability. Additionally, UiO-66 demonstrates exceptional resistance to various solvents, such as water, dimethylformamide, benzene, and acetone. Even under high mechanical pressures of up to 10,000 kg/cm², UiO-66 retains its structural integrity without compromise [9,10,11].

The UiO-66 framework is based on the linkage of zirconium oxide clusters Zr₆O₄(OH)₄ through 1,4-benzene-dicarboxylate (BDC) ligands in a three-dimensional structure with tetrahedral and octahedral cages of approximately 8 Å and 11 Å, respectively. The windows accessing the cages have apertures of around 6 Å [12,13]. This molecular configuration promotes an exceptional stability due to the high degree of coordination of Zr−O metal centers to the organic linkers [10,11]. Several authors have reported UiO-66 surface areas of 1000–1400 m²/g and pore volumes of 0.4–0.8 cm³/g [11,13,14,15,16,17,18,19,20]. The UiO-66 framework structure is shown in Figure 1, which was generated from the CIF file reported by Yang et al. [21]. It shows an axial view of a simulation box containing 2 × 2 × 2 unit cells of UiO-66.

UiO-66 has high potential as an adsorbent for applications in gas and liquid phase separations. It has been evaluated for the adsorption of hexane and xylene isomers [16,19,22,23], boron [24], arsenic [25], and diverse water contaminants [9]. UiO-66 has also been studied for gas storage [21] and, recently, it was studied for neon adsorption [26].

The potential application of UiO-66 in carbon dioxide (CO₂)/methane (CH₄) separation has been thoroughly studied. Cavka et al. [15] studied the adsorption of CH₄ and CO₂ on UiO-66, between 298 K and 343 K, and up to 80 and 30 bar, respectively. Yang et al. [18] studied the CO₂/CH₄ co-adsorption, concluding that UiO-66 can be a good candidate for pressure swing adsorption (PSA) applications, as it presents a good selectivity and working capacity; also, UiO-66 is easily regenerable.

UiO-66 functionalized materials have also been studied [12,27]. Xian et al. [28] prepared a series of polyethylenimine (PEI)-impregnated UiO-66 samples, showing that this approach can significantly increase the CO₂ adsorption capacity and CO₂/CH₄ selectivity. Jasuja and Walton [14] synthesized a dimethyl-functionalized UiO-66, enhancing the CO₂/CH₄ selectivity at high pressures.

Molecular simulation is a powerful and extensively utilized method for investigating gas adsorption equilibria on nanoporous materials, particularly in the realm of MOFs [29,30,31,32]. Classical forcefields, such as UFF, DREIDING, OPLS, and TraPPE, are commonly employed for these simulations. These forcefields have been previously employed with success to describe several MOFs, such as Cu-BTC, MOF-5, IRMOFs, and MIL materials, as evidenced by previous studies [26,29,33,34,35,36,37].

Regarding UiO-66, Yang et al. [18] used grand canonical Monte Carlo (GCMC) simulations to calculate the pure and mixed CO₂/CH₄ adsorption equilibrium. Additionally, molecular dynamics were employed to study the co-diffusion of mixtures containing both species. The study revealed that CO₂ preferentially occupies the UiO-66 tetrahedral cages, while CH₄ is predominantly confined to the octahedral cages. Remarkably, CO₂, being a slower-diffusing molecule, exhibited the intriguing effect of enhancing the mobility of the faster CH₄ within UiO-66. The same group also used GCMC to study UiO-66 for CH₄ storage. Other authors have used molecular simulation to study the adsorption of diverse species on UiO-66, including NH₃, H₂S, and CO₂, [38] xylene isomers, [22,23] and water [39,40].

There is compelling evidence that UiO-66 is a very promising candidate for the adsorptive separation of CO₂/CH₄. Therefore, in this work, we employ GCMC molecular simulation to predict the CO₂ and CH₄ adsorption equilibria on UiO-66 at three different temperatures: 298 K, 313 K, and 343 K, covering the pressure ranges of 0–30 bar and 0–70 bar, respectively. We evaluate the influence of using various forcefields in our calculations by comparing the theoretical results with the experimental data published by Cavka et al. [15]. In this work, we use the TraPPE parametrization in the description of the crystalline structure as opposed to the approach taken by Yang et al. [18], in which the authors used the DREIDING and UFF forcefields to describe the MOF framework.

2. Methods

Molecular Models and Solid Lattice

In our GCMC molecular simulations, we adopt a rigid model for the UiO-66 framework. This involves nullifying the bond, angle, torsion, and improper torsion terms, while omitting the van der Waals and Coulombic interactions between the framework atoms. This approach is employed to optimize computational efficiency and save valuable computational time. The UiO-66 crystallographic structure is that obtained by Yang and co-workers through X-ray diffraction and further refinements resulting from geometrical optimization based on density-functional theory (DFT) [18].

For the parameterization of the solid–fluid dispersive interactions, the Lennard-Jones (LJ) parameters of the organic linker were taken from those reported for similar organic molecules. Two different modelling approaches are studied. One approach is to unite each carbon of the framework and its bonded hydrogen(s) into a single interaction site; this is the united-atom approximation, here designated as UiO-66-UA (Figure 2a). Another approach is to treat each hydrogen and carbon atom as distinct interaction sites (UiO-66-EH), as shown in Figure 2b. Furthermore, the influence of considering electrical charges was evaluated for the case of CO₂ adsorption.

The TraPPE forcefield is employed in the simulations. In TraPPE, non-bonded interactions are described by pairwise-additive LJ 12-6 potentials and Coulombic interactions of partial charges with unlike LJ interactions computed through Lorentz–Berthelot combining rules.

Methane is modelled as a single chargeless LJ interaction site (Table 1). Carbon dioxide is represented as a three-site molecule. The intrinsic quadrupole moment is described by a base charge model. The CO bond length is 1.16 Å, and the bond angle O-C-O is 180°. As stated above, the non-bonded CO₂ interactions are described by the pairwise-additive LJ 12-6 potentials and Coulombic interactions of partial charges given in Table 1.

The GCMC calculations for CO₂ adsorption on UiO-66 were performed using three different forcefields, namely, the UA and chargeless MOF framework, the UA and partial charges placed on the framework atoms (UAq), and the explicit hydrogen forcefield with electrostatic charges (EHq). The use of these three approaches gives an understanding of the relevance of the hydrogen atoms and the electric charges of the MOF on the adsorption of CO₂ by UiO-66. The electric charges employed in UAq and EHq are those reported by Yang et al. [18]. These partial charges were extracted from density-functional theory (DFT) calculations and the Mulliken charge partitioning method, as reported in [41]; such calculations were based on the PW91 GGA functional combined with the double numerical basis set containing a polarization function (DNP). The LJ interactions were truncated at $r_{\mathrm{cut}}=14$ Å, and analytical tail corrections were applied. The Coulombic electrostatic interactions were computed via the Ewald summation method utilizing a constant number of 5 inverse space vectors), an electrostatic cutoff equal to $L_{\min }/2$, and a damping factor equal to $\alpha =5.6/L_{\min }$, where $L_{\min }$ is the shortest length of the simulation box.

Table 1. Lennard-Jones parameters and partial charges for the TraPPE forcefield.

Label	σ (Å)	ε/kB (K)	q (e)	Ref.
		Methane
CH4	3.73	148	-	[42]
		Carbon dioxide
C_CO2	2.80	27.0	+0.70	[43]
O_CO2	3.05	79.0	−0.35	[43]

Table 1. Lennard-Jones parameters and partial charges for the TraPPE forcefield.

Label	σ (Å)	ε/kB (K)	q (e)	Ref.
		Methane
CH4	3.73	148	-	[42]
		Carbon dioxide
C_CO2	2.80	27.0	+0.70	[43]
O_CO2	3.05	79.0	−0.35	[43]

The LJ parameterization of UiO-66 using the UA forcefield (UiO-66-UA) considers seven pseudo-atoms: Zr, C25, C13, C1, O25, O1, and O29, as illustrated in Figure 2a and indicated in

Table 2. The Lennard-Jones parameters of the UiO-66 framework for the UA forcefield and partial charges for the UAq forcefield. The electric charges are those reported by Yang et al. [18]. The lumped charge for C25 is (−0.121) + (+0.127) = +0.006 (see Table 3).

Label	σ (Å)	ε/kB (K)	Mw (g/mol)	q (e)	Ref.
C25	3.74	48.00	13.02	+0.006	[44]
C1	3.90	41.00	12.01	+0.625	[45]
C13	3.88	21.00	12.01	−0.002	[44]
O1	2.80	55.00	16.00	−0.582	[46]
O25	3.02	93.00	17.01	−1.179	[46]
O29	2.80	55.00	16.00	−0.741	[46]
Zr	2.78	34.72	91.22	+2.008	[47]
H25	0.00	0.00	1.01	+0.495	[39]

. Hydrogen H25 is lumped with the pseudo-atom C25, whereas the vdW contribution of the polar H connected to hydroxyl O25 is neglected. The LJ parameters for the aromatic carbons are taken from the TraPPE-UA parameterization of toluene; the [CH]aro and the [C]aro-CHy correspond to the C25 and the C13, respectively. The C1 parameters assigned correspond to those of a carbon atom belonging to a carboxylic acid/ester group. The parameters for the oxygens O1 and O29 were taken from the TraPPE-UA description of ethers, where the LJ values correspond to a CH_x-[O]-CH_y. The values for oxygen O25 and the hydroxyl H were taken from the TraPPE-UA parameterization of alcohols, CH_x-[O]-H and CH_x-O-[H]. The Zr parameterization was taken from the UFF forcefield.

The UiO-66-EH forcefield uses nine atoms, i.e., the seven atoms used in the UiO-66-UA forcefield, plus H1 and H25 (Figure 2b;

Table 3. The Lennard-Jones parameters of the UiO-66 framework for the EH forcefield and charges for the EHq forcefield.

Label	σ (Å)	ε/kB (K)	Mw (g/mol)	q (e)	Ref.
C25	3.60	30.70	12.01	−0.121	[48]
C1	3.90	41.00	12.01	+0.625	[45]
C13	3.88	21.00	12.01	−0.002	[44]
O1	2.80	55.00	16.00	−0.582	[46]
O25	3.02	93.00	16.00	−1.179	[46]
O29	2.80	55.00	16.00	−0.741	[46]
Zr	2.78	34.72	91.22	+2.008	[47]
H1	2.36	25.45	1.01	+0.127	[48]
H25	0.00	0.00	1.01	+0.495	[43]

). The LJ parameters for atoms Zr, C13, C1, O1, and O29 are the same as those for UiO-66-UA. The parameters of the atoms C25 and H1 are obtained through the TraPPE-EH parameterization of benzene, correspondent to a X(aro)-[C]-(aro)-X(aro) and [H]-C(aro). Due to information scarcity, the O25 and H25 parametrization employed corresponds to the TraPPE-UA of an alcohol.

3. Results and Discussion

3.1. Single-Component Adsorption Equilibria

In this study, we performed GCMC molecular simulations with different forcefields to investigate the adsorption of CO₂ and CH₄ on UiO-66 at three different temperatures: 298 K, 313 K, and 343 K. The pressure ranges examined were 0–30 bar for CO₂ and 0–70 bar for CH₄. The adsorption equilibria results can be analyzed through various interpretations, e.g., excess and absolute quantities. Here, we focus on reporting the adsorbed quantities in terms of the absolute amount adsorbed, as it is a quantity that is directly obtained from the GCMC calculations. The absolute amount adsorbed obtained from molecular simulations is calculated as

$$q=\frac{1}{m_{\mathrm{s}}}\frac{\langle n\rangle }{N_{\mathrm{A}}} ,$$

(1)

where $\langle n\rangle$ is the ensemble average number of adsorbate molecules in the simulation box, $m_{\mathrm{s}}$ is the adsorbent mass in the simulation box, and $N_{\mathrm{A}}$ is Avogadro’s constant.

Figure 3 compares the absolute experimental adsorption isotherms of CH₄ reported by Cavka et al. [14] at 298 K, 313 K, and 343 K, with the simulation results obtained using the UiO-66-UA and UiO-66-EH forcefields. Surprisingly, both forcefields demonstrate excellent agreement with the experimental CH₄ adsorption data, despite covering a wide range of temperatures and pressures. It is also seen that the UiO-66-UA forcefield gives slightly better predictions, although the differences between both forcefields are not very significant.

Looking at the lower-pressure region (0–10 bar) represented in Figure 3b,d,f for 298 K, 313 K and 343 K, respectively, it can be observed that, in this region, the predictions using the UiO-66-UA or UiO-66-EH forcefields result in minimal differences in the adsorption data calculated. As the temperature increases, larger differences between the simulation data and experimental results are observed.

Figure 4 compares the CO₂ adsorption equilibrium data measured experimentally by Cavka et al. [15] and our molecular simulation data obtained using the three different forcefields. For the sake of clarity, each experimental adsorption equilibrium isotherm is compared to the predicted GCMC data in a single graph (Figure 4a–c for 298 K, 313 K, and 343 K, respectively). The results obtained show that the UiO-66-UA forcefield, which neglects the electric charges (UiO-66-UA) in the solid, is the more accurate forcefield for the adsorption equilibrium predictions within the pressure and temperature ranges studied. It is widely acknowledged that simulated adsorption data are highly sensitive to the charge parametrization used. In cases where the partial charges are not well-tuned for a specific solid–fluid system, their impact on the results can be more detrimental than if they were entirely neglected, as seems to be the case here.

A GCMC simulation of methane condensation at its normal boiling point, $T_b=111.65\mathrm{K}$, within the porous structure of the MOF gives a maximum methane loading of ca. 0.0118 mol/g. By multiplying this value with the molar volume of methane at its normal boiling point, $v_m=37.98{\mathrm{cm}}^3/\mathrm{mol}$ converts it into an equivalent pore volume. The obtained value of 0.45 cm³/g is in agreement with the theoretical pore volume of pristine UiO-66, but larger than the experimental pore volume of 0.36 cm³/g reported by Cavka et al. [15] for their sample, as measured from argon adsorption at 87 K ($p/p_0=0.5$). This suggests that, based solely on gravimetric considerations, one might expect the simulated saturation adsorption capacity to surpass the experimental one by a factor of approximately 0.45/0.36. However, this discrepancy is not observed in our comparative analysis of the methane and carbon dioxide data.

It is important to note that UiO-66 is known to be susceptible to defects. Several authors [49,50,51] have argued that most, if not all, UiO-66 samples contain a significant number of missing linkers and/or missing cluster defects. Furthermore, research has shown that the presence of such defects can significantly impact the adsorption isotherms [40,52]. Given the likelihood that the sample of Cavka et al. [15] may also contain defects, it is somewhat surprising that the simulations, assuming a perfect crystal, align so closely with the experimental measurements.

The GCMC data obtained for both adsorbates show that explicitly considering the hydrogen atoms in the MOF modelling does not increase the accuracy of the simulated adsorption equilibria. This is observed more clearly in the calculations using the forcefield UiO-66-EHq for CO₂ (e.g., 298 K, Figure 4a,b), where the isotherm does not have the same trend as the experimental data; i.e., the calculations overestimate the experimental values at lower pressures and, at a higher pressure, the obtained simulations underestimate them. The simulation results indicate that the explicit description of the UiO-66 electrical charges does not benefit the prediction of CO₂ adsorption equilibrium. This indicates that the electrostatic parameterization obtained by Yang and co-workers [18,21] through quantum mechanics calculations is not suitable for being combined with the TraPPE forcefield.

The calculations show that, similarly to CH₄, the best forcefield to replicate the CO₂ experimental data is the UiO-66-UA. Figure 5 presents an overall comparison between the experimental data and the simulation predictions for both CO₂ and CH₄ using the UiO-66-UA forcefield for the three temperatures studied. The data in Figure 5b show that, with the increase in temperature, the CO₂ simulated data tend to slightly overestimate the experimental results.

3.2. Adsorption Energetics

The isosteric heat of adsorption, $Q_{\mathrm{st}}$, quantifies the strength of the interactions between the adsorbate molecules and the MOF solid lattice. In this work, the heat of adsorption is obtained directly from statistical thermodynamics through the GCMC simulations using the co-variance formulation of Nicholson and Parsonage [53],

$$Q_{\mathrm{st}}=\frac{\langle U\rangle \langle n\rangle -\langle Un\rangle }{\langle n^2\rangle -\langle n\rangle \langle n\rangle }+k_{\mathrm{B}}T$$

(2)

where the triangular brackets denote an ensemble average, $n$ the number of adsorbate molecules in the simulation box, and $U$ the configurational energy. The calculations rely on numerous particle insertions and deletions in the grand canonical ensemble. As a consequence, performing long simulations becomes necessary to achieve accurate statistics for calculating averages, particularly when dealing with very low loading conditions where the number of molecules is significantly reduced [54].

The isosteric heats calculated with the UiO-66-UA forcefield are presented in Figure 6 for both gases studied. Regarding the theoretical $Q_{\mathrm{st}}$ values for CH₄, they range from 16 to 18 kJ/mol, slightly decreasing up to a loading of 2 mol/kg and then increasing the $Q_{\mathrm{st}}$ up to the maximum loading values studied (~6 mol/kg). The initial decrease in the isosteric heat of adsorption with loading suggests the influence of adsorbent heterogeneity [55].

The $Q_{\mathrm{st}}$ calculated from the experimental data reported by Cavka et al. [15] compares very well with the theoretical calculations, although the initial decrease in $Q_{\mathrm{st}}$ is not observed. Figure 6 also presents the results for $Q_{\mathrm{st}}$ as a function of the CO₂ loading on UiO-66. The theoretical $Q_{\mathrm{st}}$ values obtained through simulations not considering the electrical charges are within 20 and 27 kJ/mol and show a plateau at low loadings (up to 2 mol/kg) and, for higher CO₂ loadings, an increase is observed. The isosteric heat of adsorption calculated from the experimental data presents the same trend as the theoretical values with differences smaller than 2 kJ/mol. Additionally, calculations considering the electrical charges were also performed, and the obtained results show that this approach describes more accurately the heat of adsorption derived from the experimental data at low coverage, as depicted in Figure 6. The study of the heat of adsorption corroborates the success of the forcefields employed to replicate the experimental system.

4. Conclusions

In this work, we report the successful prediction of the CH₄ and CO₂ adsorption equilibrium on UiO-66 MOF using GCMC molecular simulations with parametrizations based on the TraPPE forcefield. We cover pressures up to 70 bar for methane and up to 30 bar for carbon dioxide; the temperature range analyzed was 298–343 K for both species.

We successfully employed two models: one considering explicit hydrogen atoms (EH) and the other lumping the MOF’s hydrogen atoms to the heavier atom it is bonded to (UA). Both approaches lead to excellent agreement with the experimental CH₄ data over the wide temperature and pressure ranges studied. Regarding the adsorption of CO₂, the utilization of electrical charges was also evaluated; the results showed that, using the parametrization employed, the most effective forcefield to describe CO₂ adsorption is a chargeless UA forcefield.

The CH₄ isosteric heat of adsorption obtained from the experimental data can be described with high accuracy using the UiO-66-UA forcefield. In the case of CO₂, it is observed that the use of electrical charges enhances the prediction of the heat of adsorption, especially in the low-coverage region.