Host–guest interactions in framework materials: Insight from modeling Open Access

At a first glance, it may appear that water filtration with activated carbon, sound absorption with metal foams, and water softening using zeolites in laundry detergents have very little in common. However, all these processes rely on materials possessing tiny pores—microporous materials—with diameters in the nanometer range, which is less than the wavelength of visible light.

Naturally occurring porous materials have been used for a long time. For example, charcoal is a naturally porous substance, known to and used by the ancient Egyptians, Greeks, and Romans.¹ In the Naturalis Historia, Pliny the Elder describes the utility of charcoal in the treatment of “carbuncles.” (A carbuncle is a cluster of boils that form a connected area of infection under the skin.) Nowadays, charcoal is still used in various beauty products. In 1854, John Stenhouse invented face masks employing charcoal to filter the air, which proved to be useful in hospitals.² Since those early discoveries, numerous porous materials have been developed and manufactured artificially to address a range of laboratory and industrial needs.

Owing to their structural and compositional diversity and tunability of physical and chemical properties, metal–organic frameworks (MOFs) and covalent organic frameworks (COFs) are particularly promising candidates for various applications requiring microporous materials. MOFs are constructed from polynuclear metal-containing building units (called secondary building units, or SBUs) joined by organic linkers into various network topologies. Since demonstrating permanent porosity in 1994 by Yaghi et al.,³ MOFs have continuously attracted the attention of researchers worldwide. COFs are promising organic analogues of MOFs, first described in 2005⁴ and consisting of elements typical for organic molecules. The ability of MOFs and COFs to capture, store, and release guest molecules is, arguably, their most prominent feature, exploited in gas storage,^5–7 gas separation,^8,9 energy storage,¹⁰ gas adsorption,^11–13 catalysis,^14–17 sensing,^18–20 and biomedical applications.^21,22 The nature, behavior, and quantity of absorption depend on pore size and volume, surface area, chemical environment, and, in particular, on the host–guest intermolecular interactions.

The ultimate goal of all efforts in materials science is to make better materials, either through the discovery of new or the improvement and repurposing of existing systems. In silico modeling facilitates these efforts by elucidating the fundamental mechanisms of action and establishing the structure–property relationships. An extensive review by Mancuso et al.²³ discusses general electronic structure approaches to modeling metal–organic frameworks, while Gagliardi and co-workers discuss modeling techniques of MOFs for catalytic applications with particular focus on natural gas conversion.²⁴ In 2015, Lee et al. published an insightful overview article on studying the small-molecule interactions in MOFs by combined theoretical and experimental approaches.²⁵ More recently, Qiao and co-workers traced the transition from single simulations to high-throughput screening (HTS) and machine learning of MOF adsorbents.²⁶ In this Review, we focus more narrowly on the computational studies of the host–guest complexes between framework materials and their molecular targets, exemplifying common simulation techniques and discussing their advantages and shortcomings. The diversity of framework materials and their applications requires simulations at various time and lengths scales. If the objective is an overview of an entire system containing thousands of atoms, atomistic methods based on parametrized force fields are the preferred, albeit less precise choice. Electronic structure methods not only offer more accurate but also more costly means to identify specific adsorption sites and quantify the host–guest interactions (HGIs). These approaches, either used independently or combined in a multiscale manner, are general and can be applied to study host–guest interactions in other crystalline and porous materials, such as zeolites. Finally, up-and-coming machine learning techniques coupled with high-throughput data generation allow rapid screening of thousands of existing and hypothetical frameworks.

The playing field that is comprised by the pore space within framework materials is a battleground for guests to follow trajectories related to adsorption/desorption and diffusion. These processes, in unison, are responsible for the observed performance in guest separation and storage. For example, a membrane's separation performance is dictated by permeability, which is the product of diffusivity and concentration, in turn related to diffusion and adsorption, respectively.²⁷ The processes of diffusion and adsorption appear intertwined but they relate to distinct temporal domains: diffusion is a kinetic property, and adsorption is an equilibrium property. As a result, distinct computational strategies are often employed to examine each of these processes.^28–33 

Throughout the literature, diffusivity is usually treated with molecular dynamics (MD) simulations that sample the trajectories of guests by solving equations of motion for the particles within the system, including the forces and potential energies between particles.³⁴ We note that diffusion can also be estimated by Monte Carlo (MC) simulations and transition state theory.³⁵ Simulating the potential energy and forces for this application is achieved using classical approaches to treating guest–guest and guest–framework interactions that accelerate the simulations permitting the simulation over long time scales (>1 ns), required to accurately describe statistics of diffusion processes. Molecular diffusion in fully flexible models, where both the guest and the pore walls move, is preferred but only if a suitable force field is available. A force field is a mathematical model that describes the interactions between atoms or molecules, and it is crucial for the accuracy of MD simulations. In recent years, a number of force fields have been developed specifically for MOFs. One example of a force field that has been shown to be effective for MOFs is the ZIF-8 force field developed by Zhang et al.³⁶ Similarly, Wu et al. used a bespoke force field to simulate CO₂ transport diffusion in UiO-66 and reported a value in good agreement with experimental data.³⁷ It is important to note that these bespoke force fields are usually not applicable to other materials. Force fields are constantly improved by the community and can be derived from ab initio data.³⁸ In order to overcome this limitation, several force fields have been introduced in recent years that are more universal in their applicability. One such example is the universal force field for MOFs (UFF4MOF), which is based on a universal force field for small molecules.³⁹ Addicoat et al.⁴⁰ and Coupry et al.⁴¹ supplemented the original universal force field by adding new atom types for metal elements to better describe the coordinate bonds in MOFs. With the recent development of more universal force fields for MOFs, it is becoming possible to study diffusion in a wider range of materials using MD simulations.

Adsorption at very low pressures, for example in the Henry's law regime, can be sampled by completely random sampling to simply capture the guest–framework interactions. This random sampling is employed by Monte Carlo (MC)—a method that focuses on static properties of systems and does not involve the evolution of systems over time. In contrast to molecular dynamics (MD), where each state of the system depends on the previous state, in MC, there is no connection between the snapshots (states) of the system. Most MC algorithms are based on modifying the current snapshot by performing changes called “moves” and using an “acceptance rule” to determine whether the new state should be accepted or rejected. However, MC is rarely applied to chemically complex molecules and frameworks and has a reputation for being difficult to debug and test. The RASPA software is a general purpose classical simulation package for MC and is used throughout the literature.⁴⁶ MC can uncover the interaction strength for the guest in the framework. However, more complex MC moves correctly sampling the thermodynamics, such as grand canonical Monte Carlo (GCMC) simulations, are necessary to compute the loading of guest at a specific pressure, in order to simulate an isotherm.⁴⁷ Configurational bias GCMC can be employed to simulate the flexibility of guests.⁴⁸ In general, Monte Carlo simulations involve the sampling of millions of moves that include trial displacement, insertion, and removal of guest molecules in the framework structure. The number of guests adsorbed in framework will evolve over the course of GCMC simulation, and after a period of equilibration, it will fluctuate around an equilibrium value.^46,49 Snapshots from GCMC simulations are useful to demonstrate how guests may pack in the porous space [Fig. 1(b)]. This equilibrium imitates the experimental system, in which the adsorbed phase is at equilibrium with a gas reservoir. These approaches allow for the system configuration to sample the equilibrium case that cannot be sampled efficiently through molecular dynamics simulations.

For even more rigorous quantitative insights, traditional molecular dynamics simulations can be combined with ab initio modeling, e.g., at the density functional theory (DFT) level, to simultaneously sample configurational space, describe the electronic structure of the host–guest complexes, as well as take into account entropic and enthalpic effects.^50–53 This approach is described as ab initio molecular dynamics (AIMD) and is the most accurate method to study the dynamic behavior of guest molecules in porous materials. With this approach, one can establish a relation between the electronics of the host–guest system and its macroscopic properties. It is also well-suited to treat solid–liquid interfaces. However, AIMD is extremely computationally intensive, and this often results in limited sampling. Consequently, AIMD has not been widely applied to guest molecules in MOFs and COFs, apart from a few notable exceptions where it was employed to (i) model the uptake and release of gases, determine its temperature dependence and its effect on structural and electronic properties of MOFs;^54–57 (ii) model breathing and structural flexibility upon adsorption;^58–62 (iii) study the degradation and formation of MOFs in water;^63–67 (iv) investigate proton conduction of guests in MIL-53;⁶⁸ and (v) study ferroelectricity and how it is influenced by trapped water molecules.⁶⁹ 

The nature and strength of the interactions, established upon the guest's diffusion and adsorption in the host's pore, often determine—and can, therefore, be exploited to tune—the selectivity and efficacy of the framework material used to capture, separate, catalytically activate, and/or transport small molecules. These HGIs occasionally involve covalent bonding, but are more commonly non-covalent (intermolecular), spanning from very strong interactions (e.g., hydrogen or dative bonds to metal atoms in MOFs) to relatively weak interactions (e.g., van der Waals forces as is the case of inert gases). We note that dispersion interactions are the dominate energy contribution to most HGIs at the short range.

Accurate in silico characterization of the host–guest complexes generally relies on electronic structure theory in a “so-called” static setting, in which single-point computations are performed on either an equilibrium structure (partially or fully relaxed) or on an experimental structure. The choice of the simulation technique is often guided by the requirements of a specific application: the relative adsorption strengths of different guests are critical for separation processes, the exact location and structure of a guest are important for a catalytic mechanism, and the host–guest interaction type is central to sensing. Practical applications of framework materials vary depending on the size of the molecular guest relative to the pore space. Small guests are molecules made up of only a few atoms and occupying a relatively small portion of the pore volume. Typical examples include gaseous molecules, such as H₂, N₂, O₂, CO₂, H₂O, and light hydrocarbons, e.g., methane, ethane, ethylene, etc. For these small guests, framework materials, owing to their controllable and adjustable porosity, structure, and surface functionalities, are utilized primarily to separate mixtures (such as ethane and ethylene, methane and CO₂, O₂, and N₂, etc.), as well as for capture (particularly CO₂) and storage (e.g., H₂).^6,8 Larger guests, such as drug molecules, are typically dissolved in solvents and loaded together with the solvent. Their applications often involve MOFs for catalysis, drug delivery, and sensing.^9,24,70 Since these guest can form more and stronger interactions with the framework, a detailed analysis of these interactions becomes critical.

Two types of structural models are usually considered: periodic and finite (cluster) models. Periodic computations account for the crystalline (periodic in three dimensions) nature of the material, but are computationally demanding, and therefore, often restricted to density functional theory. In contrast, finite models (clusters), cut out from the periodic structure and capped with (most commonly) hydrogen atoms, can be treated not only at the DFT but also at correlated wavefunction theory levels, such as coupled cluster and perturbation theory. Choosing a model is a balancing act between the precision of the computational method and that of the structural model.

Periodic approaches often use density functional theory methods from the lower rungs on the Jacob's ladder, such as the local density approximation (LDA) and generalized gradient approximation (GGA). These methods by construction lack the description of dispersion, crucial in the host–guest interactions.⁷¹ To remedy this shortcoming and achieve agreement with the experimental data, either the empirical dispersion corrections, such as Grimme's D2 and D3 featuring the R⁶ term⁷² are added to the “bare” DFT results, or the so-called van der Waals (vdW) DFT functionals,⁷³ which self-consistently account for dispersion through non-local correlation, are employed. Siegel and co-workers benchmarked these methods against experimental enthalpies of CO₂ adsorption in four prototypical metal–organic frameworks: M/DOBDC (M = Mg, Ni, and Co) and 5 Cu-HKUST-1.⁷⁴ LDA methods alone were found to significantly overbind, and GGA—to underbind carbon dioxide; addition of D2 dispersion correction notably improved these results both qualitatively and quantitatively, albeit still underestimating the CO₂ adsorption enthalpies by ca. 7 kJ mol⁻¹ compared to experiment. Nonempirical vdW density functionals, particularly revPBE, afforded the best agreement (within ∼2 kJ mol⁻¹) with the reference experimental data (Fig. 2). The same authors later expanded the range of the studied metal atoms in the two MOFs, i.e., M-DOBDC with M = Mg, Ca, Sr, Sc, Ti, V, Mo, and W, and M-HKUST-1 with M = Be, Mg, Ca, Sr, and Sc, reaching very similar conclusions regarding the DFT performance.⁷⁵ 

LDA and GGA functionals also lack a proper description of the strong on-site Coulomb interaction of localized d-electrons of metal atoms. This is especially important for open metal sites (OMSs), which are formed when solvent is removed from a crystalline metal–organic framework, and which largely define the interactions with the guest molecules. To address this issue, DFT methods are appended with an additional Hubbard term and are consequently denoted as DFT + U approaches.⁷⁶ Investigating the interactions of isostructural metal–organic frameworks based on M-MOF-74 (where M = Mg, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, and Zn) with various gaseous and small hydrocarbon molecules (including H₂, CO₂, N₂, H₂O, H₂S, NH₃, CH₄, C₂H₄, etc.), Lee et al. demonstrated that the Hubbard U corrections are necessary to correctly predict not only of the ground-state multiplicity of the transition metals in MOFs but also of the binding energies in their host–guest complexes.⁷⁷ On the other hand, Kim et al. showed that adding Hubbard U-corrections to PBE-D2 binding energies between small hydrocarbons and Fe-MOF-74 significantly worsens their agreement with experiment.⁷⁸ This demonstrates the difficulties faced when investigating OMSs.

To reproduce experimental heats of adsorption beyond qualitative trends, zero-point vibrational energy and thermal energy corrections must be added to compute electronic energies of the host–guest binding. This can be achieved by numerical differentiation of the analytical first derivatives to estimate vibrational frequencies. Heats of adsorption of N₂, CO, and CO₂ on Mg-MOF-74, computed in this manner at the B3LYP-D^* level of theory, were shown to be in a reasonably good agreement with the data from variable-temperature infrared spectroscopy.⁷⁹ Delle Piane et al. demonstrated that dispersion correction is needed to predict spontaneous adsorption of ibuprofen on MOF MCM-41 walls (in accordance with experimental observations), while the lack of dispersion correction lead to positive values of the free energies of adsorption.⁸⁰ Vibrational frequencies, computed using a finite-differences approach to obtain a Hessian matrix of second derivatives of the total energy, were used by Lee et al.⁷⁷ to assess the quantum nuclear zero-point vibrational energies and finite-temperature thermal energies under the harmonic approximation at the vdW-D2 + U2 level of theory. Figure 3(a) clearly demonstrates the need to include these effects to reproduce experimental heats of adsorption.

While many periodic approaches use ab initio pseudopotentials and plane-waves (such as the projector augmented wave method),⁸¹ a fully periodic hybrid Hartree–Fock/DFT simulation of CO₂ adsorption on MIL-100(MIII) with MIII = Al, Sc, Cr, and Fe using instead atom-centered basis sets was performed by D'Amore et al.⁸² Specifically, dispersion-corrected B3LYP and MO6 computations employed a triple-zeta pob-TZVP/6–311G(d,p) basis set on the metal atoms, a 6–311G^** Pople basis set on the H, C, and O atoms of the framework, and a TZP basis set on the CO₂ molecule. A massively parallel version of the CRYSTAL code was used in these computations involving over 50 000 basis functions. The computed interaction energies, corrected for the basis set superposition error, are in reasonable agreement with experiment.

Computed binding energies can be used to assess the performance of diverse MOFs in the separation of gas mixtures. For example, Sholl and co-workers considered the M-BTC framework, where BTC is the 1,3,5-benzenetricarboxylic acid, and M = Mg, Ti, V, Cr, Mo, Mn, Fe, Ru, Co, Ni, Cu, and Zn, for the separation of the water, carbon monoxide, ethane, ethene, and ethyne [Figs. 3(b) and 3(c)], demonstrating a relationship between binding energies and separation efficiencies.⁸³ Similarly, Cunha et al. rationalized the impact of the ligand functionalization on the release of caffeine from functionalized MIL-53 via binding energies.⁸⁴ Furthermore, periodic DFT computations allow elaborating the mechanistic details of, e.g., separation of C₂H₄ and C₂H₆ on ZIF-7,⁸⁵ conversion of NO into NO₂ in Mn-MOF-74,⁸⁶ fixation of CO₂ in HbMOF1,⁸⁷ or catalytic hydrogenation of CO₂ by a Lewis pair functionalized MOF, UiO-66-P-BF₂.^88,89

Apart from structural and energetic considerations, ab initio modeling allows characterizing the electronic structure of the host–guest complexes in organic frameworks. Bader charge analysis⁹⁰ is often employed to establish the oxidation state of the metal node in MOFs, and to rationalize the extent of charge transfer between the framework and the guest.^85,91,92 For example, Canepa et al. suggested that a relationship exists between the charge of the metal ion in M-MOF-74, where M = Be, Mg, Al, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Sr, Zr, Nb, Ru, Rh, Pd, La, W, Os, Ir, and Pt, and the strength of adsorption of small molecules (H₂, CO₂, CH₄, and H₂O) in its pores.⁹³ A similar observation was also made by Koh et al. for CO₂ adsorption on M-DOBDC and M-HKUST-1, where M = Be, Mg, Ca, Sr, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Mo, W, Sn, and Pb.⁷⁵ These findings put forward the metal ion charge as a potential descriptor of adsorption energy for the high-throughput screening of MOFs (see Sec. V A below). Jensen et al. simulated the transition dipole moments and the photoluminescence spectra for a Zn-based MOF, RPM3-Zn, in the presence of nitroaromatic guests.⁹⁴ Computed photoluminescence quenching was in a good agreement with the experimentally measured one and is attributed to the shift of the lowest unoccupied molecular orbital from the framework to the guest. This selective photoluminescence quenching by nitroaromatics potentially enables the sensing of explosives by the studied MOF.

Periodic computations for MOF–guest systems generally focus on gaseous guests. On one hand, this is because applications involving gases, such as hydrogen storage or flue gas separation, are among the most common ones. On the other hand, this can be attributed to the difficulty of obtaining good starting geometries for large guests in MOFs. Starting structures can be based on chemical intuition,⁹⁵ electrostatic complementarity,⁸⁰ or modeling using force fields.^84,96,97 Although the choice of a starting geometry is crucial for correctly predicting the final adsorption site, it is often overlooked in the literature.^94,98,99 Devautour-Vinot et al.⁹⁶ examined the adsorption of caffeine on functionalized UiO-66(Zr) MOFs with periodic DFT computations. This MOF has tetrahedral and octahedral cages. The authors found that all starting configurations in the tetrahedral cage converge to the same geometry and that this adsorption site is preferred over the ones in the octahedral cage. Despite these difficulties, there are many reasons why periodic computations are important for large guests. The more space the guest occupies in the pore, the more sites are available for its interactions with the host. To describe all these interactions properly, either very large finite models or periodic computations are required. Moreover, the steric hindrance within the framework can only be captured in the periodic setting. Finally, and importantly, periodic computations are necessary to model the long-range interactions and crystal field effects.

While periodic models adequately reflect the realistic structure of the framework material, their sheer size prohibits the use of advanced ab initio methods, capable of accurately describing static and dynamic electron correlation, which can be critical in multireference (e.g., containing open-shell metal sites) and weakly bound systems (such as many host–guest complexes). To capture these effects, density functionals beyond GGA, as well as single- and multireference post-HF methods can be employed provided the studied system is small enough to make such computations feasible.¹⁰⁰ This is achieved by “cutting out” smaller, finite clusters out of the periodic structures (experimental or optimized in silico) of the guest-occluded frameworks, and “saturating” them with, e.g., hydrogen atoms to preserve charge. The cluster is typically defined based on chemical intuition and includes the metal center (in the case of MOFs), the guest molecules, and the parts of the framework interacting with or closest to it (Fig. 4); to the best of our knowledge, this process is yet to be automatized. Numerous methods and codes are available to analyze the energy, electron density, and orbitals in the finite setting. For example, the latter, straightforwardly obtained with finite computations, must instead be constructed from projections in the periodic case.

Many common exchange-correlation density functionals lack in their treatment of dynamic correlation, which is the key component of the attractive interactions between small gas molecules and organic frameworks. To overcome this limitation, in their studies of N₂/CH₄ separation in V-MOF-74 and Fe-MOF-74, Lee et al. employed a non-local Rutgers-Chalmers van der Waals density functional with the Hubbard U correction (vdW-DF2+U), several Minnesota density functionals containing the kinetic energy density term (M06-L, M06, and M11-L), and two wavefunction theory methods: local-pair natural-orbital coupled cluster theory with single and double excitations (LPNO-CCSD), and complete active space second-order perturbation theory with counterpoise corrections (CASPT2-CP).¹⁰¹ DFT computations were performed on the periodic model and clusters of two sizes: a “large” cluster containing 88 atoms, including three metal centers, and a “small” cluster composed of 19 atoms, including one metal center [Fig. 4(a)]; the wavefunction theory methods were only feasible for the finite models. Across all models and methods there is a good agreement in terms of the adsorption preference for H₂ vs CH₄ on the two studied MOFs: Fe-MOF-74 does not display any appreciable selectivity, while for V-MOF-74, N₂ binds much stronger than CH₄ (Table I). Similarly, Supronowicz et al. observed an excellent agreement (less than 2 kJ mol⁻¹ difference) between B3LYP-D3/TZVP and CCSD(T)/CBS interaction energies of H₂O with a Cu₂(HCOO)₄ finite model of the HKUST-1 MOF.¹⁰² They also showed that increasing the basis set size from TZVP to def2-QZVPP does not change the qualitative trends in the binding energies of various gases on the finite model of HKUST-1, and leads to only minor decrease in their absolute values. Koukaras et al. also showed that PBE/def2-TZVPP interaction energies of amino acids (glycine and tyrosine) in IRMOF-14 are in good agreement with the MP2 results.¹⁰³ 

Multiconfigurational wavefunction methods become indispensable when an accurate description of the open metal sites is required. Stoneburner et al. investigated interaction energies between various gases and M-catecholate MOFs for M = Mg²⁺, Sc²⁺, Ti²⁺, V²⁺, Cr²⁺, Mn²⁺, Fe²⁺, Co²⁺, Ni²⁺, Cu²⁺, and Zn²⁺ using complete and restricted active space (CAS and RAS, respectively) formalisms of the self-consistent field and second-order perturbation theory (SCF and PT2, respectively), in addition to various flavors of DFT.^104,105 Again, a good agreement was achieved for the overall trends across all tested methods [Fig. 5(a)]. Importantly, multireference methods not only provide more accurate energetics compared to the DFT but also yield physical insight into the host–guest interactions. Bernales et al. turned to multiconfigurational methods in conjunction with the DFT geometries to study the catalytic mechanism of ethylene dimerization by a cobalt(II) or a nickel(II) single-site catalyst supported on the zirconia-like nodes of the NU-1000 MOF [Fig. 4(b)].¹⁰⁶ Analysis of the complete active space molecular orbitals in the rate-limiting transition state allowed explaining superior catalytic activity of the Ni(II) active site compared to Cu.

Depending on the studied system and its application, the size of the cluster model must be adjusted. For example, a cluster consisting of only the linker was sufficient to investigate the influence of OH linker functionalization on the drug delivery performance of IRMOF-16.¹⁰⁷ Pirillo and Hijikata studied how the adsorption properties of several gases at the open metal site of a rhodium paddle-wheel unit are influenced by guests adsorbed on the unit's other side.¹⁰⁸ To consider this effect across a metal–metal bond, a cluster model containing only the metal node was sufficient. Most cluster models consist of one metal node and one or several linkers, which allows capturing (to a certain extent) the interplay between these structural units.¹⁰⁹ However, when multiple adsorption sites in a pore are likely important, the model must contain a large portion of, or even an entire pore. For instance, Supronowicz et al. used a model of an entire pore of HKUST-1 to study the interaction of various biologically important organic molecules with this MOF containing undercoordinated copper centers.¹¹⁰ Ernst and Gryn'ova probed the importance of the appropriate cluster model choice for the adsorption of 4,4′-bipyridine and 1,2-bis(4-pyridyl)ethane on GW-MOF.¹¹¹ Considering two model clusters for each guest, adsorption energies were computed for both the periodic and the finite models [Figs. 4(e) and 4(f)] at a range of DFT levels; these were subsequently refined and analyzed using energy decomposition analyses (EDAs)¹¹² to identify the physical underpinnings of the host–guest interactions (see Sec. III D below). Larger clusters including two metal centers closely reproduced the trends in the interaction energy from the periodic model, while predictions for the smaller clusters led to qualitatively wrong trends [Fig. 5(b)].

Cluster models are often employed to investigate the mechanisms of reactions inside the framework materials using density functionals beyond the generalized-gradient approximation, common in the periodic approaches. For example, the M06-L functional, which includes short-to-medium-range van der Waals interactions, has been successfully applied to study mechanisms of catalytic ethylene hydrogenation and dimerization by Ir(CO)₂ and Ir(C₂H₄)₂ complexes on UiO-66 and NU-1000 supports [Fig. 4(d)],¹¹³ catalytic propane oxidative dehydrogenation by activated Co-SIM+NU-1000 MOF,¹¹⁴ and carbonyl-ene reaction between encapsulated formaldehyde and propylene on M3(btc)₂ (M = Fe, Co, Ni, Cu and Zn).¹¹⁵ Chen et al. used B3LYP to investigate the catalytic hydrolysis of organophosphate warfare agents by MOF NU-1000 after benchmarking HF, M06, M06–2X, B3LYP, and a modified complete basis set method (CBS-QB3) on the uncatalyzed reaction.¹¹⁶ Islamoglu et al. used the M06-L functional to study the effect of amino groups in proximity of metal nodes on the organophosphate decomposition reaction.¹¹⁷ 

Complete-active space approaches, CASSCF and CASPT2, can be employed to refine the DFT energies in multireference systems, as was done in the study of CH-activation mechanism by MIL (Materials of Institute Lavoisier) metal–organic frameworks.¹¹⁸ Such multiconfigurational approaches become indispensable when modeling optically active systems, such as MOF-based luminescent sensors for small molecules. Hidalgo-Rosa et al. investigated the luminescence quenching in a Eu-MOF sensor for aniline¹¹⁹ and a Cd-MOF sensor for 4-nitroaniline¹²⁰ at a cluster model level using CASSCF and NEVPT2 methods in addition to time-dependent density functional theory, elucidating the energy transfer pathways crucial for the turn-off luminescence mechanism of these chemosensors (Fig. 6). To investigate the hydrogen bond between formaldehyde and the luminescent MOF [Zn(NH₂bdc)(bix)]_n, Yao et al. also used time-dependent density functional theory (CAM-B3LYP functional) and established that luminescence is due to the ligand-to-ligand rather than ligand-to-metal charge transfer.¹²¹ 

Instead of choosing between the periodic and finite models, combining them can yield an in-depth insight into the electronic structure whilst preserving a realistic representation of the system. Such “hybrid” approaches often employ periodic computations at the LDA and GGA level with subsequent geometry and energy refinement on a cluster model at higher-rung DFT and post-HF levels of theory. Dixit et al. demonstrated that periodic PBE computations adequately describe the hydrogen adsorption on MOF-5 (M = Li, Be, Mg, and Al) when compared to hybrid DFT and MP2 cluster modeling.¹²² A hybrid approach coined MP2:B3LYP+D^* was introduced to quantify the interaction of CO and CO₂ with the CPO-27-M metal–organic frameworks (M = Mg²⁺, Ni²⁺, Zn²⁺): the structures of adsorption complexes were generated at the B3LYP-D^* level with periodic boundary conditions, and their interaction energies were computed for representative clusters using perturbation theory.¹²³ This approach reproduces the experimental sequence of binding energies for the metal centers, and yields heats of adsorption within ca. 2 kJ mol⁻¹ of the experimental values. Similarly, a hybrid MP2:(PBE+D2)+ΔCCSD(T) method, in which PBE-D2 is used to optimize the adsorption complexes and compute their anharmonic frequencies, and the accurate electronic energies of cluster models are computed at the MP2 level with a correction for higher-order correlation effects at the CCSD(T) level, was used to estimate the adsorption energies of CO on CPO-27-Mg.¹²⁴ The computed adsorption isotherms are in close agreement with experiment, with only ∼1 kJ mol⁻¹ deviation in Gibbs free energies. Nachtigall and co-workers introduced a so-called DFT/CC approach, in which a correction term, computed for a cluster model using DFT, CCSD(T) and ΔMP2 extrapolation to complete basis set limit, is added to the periodic PBE/PAW results to approximate the CCSD(T)/CBS accuracy (Fig. 7).^125,126 When applied to model CO₂ adsorption on CuBTC MOF, this method yields a very good agreement with experiment.¹²⁷ Finally, Yu et al. developed an MP2-based QM/MM method to compute the binding energy of CO₂ on CPO-27-M (M = Mg, Mn, Fe, Co, Ni, Cu, and Zn).¹²⁸ As in the previous examples, the adsorption complex geometries were first relaxed using periodic DFT (PW91 functional), and subsequently used to cut out finite clusters. Geometries of these clusters were further refined at the B3LYP/6-31+G(d,p) level, and their energies computed at the MP2/6-31+G(d,p) level of theory. These finite-model ab initio adsorption energies were then corrected to account for the periodicity of the real system using molecular mechanics (Dreiding and Universal force fields). This approach afforded good agreement with experimental heats of adsorption (Table II).

A good tactic in football involves understanding the interaction between the players and the ball in detail. Likewise, a good tactic for designing MOFs and COFs requires detailed insight into the host–guest interaction mechanism and type. Diverse approaches to analyzing the non-covalent interactions are well established for molecular chemistry and reactivity,^129,130 biomolecular complexes,^131,132 molecular crystals,^133,134 but are somewhat less common in MOF-guest systems. Among these approaches, energy decomposition techniques discern the energetic contributions to the total interaction energy and the physical forces behind the interactions, while analysis of canonical and/or localized orbitals, as well as electron density and its derivatives allow localizing, quantifying, and visualizing the interactions.

Energy decomposition analyses (EDAs) yield quantitative partitioning of the interaction energy into physically meaningful components, e.g., electrostatic, Pauli repulsion, polarization, orbital mixing, and dispersion term. Tsivion et al. used absolutely localized molecular orbital energy decomposition analysis (ALMO-EDA) to study the interaction of hydrogen molecules with typical MOF linkers.¹³⁵ While interactions with the non-metalated linkers are driven by dispersion and weak charge transfer, they range from significant charge transfer to mostly electrostatic in metalated linkers depending on their polarity. Similarly, the ethylene dimerization by fluorinated organic ligand on palladium MOF catalysts¹⁰⁹ and the H₂ adsorption on a vanadium MOF¹³⁶ were analyzed with the ALMO-EDA approach. The Kitaura–Morokuma-Ziegler EDA was used to explain the selective adsorption of certain gas molecules on rhodium and copper paddle wheel nodes,¹⁰⁸ and to examine the interactions of H₂, CH₄, C₂H₄, C₂H₆, CO₂, and CS₂ with a Hofmann-type MOF, modeled by a gas molecule sandwiched between two pyrazine molecules.¹³⁷ Furthermore, the Morokuma–Ziegler EDA, the natural orbitals for chemical valence (NOCV) approach, and the non-covalent interactions index (NCI) were employed to investigate the charge transfer channels between MOF [Zn₂(OBA)₄(BYP)₂] and nitrobenzene, and between MOF [Cd₂(H₂L)₂(H2O)₅]_0.5H₂O explaining the turn-off mechanism of luminescence in these chemosensing MOFs.^120,138 Symmetry-adapted perturbation theory (SAPT) revealed dispersion and electrostatics as main forces behind the sorption of H₂ on the open metal sites of HKUST1-MOF, while dispersion and induction were found to be the main contributors to the selective adsorption of krypton over xenon in an ultramicroporous MOF [Ca(C₄O₄)(H₂O)] [Fig. 8(a)].^139,140

Koukaras et al. investigated IRMOF-14 and IRMOF-16 for the delivery of an anticancer agent tamoxifen.¹⁴¹ Organic linkers of the two MOFs were modified by inserting a hydroxyl group to enable hydrogen bonding with tamoxifen. Based on the spin-component scaled second-order Møller–Plesset perturbation theory [SCS(MI)-MP2] computations, OH-IRMOF-16 has a higher binding affinity toward tamoxifen than OH-IRMOF-14. To further analyze the interactions responsible for this preference, the non-covalent interaction index method was applied to the cluster model of hydroxy-functionalized IRMOF-16. NCI confirmed the strong attractive hydrogen bonding interaction between the nitrogen atom of tamoxifen and the hydroxyl unit of the MOF linker [Fig. 8(b)]. Tan et al. used DFT and NCI to explain the experimentally observed high adsorption selectivity for methanol and ethanol over other volatile organic compounds on DPPB-2 (the amorphized form of MOF DPPB-1).¹⁴² 

Ernst and Gryn'ova compared SAPT and Bickelhaupt–Baerends EDAs, as well the NCI and density overlap regions indicator (DORI) tools for two types of cluster models of the GW-MOF hosting 4,4′-bipyridine or 1,2-bis(4-pyridyl)ethane. Both EDA methods yielded the same qualitative results (Fig.9), whilst both NCI and DORI confirmed that the host–guest interactions are dominated by hydrogen bonds.¹¹¹ 

So far, this review focused on the host–guest interactions in the context of capture/adsorption, storage, transformation, and release of the guests. However, guests can also affect the macroscopic bulk properties of the MOFs. A prominent example is the breathing effect of certain MOFs possessing a dynamic framework. In their pioneering work, Loiseau et al. investigated the structural basis for the large breathing effect observed for MIL-53 upon hydration. The water molecules in the one-dimensional channels of this highly flexible MOF form hydrogen bonds with the framework, causing a substantial change in the pore size whilst preserving the topology.¹⁴³ The breathing effect influences not only the pore size but also the magnetic¹⁴⁴ and mechanical¹⁴⁵ properties of MOFs.

Dong et al. evaluated the elastic tensors of MOF-901 free of solvent and occluded with 5.7 and 11.4 wt. % methanol. According to DFT (PBE) computations of the unit-cell compression, the stability of MOF-901 increases with increasing methanol content, likely due to methanol molecules filling the pores and leading to greater structural stability.¹⁴⁶ Canepa et al. computed the elastic constants, bulk, shear, and Young's moduli of empty MOF-74-Zn and loaded with H₂, CO₂, CH₄, and H₂O using the vdW-DF functional.¹⁴⁷ The incorporation of the adsorbates generally leads to a substantial increase in the elastic constants, indicative of the loss in the MOF's flexibility. The largest increase is induced by water, which was also found to have the strongest interaction with the framework. Furthermore, the adsorbate-induced changes to the MOF's heat capacity, computed from the ab initio phonon frequencies, were found to be related to the host–guest interaction energy: the larger the energy, the larger the heat capacity [Fig. 10(a)]. Specifically, the largest heat capacity among the studied guests were found for water, in line with the interaction energies: −73.9 kJ mol⁻¹ for H₂O, −52.4 kJ mol⁻¹ for CO₂, −40.1 kJ mol⁻¹ for CH₄, and only −20.9 kJ mol⁻¹ for H₂.

Water molecules in the pores also affect the electric properties of MOFs. Using a combination of periodic DFT computations, impedance measurements, accurate X-ray diffraction measurements and charge density modelling, Scatena et al. showed that the dielectric constant of HKUST-1 differs substantially depending on whether its pores are empty, half-filled, or filled with water molecules [Fig. 10(b)].¹⁴⁸ For the same MOF, it was shown that introducing 7,7,8,8-tetracyanoquinododimethane (TCNQ) guest molecules leads to a six orders of magnitude increase in the electrical conductivity.¹⁴⁹  Ab initio calculations suggest that every TCNQ guest molecule bridges four copper paddle wheel units. This guest molecule introduces an additional empty band into the band gap of the MOF, resulting in a new charge transfer band. Similarly, guests have been shown to affect the framework's magnetic properties. For example, periodic vdW-DFT computations were employed to investigate the spin-crossover transition temperature in empty and guest-occluded {Fe(pyrazine)[Pt(CN)₄]} MOF.¹⁵⁰ Crystal orbital displacement curves were computed to analyze the stabilization or destabilization due to the host–guest interactions and the charge transfer between the host and the guest.

Luminescent MOFs find potential application in light-emitting diodes, lasers, and as chemical sensors.⁹⁹ When used as the latter, their fluorescence (and phosphorescence) is selectively turned on or off by guest molecules based on their type and/or quantity and tuned by the host–guest interactions. Yang and Yan studied the fluorescence of 4-(dicyanomethylene)-2-methyl-6-(4-dimethylaminostyryl)-4H-pyran (DCM) incorporated into stilbene-MOF and IRMOF-8. DCM is a laser dye with red fluorescence, while these MOFs show blue emission under UV irradiation. The incorporation of DCM into them resulted in a blue/red two-color emission. Periodic DFT computations demonstrated that an energy transfer from the MOF to the DCM molecule occurs during the excitation process. Furthermore, the authors probed the luminescence response of the host-guest systems when volatile organic solvents were added. This affected the emission wavelength and the intensity ratio of the blue to red emission, which demonstrates that these MOF-DCM systems could be useful ratiometric luminescence sensors. The phosphorescence of two Zn-terephthalate (TPA) MOFs, 1-DMF and MOF-5, upon treatment with a pyridine solution was also investigated.¹⁵¹ Frontier orbital analysis revealed that in the empty MOF the HOMO is mainly located on Zn₂O₈ metal nanoclusters, while the LUMO is primarily located on the TPA ligands. Pyridine guest introduces intermediate energy levels participating in the photoemission process and shifting the wavelength of the emission. Hidalgo-Rosa et al. rationalized the selective fluorescence quenching in [Zn₂(OBA)₄(BPY)₂] MOF by nitrobenzene but not by toluene using time-dependent DFT and finite models.¹¹⁹ They showed that the fluorescence emission is due to a linker-to-linker charge transfer involving π-type orbitals. DFT geometry optimizations revealed that the absorption of nitrobenzene leads to the formation of host-guest hydrogen bonds. In this system, the LUMO is now located on nitrobenzene, from which the energy dissipates through the nonradiative transition resulting in the turn-off of the fluorescence. Instead, toluene does not introduce a LUMO between the HOMO and LUMO of the empty MOF.

Given the enormous number of possible metal–organic frameworks, selecting the structure most suited to a specific application (i.e., a specific guest molecule) necessitates efficient and rapid techniques to explore this vast chemical space. High-throughput screening (HTS) approaches benefitting from continuously increasing computer power and improving simulation algorithms enable searching thousands of structures according to target criteria. These approaches typically involve either constructing a new or utilizing an existing database of MOFs, such as the Computation-Ready, Experimental (CoRE) metal–organic framework database,¹⁵² computing their relevant physico-chemical properties (descriptors) using an appropriate, i.e., simultaneously sufficiently accurate and computationally feasible method, and ranking the structures so as to identify the most promising candidates.^153,154

In 2012, the Snurr research group constructed an extensive database of hypothetical MOFs (ca. 138 000 structures) from the building blocks of existing frameworks;¹⁵⁵ this hMOF database is available at http://hmofs.northwestern.edu. Since then, it has been repeatedly employed to search for MOFs with desirable adsorption properties, and to identify the relevant structure–property relationships. In the seminal publication,¹⁵⁵ the hMOF database was screened for methane-storage capacity. The most promising candidate, NOTT-107, was then synthesized and shown to outperform the best material for CH₄ storage at the time, PCN-14. Improved materials were discovered in the hMOF database for Xe/Kr separation,¹⁵⁶ and for many other applications by identifying the relevant screening criteria. For example, a clear correlation was established between purely structural (e.g., pore size, surface area, and pore volume) and chemical features (such as functional groups) of MOFs and their gas adsorption properties (Fig. 11).¹⁵⁷ To rapidly identify MOFs with superior hydrogen storage capacity, a simple metric—the binding fraction, i.e., the fraction of the unit cell volume within a given distance of the framework—was introduced and shown to reliably reflect the gas storage capacity of the framework.¹⁵⁸ The ratio of Henry's law constants between CO₂ and H₂O was used to identify MOFs with high CO₂ selectivity under high humidity conditions.¹⁵⁹ 

The Grand Canonical Monte Carlo method is among the most frequently used simulation techniques in high-throughput workflows due to its relatively low computational cost. GCMC high-throughput screening allows estimating typical adsorbent selection criteria, such as adsorption selectivity, adsorbent performance score (APS), sorbent selection parameter, working capacity, and regenerability of MOFs. Candidate structures are generally ranked according to a combination of these metrics. This approach was used to discover frameworks with high CO₂/CH₄ selectivity and total loading above 0.5 mol kg⁻¹ among ca. 3000 existing MOFs,¹⁶⁰ candidates for separating the dibranched hexane isomers from their linear and monobranched counterparts out of 12 351 existing MOFs,¹⁶¹ and frameworks for CH₄/H₂ separation among almost 55 000 structures from the Cambridge Structural Database.¹⁶² In the latter example, a simple mathematical model was also proposed to predict the CH₄/H₂ selectivity.

Density functional theory computations, either in a periodic or in a finite fashion, can be included at the later stages of a high-throughput screening workflow to refine the predictions from GCMC or other semi-empirical models for selected candidates. Demir et al. screened the CoRE MOFs database for O₂ and N₂ binding using GCMC for the initial selection, and periodic and cluster DFT to compute accurate binding energies.¹⁶³ To find the best candidates from the same database for propane/propene separation by diffusion, Pramudya et al. combined molecular dynamics, nudged elastic band (NEB) with the classical UFF-FM force field, and DFT computations in a high-throughput workflow.¹⁶⁴ DFT was particularly necessary when modeling MOFs with long flexible linkers, as the NEB UFF-FM method could not provide reliable energy barrier estimates for these systems. Periodic DFT simulations were also used to obtain accurate geometries of 936 anion-pillared MOFs prior to GCMC simulations of Xe and Kr adsorption in the high-throughput screening of MOFs for Xe/Kr separation.¹⁶⁵ 

The entire screening procedure can be performed at the ab initio level instead of classical simulations, typically using the LDA or GGA density functionals. For example, Canepa et al. investigated the adsorption of H₂, CO₂, CH₄, and H₂O in MOF-74-M (M = Be, Mg, Al, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Sr, Zr, Nb, Ru, Rh, Pd, La, W, Os, Ir, and Pt) using periodic vdW-DF computations, and discovered that only a few noble metals afford selective adsorption of CO₂ over H₂O [Fig. 12(a)].⁹³ Rosen et al. developed a fully automated, high-throughput periodic density functional theory workflow and applied it to screen the CoRE database for MOF catalysts for methane activation.¹⁶⁶ Periodic PBE-D3(BJ) computations revealed that the oxidation of the open metal site, and not the C–H bond activation in CH₄, is the rate determining step, suggesting that improved MOF heterogeneous catalysts can be achieved with low-valence, redox-active open metal sites. Extending this work, the authors screened 60 MOFs for the catalytic activation of methane and demonstrated that the active site formation energy is a reliable descriptor of both the thermodynamic stability of the metal-oxo active site and the framework's ability to activate the C–H bonds [Fig. 12(b)].¹⁶⁷ Furthermore, Rosen et al. screened a series of MOFs comprised of 6 MOF families (MOF-74, MOF-74-S, MIL-88B, MIL-88B-OH, and MAF-X, where MAF = metal-azolate framework, X = bridging anion, e.g., μ-Br–, μ-Cl–, μ-F–, μ-SH–, or μ-OH–) and diverse metals (V²⁺, Cr²⁺, Mn²⁺, Fe²⁺, Co²⁺, Ni²⁺, Cu²⁺, and Zn²⁺ cations for MOFs with M²⁺ binding sites, and Sc³⁺, Ti³⁺, V³⁺, Cr³⁺, Mn³⁺, Fe³⁺, Co³⁺, and Ni³⁺ cations for MOFs with M³⁺ binding sites) for O₂ and N₂ adsorption.¹⁶⁸ The MOF unit cell shapes and volumes were optimized at the periodic PBE-(D3)BJ level of theory with Hubbard U correction, and the atomic positions and adsorption energies were obtained at the M06-L level. Gas molecules were systematically placed inside the pores using the MOF adsorbate initializer program,¹⁶⁶ and the adsorption modes were analyzed using the CrystalNN bonding topology algorithm.¹⁶⁹ This study elucidated the effect of the bridging anion on the gas adsorption strength and demonstrated how this effect can be utilized to develop more selective MOF adsorbents. Finally, Fumanal et al. based their high-throughput screening of MOF photocatalysts on two energy-based descriptors: the UV-vis light absorption capability and the band energy alignment.¹⁷⁰ Both descriptors were derived from the electronic structure features (bandgap, ionization potential, etc.) assessed at the periodic PBE and (selectively) PBE0 levels of theory.

Apart from physical and chemical descriptors, the topology of the framework itself can be used in the high-throughput screening, as was demonstrated by First et al. for the separation of gases (O₂/N₂, CO₂/N₂, CO₂/CH₄, and CO₂/H₂) and chemicals (propane/propylene, ethane/ethylene, styrene/ethylbenzene, and xylenes) in MOFs.¹⁷¹ 1690 frameworks were assessed based on the shape selectivity criterion, which considers the minimum-energy pathway of a guest molecule through the material and is evaluated using the MOFomics code.¹⁷² 

High-throughput screening generates extensive sets of structural and computed data; thus, tailored approaches are needed to analyze it, to discover promising candidates, and to arrive at practical material design guidelines. Moghadam et al. developed an interactive visualization tool, which assists in analyzing the multidimensional structure–property plots, obtained in the HTS, and identifying the relationships between structural properties and adsorbent performance.¹⁷³ Applying this technique in combination with the GCMC and DFT computations to a database of 2932 existing MOFs, the authors identified UMCM-152 MOF as a promising candidate for O₂ storage (Fig. 13); this prediction was confirmed experimentally, resulting in a world-record performance, 22.5% over the best performing system at the time. The interactive 5D visualization tool for comprehensive data mining is accessible at http://aam.ceb.cam.ac.uk/mof-explorer.

While high-throughput screening certainly offers a relatively rapid means to map the large structure–property space of the framework materials, the computational and time costs of the commonly employed GCMC and DFT methods remain a major bottleneck when applied to databases featuring more than 100 000 materials. Moreover, large amounts of data, generated in high-throughput workflows, are, on one hand, challenging to analyze, but, on the other hand, attractive for big-data approaches. This is where machine learning (ML) can accelerate materials discovery by as much as 2–3 orders of magnitude, offering efficient ways to predict the material properties, analyze the structure–property relationships, and even reverse-design new material architectures. Applying machine learning to materials discovery generally entails (i) preparing a sufficiently large and representative training set containing material structures and computed (e.g., through high-throughput screening) or measured properties, (ii) selecting an appropriate machine-readable representation of a material, (iii) choosing an efficient machine learning model (architecture), and (iv) visualizing and analyzing the results of the ML prediction. In the past 4 years, several review articles have been published on this topic, both focusing on specific applications in gas storage and separation,^174–176 and providing a more general overview of ML methods in this field of research.^177–179 Below we highlight recent and innovative ML approaches to the prediction and design of framework materials hosting molecular guests.

In addition to the sets of existing and hypothetical MOF structures, for which various properties have been computed with the grand-canonical Monte Carlo and/or density functional theory in a high-throughput manner (see Sec. V A), a set containing experimental and DFT-computed data for 15 713 MOF, entitled QMOF, has been developed by the Snurr group with the specific goals to train and test machine learning models, and to facilitate the discovery of new materials.¹⁸⁰ Using the SOAP similarity kernel for unsupervised dimensionality reduction with a uniform manifold approximation and projection (UMAP) algorithm, the authors constructed an insightful structure (experimental crystal structures)—property [PBE-D3(BJ) computed band gaps] map of QMOF (Fig. 14).

Choosing the right descriptor of the material structure and chemistry is crucial to obtain reliable ML predictions for a reasonably sized training set. In 2013, the Woo research group developed an atomic property weighted radial distribution function (AP-RDF) descriptor¹⁸¹ specifically tailored to gas adsorption in MOFs. The AP-RDF score employs Gasteiger's radial distribution functions, popular in chemoinformatics, and weights them according to atomic electronegativity, polarizability, and van der Waals volume, thus capturing geometric and chemical features simultaneously. Approximately 58 000 structures from the hMOF database were evaluated for their CH₂, N₂, and CO₂ uptake capacity. Both the principal component analysis (PCA) and the quantitative structure–property relationship (QSPR) method using AP-RDF scores were able to correctly identify the high-performing MOFs, with R² values for the predicted uptake capacity above 0.7. A webserver implementation of the AP-RDF scores using multilinear regression (MLR) and support vector machine (SVM) learning models for predicting the uptake capacities called MOFIA (MOF Informatics Analysis) was also developed. This tool was subsequently used to discover MOFs with high methane storage capacities based on the interplay of two structural parameters—dominant pore size and void fraction—using a nonlinear support vector machine model [Fig. 15(a)].¹⁸² Similarly, Halder and Singh found the void fraction to be the most significant factor defining the C₂H₆/C₂H₄ selectivity when screening the hMOF database,¹⁸³ while Yang and co-workers identified the pore limiting diameter, followed by the largest cavity diameter and porosity, as most important for the gas separation efficiency of existing MOFs using principal component analysis.¹⁸⁴ For the selective adsorption of p-xylene over o- and m-isomers in MOFs, the pore limiting diameter and largest cavity diameter were again found to me the most important performance descriptors.¹⁸⁵ AP-RDF descriptors were also utilized to screen MOFs for CO₂ capture in conjunction with a classifier machine learning model.¹⁸⁶ Specifically, instead of training ML to predict the absolute uptake capacities, an SVM with the AP-RDF scores was trained to classify whether a given MOF has low or high CO₂ uptake capacity. Out of 292 050 tested MOFs, only ca. 10% were selected by the classifier for subsequent high-throughput screening, thus greatly diminishing the associated computational effort. In another study, three AP-RDF descriptors were combined with six geometric descriptors in a gradient boosted trees regression (GBTR) to predict CO₂ working capacities and CO₂/H₂ selectivities of 358 400 hypothetical MOFs with R² values of 0.944 and 0.872, respectively [Fig. 15(b)].¹⁸⁷ This approach allowed identifying the 1000 best performing MOFs for CO₂/H₂ separation.

In 2017, Pardakhti et al. appended typical structural descriptors of MOFs (such as density, maximum and dominant pore diameter, void fraction, and gravimetric surface area), computed using Monte Carlo approaches, with various descriptors of the framework's atomic composition, extracted from crystallographic data.¹⁸⁸ These include the type and number of each atom, degree of unsaturation, ratios of the quantities of different atoms, electronegativity, and metallic percentage. Across all tested ML models, a combination of structural and chemical descriptors leads to better prediction of the methane uptake compared to models using descriptors of only one class (Table III). The best results—an R² of 0.98 and a mean absolute percent error of about 7%—were achieved with the random forest algorithm, which took approximately 2 h on a single personal computer to train the model and predict adsorption capacity of 130 398 MOFs. The Froudakis group employed a similar concept of ML descriptors for MOFs, based on the types of atoms in the material rather than its building blocks, and also achieved the best performance in predicting methane and carbon dioxide adsorption capacities of 137 953 hypothetical MOFs with a combination of geometric and atom type descriptors.¹⁸⁹ Moreover, a model, trained on MOFs only and using the atom types in addition to the structural descriptors, was able to accurately predict the methane adsorption capacity in 69 839 covalent organic frameworks (Table IV).

Machine learning models come in many flavors.¹⁹⁰ The tree-based methods, such as decision trees (DTs), regression trees (RTs), and random forests (RFs), recursively split the training sample according to simple decision rules. Kernel-based methods, e.g., support vector machines (SVMs), kernel ridge regression (KRR), and Gaussian process regression (GPR), non-linearly transform the inputs into a higher dimensional space. Finally, artificial neural networks (NNs) are repeated compositions of simple functions, typically arranged in consecutive layers. To identify the best ML model for a specific type of data and chemical problem, several models are often compared in terms of their learning speeds and prediction accuracies. Yang et al. evaluated four ML algorithms—DT, RF, SVM, and backpropagation neural network (BPNN)—for predicting the gas separation performance of 6013 computation-ready, experimental metal–organic framework membranes (CoRE-MOFMs) using fivefold cross-validation.¹⁸⁴ Wu and co-workers compared the predictive power of support vector machine, random forest regression, and gradient boosting regression tree models for methane adsorption in 130 397 hypothetical MOFs.¹⁹¹ In both studies, the random forest algorithm outperformed all other considered ML models. For the training set based on experimental results of CO₂ fixture with MOF catalysts from approximately 100 published papers, SVM, stochastic gradient descent (SGD), and NN models performed better than decision trees and K-nearest neighbor classification models (Fig. 16).¹⁹² Finally, the backpropagation neural network was shown to outperform the partial least-square (PLS) method in predicting the adsorption capacity of organic sulfur gases in 4764 CORE-MOFs after being trained on the hMOF database.¹⁹³ 

In 2020, the Froudakis group addressed the problem of the model choice by adapting an Automated Machine Learning (AutoML) architecture to train statistical and ML models and estimate their predictive performance for MOFs' chemical properties.¹⁹⁴ This approach allows tuning the hyper-parameters of the ML models and avoids the common issue of overfitting, accurately assessing predictive performance even for relatively small (less than 100 datapoints) sets of experimental data. The pipeline is implemented in the Just Add Data Bio, or JADBio tool, and an illustrative prediction of CO₂ and CH₂ adsorption under various thermodynamic conditions is available at https://app.jadbio.com/share/86477fd7-d467-464d-ac41-fcbb0475444b.

Very recently, Van Speybroeck and co-workers developed an approach for constructing machine learning potentials (MLPs) for MOFs.¹⁹⁵ In contrast to other MLP construction techniques that require costly molecular dynamics trajectories at the density functional theory level as training data, this new incremental learning scheme employs machine learning potentials to generate the training data, metadynamics to ensure broad sampling of the phase space, and a message passing neural network for efficient on-the-fly learning. In this manner, accurate and transferrable ML potentials for MOFs can be generated from as little as a few hundred single-point DFT computations.

Automated workflows combining high-throughput screening with machine learning prediction offer unparallel speed of materials discovery. For example, Wand et al. combined hierarchical screening of the CORE-MOF database based using GCMC simulations with the crystal graph convolutional neural networks algorithm to develop a predictive model of methane adsorption (Fig. 17).¹⁹⁶ This model was transferred to randomly selected covalent organic frameworks and zeolitic imidazolate framework materials, as well as to the hMOF database, in which a material with the maximum working capacity close to the Department of Energy's 2015 target was discovered. Reverse design of a metal–organic framework for efficient CO₂ capture in humid conditions was achieved by Zhang and co-workers.¹⁹⁷ 27 combinations of metal nodes and framework topologies were screened using a combination of Monte Carlo tree search (MCST) and recurrent neural network (RNN) using three performance descriptors: high adsorption performance, experimental accessibility, and good hydrophobicity. Several newly designed systems were found to display a good trade-off between adsorption selectivity and uptake capacity whilst being sufficiently hydrophobic and synthetically viable. Very recently, Chen and co-workers screened several experimental and hypothetical databases in search of MOFs with open Cu sites for efficient isobutene/isobutane separation.¹⁹⁸ Configuration-bias Monte Carlo was used to simulate adsorption, and molecular dynamics—diffusion of the hydrocarbons in MOFs. The generated data were then used to train five ML models in conjunction with seven physical descriptors to predict the separation performance of MOFs. Finally, a materials-genomics strategy was applied to identify the best-performing genes (nodes and linkers) and cross-assemble them into novel MOFs. On this basis, five MOFs—four from the screened databases and one newly assembled—were found to have high thermal stability, isobutene uptake, and isobutene/isobutane selectivity. Stable adsorption configurations of these candidates were verified using GGA and hybrid DFT computations.

Covalent organic frameworks are relatively recent, metal-free analogues of MOFs, offering a number of advantages over their famous ancestors while retaining the large internal surface area, ordered pores, structural diversity, tunable functionalization, and the ability to absorb, store, and release guest molecules. Compared to MOFs, COFs are more lightweight, environmentally sustainable, and thermally stable. Whereas MOFs can form relatively strong interactions with guest molecules via their empty metal sites, the host–guest interactions in COFs are exclusively non-covalent. Many COFs consist of 2D layers of covalently connected nodes and linkers, which are assembled in the third dimension via non-covalent interactions. The stacking of the 2D layers is crucial for the properties of the COF¹⁹⁹ yet is difficult to determine.^200,201 While the x-ray diffraction is the predominant method to structurally characterize COF crystals, the diffraction peaks represent an average over the measured portion of the crystal and, therefore, random layer offsets are difficult to determine and to distinguish from eclipsed stacking. The structural characterization of COFs is further complicated by the difficulty of growing single crystals suitable for single crystal x-ray diffraction. Thus, they are often characterized with powder XRD. Even more challenging is the structure determination of COFs encapsulating guest molecules. Consequently, even though a few experimentally determined host-guest structures exist for MOFs, to the best of our knowledge, none have been reported for COFs.

Methodological approaches to the modeling of host–guest interactions in COFs are fairly similar to those used for MOFs, e.g., density functional theory and other ab initio methods, molecular dynamics and Monte Carlo simulations, as well as their multiscale combinations.²⁰² The absence of the metal atoms has several important implications for the modeling of COFs. On the one hand, it alleviates the need to employ relativistic and Hubbard U-corrections. On the other hand, the choice of the cluster model is less intuitive than in MOFs and requires cleaving covalent single and even double bonds.^203,204

Classical molecular dynamics simulations, often combined with periodic DFT modeling, were employed to investigate the interactions of COFs with anticancer drugs 5-fluorouracil²⁰⁵ and quercetin,²⁰⁶ pollutant bisphenol A [Fig. 18(a)],²⁰⁷ to rationalize the separation of CO₂/N₂ mixtures in a 2D-COF,²⁰⁸ and to explain the ultrahigh adsorption affinity for dispersed anionic pollutants in a new three-dimensional framework, 3DCOF-g-VBPPh3Cl [Fig. 18(b)].²⁰⁹  Ab initio molecular dynamics simulations are often evoked to elucidate the interaction sites and to estimate the temperature- and pressure-dependent loading capacities of COFs, particularly for hydrogen storage applications [Fig. 18(c)].^210–212 

Grand canonical Monte Carlo simulations remain the method of choice for modeling the dynamical behavior of COFs upon guest capture and/or release. Tong et al. combined GCMC and configurational bias MC with cluster and periodic DFT computations to investigate the structure–property relationships across a diverse set of 46 COFs for CO₂ capture.²⁰⁸ Large “adsorbility” of adsorbates was put forward as a promising descriptor of the framework's separation performance. Using a similar multiscale approach, Das and Mandal explained the selectivity of a triazine-based benz-bis(imidazole)-bridged COF (TBICOF) toward CO₂ over N₂ and CH₄ [Fig. 19(a)], and toward benzene over cyclohexane.²¹³ Chen and co-workers supplemented periodic DFT and GCMC simulations with the extended transition state—natural orbitals for chemical valence (ETS-NOCV) analysis to study the selectivity of a diamond-topology covalent organic framework, COF-300 (dia-c5), toward various C4 hydrocarbons [Fig. 19(b)].²¹⁴ Diverse topological parameters of the framework were found to be useful in tuning the selectivity and predicting the material's synthesizability. Metropolis Monte Carlo together with periodic DFT were used to investigate the adsorption of environmental pollutants—polybrominated diphenyl ethers—in a reticulated framework, TAPT-DMTA-COF.²¹⁵ 

Several databases of COF structures have been published to enable the high-throughput screening of these materials for specific applications. In 2017, the Zhong group introduced²¹⁶ and in 2018 updated²¹⁷ a computation-ready, experimental covalent organic framework (CoRE COF) database, which includes almost all existing experimental structures. A year later, the database was updated with new structures. The solvent-free and disorder-free structure files of over 280 COFs in 12 topologies are publicly available online at https://core-cof.github.io/CoRE-COF-Database/. Evaluating the structures in the CoRE COF database for CH₄ delivery, the authors showed that the top performing systems have large volumetric surfaces and the pore channels are passable in three dimensions.²¹⁷ Aksu et al. screened CoRE-COF for CO₂/H₂ separation by evaluating the CO₂ and H₂ permeabilities and selectivities of COF membranes using a combination of MD and GCMC.²¹⁸ Materials with narrow pores and low porosities were found to be best for separating CO₂ from H₂, while COFs with the opposite characteristics efficiently separate H₂ from CO₂. Using a similar computational setup, this research group also demonstrated that COFs outperform conventional adsorbents (zeolites and activated carbons) in terms of selectivity of CO₂ separation from flue gas and offer better working capacities than MOFs.²¹⁹ In 2019, the Smit group went beyond the “CoRE” approach, introducing a set of over 300 CURATED (Clean, Uniform, and Refined with Automatic Tracking from Experimental Database) COFs together with DFT geometries and point charges, available at https://www.materialscloud.org/discover/curated-cofs.²²⁰ A workflow for screening this database for promising CO₂ adsorbents was also developed and encoded in the Automated Interactive Infrastructure and Database for Computational Science (AiiDA).²²¹ 

In 2018, the Smit group presented a database of 69 840 hypothetical covalent organic frameworks built from 666 organic linkers and four established synthetic routes.²²² The database was screened for methane storage capacities with GCMC, resulting in a top performer, composed of carbon-carbon bonded triazine linkers in the tbd topology, with a deliverable capacity beating the best existing methane storage materials. Optimal values for several topological features of the frameworks, including density and surface area, were proposed to achieve the best adsorbed natural gas storage performance (Fig. 20). This database was also screened by the Smit group for carbon capture using parasitic energy as a key performance descriptor.²²³ Out of over 69 000 hypothetical COFs, more than 70 outperformed the best experimental COFs, and several performed similarly to Mg-MOF-74. Finally, Lan and co-workers presented a methodology for a high-throughput construction of covalent organic frameworks based on the materials genomics strategy mimicking their natural growth processes.²²⁴ From a library of 130, a database of ca. 470 000 COFs with known and ten new topologies was constructed. Two three-dimensional COFs with a novel ffc topology were successfully synthesized.

In 2023, a Ready-to-use and Diverse Database of Covalent Organic Frameworks with Force field based Energy Evaluation (ReDD-COFFEE) was published and made freely available via Materials Cloud (https://doi.org/10.24435/materialscloud:nw-3j).²²⁵ It contains 268 687 DFT-optimized hypothetical COF structures together with bespoke force field parameters for each framework. The structures were generated from 279 diverse secondary building blocks, arranged in 1114 distinct two- and three-dimensional topologies; only structurally stable and likely synthesizable systems were included in ReDD-COFFEE. System-specific periodic force fields were obtained from cluster force fields, derived in turn for each SBU by fitting to ab initio data, and were shown to outperform a popular generic UFF in describing the crystal structures of several experimentally characterized COFs. Finally, a subset of 10 000 structures from ReDD-COFFEE was screened using GCMC for vehicular methane storage, yielding candidates with outstanding volumetric deliverable capacities.

In silico modeling of the host–guest complexes between organic framework materials and their molecular guests, such as gases, small hydrocarbons, and therapeutics, offers valuable insight into the dynamic and static, structural and electronic aspects of the host–guest interactions, which can be used to rationalize experimental observations and identify the best performing systems (Fig. 21). Molecular dynamics and Monte Carlo methods allow modeling guest diffusion through and adsorption in the framework's pores, simulating the adsorption isotherms and evaluating the loading capacities. Periodic density functional theory with dispersion, Hubbard U (for MOFs), and thermal corrections provides refined structures and heats of adsorption, while cluster models treated with more accurate DFT, post-HF, and multireference methods yield information about the precise physical nature of the interactions, reaction mechanisms, and response properties of the framework–guest complexes. Interaction analysis and partitioning schemes deliver a fundamental understanding of the role of the framework's building blocks and functionalization. Employing some of these modeling techniques in a multiscale manner, high-throughput screening and machine learning are performed on the databases of existing and hypothetical framework materials to not only rank the candidates according to their performance in a given application, but also to elucidate the structural parameters determining performance, thus opening the doors toward tuning of the known and reverse-designing the new MOFs and COFs.

If computational modeling is so powerful, then why have not the best (most stable, lightweight, sustainable, high-capacity, responsive) organic frameworks been already developed for applications in molecular capture, separation, transport, storage, catalysis, etc.? For one, this area of research is plagued by the common dilemma of in silico modeling: realistic representation of the system vs accuracy of the theoretical approach. No method is capable of describing large, periodic frameworks hosting numerous guest molecules under relevant conditions (solvent, pH, pressure, temperature, irradiation, etc.) across practically relevant timescales whilst fully taking into account static and dynamic correlation, non-adiabaticity, relativistic effects, and many other challenging aspects of the electronic structure. Whether all these factors are important, and to what extent, strongly depends on the investigated host-guest system. Cluster models address some of these methodological problems but sacrifice realistic system representation and necessitate cluster selection—a procedure so far lacking well defined rules and automatized implementations. With a few notable exceptions,^226–230 energy decomposition and density and wavefunction partitioning schemes, helpful in understanding the intricacies of the host–guest interactions, remain limited to these finite and relatively small cluster models. Similarly, predictions of the influence of guest occlusion on the solid-state properties of the frameworks, such as their band structures and densities of states, remain scarce. Large realistic models are treatable with molecular mechanics force fields and GGA and LDA DFT methods; however, using these “cheaper” approaches often leads to compromised numerical accuracy. Since it is rarely computationally feasible to perform accurate high-level electronic structure theory computations on large, realistic models, the reliability of the various approximations can only be verified against experimental data. The latter is generally limited, preventing across-the-board validation of these methods. For the same reason, high-throughput screening and machine learning predictions utilize computational rather than experimental data and are, therefore, subject to the same error bars as the underlying in silico methods. Moreover, databases of experimental structures, fed into HTS and ML, are often biased toward high-performance materials. Insufficient sampling in the regions of chemical space containing low-performing and experimentally overlooked organic frameworks compromises the reliability of machine learning when it is used to predict the properties of novel and unconventional materials. Finally, the disconnect between theory and experiment is perhaps most pronounced in the practical viability of in silico predicted materials, which are often “good on paper” (in a computer), but are unstable, hard or impossible to manufacture, especially at scale, etc.

In the light (or, more appropriately, the gloom) of the previous paragraph, is there no hope for the computationally designed breakthrough for organic frameworks supporting tailored interactions with molecular guests? In this Review, we have collected illustrative, successful recipes from the literature for addressing the aforementioned challenges. Individual simulation techniques, restricted to specific system sizes and timescales at a given accuracy, are combined in multilevel and multiscale approaches balancing large system sizes with reliable predictions. Computations on thousands of framework materials using multiple codes are facilitated by automated workflows, and their results are shared freely and openly via community-ran repositories. Together with unbiased databases of hypothetical MOF and COF structures, published in open access, they are actively used in training machine learning models, which can predict not only the interaction features and performance metrics, but even the synthesizability of the frameworks.^231–234 With accurate quantum-chemical insight on one side and big-data analysis on the other, precise design guidelines based on the compositional and topological parameters of the frameworks are elucidated. Ultimately, they yield predictions of candidate materials with outstanding performances, at times overtaking the best systems in use at the time. Occasionally, these predictions are subsequently validated by experiment.

Which challenges should be tackled next? In our opinion, automated cluster selection will substantially facilitate the computational analysis of the host–guest interactions in organic frameworks and is likely to benefit greatly from cross-disciplinary approaches. Furthermore, community standards for in silico generated data, including metadata records, will enable post-processing and reuse of these data in future high-throughput and machine learning studies, significantly saving computational resources and time whilst simultaneously expanding the training sets in size and heterogeneity. Similarly, more experimental data, especially measured on a range of systems under identical conditions, as well as reporting of the “failed” experiments, will allow broader benchmarking of the methods and development of more universally applicable predictive models. Ultimately, the key to the continuing success in the rapidly growing field of metal–organic and covalent organic framework adsorbents, filters, catalysts, and delivery agents is in communication between the main players—method developers, computational materials scientists, and experimentalists.

M.E. and G.G. acknowledge the support from the Klaus Tschira Foundation and funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement No. 101042290 “PATTERNCHEM”). J.D.E. is the recipient of an Australian Research Council Discovery Early Career Award (Project No. DE220100163) funded by the Australian Government.

The authors have no conflicts to disclose.

Michelle Ernst: Conceptualization (equal); Visualization (supporting); Writing – original draft (lead); Writing – review & editing (equal). Jack Evans: Writing – review & editing (equal). Ganna Gryn'ova: Conceptualization (equal); Funding acquisition (lead); Supervision (lead); Visualization (lead); Writing – original draft (supporting); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.