Polarizable Force Fields for CO2 and CH4 Adsorption in M-MOF-74.

The family of M-MOF-74, with M = Co, Cr, Cu, Fe, Mg, Mn, Ni, Ti, V, and Zn, provides opportunities for numerous energy related gas separation applications. The pore structure of M-MOF-74 exhibits a high internal surface area and an exceptionally large adsorption capacity. The chemical environment of the adsorbate molecule in M-MOF-74 can be tuned by exchanging the metal ion incorporated in the structure. To optimize materials for a given separation process, insights into how the choice of the metal ion affects the interaction strength with adsorbate molecules and how to model these interactions are essential. Here, we quantitatively highlight the importance of polarization by comparing the proposed polarizable force field to orbital interaction energies from DFT calculations. Adsorption isotherms and heats of adsorption are computed for CO2, CH4, and their mixtures in M-MOF-74 with all 10 metal ions. The results are compared to experimental data, and to previous simulation results using nonpolarizable force fields derived from quantum mechanics. To the best of our knowledge, the developed polarizable force field is the only one so far trying to cover such a large set of possible metal ions. For the majority of metal ions, our simulations are in good agreement with experiments, demonstrating the effectiveness of our polarizable potential and the transferability of the adopted approach.

Introduction

The society’s demand for energy and how it is currently satisfied interweaves strongly with anthropogenic CO2 emissions and hence to the changing climate.[1−3] It is evident that, to maintain present living standards, the energy sector needs to be altered drastically.[4] New environmentally friendly ways of transforming energy have to be implemented on a large scale.[5] This significant change of the energy sector is, however, still years from being fulfilled.[6] New technologies need to be developed and further improved.[7] Near-term measures include the considerable reduction of CO2 emitted by conventional power plants.[8] To reduce CO2 emissions of current power plants, CO2 needs to be separated from, e.g., the flue gas.[8,9] Besides carbon capture, CO2 removal is also crucial for other technologies, e.g., the purification of natural gas.[10−12] A promising technology for the efficient separation of large quantities of CO2 is the separation via solid adsorbents.[13,14] In this context, metal–organic frameworks (MOFs) have received substantial attention.[6,15−17] MOFs are a relatively new and versatile type of material with various possible application areas such as in gas separation,[18−22] gas storage,[23−26] gas and liquid separation,[7,27−31] catalysis,[32,33] sensing,[34] drug delivery,[35,36] microelectronics,[37,38] and biotechnology.[39−41] MOFs are constructed of metal ions or clusters connected by organic linkers.[42] In recent years, a tremendous number of new MOFs has been synthesized[43] and an almost infinitive number seems to be theoretically possible.[44] By adjusting the combination of metal ions and the organic parts, the properties of MOFs are widely tunable and materials with exceptionally large surface areas can be created.[45] The pore geometry can be customized to enhance the separation of molecules due to the topology.[46,47] In addition to the geometry, the chemical composition can be tuned to further improve separation performance.[46,48,49] For instance, coordinatively unsaturated metal ions, so-called open-metal sites, can be embedded on the surface of the pore structure.[50] These metal ions are accessible for guest molecules and therefore interact strongly with certain adsorbate molecules.[51] Therefore, the uptake of some adsorbate molecules can be increased significantly.[52] Understanding and predicting the interaction of open-metal sites with adsorbates is crucial for the design of new customized adsorbent materials.[53] A challenge that is inherent with the enormous number of possibilities in the synthesis of MOFs is the selection of the best one for a particular application.[47] Experimental screening of hundreds of thousands of MOFs is impractical. A large effort has been made on developing computational screening approaches to facilitate the selection.[13,44,54−56] Today, it is possible to predict adsorption properties for large sets of existing and hypothetical MOFs based on molecular simulations.[57,58] A prerequirement for this kind of computational screening is a force field that represents the molecular interactions reasonably well for all materials under investigation. Unfortunately, the existing generic force fields do not fulfill this prerequirement for all MOFs.[52,59−64] Especially, the promising class of MOFs with open-metal sites has been shown to be poorly described by generic force fields and research has been focused on developing improved force fields for these materials.[62,65−68] To illustrate the failure of generic force fields, Figure compares experimental measurements from Herm et al.[18] to the CO2 uptake in Mg-MOF-74 calculated from grand-canonical Monte Carlo simulations applying the generic UFF force field[69] for Mg-MOF-74 and the TraPPE force field[70] for CO2 (i.e., standard generic force fields commonly used for porous materials[64,71−73]).

Figure 1

Comparison between the experimental adsorption isotherm of CO2 in Mg-MOF-74 from Herm et al.[18] and the simulated one using the UFF force field[69] for Mg-MOF-74 and the TraPPE force field[70] for CO2 at 313 K.

The CO2 uptake of Mg-MOF-74 predicted from molecular simulation for low fugacities with the UFF force field is considerably lower than the experimentally determined one. In particular, this region is highly relevant to carbon capture. The distinct inflection in the adsorption isotherm is also not depicted, which suggests that the strong affinity of CO2 close to the open-metal sites is not modeled correctly with the UFF force field. To obtain accurate force fields, several studies have been conducted in which force fields for individual MOFs are matched to interaction energies computed with quantum mechanical methods.[49,62−64,66,74−76] In some of these studies, the applicability of the customized force fields was also investigated for MOFs with very similar topology and composition.[49,63,64] Borycz et al.[71] used this methodology to investigate the influence of the exchange of the metal ion for a MOF without open-metal sites. Addicoat et al.[77] designed an extension to the UFF force field to capture the structure of MOFs. Moreover, Vanduyfhuys et al.[78] developed a software package called QuickFF to automatically derive force fields for MOFs from ab initio input. The iso-structural M-MOF-74 has been pointed out to be a well suited study case to investigate the influence of different metal ions on the adsorption properties of small molecules[51,64,79,80] and thereby further improve force fields. The pore structure of M-MOF-74 is only slightly influenced by the exchange of the metal ion, whereas the adsorption properties can change considerably. Understanding and describing the underlying interactions of open-metal sites with guest molecules is of fundamental interest[64] and can help to find trends and design even better adsorbent materials.[71] M-MOF-74 is built of one-dimensional hexagonal pores with a diameter around 11 Å[81] and exhibits a particular high density of open-metal sites.[82]Figure shows an extract from the periodic structure of Mg-MOF-74.

Figure 2

Extract from the periodic structure of Mg-MOF-74. Mg, C, O, and H atoms are represented in green, gray, red, and white, respectively.

In recent studies, Mg-MOF-74 has been shown to be a promising candidate for carbon capture due to its high CO2 uptake capacity at low partial pressures[18,19,52,81,83,84] and for natural gas sweetening.[64] Subsequently, M-MOF-74 has been extensively investigated for various gas separations.[18,52,63,75,81,85] Among others, experimental studies include adsorption measurements of CO2,[18,19,46,51,52,81,83,84] CO,[51,86] CH4,[51,79,87,88] C2H6,[51,87−90] C2H4,[87−90] C2H2,[87,88] C3H8,[87−90] C3H6,[87−90] Ar,[51] O2,[20] and N2.[20,51] Adsorption sites have been investigated via neutron and X-ray powder diffraction to determine the binding geometry.[46,79,91−93] As a complement to experiments, various quantum mechanical studies have been conducted to theoretically investigate adsorption sites,[46,91] adsorption energies,[79,80] and the underlying contributions and mechanisms[48,53] for a large number of guest molecules. Moreover, the mechanism of competitive adsorption has been studied by Tan et al.[94] These authors found that kinetic effects can play a significant role in the replacement of adsorbate molecules close to the open-metal sites. Molecular simulations have been used to investigate the adsorption behavior at uptakes larger than one guest molecule per open-metal site[49,62,63,66,74−76] and the hopping of guest molecules between open-metal sites.[95,96] Despite the significant progress, it is still a major challenge to accurately capture the change of interaction strength with varying metal ions in M-MOF-74 in molecular simulations. Several simulation studies have been conducted to reproduce the adsorption behavior of some of the M-MOF-74 structures.[49,62,63,66,74−76] In these studies, standard interaction potentials are reparametrized to reproduce guest–host interactions from quantum mechanical calculations. In some quantum mechanical studies, it is suggested that guest molecules are polarized in the vicinity of the open-metal sites in M-MOF-74[46,48,49,51,53,91] and that this interaction contributes to the enhanced CO2 affinity. Standard force fields do not include this effect directly and therefore separate adjustments of the force field parameters may be necessary for every new structure. In this context, polarizable force fields for porous materials have been suggested,[49,97−100] due to the observation of polarization in several other MOFs.[101−103] However, we underline the potential of polarizable force fields for the study of adsorption in grand-canonical Monte Carlo simulations of M-MOF-74. We anticipate that considering polarization explicitly can help to create force fields that overcome the shortcomings of current generic force fields.

In this manuscript, we evaluate the potential of explicit polarization to improve the issue of limited force field transferability using MOFs with open-metal sites. In particular, we study the adsorption of CO2 and CH4 in M-MOF-74 with M = Co, Cr, Cu, Fe, Mg, Mn, Ni, Ti, V, and Zn. We extend the previously developed polarizable force field for CO2 in Mg-MOF-74[118] to structures based on nine more metal ions and CH4 without additional fitting parameters. Subsequently, we conduct grand-canonical Monte Carlo simulations, and compare our results to results using other force fields[64,69] and experiments.[19,46,64,79,81,104] Thereby, it is shown that polarizable force fields have the potential to improve the transferability of force fields describing porous materials.

Methodology

Force fields describing intermolecular interactions are the foundation of molecular simulations.[59,105] By definition, molecular simulations represent the behavior of a system for a given force field. The capability of the force field to describe the true molecular interactions determines its applicability.[106] An additional desirable characteristic of force fields is transferability.[59,77,78] Ideally, a force field should be able to describe the experimentally observed behavior for a preferably large set of systems. Generic force fields like UFF,[69] DREIDING,[107] and OPLS[108] have been designed for organic, biological, and inorganic materials.[59] However, if the conditions of the system under investigation vary from the ones the force field was developed for, the resemblance of the real system behavior may be poor.[59] Many studies focus on regaging force field parameters to capture experimental behavior.[49,62−64,66,74−76,109] This approach works well for individual systems.[64] Nevertheless, the transferability of the created force fields is likely to be limited to structures with very similar topology and chemical composition.[64] In addition, for each pairwise interaction and all point charges, new parameters are required which are all mutually dependent.[64] Repeatedly readjusting force field parameters for every new system is cumbersome. The reason for the limited transferability of this approach could be attributed to the implicit consideration of the interactions that are exceptional for a particular system. Another disadvantage is that the predictive potential of molecular simulations is largely lost when force fields become completely empirical and need to be readjusted for every new system. A more sustainable approach is to develop force fields with a broader applicability due to a physically motivated extension which considers these exceptional interactions. Several studies focus on force field improvement of gas adsorption in MOFs.[59,77,78] Unfortunately, a force field with general applicability for the adsorption of small molecules in MOFs does not exist.[59] Especially, the modeling of MOFs with open-metal sites represents a challenge.[62,66,67] Quantum mechanical calculations of CO2 adsorption in Mg-MOF-74 suggest that polarization of CO2 in the vicinity of Mg ions is important and significantly contributes to the interaction energy.[46,48,49,51,53,91] In contrast, charge transfer between CO2 molecules and the MOF framework seems to be negligible.[46,53,104] Several methods have been proposed for considering polarization in molecular simulations, i.e., the induced dipole method, the fluctuating charge method, and the shell method (also known as Drude oscillator and charge-on-a-spring model).[110−115] For molecular dynamics simulations, these methods are well established.[100,110,116,117] However, the many-body nature of polarization makes these algorithms more suitable for molecular dynamics simulation in which all molecules are moved in every simulation step. This is in contrast to Monte Carlo simulations in which usually only one molecule is moved.[117] Hence, in Monte Carlo simulations, more steps are required to create independent configurations of the system. Normally, this is unproblematic, since the interactions need to be computed only for the moved molecule. However, when considering polarization, the interactions between all molecules change and have to be recomputed due to the many-body nature of polarization for every step. This leads to a less frequent consideration of polarization in Monte Carlo simulations. As described in our initial study,[118] we use the procedure developed by Lachet et al.[119] to mitigate this limitation. The procedure uses the induced dipole method in which the induction energy Uind is expressed aswhere μ is the induced dipole, E0 is the permanent electric field created by the static partial charges at interaction site i, and N is the total number of interaction sites in the system. The energy contribution of the induction energy has to be computed in every Monte Carlo step. In this way, the difference in induction energy to the previous configuration can simply be added as another energy term in the acceptance rule of the Monte Carlo algorithm. Higher order induced multipoles are not explicitly incorporated in the induced dipole method. In a similar system, Lachet et al.[119] estimated the related error to be less than 5% of the total induction energy. Special for the approach of Lachet et al.[119] is that it accounts solely for polarization between the framework and adsorbate molecules and that it neglects polarization caused by induced dipoles, so-called back-polarization. Using these assumptions, eq can be rearranged towhere α is the atomic polarizability of interaction site i and n is the number of interaction sites of the moved molecules. Thereby, an iterative scheme is avoided and the computational costs of the method are drastically reduced. In fact, the computational costs can be similar to simulations without considering explicit polarization. In the case of, e.g., a translation move of a single molecule only, the n interaction sites of this molecule have to be evaluated to determine the change in the induction energy. Lachet et al.[119] showed that the error in energy introduced by this assumption is around 6% in a xylene NaY zeolite system. To verify the contribution of back-polarization in Mg-MOF-74, in Figure , the total interaction energy of a CO2 molecule approaching the Mg ion with and without consideration of back-polarization for the developed polarizable force field is compared.

Figure 3

Total energy of a single CO2 molecule in Mg-MOF-74 calculated using the developed polarizable force field as a function of the distance to the open-metal site. Comparison between interactions with and without back-polarization.

The influence of back-polarization increases with decreasing distance between the CO2 molecule and the metal ion. For the most favorable position at approximately 2.4 Å, the difference in total energy is approximatively 7%. This deviation seems to be acceptable in comparison with the considerable speedup of the simulations. Besides polarization, repulsion and dispersion interactions are considered via a standard Lennard-Jones potential and static charge distributions are modeled via point charges. When explicitly accounting for polarization, one has to ensure that the force field parameters describing the remaining interactions do not include an implicit polarization contribution which would have to be removed. Otherwise, the contribution of polarization would be double counted, once implicitly and once explicitly. The removal of implicit polarization is necessary if a standard force field is used as the starting point for the development of a polarizable force field, because current force fields are likely to be calibrated to reproduce certain experimentally observed properties. For example, the TraPPE force fields for CO2 and N2 are fitted to reproduce experimental vapor–liquid equilibria of the pure components and their mixtures with alkanes without explicitly considering polarization.[70] Hence, in the fitting of these force fields, all present interactions are indirectly considered and the resulting potential parameters are effective parameters. As the starting point for our polarizable force field, we use the UFF[69] and TraPPE force fields.[70] These are standard force fields frequently used for molecular simulations of porous materials.[64,67,71,120−123] To remove the contribution of implicitly considered polarization to the interaction potential, a global scaling parameter λ is applied to all Lennard-Jones energy parameters developed without explicit polarization. A simple procedure is chosen to verify the applicability of polarizable force fields rather than attempt to perfectly reproduce experimental results. Here, we reduce the Lennard-Jones energy parameters ε taken from the UFF and TraPPE force fields with respect to their atomic polarizabilities viawhere α and αmax are the atomic polarizabilies of interaction site i and the largest atomic polarizability, respectively. The scaling parameter λ can vary between 0 and 1. Thereby, it is assured that nonpolarizable interaction sites (α = 0) are unchanged and that the potential energy parameters of the atoms with the largest polarizability are reduced the most. A more detailed description of the derivation and the algorithm of the presented method can be found in our previous publication.[118] The required atomic polarizabilities α are taken from the literature.[124,125] Many different values for atomic polarizabilities can be found for every atom.[119,124−128] Their values can differ significantly depending on the experimental procedure or the theoretical assumptions made.[129] Hence, a global scaling factor ζ is used to adjust the magnitude of the atomic polarizabilities taken from the literature αlit with respect to the chosen interaction potential according toThereby, the ratio between the individual atomic polarizabilities is not affected to ensure a reasonable relative contribution of polarization between the atoms. This kind of scaling procedure for atomic polarizabilities is frequently used in the literature,[130,131] and the scaled polarizabilities adopted in this study have comparable magnitudes to previous molecular simulation studies.[99,119,132] In this manuscript, the values of ζ and λ are adjusted to reproduce the experimental adsorption isotherm for CO2 in Mg-MOF-74. In a first step, the low fugacity region of the simulated adsorption isotherm and the heat of adsorption for CO2 in Mg-MOF-74 are tuned by scaling all atomic polarizabilities with ζ. In the low fugacity region, CO2 molecules adsorb close to the open-metal sites where polarization interactions are of particular importance. Subsequently, the scaling parameter λ is adapted to remove the implicit contribution of polarization from the Lennard-Jones potential. Therefore, the value of λ is lowered to match the high fugacity region of the experimentally determined adsorption isotherm. In this region, the centers of the channels of Mg-MOF-74 are filled with CO2 molecules. The locations in the centers of the channels are further away from the open-metal sites, and therefore, polarization is less important. By applying this two-step procedure, we divide the interaction energy into the underlying physical contributions without using an elaborated approach. For the remaining M-MOF-74 structures and for CH4, the scaling factors determined for the Mg structure with CO2 are used. Thereby, the transferability of the approach is investigated. The procedure is chosen to verify if the polarizable force field has the potential to describe the difference between the different metal ions embedded in M-MOF-74.

Simulation Details

Grand-canonical Monte Carlo simulations implemented in the RASPA software package[133,134] are conducted to compute the uptake and heats of adsorption of CH4 and CO2 in the different structures of the M-MOF-74 (M = Co, Cr, Cu, Fe, Mg, Mn, Ni, Ti, V, and Zn) family. The uptakes are computed for varying fugacities, for pure components and mixtures at 298 K. DFT-optimized, all atomic MOF structures with atomic charges are taken from Lee et al.[80] In the simulations, the structures are considered to be rigid. Lennard-Jones parameters for CH4 and CO2 are taken from the TraPPE force field.[70] Interactions between guest molecules are not modified and computed according to the TraPPE force field. The UFF force field[69] is used for the atoms of M-MOF-74. Cross-interactions are calculated via the Lorentz–Berthelot mixing rule from atomic parameters.[135] All Lennard-Jones energy parameters εscale used in the simulations are adjusted according to eq with λ = 0.7. Thereby, we account for previously implicitly considered polarization. The Lennard-Jones potential is truncated at a cutoff distance of 12.8 Å without tail corrections. To mimic the behavior of the continuous system, i.e., a repetition of identical one-dimensional pores, periodic boundary conditions are applied in all directions (see Figure ). The simulated system is composed of multiple unit cells to ensure a minimum distance of more than twice the cutoff radius between periodic images. The Ewald summation technique with a relative precision of 10–6 is used to calculate electrostatic interactions between static point charges.[136] Explicit polarization is considered via the induced dipole method.[117] Polarization is exclusively considered between the framework and adsorbate molecules. Additionally, back-polarization is neglected to achieve reasonable simulation times. The required atomic polarizabilties α are taken from Shannon[124] and van Duijnen and Swart[125] and are scaled with a ζ value of 0.09. All force field parameters are summarized in Tables S1–S11 (Supporting Information). For the comparison with experimental results and simulation results of others, we use the Peng–Robinson equation of state to convert pressures to fugacities.[137] The DFT calculations to determine the orbital interaction energy (as explained below) are performed with the Amsterdam Density Functional (ADF) package.[138,139] The B3LYP-D3 exchange-correlation functional[140−144] is used with a TZP-STO basis set and a large frozen core. A fragment analysis is performed between CO2 and Mg-MOF-74 to assess the net interaction between these two fragments. Using the energy decomposition analysis scheme by Ziegler and Rauk,[145−147] the interaction energy Δint between the two fragments is decomposed intoΔVelstat comprises classical electrostatic interactions between unperturbed charge distributions of the deformed fragments. ΔPauli describes the Pauli repulsion energy and corresponds to the destabilizing interactions between occupied orbitals. The Pauli repulsion energy is responsible for steric repulsion. Δoi represents the orbital interaction energy which accounts for charge transfer and polarization.[138]

Results and Discussion

As an initial step in the evaluation of the developed polarizable force field, it is important to investigate the role of polarization in the adsorption of guest molecules in M-MOF-74. In Figure , we compare the polarization energy of a CO2 molecule approaching the Mg ion of Mg-MOF-74 estimated with the developed polarizable force field to the orbital interaction energy calculated from ADF.

Figure 4

Comparison of the polarization energy computed with the developed polarizable force field without considering back-polarization and the orbital interaction energy from DFT calculations as a function of the distance between a CO2 molecule and the Mg ion.

The orbital interaction energy should be a good approximation for the polarization energy, since no reaction is taking place and considerable charge transfer is not expected for very similar configurations of CO2 inside Mg-MOF-74.[46,53,104] For relevant distances, both methods show a comparable trend for the energy contributions. The most relevant distance between the CO2 molecule and the Mg ion is where the total energy is the lowest (i.e., 2.3–2.5 Å, as shown in Figure ). At this distance, the polarizable force field predicts that the polarization energy of a single CO2 molecule in Mg-MOF-74 has a significant contribution of around 30% to the total energy. For larger distances, the contribution of polarization decreases rapidly. In our previous study,[118] we investigated the quality of the developed polarizable force field by comparing the total energy of random CO2 positions inside Mg-MOF-74 with detailed DFT calculations from Lin et al.[75] Thereby, also less favorable positions further away from the Mg ions are probed. These positions are occupied after all open-metal sites are saturated. In general, the polarizable force field describes most positions considerably better than the UFF force field and with a quality comparable to a nonpolarizable force field that has been developed by readjusting the majority of force field parameters.[64] The scheme applied here is considerably simpler.[118] The resulting adsorption isotherms for CO2 in Mg-MOF-74 in comparison to experimental measurements, the UFF force field, and the DFT-derived nonpolarizable force field of Mercado et al.[64] are shown in Figure a.

Figure 5

Comparison between the experimental results of Herm et al.[18] (open), Queen et al.[46] (yellow), Yu et al.[104] (orange), and Dietzel et al.[81] (brown) and simulation results using the developed polarizable force field (black), the UFF force field (blue) and the DFT-derived nonpolarizable force field of Mercado et al.[64] (green) for CO2 in Mg-MOF-74. (a) Adsorption isotherm at 298 K (Herm et al.,[18] 313 K); (b) heat of adsorption as a function of uptake.

The simulation results with the polarizable force field clearly display the inflection of the experimental adsorption isotherm. The predicted behavior is significantly better than that with the UFF force field. This is expected, because the scaling factors are adjusted to reproduce the experimental data. The overall agreement with the experimental measurements is comparable with the DFT-derived nonpolarizable force field of Mercado et al.[64] Both force fields can predict the low fugacity region which is particularly important for carbon capture. For higher fugacities, simulations with all compared force fields predict higher CO2 uptakes in comparison to the experiments. As pointed out earlier, this can be attributed to the fact that a certain degree of inaccessibility due to diffusion limitation or defects in the crystal structure is inherent with experimental structures.[49,81] In the limit of very high CO2 uptakes, the guest–host interactions become less important and the adsorption is dominated by the accessible volume for CO2.[71] In the development of the used TraPPE force field, the CO2–CO2 interactions were adjusted to reproduce the vapor–liquid equilibria and it describes the density per void volume well. Therefore, the uptake of CO2 predicted using the polarizable and UFF force field converges in the high fugacity region. It should be noted that Mercado et al.[64] scaled the calculated CO2 uptakes with 0.85 to account for inaccessibility of open-metal sites. This scaling procedure mainly improves the agreement between experiments and computations for the high fugacity region. Figure b shows the heat of adsorption as a function of CO2 molecules per metal ion. The distinct inflection of the adsorption isotherm caused by the strong affinity of the CO2 molecule toward the metal ions is reflected by the change of the heat of adsorption with increasing gas uptake. The calculated heat of adsorption has an inflection at exactly one CO2 per metal ion. Before and after the rapid decrease at one CO2 molecule per metal ion, the heat of adsorption increases slightly. This increase can be related to a rise in the total number of adsorbed CO2 and therefore a larger contribution of the CO2–CO2 interactions to the total energy. Similarly, the experimental heats of adsorption increase initially. In general, the experimental heats of adsorption have to be regarded with wariness. The heat of adsorption is not measured directly but calculated according to −qst/R = ∂(ln p)/∂(1/T) at constant loading and averaged over adsorption isotherms at different temperatures.[104] The temperatures considered vary for all experimental studies. The experimental curves consistently show a drop in the heat of adsorption for lower ratios of the number of guest molecules and metal ions than the simulation results. Different sets of experimental adsorption isotherms show inflections at different uptakes of CO2. This is another indication for defects in the crystal structure and the blocking of some of the metal ions of the experimental structures. Haldoupis et al.[49] further investigated the effect of blocking for the Co, Cu, Mn, and Ni based structures. These authors illustrate that varying levels of pore accessibilities can explain the discrepancy between different experimental studies. Especially, the unavailability of open-metal sites can explain the drop in the heat of adsorption prior to the complete saturation of these sites. This can be caused by residual solvent molecules binding to the open-metal sites. According to Haldoupis et al.,[49] these residual solvent molecules could reduce the number of accessible open-metal sites by 20–30%, while only slightly affecting the accessible surface area and the accessible volume. Previously, the good agreement between the experimental BET surface area and pore volume and the theoretical void space in the empty crystal structure made Dietzel et al.[81] suggest that their MOF was fully activated. In contrast, Haldoupis et al.[49] concluded that a combined effect of crystalline defects and residual solvent molecules is most likely to cause the difference between simulations and experiments. The focus of this study is to evaluate the applicability of a polarizability force field for describing the interactions of guest molecules with different metal ions. In this regard, the consideration of residual solvent molecules and defects does not seem to be crucial. For the remaining M-MOF-74 structures, the parameters adjusted for Mg-MOF-74 and CO2 are used. Thus, the calculated values are predictions based on the two global scaling parameters λ and ζ adjusted for Mg-MOF-74. To obtain an overview of the results, we divided the predictions for CO2 in the M-MOF-74 structures into three groups. The first group consists of Co, Cu, Ni, and Zn, the second of Cr, Ti, and V, and the third of Fe and Mn. In Figure , the computational results for the first group are presented and compared to experimental measurements of Queen et al.[46] and the UFF force field.

Figure 6

Comparison between the experimental results of Queen et al.[46] (yellow) and simulation results using the developed polarizable force field (black) and the UFF force field (blue) for CO2 in the Co (■), Cu (●), Ni (▶), and Zn (⧫) based structures. (a and c) Adsorption isotherms at 298 K; (b and d) heats of adsorption as a function of uptake.

For these structures, the developed polarizable force field is able to describe the experimental measurements well. This is most striking in comparison to the UFF force field which is not able to model the differences between the metal ions. For the polarizable force field, the largest deviations in the adsorption isotherms can be observed for Cu and Co based structures. In agreement with the experimental data, the simulation results for Co and Ni based MOFs show a less distinct inflection for the CO2 adsorption isotherm than for Mg-MOF-74. The experimental adsorption isotherms for Cu and Zn based structures do not show an inflection for CO2, which is also accurately predicted in the simulations applying the polarizable force field. The comparison to the experimental results of Yu et al.[104] and the simulation results of Mercado et al.[64] can be found in Figures S22–S11 (Supporting Information). The results of Mercado et al.[64] match the experimental adsorption isotherms well. It is worth mentioning that Mercado et al.[64] did not perform simulations for Cu-MOF-74. These authors did not develop a force field for Cu-MOF-74, because of an elongation of the unit cell in the c-direction[80] in comparison to the other M-MOF-74 structures. Although our results for the Cr based structure deviate from the experimental results (compare Figure a), the elongation does not seem to be problematic for the general applicability of our approach. The calculated heats of adsorption shown in Figure b and d have a similar quality as that for the Mg based structure. The largest discrepancy between simulations and experiments can be observed for Zn-MOF-74. This is surprising, because the calculated adsorption isotherm agrees very well with experiments and is very similar to the calculated adsorption isotherm for Co-MOF-74 with a similar heat of adsorption. Similar to Mg-MOF-74, the heat of adsorption derived from experiments shows an inflection significantly before an uptake of one CO2 molecule per metal ion. As mentioned previously, residual solvent molecules are likely to cause this shift in the heat of adsorption,[49] since less open-metal sites are accessible. The simulations predict a similar behavior for all structures after all open-metal sites are saturated with CO2. In this region, the CO2 molecules start to accumulate in the centers of the channels. The geometry of the channels is almost identical for all types of M-MOF-74, and the CO2 molecules are sufficiently far away from the metal ions to not be significantly affected by polarization. Overall, the polarizable force field seems to have the potential to capture the different degrees of polarizations related to the different metal ions for these four structures. Moreover, in contrast to the UFF force field, the polarizable force field is able to predict the correct order of adsorption strength for the Co, Cu, Ni, and Zn based structures. The computationally predicted adsorption isotherms and heats of adsorption for the second group are compared to the computational results with the UFF force field in Figure .

Figure 7

Comparison between the simulation results using the developed polarizable force field (black) and the UFF force field (blue) for CO2 in the Ti (■), V (●), and Cr (▶) based structures. (a) Adsorption isotherms at 298 K; (b) heats of adsorption as a function of uptake.

For these structures, no experimental adsorption measurements are available. To the best of our knowledge, these structures belong to the group for which the experimental syntheses are still challenging.[148] The simulations predict a very distinct inflection for the adsorption isotherms of Ti- and V-MOF-74, while the one for Cr-MOF-74 does not show an inflection. The predictions for the Ti and V based structures agree with theoretical predictions of Park et al.[42] These authors expect the structures to have even stronger interactions with CO2 than Mg-MOF-74 which is supported by our simulations. The trend of the adsorption isotherms for the three MOFs is reflected in the heats of adsorption. The possibility to predict large differences between adsorption behavior shows further that polarizable force fields have the potential to describe such significant differences in the adsorption behavior. Again, the UFF force field predicts a totally different adsorption behavior and a smaller difference between the metal ions (compare Figure ). The remaining M-MOF-74 structures are based on Fe and Mn. A notably large discrepancy between simulations and experiments is found for these structures, as shown in Figure .

Figure 8

Comparison between the experimental results of Queen et al.[46] (yellow) and simulation results using the developed polarizable force field (black) and the UFF force field (blue) for CO2 in the Fe (■) and Mn (●) based structures. (a) Adsorption isotherms at 298 K; (b) heats of adsorption as a function of uptake.

The experimental results for the Mn and Fe based structures are very similar. Both structures show weaker interactions between the metal ions and the CO2 molecules than for the Mg based structure. As can be seen, the developed polarizable force field significantly overestimates these interactions. The UFF force field is able to capture the adsorption behavior better. Several reasons for the overestimation are possible, and we feel a combination of different effects is most likely. Interestingly, Mercado et al.[64] also failed to obtain a reasonable force field for Mn-MOF-74 based on fitting the interaction potential to quantum mechanical energies. This further suggests that the explanation for the failure could be rather complicated. For example, ferromagnetic effects which are not considered in the polarizable force field may play a more important role for Mn and Fe than for the other structures. Additionally, the initial force field parameters taken from the UFF force field for the two metal atoms could be of particularly bad quality. A comparison between the values for Fe from the DREIDING (ε/kB = 27.677 K, σ = 4.045 Å) and the UFF force field (ε/kB = 6.542 K, σ = 2.59 Å) shows the huge difference. Simulations based on the DREIDING parameters result in a totally different prediction of the adsorption behavior. This is illustrated by the Henry coefficients of CO2 in Fe-MOF-74 in the limit of infinite dilution condition[149] we computed for both sets of force field parameters using Widom test particle insertions[149] (Supporting Information). The sensitivity of the system might also play an important role. In addition, the quality of the selected level of theory for the structure optimization could be better for some of the metal ions than for others.[64] The failure of the polarizable force field to predict the behavior of these two structures needs to be further investigated. Nevertheless, it does not diminish the potential of polarizable force fields for the description of MOFs with open-metal sites. To further verify the applicability of polarizable force fields, grand-canonical Monte Carlo simulations are performed for CH4 in the M-MOF-74 series. The separation of CO2 and CH4 is industrially relevant.[12] Additionally, CH4 is explicitly chosen to examine the suitability of the polarizable force field to capture the varying influence of different metal ions in M-MOF-74. Both CO2 and CH4 have a similar polarizability[125] but show a totally different adsorption behavior in the series of M-MOF-74. Previous studies explain the difference with a combination of electrostatic interactions caused by the permanent quadrupole of CO2 and polarization.[51] As an example, in Figure , the predicted adsorption isotherms for CH4 in the Mg, Co, Ni, and Zn structures are compared to experimental measurements, simulations of Mercado et al.,[64] and the UFF force field.

Figure 9

Comparison between the experimental results of Wu et al.[79] (violet), Mason et al.[19] (gray), Dietzel et al.[81] (brown), and Mercado et al.[64] (cyan) and simulation results using the developed polarizable force field (black), the UFF force field (blue), and the DFT-derived nonpolarizable force field of Mercado et al.[64] (green) for CH4. (a) Mg-MOF-74, (b) Co-MOF-74, (c) Ni-MOF-74, (d) Zn-MOF-74.

The data for the remaining structures and the heats of adsorption are provided in Figures S12–S21 (Supporting Information). For low fugacities, simulations with the Mn based structure show an unphysical behavior (compare Figure S17). The DFT-optimized structure does not seem to be perfectly symmetrical. Four adsorption sites close to Mn ions are much stronger for CH4 than the remaining adsorption sites. These adsorption sites are occupied with CH4 for all fugacities. To investigate this unusual behavior, we conducted Widom test particle insertions to compute Henry coefficients of CH4 in Mn-MOF-74 in the limit of infinite dilution condition[149] with and without blocking of these four adsorption sites. In addition, simulations were performed with Lennard-Jones parameters for Mn from the UFF force field and with parameters for Zn from the DREIDING force field. The results are presented in the Supporting Information and show that the four adsorption sites are responsible for the behavior. For the DREIDING parameters, the problem does not occur. The locations of the CH4 molecules are different due to the different Lennard-Jones parameters and hence interactions are much smaller. Overall, the UFF force field performs better than the developed polarizable force field for CH4. The polarizable force field underpredicts the uptake for CH4. This is due to the very simple and crude procedure to determine the force field parameters. Actually, not a single parameter was adjusted for CH4. CO2 and CH4 have a similar polarizability,[125] but CO2 is modeled with three interaction sides and CH4 with only one interaction side. Hence, in our force field, CH4 has a larger assigned polarizability than CO2 and the Lennard-Jones energy parameter for CH4 is more reduced. Besides, CH4 is modeled without point charges. It is very reassuring that the adsorption isotherms of CH4 computed with the polarizable force field do not show a distinct inflection which is in agreement with experimental results. As expected by Mishra et al.[51] and observed in this study, the strong interactions in the case of CO2 are caused by a superposition of static polarity and polarization. Hence, the adsorption of the uncharged CH4 molecule is far less affected by the different metal ions and a similar behavior is observed for all M-MOF-74 structures. This is correctly captured by the polarizable force field. In addition, the large CH4 molecules have a larger distance to the metal ions and hence almost no dipole is induced. It would be straightforward to improve the performance of the polarizable force field by introducing a molecule dependent scaling parameter to account for the unequal number of interaction sides. The simulation results of Mercado et al.[64] are better than the ones with the UFF force field and the polarizable force field. The experimental adsorption isotherms are well reproduced. However, the same procedure as that for CO2 was performed to fit the force field parameters for all metal ions separately and the CH4 uptake is scaled with a factor of 0.85 to account for inaccessible open-metal sites and blocked pores. In comparison to experimental measurements, the UFF, and the polarizable force fields, the curvature of the computed adsorption isotherms for CH4 of Mercado et al.[64] seems to be systematically different. The slope of the adsorption isotherm decreases at much lower fugacities. The reason for this behavior should be investigated to further improve the approach. Finally, it is possible to make predictions for the adsorption of mixtures of CH4 and CO2 with the developed polarizable force field. Measuring the gas uptake of mixtures in MOFs is more complicated experimentally and is not often done. Therefore, computational prediction of mixtures is very useful to predict the capability of a material to separate gases. In Figure , the predictions of an equimolar mixture of CH4 and CO2 are shown for the polarizable force field and compared to predictions based on the UFF force field for Mg-MOF-74 and Zn-MOF-74.

Figure 10

Comparison between the simulation results using the developed polarizable force field (black) and the UFF force field (blue) for an equimolar mixture of CO2 and CH4 in the Mg (up) and Zn (down) based structures. (a and c) Adsorption isotherms at 298 K and (b and d) heats of adsorption as a function of uptake in Mg-MOF-74 and Zn-MOF-74, respectively.

For these structures, the single component adsorption isotherms are reproduced reasonably well and therefore the mixture adsorption isotherms are expected to be predicted reasonably. Predictions of mixtures being adsorbed in the other M-MOF-74 structures can be found in Figures S22–S31 (Supporting Information). Close to the open-metal sites, the interactions between individual guest molecules with the framework are predominant. Inside the channels, the interactions between CO2 and CH4 molecules should be dominant which are reproduced well by the TraPPE force field. Due to the strong interactions for CO2 with the framework, a significantly lower uptake of CH4 is predicted for Mg-MOF-74 with the polarizable force field. Especially, for low fugacities, mainly CO2 is adsorbed. Only after most of the Mg sites are saturated CH4 is taken up. In contrast, the UFF force field predicts an uptake of CH4 even for very low fugacities. The depicted heats of adsorption show the released heat of adsorption for the mixture as a function of CO2 uptake per metal ion. For Mg-MOF-74, the initial adsorption behavior is very similar to the one for pure CO2, because initially almost exclusively CO2 is adsorbed. In contrast, in Zn-MOF-74, CH4 is already adsorbed for low fugacities. The trend of the heat of adsorption for the mixture deviates from the ones for the pure components. The Zn based structure seems to be less suitable to separate CO2 from CH4. The UFF force field predicts an initially larger uptake of CH4 in comparison to the polarizable force field. This can partially be attributed to smaller interactions between CH4 and the framework for Zn-MOF-74 with the polarizable force field.

Conclusions

The simulations using the developed polarizable force field agree reasonably well with experimental measurements for most of the investigated structures of the M-MOF-74 family. The quality of the predictions for CO2 is significantly better than with the UFF force field and for most cases comparable to structure specific force fields developed with more elaborated schemes. The polarization energy computed with the polarizable force field shows a behavior similar to the orbital interaction energy determined from DFT calculations. In principle, these energy contributions should be similar if no reaction and no charge transfer takes place. The conducted procedure of first scaling atomic polarizabilities and subsequently adjusting the Lennard-Jones interaction parameters including implicit polarization is simple and requires relatively little effort. The two global scaling factors used here are exclusively tuned for the Mg based structure and CO2. Hence, the results for the other structures and CH4 are predictions. For CH4, no inflection is observed although the largest polarizability is assigned to CH4. The predictions for CH4 adsorption could be significantly improved by a molecule specific adjustment of the polarizability. The concept of only considering explicit polarization between guest molecules and the framework and neglecting back-polarization seems to be a well suitable approach to study adsorption phenomena in porous materials. The assumptions considerably enhance the computational performance of Monte Carlo simulations while using polarizable force fields. Actually, the computational time can be similar to Monte Carlo simulations without a polarizable force field. This is an important assessment, because Monte Carlo simulations are the method of choice for the prediction of adsorption properties in porous materials. Future work will focus on further developing polarizable force fields and deriving a consistent set of parameters from quantum mechanical calculations to avoid fitting to experimental data. We believe that this can lead to force fields with better physical justification and improved transferability. This is crucial for the usage of Monte Carlo simulations for material screening and to make meaningful predictions. Polarizable force fields for Monte Carlo simulations are also promising for other systems with a significant polarization contribution,[116,150] i.e., water,[117] systems including ions,[116,132] or xylenes.[119]