Quantum chemical modeling of hydrogen binding in metal–organic frameworks: validation, insight, predictions and challenges

Abstract

A detailed chemical understanding of H₂ interactions with binding sites in the nanoporous crystalline structure of metal–organic frameworks (MOFs) can lay a sound basis for the design of new sorbent materials. Computational quantum chemical calculations can aid in this quest. To set the stage, we review general thermodynamic considerations that control the usable storage capacity of a sorbent. We then discuss cluster modeling of H₂ ligation at MOF binding sites using state-of-the-art density functional theory (DFT) calculations, and how the binding can be understood using energy decomposition analysis (EDA). Employing these tools, we illustrate the connections between the character of the MOF binding site and the associated adsorption thermodynamics using four experimentally characterized MOFs, highlighting the role of open metal sites (OMSs) in accessing binding strengths relevant to room temperature storage. The sorbents are MOF-5, with no open metal sites, Ni₂(m-dobdc), containing Lewis acidic Ni(II) sites, Cu(I)-MFU-4l, containing π basic Cu(I) sites and V₂Cl_2.8(btdd), also containing π-basic V(II) sites. We next explore the potential for binding multiple H₂ molecules at a single metal site, with thermodynamics useful for storage at ambient temperature; a materials design goal which has not yet been experimentally demonstrated. Computations on Ca²⁺ or Mg²⁺ bound to catecholate or Ca²⁺ bound to porphyrin show the potential for binding up to 4 H₂; there is precedent for the inclusion of both catecholate and porphyrin motifs in MOFs. Turning to transition metals, we discuss the prediction that two H₂ molecules can bind at V(II)-MFU-4l, a material that has been synthesized with solvent coordinated to the V(II) site. Additional calculations demonstrate binding three equivalents of hydrogen per OMS in Sc(I) or Ti(I)-exchanged MFU-4l. Overall, the results suggest promising prospects for experimentally realizing higher capacity hydrogen storage MOFs, if nontrivial synthetic and desolvation challenges can be overcome. Coupled with the unbounded chemical diversity of MOFs, there is ample scope for additional exploration and discovery.

---

1 Introduction

Touted as the Swiss Army Knife for decarbonization, hydrogen provides a versatile alternative to fossil fuels applicable across multiple economic sectors, including transportation, power generation, chemical production (e.g. hydrocarbons, ammonia), and industrial manufacturing (e.g. steel).^1,2 The future hydrogen economy^3–5 will be driven by clean hydrogen, using the reactive H–H chemical bond as the medium for energy storage. Fuel cell-based energy release from H₂ consumption yields only water, which can, in turn, be regenerated into H₂via the reverse processes of water oxidation and proton reduction, resulting in minimal CO₂ generation if green electricity is used. Four serious barriers stand in the way of this vision: (i) obtaining H₂ from green sources at a viable cost, (ii) the storage challenge (which is the motivating issue for this work), (iii) efficient and economical fuel cells, and (iv) implementation of large-scale infrastructure.

The H₂ storage challenge is well-recognized^6–17 with potential solutions ranging from storage at cryogenic temperatures, or under high pressures,^18,19 or in materials like metal hydrides^20,21 and nanoporous frameworks.^22–27 The low volumetric energy density of hydrogen when compared to gasoline²² leads to storage at pressures above 350 bar or at liquid hydrogen temperatures, each with very significant handling costs.^28,29 Metal–organic frameworks (MOFs) potentially solve this problem since they densify hydrogen in the solid state, allowing its storage and release at more modest pressures.^22,23

The original goal of hydrogen storage research was to address light vehicle needs. In pursuit of H₂ powered vehicles, the U.S. Department of Energy (DOE) has set a volumetric capacity target of 40 g L⁻¹ for the year 2025 at temperatures in the range of −40 to 60 °C (233 K to 333 K), and a delivery pressure in the range of 5 to 12 bar.^15,30 To meet the system level targets, the storage capacity of material-based solutions likely needs to match the density of liquid H₂ (71 g L⁻¹) due to the limited tank space in light-duty vehicles. Consequently, the US DOE has been exploring the potential use of H₂ for heavy-duty applications (e.g., trucks, rail, maritime) as well as power sector applications (e.g., backup and stationary power).³¹ Indeed, recent techno-economic analysis suggests that the cost of hydrogen storage in MOFs for backup power could be lower than the cost of a liquid hydrogen storage system and comparable to the cost of compressed H₂ at 350 bar,²⁹ although significant improvements in materials design and manufacturing are necessary for the practical realization of a MOF-based technology. Use of MOFs for H₂ storage stands in the context of the explosion of exciting MOF-based science^32,33 that makes this area an active frontier of chemistry and materials science,³⁴ across areas ranging from gas separation³⁵ to water purification³⁶ to catalysis.^37,38

Towards the goal of rational materials design for H₂ storage in MOFs, it is essential to understand the interactions between hydrogen and MOFs at a fundamental level. As we shall see, it is particularly interesting to explore the interactions between hydrogen and open metal sites in MOFs. In a general sense, these interactions originate from the normal driving forces of host–guest binding: permanent and induced electrostatics, dispersion, charge-transfer, and Pauli repulsions.^39–41 However, the details of site-specific binding can be quite complex and multi-faceted, encompassing a range of weak dispersion interactions as well as stronger associations with local charges or an accessible metal site. MOFs, with their high surface-to-volume ratio, are ideally suited for van der Waals interactions with H₂,^42,43 but these weak interactions alone do not suffice for storage under ambient or near-ambient conditions (down to −50 °C).^22,23 Put simply, molecular H₂ is relatively inert in its interactions with the stable organic linkers.

Fortunately, the modular nature of MOFs allows for the incorporation of more diverse binding motifs, either at open (coordinately unsaturated) metal sites (OMSs)^44,45 or on linkers⁴⁶ which may themselves be metal-decorated.⁴⁷ While outside our present scope, other porous materials can also support OMSs, such as zeolites where cation exchange reactions can replace protons with metal ions, such as alkalis, alkaline earths or transition metals.⁴⁸ The flexibility of MOFs makes them unique in this arena and enables tuning of hydrogen interactions with the MOF surface through judicious chemical modifications, albeit subject to the constraints of synthetic feasibility and suitable mechanical stability of the resulting MOF. A great amount of experimental effort over the past two decades has resulted in considerable progress toward this objective.^{22,23,44–46} Our purpose here is not to review that ground-breaking work, but rather to present computational explorations that illustrate some key aspects of the prospects for going further, particularly with regard to potential binding site design paradigms. The first main topic is discussing optimal binding site characteristics based on thermodynamic considerations that maximize the usable storage capacity of H₂ in MOFs at ambient to near-ambient temperatures under the constraints of a fixed pressure swing between loading and unloading conditions.^49–51 This analysis leads to an optimal value of the binding free energy for chosen working conditions. In turn, this determines a set of optimal enthalpy–entropy tradeoffs that are coupled to optimizing the density of accessible sites.

The exploration of hydrogen-binding interactions in MOFs that have been prepared or are possible candidates for preparation is possible with accurate modern density functional theory (DFT).^52,53 As we discuss in the second section, such calculations yield good numerical predictions on the one hand, and, on the other hand, insights into the origin of binding strength via the tools of energy decomposition analysis (EDA).^54,55 The fidelity of the calculations depends upon the model of the MOF: we use cluster calculations that aim to capture all relevant structural and electronic features close to a hydrogen binding site; cluster modeling is also reviewed.

The third section summarizes DFT calculations on H₂ binding in four separate MOFs, each of which has already been experimentally characterized. These examples span a range of H₂ binding energies spanning the physisorption regime through to relatively strong chemisorption. From weakest binding to strongest binding, we begin with MOF-5 (Zn₄O(bdc)₃, bdc²⁻ = 1,4-benzenedicarboxylate),⁵⁶ which lacks open metal sites, followed by Ni₂(m-dobdc) (m-dobdc⁴⁻ = 4,6-dioxido-1,3-benzenedicarboxylate) which has Ni(II) sites,^25,57 and then the exceptionally strong H₂ binding in Cu(I)-MFU-4l.^58,59 The role of open transition metal sites in MOFs is one pathway to go beyond physisorption binding strengths,⁴⁴ and even in the case of Cu(I)-MFU-4l, to a binding that is stronger than optimal. The final example of this section is the near-optimal binding in V₂Cl_2.8(btdd) (H₂btdd = bis(1H-1,2,3-triazolo[4,5-b],[4′,5′-i])dibenzo[1,4]dioxin).²⁷ These MOFs are each illustrated in Fig. 1.

Fig. 1 Metal−organic frameworks (MOFs) are promising materials for H₂ storage because they combine a porous structure with an enormous scope for chemical tailoring of H₂ binding sites, as illustrated by the crystal structures shown here. (a) MOF-5, (b) Ni₂(m-dobdc), (c) Cu(I)-MFU-4l and (d) V₂Cl_2.8(btdd). Grey, blue, red, and light green represent carbon, nitrogen, oxygen, and chlorine, respectively. Zinc sites are shaded blue-grey in (a), nickel rendered green in (b), copper rendered orange in (c), and vanadium grey in (d).

The final topic we address is the question of binding multiple H₂ molecules to a single site with appropriate binding energy. Whilst well-recognized as a goal that could be a critical multiplier to the site density in terms of usable storage capacity, it has not yet been realized experimentally, notwithstanding exciting progress⁶⁰ and a range of predictions.^61–70 In the fourth section, we review two promising existing suggestions, and also present some new computational results. We first consider the use of main group ions, Ca²⁺ and Mg²⁺, binding to catecholate functional groups,⁶⁸ which are derivatives of MOF linkers, and therefore represent a pathway to post-synthetic modification of MOFs. Additionally, we present new calculations on calcium porphyrin. We then turn to explore other open transition metal sites, beginning with our previous DFT calculations that demonstrated the feasibility of two hydrogen molecules coordinating at each V(II) site of V(II) exchanged MFU-4l.⁷⁰ Exploring the rare +1 oxidation state, we predict the coordination of three hydrogen molecules is possible at each Sc(I) or Ti(I) metal site in metal exchanged MFU-4l. We conclude with some discussion on the prospects for future advances in MOF-based storage materials, as well as identifying areas where progress in computational modeling would be desirable.

2 Optimal binding for hydrogen storage

Bhatia and Myers pioneered the understanding of optimal hydrogen storage conditions in porous adsorbents.⁴⁹ They proposed an optimal enthalpy change (ΔH) of −15.1 kJ mol⁻¹ for hydrogen storage at 298 K and a pressure swing of 1.5–30 bar, presuming a standard entropy change ΔS° ≈ −8R. Their standard entropy change is derived from a Langmuirian analysis of H₂ adsorption in cylindrical pores of silica and slit pores of carbon. However, Bae and Snurr have shown that Bhatia and Myers' entropy change assumption might not hold for materials with strong H₂ binding sites. By utilizing Grand Canonical Monte Carlo (GCMC) simulations, Bae and Snurr deduced an optimal isosteric heat −Q_st = ΔH_ads for H₂ binding at around −20 kJ mol⁻¹ under a slightly different pressure swing of 1.5–120 bar at 298 K.⁷¹ Experimental studies on zeolites have estimated the optimal enthalpy to be in the range of −22 to −25 kJ mol⁻¹ for 1.5–30 bar pressure swing at 298 K.⁷²

Despite its limitations, the Langmuir model still offers valuable insights into the thermodynamics of H₂ adsorption at specific binding sites within a MOF. This model, rooted in the balance between the energy lowering upon H₂ binding and its entropic drive to gain translational freedom, suggests that the surface coverage (θ) at thermal equilibrium can be defined as . Here, p is the relative pressure to a reference pressure P₀ of 1 bar. The binding constant K(T), which reflects the ratio of adsorbed vs free adsorbate molecules, can be determined by the standard free energy change ΔG°(T) = −RT ln K on binding at the reference pressure.

Examples of hypothetical binding isotherms are shown in Fig. 2(a). The usable coverage, represented by U, is calculated as the difference in site coverage between loading and unloading pressures (p_L and p_UL) expressed as U(T,p_L,p_UL) = θ(T,p_L) − θ(T,p_UL). This dimensionless usable coverage is multiplied by the density of binding sites in the material to obtain the volumetric capacity, which is commonly expressed in g L⁻¹. The usable coverage describes the amount of hydrogen desorbed at p_UL after being loaded at p_L, for a specific temperature. Our problem is to maximize usable capacity by tailoring the material properties of the MOF. The amount of gas released is greater when the pressure swing (from p_L to p_UL) has surface coverage approaching the saturation limit at p_L and as close to zero as possible at p_UL. Given that p_L, p_UL and T should be viewed as specified device parameters, the development of an optimal MOF, therefore, amounts to optimizing the binding constant K(T) to maximize usable coverage for those constraints. Separately, the site density should also be optimized to maximize usable capacity. As illustrated in Fig. 2(a), it is important to avoid steep uptake at very low pressures (i.e. at p_UL), which occurs with strong binding, as well as small uptake at low to moderate pressures (i.e. at p_L), which occurs with weak physisorption. The optimal standard binding free energies, ΔG°*(T), which maximize usable capacity, can be found by variation.⁴⁹

Fig. 2 An analysis of optimal binding interactions for high-density hydrogen storage on MOFs. The target H₂ standard binding free energy of a site optimizes the usable capacity (U) of a storage material, which (for given temperature, site density and either 0 or 1 H₂ per site) is proportional to the change in binding site occupancy (coverage, θ) between a high loading pressure (e.g. p_L = 100 bar in (a)) and a low unloading pressure (e.g. p_UL = 5 bar in (a)) (a) Binding isotherms for the Langmuir model for 3 different binding free energies at room temperature (T = 298 K), illustrating the loss of usable capacity with too strong binding and too weak binding. (b) The dimensionless usable capacity U(T,p_L,p_UL) as a function of the standard free energy of binding at room temperature, for p_L = 100 bar, p_UL = 5 bar. (c) Optimal binding free energy ΔG°*(T) that correspond to U_max as a function of temperature for 6 different choices of the loading and unloading pressures p_L, p_UL.

Fig. 2(b) shows that U(T) is maximal for a binding free energy change of +7.7 kJ mol⁻¹ at 298 K for a pressure swing of 5–100 bar. This optimal binding isotherm for single-site Langmuir adsorption corresponds to the blue curve at the centre of Fig. 2(a), and can utilize 63% of the available binding sites to release hydrogen from the material.

ΔG°* depends on both temperature and the chosen pressure swing. Fig. 2(c) charts ΔG°* from 223 to 323 K (−50 to 50 °C) for loading pressures p_L = 100 and 170 bar and unloading pressures p_UL = 2, 5 and 12 bar. This data shows that the optimal standard free energy change ΔG°*(T) that enables maximal recovery of hydrogen from the sorbent lies in the range of 5–10 kJ mol⁻¹ under ambient to moderately cooled conditions and these chosen pressure swings.

While we quantify optimality conditions in terms of the standard free energy change on H₂ binding based on the Langmuir model, this is admittedly an approximation, and there may be coupling between binding sites if they are close. While simplified enthalpy–entropy correlation functions for H₂ binding in nanoporous frameworks exist,⁷³ they are not readily generalizable across different framework topologies and binding sites, particularly when the binding arises from the synergy between the molecular orbitals of the adsorbate and the host. Therefore, optimal binding interactions are best described in terms of the binding free energy without loss of generality. Clearly, a range of different enthalpy–entropy tradeoffs can achieve a particular target free energy of binding. All else being equal, the most desirable tradeoff, from a heat management perspective in H₂ storage, is the least loss of entropy and least gain in enthalpy upon binding.

3 Models

With some idea of the target binding free energy and binding enthalpy in mind from the previous section, we consider how to computationally predict and understand the H₂ binding energy (Fig. 3). As illustrated in Fig. 1, H₂ binding sites in MOFs exist inside a periodic three-dimensional material, and upon the approach of a guest molecule to a metal site, there may be strong binding contributions reflecting local chemical interactions at the site, supplemented by weaker long-range dispersion interactions. Three approximations are needed to (approximately) predict binding energies. First, the binding site is modelled, either using a periodic boundary condition (PBC) approximation that neglects disorder or a cluster model that neglects long-range interactions. Second, the exact equations of quantum mechanics are replaced by a density functional theory (DFT) model, which reduces the formal exponential complexity to polynomial (roughly cubic scaling in the number of atoms). Third, the DFT equations are solved using a finite basis set to make the prefactor of polynomial scaling manageable. We shall elaborate briefly on all three of these aspects in the remainder of this section, as well as discuss a fourth issue, which is how to obtain some physical insight from the calculations beyond predicted observables.

The choice between a PBC model and a cluster model is made for MOF binding site modeling based on two considerations. First is physical appropriateness. The PBC approach is most appropriate when long-range interactions are important. However, using the most accurate density functionals such as ωB97M-V is not presently computationally tractable with PBC approaches. Hence, the cluster approach is often favoured when an accurate description of local binding is of primary importance. PBC models are advantageous for examining multiple site adsorption, including adsorption on linkers and interactions over larger scales of the framework, which complements the localized interaction studies provided by cluster models presented here.

Fig. 3 Workflow for DFT cluster calculations of H₂ binding in MOFs. The top panel illustrates the selection of an appropriate cluster model from crystallographic data. The cluster should be centered on the binding site, with all chemically important functional groups around that site included, followed by addition of suitable passivating groups to terminate all dangling bonds at the cluster periphery. The cluster is then DFT-optimized to obtain equilibrium geometries with and without H₂ ligation. The lower panel (going right to left) is a schematic illustration of the energy decomposition analysis (EDA) used to understand the main physical contributions to binding, ΔE_ads, typically using a larger basis set that is augmented on atoms near the binding site. The boxes schematically illustrate EDA constraints that isolate those physical contributions by controlling the atomic orbital to molecular orbital transformations on each fragment (MOF and H₂) at each stage of the EDA sequence.⁵⁵

All subsequent figures illustrating MOF modeling show the cluster chosen to replace the extended framework. Computational tractability versus physical reality likewise governs the choice of the cluster because a factor of two increase in the number of atoms leads to almost an order of magnitude increase in compute costs. On the other hand, the ideal cluster model should accurately capture metal–ligand interactions at short distances while adequately accounting for the dominant dispersion forces at comparatively longer ranges that govern binding in MOF pores. As an example of the considerations involved, we can examine the binding of H₂ at the Cu(I) site in MFU-4l. This scenario can be depicted through a cluster model of the node truncated at its benzotriazolate extremities such that the arene is capped with hydrogen atoms (Fig. 3). These models allow us to validate binding energies and geometries against experimental crystal structure data and measured isosteric heats.^58,59 These validations, in turn, enable the application of a smaller cluster model for computations of higher computational demand, such as frequency calculations. Complete ab initio relaxation of the cluster was performed in response to H₂ loading for MOF topologies that can represent the nodal unit as a closed geometric shape, like in MOFs with a cubic crystalline structure such as MOF-5 and MFU-4l. The MOF V₂Cl_2.8(btdd), which features 1D dimensional chains truncated to the two nearest neighbours, was also subjected to a full cluster relaxation based on its extensive experimental validation.²⁷ On the other hand, Ni₂(m-dobdc), a MOF exhibiting similar 1D chains, the optimization of the cluster was achieved by fixing the bare node's coordinates using values from neutron diffraction data, an approach that has been previously demonstrated to be effective.⁵⁷ A comprehensive understanding of the cluster models utilized in this work can be found in the original literature.^{27,57,59,70,74}

4 Methods

Turning to density functional theory, modern DFT is based on an in-principle exact mapping of the 3n-dimensional many electron integro-differential Schrödinger equation onto a set of n coupled 3-dimensional equations, the Kohn–Sham equations.⁷⁵ However, the KS equations depend on an unknown universal functional that describes exchange and correlation (XC), which in practice is modelled. The resulting approximate functionals are commonly grouped onto five rungs of a metaphorical Jacob's ladder that ascends from the primordial ooze of the Hartree world, where electron correlation is not described at all, to the heaven of chemical accuracy.⁷⁶ In practice, functionals on rungs 1–3 (the local spin density approximation (LSDA), generalized gradient approximations (GGAs) and meta-GGAs) can be tractably evaluated with PBC models. Hybrid density functionals on rung 4, which include a portion of exact exchange,⁷⁷ are computationally more demanding than local or semi-local functionals. This can make calculations with PBC models using hybrid functionals particularly resource-intensive, especially for large systems. However, the accuracy achievable at the 4th rung is known to be substantially higher across a range of chemical energy differences^76,78 that include the noncovalent interactions critical to H₂ binding. The comparative performance is illustrated for 275 H₂ binding energies⁵² (ranging from very weak to very strong) in Fig. 4. It is to gain the benefit of this higher accuracy that we employ cluster modeling in this work.

Fig. 4 Regularized mean absolute percent errors (RegMAPE) (in percent) for a wide range of density functionals assessed against high accuracy benchmarks for 275 H₂ binding energies (defining the H2bind275 dataset). The results are arranged according to the rungs of Jacobs ladder. The error metric, RegMAPE, is constructed to give most weight to binding energies in the key interaction range of 15–25 kJ mol⁻¹ (see ref. 52 for definition). Reproduced from Veccham et al., J. Chem. Theory Comput., 2020, 16(8), 4963–4982; Copyright 2020 American Chemical Society.

Specifically, we employ the ωB97M-V functional,⁷⁹ which is an accurate range-separated hybrid meta-GGA that includes non-local VV10 dispersion,⁸⁰ designed by the combinatorial “survival of the most transferable” protocol.⁷⁶ The VV10 non-local correlation functional provides a proper description of long-range electron correlation responsible for non-covalent interactions.⁸⁰ Metal–ligand interactions due to charge-induced dipoles and orbital overlap entail accurate prediction of main group-transition metal chemistry. ωB97M-V is the top-performing hybrid density functional in several large assessments, including the MGC84 database,⁷⁶ the large and diverse GMTKN55 benchmark dataset,^78,81,82 and the TMC151 transition metal database.⁸³ These findings are buttressed by an extensive benchmark^52,53 of density functionals for hydrogen storage.

Once a density functional is selected, the next step in cluster modeling is to choose an appropriate atomic orbital (AO) basis set, which is used to represent the unknowns, the KS molecular orbitals (MOs). The compute cost of rung 4 DFT calculations rises as the 4th power of the number of AOs per atom,⁸⁴ so again, there is a steep computational penalty to using very large basis sets. At the same time, not every property is equally sensitive to the quality (i.e. size) of the AO basis. In particular, optimized geometries offer the possibility of cancelling the basis set incompleteness error from one set of coordinates to another nearby set of coordinates with virtually identical chemistry and bonding. Therefore geometries are far less basis set sensitive (and are also less functional sensitive) than absolute binding energies since the latter is the difference between two very different configurations.

With these considerations in mind, geometries were optimized in this work using the relatively small def2-SVP basis,⁸⁵ which contains f-polarization at the transition metal site, and p-polarization on H to account for polarization adequately. The ωB97M-V functional was used (with a few exceptions discussed where relevant when B3LYP-D2 was used). Geometries were converged to 3 × 10⁻⁴ kJ mol⁻¹ in energy, and 3 × 10⁻⁵ a.u. in the maximum gradient component. Single point calculations for electronic binding energy, ΔE_ads, utilized ωB97M-V with the augmented def2-TZVPPD basis⁸⁵ at the metal site and dihydrogen, which includes 2f and 1g polarization functions at the metal site and 3p and 1d polarization functions for hydrogen. All single point energy calculations were counterpoise corrected for basis set superposition error, following best practices with basis sets of this size.⁷⁶

Thermochemistry was determined using the B3LYP functional with Grimme's empirical dispersion corrections⁸⁶ in the def2-SVP basis with the same convergence criteria mentioned above, followed by Hessian evaluations at the minimum energy configuration. The total enthalpy change (ΔH(T)) was calculated as ΔH(T) = ΔE_ads + ΔZPVE + ΔH_vib(T) + ΔnRT, with the free energy of binding estimated as ΔH(T) − TΔS(T), with ΔS(T) and ΔH(T) computed using the rigid–rotor–harmonic oscillator (RRHO) approximation. Low-frequency normal modes (≤100 cm⁻¹) were set to 100 cm⁻¹ to account for soft vibrations.⁸⁷ All computations were conducted using QChem 5.3 software⁸⁸ and a development version of Q-Chem 6.1.

If the approximations described above were well chosen, the resulting DFT calculations should be a reliable numerical experiment that yields the thermodynamic and structural information needed to characterize H₂ binding in MOFs. To provide additional insight into why those properties take the values they do, we employ the absolutely localized molecular orbital (ALMO)-EDA to reveal the interaction energy as being due to distinct physically interpretable components,^{54,55,89–91} relative to the noninteracting MOF and H₂. The electronic component of the adsorption enthalpy is expressed as:

ΔE_ads = ΔPrep + ΔFrz^DF + ΔDisp + ΔPol + ΔCT (2)

The preparation energy, ΔPrep quantifies the energy required to deform the unbound MOF and H₂ fragments into their bound states.⁹¹ Using frozen MOs from each fragment, the energy change associated with that frozen density is ΔFrz. Dispersion interactions (ΔDisp) are separated from ΔFrz using the difference between the target density functional (e.g. ωB97M-V⁷⁹ in this work), accounting for non-local van der Waals interactions,⁸⁰ and a suitable dispersion-free functional (in our case, the Hartree–Fock approximation).⁹² This separation leaves a dispersion-free frozen interaction, ΔFrz^DF, associated with permanent electrostatics and Pauli repulsions. The polarization term (ΔPol) reflects the energy reduction resulting from the response of each fragment to the electric field of its neighbour, computed using fragment electrical response functions (FERFs).⁹³ Charge transfer (ΔCT) between the MOF and H₂ is estimated through variational minimization of the system density and complete orbital relaxation, encompassing forward donation from hydrogen σ to the metal d_σ and back-donation from the metal d_π to the hydrogen σ*. The CT component can also be exactly decomposed into “complementary occupied-virtual orbital pair” (COVP) contributions.⁸⁹ Often only one COVP dominates forward and back donation and that most important COVP(s) can then be visualized to yield an orbital picture of the CT process.

5 Characterizing open metal sites for H₂ storage

MOF-5

MOF-5 has been extensively studied due to its potential for hydrogen storage and is an archetypal metal–organic framework with broad implications for gas storage and separation applications.^{32,34,56,94–97} A MOF with no open metal sites, MOF-5 is composed of Zn₄O(bdc)₃ (bdc²⁻ = 1,4-benzenedicarboxylate) units forming a cubic three-dimensional extended porous structure with the unit cell pictured in Fig. 5a. H₂ binds to MOF-5 in an array of binding sites classified as α, β, γ and L depending on its relative positioning with respect to the interstitial oxide at the centre of the Zn₄O tetrahedra, and the linker.^98,99 The strongest binding α site in MOF-5 depicted with a cluster model in Fig. 5b binds H₂ with an experimentally measured enthalpy of −7 kJ mol⁻¹ between 60 to 70 K¹⁰⁰ (lower values of ∼3.5 kJ mol⁻¹ have also been reported¹⁰¹). Our calculations suggest that hydrogen binds at the cup site α in MOF-5 at a distance of 3.43 Å from the interstitial oxide at the centre of the node with ΔE = −8.6 kJ mol⁻¹, which is in qualitative agreement (slightly stronger) than the experiment, and previous computations of H₂ binding in the material at the MP2 level of theory.⁹⁹ This binding is driven primarily by dispersion, as seen by the large ΔDisp of −16.3 kJ mol⁻¹ and weak ΔPol and ΔCT terms (see Fig. 9). A secondary charge transfer energy lowering of −3.4 kJ mol⁻¹ dominated by back-donation from lone pairs on polarized oxygens of the bdc²⁻ linkers to the hydrogen σ* is visualized in Fig. 5c.

Fig. 5 Illustration of hydrogen binding at the cup site in MOF-5, a framework without open metal sites. (a) Cubic unit cell of MOF-5. (b) Cluster model for MOF-5 depicting relatively weak H₂ binding at the α (cup) site, as implied by the long coordination distance. Dative contributions to H₂ binding are very weak in this dispersion-dominated interaction, consistent with the delocalized character of the framework acceptor and donor COVP orbitals shown in (c) and (d) together with the corresponding H₂ σ donor and σ* acceptor orbitals. Distances are labelled in Å with Zn(II) rendered in blue-grey, O in red, H in light grey and C in dark grey. The orbitals were plotted with an isosurface of ±0.07 Å. The smaller formate model was used to visualize charge transfer COVP orbitals.

We must also consider the intricacies of zero-point vibrational energy change in our computational models of MOF-5. Earlier simulations by Sillar, Hoffman, and Sauer have demonstrated that hydrogen molecules exhibit free rotational movement when bound to MOF-5's cup sites.⁹⁹ Using MP2 simulations in conjunction with the multi-site Langmuir model, they were able to replicate experimental binding curves, assuming H₂ retains its rotational freedom in the adsorbed state. In alignment with this, our calculations employ a simplified formate model for the cup site described in previous work^74,99 for thermochemistry calculations and assume unhindered rotation of H₂ in its bound state. This methodological choice leads to predicted binding enthalpies of −10.0 kJ mol⁻¹ at 77 K, which, although stronger than the experimentally observed −7.7 kJ mol⁻¹, serves to illuminate future avenues for refining vibrational ground state predictions. The limitations of the RRHO approximation become particularly pronounced for light elements like hydrogen, especially when interacting anharmonically with MOF-5's delocalized soft modes. A comprehensive solution would necessitate computations of the vibrational ground state allowing full nodal relaxation in response to H₂ uptake and anharmonic coupling between delocalized modes of the framework, posing a formidable computational challenge and presents an avenue for future methodological development.

Contributors to binding due to frozen (which includes the effect of electrostatic or Coulomb interactions and Pauli repulsion), dispersion, polarization, and charge transfer are summarized in Fig. 9.

Ni₂(m-dobdc)

A MOF with Ni(II) sites, Ni₂(m-dobdc) has surpassed the volumetric usable capacity of MOF-5 under ambient conditions.²⁵ As the leading MOF for ambient temperature hydrogen storage, Ni₂(m-dobdc)'s volumetric usable capacity has been measured at 11.0 g L⁻¹ at 298 K at a pressure swing of 5–100 bar (and 23.0 g L⁻¹ with a temperature swing between 198 K and 298 K). This MOF is characterized by the m-dobdc linker, and derives from MOF-74 family of metal–organic frameworks^{25,57,102–109} which feature hexagonally packed cylindrical channels with cross-sections that are lined with Ni(II) sites at a high volumetric density (see Fig. 6a).

Fig. 6 Illustration of H₂ binding at the open metal site in Ni₂(m-dobdc). Panel (a) displays the hexagonal unit cell of Ni₂(m-dobdc), while panel (b) shows the cluster model, with 2 Ni sites, with one on the right binding H₂. Dative interactions are as important as dispersion in the binding, and are dominated by forward donation. For the right-hand Ni site, panel (c) shows the orbitals engaged in forward donation (donor bold, acceptor pale); showing the central role of an empty Ni d orbital as acceptor. Back donation is relatively unimportant, consistent with panel (d) which shows a very diffuse MOF donor orbital (bold) weakly coupling to the H₂ σ* acceptor orbitals(pale). Distances are labelled in Å. The orbitals were plotted with an isosurface of ±0.07 Å.

Given its significance for storage and sequestration applications, MOF-74 and its derivatives have been simulated extensively with electronic structure methods,^57,108,110 and with force fields.^111,112Fig. 6b shows the DFT optimized geometry for a two-site model of the Ni₂(m-dobdc) node where H₂ binds at the distorted octahedral Ni(II) sites (shown in green) at centre of mass (COM) distances of 2.16 Å away from the Ni(II) site (powder neutron diffraction indicates a COM distance of 2.18 Å).⁵⁷Fig. 9, Column II shows binding energy estimates for H₂ binding to the open metal site. The net electronic energy lowering in Ni₂(m-dobdc) (ΔE = −14.5 kJ mol⁻¹) is stronger than MOF-5 and driven by a predilection for polarization (ΔPol = −15.9 kJ mol⁻¹) and charge transfer (ΔCT = −23.8 kJ mol⁻¹). Charge-flow is dominated by the forward donation from H₂ σ to the metal d manifold visualized in Fig. 6(c). Our ΔH(T) of −14.3 kJ mol⁻¹ under ambient conditions aligns well with the experimentally measured value of 13.7 kJ mol⁻¹.²⁵ Though the polarization of H₂ by the local charge at the Ni(II) centres makes H₂ binding enthalpy in this framework stronger than MOF-5, it still falls short of the optimal range sought for attaining maximum usable coverage in sorbent materials.

Simulations of the 2-site cluster model of the node were performed with MOF coordinates fixed, a constraint encapsulated in ΔPrep, the energy required to achieve the bound configurations for both the node and H₂. Geometry optimizations with the MOF cluster fixed led to SCF convergence problems with our preferred functional, ωB97M-V, which we circumvented by employing B3LYP-D2 for geometry optimizations and frequency analysis. Single-point vertical EDA computations were conducted using ωB97M-V. Thermochemical calculations accounted for partial Hessian evaluations of the bound H₂ using the RRHO approximation. It's worth noting that, akin to the case of MOF-5, we assumed free rotation of H₂ in its bound state. This assumption is in line with prior literature and serves as a simplification that achieves better agreement with experimental observables. However, it should be acknowledged that this approximation reflects limitations in vibrational modeling, and future methodological improvements would be desirable.

Cu(I)-MFU-4l

With an isosteric heat of −32 kJ mol⁻¹, one of the strongest for H₂ binding in MOFs, Cu(I)-MFU-4l has been of interest with regard to H₂ binding.^{58,59,113,114} This material comprises cubic unit cells (Fig. 7a) with Kuratowski Secondary Binding Units (SBUs)¹¹⁵ (Fig. 7b), which feature coordinatively unsaturated and initially trigonal pyramidal Cu(I) sites which become pseudo-tetrahedral upon H₂ coordination. The binding enthalpy of H₂ to the Cu(I) site is driven by charge transfer to and from the metal site (see Column III in Fig. 9). H₂ approach to the binding site results in a pyramidalization at the binding Cu(I) site reflected by a large ΔPrep or geometric distortion term, resulting in a destabilization of the d_xz, d_yz orbitals that can back-donate effectively to the hydrogen σ*.⁵⁹ However, the binding in Cu(I)-MFU-4l is excessively strong, with the standard change in binding free energy predicted to be 2.39 kJ mol⁻¹ at a temperature of 298 K, which falls outside the desirable range of between 5–10 kJ mol⁻¹. Additionally, the density of open metal sites within MFU-4l is relatively low, necessitating refinement of both the local composition at the binding site and the global framework topology to establish this material as a viable option for high-capacity hydrogen storage. Within the context of strong Cu(I)-H₂ interactions in a metal–organic framework, there is a recent report that 1% of Cu sites in NU-2100 exhibited high isosteric heat of adsorption (32 kJ mol⁻¹).¹¹⁶ However, computations of the material with periodic DFT at the PBE-D3BJ level of theory predict weaker binding with adsorption energies between −9 to −11 kJ mol⁻¹.

Fig. 7 Illustration of strong H₂ binding at the Cu(I) site in Cu(I)-MFU-4l. Panel (a) shows the cubic unit cell of the material, and (b) shows the cluster model representation H₂ binding in Cu(I)-MFU-4l; note the short 1.7 Å Cu–H₂ distance. Tetrahedral Cu(I) binding site with significant charge transfer orbitals for forward- and back-donation are shown in (c) and (d). Cu(I)-H₂ back-bonding interactions are strong at this binding site (see also Fig. 9). Distances are labelled in Å with Zn(II) rendered in blue-grey, Cu(I) in bronze, N in blue, H in light grey, C in dark grey, and Cl in light-green. The orbitals were plotted with an isosurface of ±0.07 Å. A smaller cluster model truncated at triazolate extremities was used for the visualization of charge transfer orbitals.

V₂Cl_2.8(btdd)

Gagliardi et al. initially put forth V(II)-MOF-74 as a candidate for the separation of N₂ and CH₄.¹⁰⁸ Subsequently, Jaramillo et al. synthesized coordinatively unsaturated dicationic V(II) substitution in V₂Cl_2.8(btdd) (depicted in Fig. 8),¹¹⁷ thereby demonstrating the viability of this material for chemical handle-based separation through the preferential activation of π* orbitals in guest molecules. Further research has revealed the material's efficacy for H₂ storage under ambient conditions.²⁷ This MOF has a mixture of V(II) and V(III) sites that are formed by oxidation of V(II) during synthesis and, therefore, can be effectively modeled as a two-site cluster with a halide capping the V(III) site adjacent to the V(II) OMS where H₂ can coordinate (see Fig. 8(b)). A mechanism similar to that in Ni₂(m-dobdc) prevails here where ΔPol = −6.1 kJ mol⁻¹ polarizes the H₂ in preparation for charge transfer. ΔCT, however is much stronger than Ni₂(m-dobdc) for H₂ binding to the π-basic V(II) site (see column IV of Fig. 9) with back-bonding attenuated when compared to Cu(I) in MFU-4l. The acceptor orbital for forward donation from H₂ has d_z² character on V. Back donation occurs from the d_π orbital to the H₂ σ* orbital.

Fig. 8 Illustration of H₂ binding at the V(II) site in V₂Cl_2.8(btdd). Panel (a) shows the unit cell and (b) shows the DFT optimized geometry for a two-site cluster model representation of H₂ binding to the OMS in the material. Significant charge transfer orbitals for forward and back-donation are depicted in (c) and (d). V(II)–H₂ backbonding at this site is attenuated in comparison to Cu(I)–H₂ interactions in Cu(I)-MFU-4l. Distances are labelled in Å with V(II) rendered in grey, N in blue, Cl in green, H in light grey and C in black. The orbitals were plotted with an isosurface of ±0.07 Å.

Fig. 9 A comparative summary of energy decomposition analysis (EDA) (see eqn (2) for definition of terms) of H₂ binding to four MOFs. ΔCT contains forward (H₂ → MOF) and back-donation (MOF → H₂) contributions. The MOFs are MOF-5, a material with no open metal sites (see Fig. 5), the Lewis acidic Ni(II) site in Ni₂(m-dobdc) (see Fig. 6), and π-basic Cu(I) and V(II) sites in Cu(I)-MFU-4l (see Fig. 7) and V₂Cl_2.8(btdd) (see Fig. 8) respectively. ΔH°(T) values are reported at 298 K and 1 bar, and compared to experimental values.

Comparative analysis

This EDA provides insight into the chemical factors governing optimal interactions for high-capacity hydrogen storage, underscoring the significance of precise modulation of open metal site interactions as a primary design principle. Despite considering 4 quite distinct binding sites, the range of dispersion stabilization for H₂ binding spans a fairly compact range from −16.0 to −24.0 kJ mol⁻¹. Unique characteristics of the open metal sites in each framework lead to divergent adsorption behaviour. In MOF-5, dispersion dominates binding at the cup site, while the polarization of H₂ by unscreened charges at open metal sites facilitates stronger binding in Ni₂(m-dobdc). The Lewis acidic Ni(II) in Ni₂(m-dobdc) exhibits forward-dominant charge transfer, whereas Cu(I) in MFU-4l and V(II) in V₂Cl_2.8(btdd) act as π bases. Excessive charge transfer at the π-basic Cu(I) site in MFU-4l leads to chemisorption and a strongly repulsive ΔFrz arising from increased occupation pressure. In contrast, attenuated back-bonding at the more diffuse and electropositive V(II) site results in near-optimal binding enthalpy for H₂ storage under ambient conditions.

Fig. 10 delineates the predicted ΔG°(T) values for H₂ binding across the four investigated MOFs, spanning temperatures from 223 K (−50 °C) to 323 K (50 °C) under the RRHO approximation. MOF-5, devoid of open metal sites, exhibits binding enthalpies that are insufficient for maximizing room-temperature capacities. In contrast, Ni₂(m-dobdc) offers stronger binding but still fails to meet the requisite threshold for ambient storage. Intriguingly, Cu(I)-MFU-4l, characterized by its strongly π-basic sites, binds H₂ too strongly for practical ambient storage. Among the candidates, V₂Cl_2.8(btdd) stands out for its near-optimal enthalpy–entropy trade-off, boasting an H₂ binding enthalpy of ΔH(298 K, 1 bar) = −21.9kJ mol⁻¹. This enthalpy, in conjunction with a predicted binding free energy change of 10.0 kJ mol⁻¹ at 298 K and 1 bar, places it within the optimal range for maximizing deliverable H₂ capacity under ambient conditions as evidenced by its proximity to the dashed grey curve for optimal free energy of binding ΔG°*(T) under a pressure swing of 5–100 bar. The optimal binding free energy change of 7.7 kJ mol⁻¹ at 298 K (25 °C), and 5.8 kJ mol⁻¹ at 223 K (−50 °C) with this pressure swing are shown by black and blue encircled points. Ni₂(m-dobdc) is closer to optimality conditions for hydrogen storage under moderately chilled conditions. Finally, it is worth recalling that the density of strong binding sites is critical to usable storage capacity. In fact, V₂Cl_2.8(btdd) does not have a high enough site density to compete with the usable capacity of Ni₂(m-dobdc) at ambient temperature despite its superior binding characteristics.

Fig. 10 A comparison of free energy change ΔG°(T) for H₂ binding for the four MOFs under ambient to moderately chilled conditions for T in the range of 223 K (−50 °C) to 323 K (50 °C) and a reference pressure of 1 bar. The dashed grey curve denotes the optimal standard free energy (ΔG°*(T)) for H₂ binding for a pressure swing of 5–100 bar. Optimal free energy changes at 223 K (5.8 kJ mol⁻¹) and 298 K (7.7 kJ mol⁻¹) with this pressure swing are marked with light-blue and black circles.

6 Prospects for multi-H₂ coordination at open metal sites

The case studies in the previous section illustrated several main points. Perhaps the most important point, now generally accepted, is that the chemical diversity of MOFs permits the realization of binding sites that can closely approach the ideal hydrogen binding free energy. This maximizes the usable capacity for room-temperature H₂ storage given a specific site density. A second important point is that quantum chemical cluster calculations using high-quality density functionals can accurately evaluate the binding energy, as evidenced by good agreement with binding enthalpies. Contributions associated with vibrational motion, such as the binding entropy and the ZPE contribution to binding, are the most challenging, and there is room for further methodology improvements.

Relative to the binding of a single H₂ per site, the volumetric and gravimetric storage capacity of a MOF-based material could be greatly enhanced if a given open metal site could bind multiple hydrogen molecules. While fully exposed metal ions can bind many H₂ molecules, the challenge for a viable material is having partially exposed metal sites that are still stable. Indeed, many calculations on specific material designs allude to this possibility.^61–68,70 There is a reported instance of multiple H₂ molecules ligating at Mn in MOF, albeit with a weak interaction with the second H₂ interactions.⁶⁰ Nonetheless, the experimental demonstration of reversible binding of multiple hydrogens in the first coordination sphere of an open metal site with standard binding enthalpies of −20 kJ mol⁻¹ or higher at room temperature remains elusive. There are substantial experimental hurdles to achieving this goal, including appropriate materials design, successful synthesis, and full desolvation of the resulting open metal sites. Exploratory quantum chemical calculations can provide insight into binding sites capable of accommodating multiple H₂ molecules in the desired range, and the purpose of this section is to review and extend some existing work towards this goal.

Catecholates are the doubly deprotonated forms of catechols (for which the parent compound is 1,2-benzenediol, or 2-hydoxyphenol). The resulting dianion can form strong chelating bonds with divalent metal ions of appropriate size,¹¹⁸ where the doubly coordinated metal is significantly exposed. In the context of hydrogen storage, the idea of a quite open metal site on a catecholate, with only 2-fold coordination, is potentially exciting because of the prospect that strong yet reversible binding of more than one H₂ could be achieved, thereby allowing for higher-density storage. Furthermore, incorporating the catechol functionality on a MOF linker¹¹⁹ potentially opens the door for post-synthetic modifications to incorporate the open metal sites. Indeed, there are successful synthetic precedents for incorporating catechol-based linkers into MOFs,^120,121 and demonstrating postsynthetic metalation.¹¹⁹ That progress has synergized with and inspired a range of computational studies of metal-catecholates,^122–124 including for storage of multiple hydrogens at a single site.^66,68

In particular, magnesium and calcium catechols have demonstrated the capability of binding multiple equivalents of H₂, where the electronic binding energy ΔE_ads of each hydrogen may exceed −20 kJ mol⁻¹.⁶⁸Fig. 11 shows the binding of four equivalents of H₂ in calcium catecholate and a summary of the binding energies of multiple H₂ equivalents in catecholate-decorated MOFs. It is promising that the binding energies of more than one H₂ at a single Mg and Ca catecholate site are in the ideal range for usable capacity with a 5–100 bar pressure swing. However, the highly exposed alkaline earth metal binding site will strongly coordinate with the solvent as well, and therefore, part of the synthetic challenge will be the problem of desolvation. If even partial desolvation is possible, then there may still be potential for a usable material because computations have demonstrated that Ca-catecholate can bind two hydrogen equivalents in the ideal range even when there is a residual solvent molecule (acetonitrile) still coordinated at the binding site.⁶⁸ On the other hand, in the corresponding catMg material, the H₂ binding energies were significantly reduced by the presence of a coordinating CH₃CN molecule.

Fig. 11 The binding of multiple H₂ molecules to a single site is possible if the metal is exposed or under-coordinated. Two such examples where the adsorption energy per H₂ is remarkably invariant to the number of H₂'s bound are (a) Ca-catecholate (catCa) and (b) Ca-porphyrin (Ca-por). Calculations predict that both can bind four equivalents of H₂. (c) Differential adsorption energies ΔE_ads for catCa and Ca-por, as well as Mg-catecholate (catMg), and Mg-thiocatecholate, which show the more typical fall-off of binding with number of coordinated H₂'s. The bottom left corner is left empty since catMg can only bind three H₂ molecules.

Multiple hydrogen ligation at an open metal site can also be realized through metal porphyrins, which can be incorporated as linkers in porous materials.^125–130 PCN-224¹³¹ and PCN-221¹³² are examples of MOFs incorporating metallated porphyrins. The study of dihydrogen adduct formation in calcium-containing porphyrin systems is significant due to calcium's abundance and the capacity of Ca-porphyrin to support up to four dihydrogen adducts.¹³³ Past endeavours at quantifying the binding of these H₂ adducts in Ca-porphyrin have used the local density approximation (LDA), and the generalized gradient approximation (GGA) with periodic DFT, giving binding estimates of −24.1 kJ mol⁻¹, and −9.6 kJ mol⁻¹ per hydrogen respectively. In light of the significant difference in predicted binding energies between the LDA and GGA approximations, their limited treatment of dispersion effects, and the improvements that are possible with high-level functionals,^52,53 we employ one such functional here. This ωB97M-V functional⁷⁹ employs the non-local and highly effective VV10 form⁸⁰ for capturing dispersion interactions and is well-suited for quantifying binding at these sites.^52,53Fig. 11 compares differential binding energies for the sequential ligation of multiple equivalents of H₂ in Ca-porphyrin and compares them with calcium and magnesium catecholate functionalized materials that have been investigated earlier.⁶⁸ It is evident that the Ca-porphyrin binding energies are only a little weaker than ideal (∼−12 kJ mol⁻¹) and, furthermore, are remarkably invariant through the binding of 4 H₂ molecules.

Fig. 12 shows Ca-porphyrin binding one, two, three, and four H₂ molecules. The Ca²⁺ ion (shown in green in Fig. 12) in Ca-porphyrin is raised from the plane of the porphyrin ring, and this pyramidalization is supported by experimental evidence on the crystallized Ca-porphyrin.¹²⁹ The hydrogen molecules are bound close to the Ca²⁺ ion, assuming positions within the six-membered pockets that interlace the modified pyrrole subunits at distances ranging from 2.64 Å to 2.75 Å for the sequential ligation of four equivalents. The bound H₂ molecules are little perturbed from the gas phase, as ascertained by a 0.76 Å H–H distance (compared to 0.74 Å in the gas phase) that remains consistent across the board.¹³⁴

Fig. 12 Front-on and sideways perspectives for 4 H₂ molecules binding to Ca-porphyrin with distances from H₂ centre of mass to Ca site (green) labelled in Å (a)–(d), showing a striking invariance to the number of coordinated H₂ molecules. Charge transfer orbitals engaged in forward- and back-donation to and from the metal site for the case of 4 coordinated H₂ (e) and (f). EDA (see eqn (2) for definitions) of the differential ΔE_ads for the ligation of 4 H₂ molecules in Ca-porphyrin (g), showing a key role for dispersion (ΔDisp) and a secondary role for charge transfer (ΔCT). The orbitals were plotted with an isosurface of ±0.07 Å.

Fig. 12g presents the electronic interaction energy ΔE and energy decomposition analysis for the sequential binding of four H₂ molecules near Ca-porphyrin. The binding strength for each hydrogen is stronger than −10 kJ mol⁻¹, with dispersion forces contributing significantly to the overall binding, providing >10 kJ mol⁻¹ of stabilization in each instance. Short-range interactions such as polarization and charge transfer also contribute to the binding, with charge transfer playing a role in the net binding energy albeit to a lesser extent when compared to the binding of H₂ to Ni(II), Cu(I), and V(II) OMS in MOFs studied earlier. Fig. 12e and f shows charge transfer orbitals for forward and back-donation. While the Ca²⁺ ion does accept charge from the ligating hydrogens (Fig. 12e), back-donation occurs primarily from the nitrogen atoms on the porphyrin ring (Fig. 12f), which indicates that the ring itself an active participant in the binding. In summary, our findings corroborate that the binding mechanism of H₂ to Ca-porphyrin is slightly stronger than −10 kJ mol⁻¹ and is largely controlled by dispersion effects augmented by charge transfer. This is in contrast to binding in metallated porphyrins containing transition metals like Ti, where the Kubas interaction plays a significant role in dihydrogen binding.¹³³

Therefore, a broader potential for multiple H₂ coordination lies within the open transition metal sites of MOFs (or, possibly, decoration of linkers with transition metals, although that can also readily lead to overly strong binding¹³⁵). Owing to their unique electronic configurations, some open transition metal sites offer the potential for binding multiple hydrogen molecules, thereby setting the stage for improved usable hydrogen storage capacities. Our recent investigation⁷⁰ focused on the open metal site of V(II)-exchanged MFU-4l, where the parent scaffold supports tetrahedral V(II) sites.¹³⁶ Utilizing state-of-the-art density functional theory calculations, we predicted the binding of two hydrogen molecules at V(II), where the V(II) site starts off with four-coordinate pseudo-tetrahedral coordination.⁷⁰ This binding site also exhibited the potential for tuning the H₂ binding strength by varying the adjacent first coordination sphere halide counterion. Fig. 13(a)–(d) shows V(II) sites in MFU-4l binding 2 H₂ molecules and summarizes the quite wide range of associated adsorption energies with EDA components for each (Fig. 13(e)). While the strongest binding was predicted for F⁻ substituted nodes, the heavier halides exhibited stronger ΔE_ads for binding the second equivalent of H₂ in V(II)-MFU-4l when compared to the first. The orbital synergy between the ligating H₂ molecules and V(II) underpinned the binding, characterized by robust forward donation and uniformly weaker back-bonding stabilization. In particular, the binding strengths for each hydrogen for fluoride substituted nodes in V(II)-MFU-4l were predicted to be around −20 kJ mol⁻¹, which makes it a viable candidate for H₂ storage under ambient or moderately chilled conditions.

Fig. 13 Multiple H₂ ligations in V(II)-X-MFU-4l, (X = F, Cl, Br, I) (a)–(d) and EDA of the quite wide range of adsorption behaviour as a function of the halide counterion (e). Dative effects, particularly forward donation, dominate the binding energy. Distances are reported in Å, and energy contributions in kJ mol⁻¹ in (e).

Intrigued by the prospect that a monovalent metal without the halide at the peripheral site of the MFU-4l node could accommodate three equivalents of H₂ on the completion of its octahedral sphere, we performed a computational search of metal-substituted nodes with monovalent first-row transition metals in their high-spin configuration. In pursuit of identifying monovalent first-row transition metals capable of accommodating three H₂ molecules within the MFU-4l node, we employed a hierarchical computational strategy. Initial screening relied on a truncated small cluster model, focusing on the triazolate-terminated node. This exploratory phase singled out Sc(I) and Ti(I) within the MFU-4l framework as the most promising candidates for binding 3 H₂ molecules at a single site, presumably enabled by the relatively large size of these ions relative to later members of the first transition series. Despite the formidable challenge of synthesizing these unusual metal substitutions in MFU-4l, there are synthetic precedents for these metals in their +1 oxidation states in the literature.^138,139

For more accurate modeling, we then switched to larger benzotriazolate-terminated nodes near the binding sites to refine the Sc(I) and Ti(I) results. The OMSs for Sc(I) and Ti(I) were modelled in their high-spin ground states with S = 1, and respectively.

The results reported were based on geometry optimizations and thermochemistry calculations (vibrational corrections) at the B3LYP-D2 level of theory, while the electronic energies and the associated EDA analysis were evaluated with the ωB97M-V functional. The more computationally expensive frequency calculations for these systems were performed with the model of the node truncated at triazolate extremities.

Fig. 14 illustrates the coordination of 1, 2 and then 3-H₂ molecules to the Sc(I)-MFU-4l site. All three coordination numbers show similar centre of mass distances to the OMS (Fig. 14). Distances to the closest contact at the node for the step-wise attachment of three hydrogens is shown in purple. It is particularly striking that many of these contact distances (each to another hydrogen) are inside the sum of van der Waals radii of two hydrogens (∼2.4 Å). Sc(I) exhibits a ΔE_ads range of −30 to −35 kJ mol⁻¹, slightly weaker than Ti(I), while its ΔH(T) at standard temperature and pressure ranges between −24.7 and −26.9 kJ mol⁻¹ (see Table S2, ESI†), also attenuated compared to Ti(I)-substituted MFU-4l. To further characterize these results, we will discuss EDA, thermochemistry, and binding curves to predict the usable capacities of transmetalated MFU-4l and MFU-4 that has these sites embedded within the material.

Fig. 14 Summary of computational predictions of the ligation of three equivalents of H₂ at the Sc(I) core in Sc(I)-MFU-4l. Distances to the H₂ centre of mass are labeled in black. The quite striking invariance of this distance to the number of ligated H₂ corresponds to very similar binding energies (see Fig. 15 for details), and is a very desirable feature of this binding site. Distances to the closest contact at the node for the bound hydrogens are shown in purple. All distances are in Å.

Fig. 15 reports EDA for binding each of the 3 hydrogens at Sc(I) and Ti(I) sites, respectively. As is the case with our analysis of the four extensively characterized MOFs discussed in the previous section, the dispersion stabilization contributions to the binding energy stay consistent across the board.

Fig. 15 Binding energies and their EDA (see eqn (2)) for the binding of three equivalents of H₂ at Sc(I) and Ti(I) sites in MFU-4l. The Sc(I) sites demonstrate a weak dependence of binding on the number of bound H₂, which is desirable for H₂ storage purposes. By contrast, the Ti(I) sites show fall-off in binding with number of bound H₂. See text and Fig. 16 for the origin of the remarkable change in the ΔFrz^DF and ΔPol EDA components for the Sc(I) site with 2 bound H₂.

The electronic binding energy, ΔE_ads, is comparable for all three hydrogens. However, the EDA results indicate that when the second hydrogen binds to Sc(I), the Pauli repulsion component in ΔFrz is remarkably large in magnitude. Moreover, its ΔPol counterpart is also remarkably large and attractive compared to the first and third hydrogens. We turned to the recently developed polarization analysis¹³⁷ to seek an explanation, resulting in the POL COVPs shown in Fig. 16. Fig. 16(b) reveals that Sc d-orbital occupation changes to relieve Pauli repulsions (i.e. nearly an entire electron promoted from occupied to virtual) for the binding of the second H₂.

Fig. 16 The polarization COVPs (see ref. 137 for theory) associated with relaxing the Sc(I) site upon (a) binding the first H₂, (b) binding the second H₂, and (c) binding the third H₂. Each COVP consists of a donor orbital shown as a solid surface and an acceptor orbital shown as a triangulated surface. The number of electrons rearranged is shown in millielectrons, and it is evident that nearly an entire Sc d electron is rearranged upon the binding of the 2nd H₂; this is the origin of the remarkable EDA results shown in Fig. 15. The orbitals were plotted with an isosurface of ±0.07 Å.

The fact that the largest value of ΔPol appears to reflect relief of Pauli repulsion rather than conventional electrostatic polarization is not unprecedented: we also recently observed this in computations of the binding of water to Cu(I)-MFU-4l.¹⁴⁰ Both the structural factors discussed above and the EDA results indicate how effectively the Sc(I) site serves to pack together the 3 H₂ molecules that it binds.

This sequential binding of three H₂ molecules at a single binding site can be described by the reactions . This process can be represented by the binding polynomial, Q(p)

Q = 1 + K₁p + K₁K₂p² + K₁K₂K₃p³, (3)

where K₁, K₂, and K₃ are the respective binding constants for the sequential ligation of three H₂ molecules.

The average ligand occupancy v(p,T) at a given temperature and pressure is determined by , enabling the derivation of the binding isotherm, which in turn predicts the material's usable capacity. Distinct from a conventional multi-site Langmuir model, which assumes site independence, the stoichiometric sequential multi-site model incorporates higher-order polynomial terms to account for cooperative ligand binding effects.^99,141 This leads to a sigmoidal binding isotherm in some cases (e.g., Sc(I) in Fig. 17a).

Fig. 17 Calculated room temperature binding curves (a), where site coverage on the y-axis measures the number of H₂ bound at a given site. Panel (b) presents the predicted ambient temperature usable capacities for Cu(I)-, Sc(I)- and Ti(I)-MFU-4l and the corresponding MFU-4 frameworks, for the case of a pressure swing between 5 bar and 100 bar. The Sc(I)-MFU-4l data shows the value of coordinating multiple H₂ per site, if the binding free energy is close to the ideal value for room temperature, and this binding free energy is not strongly dependent upon the number of H₂ adsorbed at a given site. By contrast, the single H₂ binding in Cu(I)-MFU-4l is too strong for these conditions, and so is the H₂ binding in Cu(I)-MFU-4l. The higher density of sites possible in the MFU-4 frameworks versus the MFU-4l family is a critical factor that directly scales the usable capacity by nearly a factor of three.

The binding of three hydrogen molecules at the Sc(I) site in MFU-4l results in a binding isotherm that affords steep uptake at room temperature and at a pressure swing between 5–100 bar. Since the binding in Ti(I) in MFU-4l is considerably stronger, the corresponding binding curve shows strong uptake at pressures below 5 bar. These binding curves are compared with that of single H₂ binding in Cu(I)-MFU-4l to provide perspective on the enhancement of site coverage facilitated by multiple H₂ binding at a single open metal site.

The surface coverage can be translated to capacity given the density of open metal sites, for instance, from crystallographic data. Considering the relatively low density of open metal sites in the MFU-4l, we anticipate that the adsorption capacity in g L⁻¹ will be higher for MFU-4, a material with a similar topology but a higher density of binding sites. The MFU-4l unit cell, cubic with an edge of 30.91 Å, comprises eight nodes, each with four tetrahedral sites. The hydrogen adsorption amount is 3.63 g L⁻¹, if we assume that each site is binding a single hydrogen molecule. The smaller unit cell in MFU-4, 21.63 Å, leads to denser metal site packing and a maximum uptake of 10.59 g L⁻¹ if each site can bind only a single molecule, and presuming consistent unit cell dimensions after metal exchange. These estimates focus on hydrogen binding at the OMS, excluding additional adsorption within pores. For comprehensive uptake predictions, Grand Canonical Monte Carlo (GCMC) models, which factor in MOF surface area and free volume, can be utilized.^71,142

Room temperature binding curves parameterized from the free energies of binding are shown in Fig. 17a along with estimates of usable capacity for the MFU-4 and MFU-4l families of materials (see Fig. 17b) for a pressure swing of 5–100 bar. Our results indicate that the multi-ligation of H₂ in Sc(I)-MFU-4 yields impressive usable capacities of 25.65 g L⁻¹ under these conditions, which results from the synergy between multiple H₂ binding, and the nearly 3-fold increase in OMS density from using MFU-4 rather than MFU-4l. These estimates for usable capacities are conservative lower bounds since they do not account for dispersion stabilized H₂ within the MOF pores.

7 Conclusions and outlook

Decarbonization of the world economy is an urgent challenge to address the increasingly alarming rate of climate change arising from greenhouse gas emission from legacy hydrocarbon-based transportation, heating and electricity generation applications. Hydrogen has tremendous promise as a clean fuel, but this promise is inhibited by the challenge of economically transporting it from the point of generation to the point of application, dispensing it, and even storing it from time of generation to time of consumption: this is the hydrogen storage problem. One approach to the as-yet not practically solved H₂ storage problem is the use of sorbents, such as MOFs.

In this perspective article, we have discussed computational quantum chemistry modeling of H₂ binding in MOFs, across topics ranging from thermodynamic generalities to specific strengths and weaknesses of our chosen simulation methods to a range of comparisons between theory and experiment and a range of predictions. To summarize the generalities first, MOF materials design focuses on achieving maximal density of H₂ binding sites that individually possess near-optimal thermodynamic characteristics for ambient temperature storage. Of course that is subject to suitable materials robustness, cost constraints, and synthetic realizability, which are topics we do not dwell upon in this computational perspective. In principle, optimal binding at 300 K should correspond to ΔG ∼ +8 kJ mol⁻¹ to optimize usable capacity for a 5–100 bar pressure swing, or ΔG ∼ +9 kJ mol⁻¹ to optimize usable capacity for a 5–170 bar swing. Given at least partial enthalpy–entropy compensation, an optimal range for binding enthalpies emerges that is between roughly −15 and −25 kJ mol⁻¹.

Computational quantum chemistry permits, in principle, near-exact solution of the potential energy surface based on well-defined approximations to quantum mechanics. In practice, of course, there are a range of tradeoffs between the accuracy of results and the feasibility of calculations for a given system size. We have summarized key advances in the accuracy associated with DFT calculations on cluster models of H₂ binding sites of MOFs. From the perspective of H₂ binding, the best modern density functionals on rungs 3 (meta-GGAs) and 4 (hybrids) of the Jacob's ladder classification are capable of results that are typically within a few kJ mol⁻¹ of the true result for the electronic binding energy. It is worth remembering that the worst are much poorer! We also reviewed and employed the latest EDA to unravel the main driving forces behind H₂ binding. EDA permits the separation of the electronic binding energy into a frozen contribution containing Pauli repulsions and permanent electrostatics, a dispersion contribution, polarization (induced electrostatics/relief of Pauli repulsions), and charge-transfer (donor–acceptor) contributions. The latter three contributions are purely attractive.

About simulation accuracy, the situation is considerably less satisfactory for the corrections due to the vibrational motion that is necessary to obtain thermodynamic functions such as the enthalpy, entropy and free energy changes associated with binding. The standard rigid–rotor harmonic oscillator (RRHO) approach can be ineffective, and it is necessary to carefully consider that bound H₂ retains motional freedom beyond that captured in the RRHO. For instance, in line with earlier work by Sauer and colleagues,⁹⁹ we assumed free rotational movement of H₂ when bound to the cup site in MOF-5, and to Ni(II) in Ni₂(m-dobdc). The RRHO model also has known limitations when dealing with delocalized soft modes within the framework. These complexities remind us that each framework possesses unique vibrational characteristics that require customized computational approaches. They also point to a compelling avenue for future methodological refinements.

Improved computational predictions of free energies of hydrogen binding in metal–organic frameworks are challenging due to the anharmonic vibrations and coupling between multiple vibrational modes of the material; they are known to play a role in MOFs.^143,144 An in-principle exact approach comes from the path integral formulation of quantum mechanics,¹⁴⁵ which leads to well-developed path-integral methods¹⁴⁶ that require either molecular dynamics or Monte Carlo sampling. While path integrals are the right ultimate goal, their expense is prohibitive, and more approximate methods are necessary, until automatic development of cheap force fields that accurately approximate our expensive quantum mechanical calculations becomes possible, or computing power advances a few more orders of magnitude. Beyond the RRHO model, and the simple corrections we have employed are a range of approaches that should offer higher accuracy than RRHO with cost still far less than path integrals.¹⁴⁷ More work in this area would be beneficial to enable routine treatment of the beyond-harmonic quantum aspects of nuclear motion such as anharmonic couplings between MOF phonons and bound H₂ modes, which is crucial for precise predictions of free energies of hydrogen binding.

As validation examples, as well as specific cases of H₂ binding in MOFs to be analyzed and understood, we report calculations examining MOF-5, Ni₂(m-dobdc), Cu(I)-MFU-4l, and V₂Cl_2.8(btdd). MOF-5 represents the extreme of overly weak binding, as it displays exceptional hydrogen storage capabilities only at low temperatures because H₂ molecules are attracted to the framework via van der Waals forces. The other examples involve open metal sites, demonstrating the importance of orbital synergy in H₂ binding within these frameworks. Specifically, V₂Cl_2.8(btdd) exhibits slightly stronger charge transfer interactions than Ni₂(m-dobdc), and slightly attenuated back-bonding compared to Cu(I)-MFU-4l, resulting in near-optimal binding for storage at ambient conditions. V₂Cl_2.8(btdd) is, however, less promising for deployment than Ni₂(m-dobdc) because of lower usable capacity due to its lower site density in addition to air sensitivity.

Densification of H₂ with MOFs could be greatly enhanced by the binding of more than one H₂ at a specific open metal site. We have reviewed and extended existing calculations that probe the candidate sites for binding multiple H₂ molecules, starting with metal catecholates, which present very exposed metal sites which are capable of binding either 2 (Mg-catecholate) or 4 (Ca-catecholate) H₂ with binding affinities that are stronger than physisorption though likely weaker than optimal. Turning to metals, the MFU-4l and MFU-4 family of materials offer a unique tetrahedral coordination scaffold stabilising transition metal sites. If a suitable metal ion replaces Zn, our calculations show that this allows for the binding of two equivalents of H₂ in the case of divalent V(II) and three equivalents in monovalent Sc(I) or Ti(I).

As already discussed, the Sc(I) and Ti(I) oxidation states are exotic, and it may therefore not be possible to achieve post-synthetic modification of the MFU-4l and MFU-4 frameworks to create these novel low-valent sites. We do note that MFU-4l is known to support Ti(IV) and Ti(III),¹⁴⁸ which can serve as a high-valent starting point for attempting chemical reduction for instance using sodium naphthalene. Regardless, the binding of three H₂ at, e.g. Sc(I)-MFU-4l (or MFU-4) is striking for the density with which the H₂ molecules pack with near-optimal binding enthalpies. Our optimized structures reveal many contacts inside the nominal van der Waals contact distance in the ligated structures, as well as an estimate of usable capacity of ∼25 g L⁻¹ at ambient conditions. Both these aspects point to exciting potential for future discoveries which yield similar binding features and packing densities in synthetically realizable MOFs.

Author contributions

RC and MHG conceptualized the manuscript and formulated the methodology for this work. RC curated the data and wrote the original draft. RC, JT, HS and KMC visualized the data. HS implemented and provided some of the software code for the simulations. HZHJ provided insight into optimal binding conditions for hydrogen storage. MHG substantially edited the initial drafts to create the article that was then reviewed by JT, HS, YY, KMC, HF, and JRL.

Conflicts of interest

MHG is a part-owner of Q-Chem Inc., whose software was used for the calculations reported here.

Acknowledgements

The authors gratefully acknowledge research support from the Hydrogen Materials Advanced Research Consortium (HyMARC), established as part of the Energy Materials Network under the U.S. DOE, Office of Energy Efficiency and Renewable Energy, under Contract No. DE-AC02-05CH11231. Work on energy decomposition analysis was supported by the U.S. National Science Foundation through Grants No. CHE-1955643 and CHE-2313791. Work on vibrational contributions to thermodynamic functions was supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. DOE under contract no. DE-AC02-05CH11231, through the Gas Phase Chemical Physics Program. KMC acknowledges the Arnold O. Beckman Foundation for financial support.

Notes and references

1. E. L. Miller, S. T. Thompson, K. Randolph, Z. Hulvey, N. Rustagi and S. Satyapal, MRS Bull., 2020, 45, 57–64.
2. A. M. Oliveira, R. R. Beswick and Y. Yan, Curr. Opin. Chem. Eng., 2021, 33, 100701.
3. G. W. Crabtree, M. S. Dresselhaus and M. V. Buchanan, Phys. Today, 2004, 57, 39–44.
4. B. C. Tashie-Lewis and S. G. Nnabuife, Chem. Eng. J. Adv., 2021, 8, 100172.
5. S. K. Kar, S. Harichandan and B. Roy, Int. J. Hydrogen Energy, 2022, 47, 10803–10824.
6. H.-M. Cheng, Q.-H. Yang and C. Liu, Carbon, 2001, 39, 1447–1454.
7. I. Jain, C. Lal and A. Jain, Int. J. Hydrogen Energy, 2010, 35, 5133–5144.
8. D. P. Broom, Hydrogen storage materials: the characterisation of their storage properties, Springer, 2011, vol. 1.
9. Y.-h Zhang, Z.-c Jia, Z.-m Yuan, T. Yang, Y. Qi and D.-l Zhao, J. Iron Steel Res. Int., 2015, 22, 757–770.
10. Q. Lai, M. Paskevicius, D. A. Sheppard, C. E. Buckley, A. W. Thornton, M. R. Hill, Q. Gu, J. Mao, Z. Huang, H. K. Liu, Z. Guo, A. Banerjee, S. Chakraborty, R. Ahuja and K. F. Aguey-Zinsou, ChemSusChem, 2015, 8, 2789–2825.
11. E. Rivard, M. Trudeau and K. Zaghib, Materials, 2019, 12, 1973.
12. J. O. Abe, A. Popoola, E. Ajenifuja and O. M. Popoola, Int. J. Hydrogen Energy, 2019, 44, 15072–15086.
13. M. R. Usman, Renewable Sustainable Energy Rev., 2022, 167, 112743.
14. D. Wei, X. Shi, R. Qu, K. Junge, H. Junge and M. Beller, ACS Energy Lett., 2022, 7, 3734–3752.
15. M. D. Allendorf, Z. Hulvey, T. Gennett, A. Ahmed, T. Autrey, J. Camp, E. S. Cho, H. Furukawa, M. Haranczyk and M. Head-Gordon, et al. , Energy Environ. Sci., 2018, 11, 2784–2812.
16. M. D. Allendorf, V. Stavila, J. L. Snider, M. Witman, M. E. Bowden, K. Brooks, B. L. Tran and T. Autrey, Nat. Chem., 2022, 14, 1214–1223.
17. M. Wieliczko and N. Stetson, MRS Energy Sustainability, 2020, 7, E41.
18. H. Barthélémy, Int. J. Hydrogen Energy, 2012, 37, 17364–17372.
19. C. Chilev and F. D. Lamari, Int. J. Hydrogen Energy, 2016, 41, 1744–1758.
20. A. Schneemann, J. L. White, S. Kang, S. Jeong, L. F. Wan, E. S. Cho, T. W. Heo, D. Prendergast, J. J. Urban, B. C. Wood, M. D. Allendorf and V. Stavila, Chem. Rev., 2018, 118, 10775–10839.
21. A. Schneemann, L. F. Wan, A. S. Lipton, Y. S. Liu, J. L. Snider, A. A. Baker, J. D. Sugar, C. D. Spataru, J. Guo, T. S. Autrey, M. Jørgensen, T. R. Jensen, B. C. Wood, M. D. Allendorf and V. Stavila, ACS Nano, 2020, 14, 10294–10304.
22. M. P. Suh, H. J. Park, T. K. Prasad and D.-W. Lim, Chem. Rev., 2012, 112, 782–835.
23. L. J. Murray, M. Dincă and J. R. Long, Chem. Soc. Rev., 2009, 38, 1294–1314.
24. A. G. Wong-Foy, A. J. Matzger and O. M. Yaghi, J. Am. Chem. Soc., 2006, 128, 3494–3495.
25. M. T. Kapelewski, T. Runčevski, J. D. Tarver, H. Z. Jiang, K. E. Hurst, P. A. Parilla, A. Ayala, T. Gennett, S. A. Fitzgerald, C. M. Brown and J. R. Long, Chem. Mater., 2018, 30, 8179–8189.
26. J. Purewal, M. Veenstra, D. Tamburello, A. Ahmed, A. J. Matzger, A. G. Wong-Foy, S. Seth, Y. Liu and D. J. Siegel, Int. J. Hydrogen Energy, 2019, 44, 15135–15145.
27. D. E. Jaramillo, H. Z. Jiang, H. A. Evans, R. Chakraborty, H. Furukawa, C. M. Brown, M. Head-Gordon and J. R. Long, J. Am. Chem. Soc., 2021, 143, 6248–6256.
28. H. Nazir, N. Muthuswamy, C. Louis, S. Jose, J. Prakash, M. E. Buan, C. Flox, S. Chavan, X. Shi and P. Kauranen, et al. , Int. J. Hydrogen Energy, 2020, 45, 20693–20708.
29. P. Peng, A. Anastasopoulou, K. Brooks, H. Furukawa, M. E. Bowden, J. R. Long, T. Autrey and H. Breunig, Nat. Energy, 2022, 7, 448–458.
30. U. D. of Energy, DOE Technical Targets for Onboard Hydrogen Storage for Light-Duty Vehicles, https://www.energy.gov/eere/fuelcells/doe-technica-targets-onboard-hydrogen-storage-light-duty-vehicles, 2017, Accessed on: October 2023.
31. U.S. Department of Energy under the Biden-Harris Administration, US releases National Clean Hydrogen Strategy and Roadmap, MRS Bull., 2023, 48, 707–708 DOI:10.1557/s43577-023-00574-9.
32. H. Furukawa, K. E. Cordova, M. OKeeffe and O. M. Yaghi, Science, 2013, 341, 1230444.
33. H. B. Wu and X. W. Lou, Sci. Adv., 2017, 3, eaap9252.
34. O. M. Yaghi, M. O'Keeffe, N. W. Ockwig, H. K. Chae, M. Eddaoudi and J. Kim, Nature, 2003, 423, 705–714.
35. R.-B. Lin, S. Xiang, H. Xing, W. Zhou and B. Chen, Coord. Chem. Rev., 2019, 378, 87–103.
36. S. Rojas and P. Horcajada, Chem. Rev., 2020, 120, 8378–8415.
37. Q. Wang and D. Astruc, Chem. Rev., 2020, 120, 1438–1511.
38. A. Bavykina, N. Kolobov, I. S. Khan, J. A. Bau, A. Ramirez and J. Gascon, Chem. Rev., 2020, 120, 8468–8535.
39. J. Y. Saillard and R. Hoffmann, J. Am. Chem. Soc., 1984, 106, 2006–2026.
40. R. C. Lochan and M. Head-Gordon, Phys. Chem. Chem. Phys., 2006, 8, 1357–1370.
41. P. R. Kemper, P. Weis, M. T. Bowers and P. Matre, J. Am. Chem. Soc., 1998, 120, 13494–13502.
42. E. Tsivion, J. R. Long and M. Head-Gordon, J. Am. Chem. Soc., 2014, 136, 17827–17835.
43. D. Zhao, X. Wang, L. Yue, Y. He and B. Chen, Chem. Commun., 2022, 58, 11059–11078.
44. M. Dincă and J. R. Long, Angew. Chem., Int. Ed., 2008, 47, 6766–6779.
45. Ü. Kökçam-Demir, A. Goldman, L. Esrafili, M. Gharib, A. Morsali, O. Weingart and C. Janiak, Chem. Soc. Rev., 2020, 49, 2751–2798.
46. S. A. A. Razavi and A. Morsali, Coord. Chem. Rev., 2019, 399, 213023.
47. X. Zou, M.-H. Cha, S. Kim, M. C. Nguyen, G. Zhou, W. Duan and J. Ihm, Int. J. Hydrogen Energy, 2010, 35, 198–203.
48. M. Babucci, A. Guntida and B. C. Gates, Chem. Rev., 2020, 120, 11956–11985.
49. S. K. Bhatia and A. L. Myers, Langmuir, 2006, 22, 1688–1700.
50. Y.-S. Bae and R. Q. Snurr, Microporous Mesoporous Mater., 2010, 132, 300–303.
51. C. M. Simon, J. Kim, L.-C. Lin, R. L. Martin, M. Haranczyk and B. Smit, Phys. Chem. Chem. Phys., 2014, 16, 5499–5513.
52. S. P. Veccham and M. Head-Gordon, J. Chem. Theory Comput., 2020, 16, 4963–4982.
53. S. P. Veccham and M. Head-Gordon, J. Phys. Chem. A, 2021, 125, 4245–4257.
54. P. R. Horn, Y. Mao and M. Head-Gordon, Phys. Chem. Chem. Phys., 2016, 18, 23067–23079.
55. Y. Mao, M. Loipersberger, P. R. Horn, A. Das, O. Demerdash, D. S. Levine, S. Prasad Veccham, T. Head-Gordon and M. Head-Gordon, Ann. Rev. Phys. Chem., 2021, 72, 641–666.
56. N. L. Rosi, J. Eckert, M. Eddaoudi, D. T. Vodak, J. Kim, M. O'Keeffe and O. M. Yaghi, Science, 2003, 300, 1127–1129.
57. M. T. Kapelewski, S. J. Geier, M. R. Hudson, D. Stück, J. A. Mason, J. N. Nelson, D. J. Xiao, Z. Hulvey, E. Gilmour, S. A. Fitzgerald, M. Head-Gordon, C. M. Brown and J. R. Long, J. Am. Chem. Soc., 2014, 136, 12119–12129.
58. D. Denysenko, M. Grzywa, J. Jelic, K. Reuter and D. Volkmer, Angew. Chem., Int. Ed., 2014, 53, 5832–5836.
59. B. R. Barnett, H. A. Evans, G. M. Su, H. Z. H. Jiang, R. Chakraborty, D. Banyeretse, T. J. Hartman, M. B. Martinez, B. A. Trump, J. D. Tarver, M. N. Dods, L. M. Funke, J. Börgel, J. A. Reimer, W. S. Drisdell, K. E. Hurst, T. Gennett, S. A. FitzGerald, C. M. Brown, M. Head-Gordon and J. R. Long, J. Am. Chem. Soc., 2021, 143, 14884–14894.
60. T. Runčevski, M. T. Kapelewski, R. M. Torres-Gavosto, J. D. Tarver, C. M. Brown and J. R. Long, ChemComm, 2016, 52, 8251–8254.
61. T. Sagara, J. Klassen, J. Ortony and E. Ganz, J. Chem. Phys., 2005, 123, 014701.
62. E. Klontzas, A. Mavrandonakis, E. Tylianakis and G. E. Froudakis, Nano Lett., 2008, 8, 1572–1576.
63. A. Mavrandonakis, E. Klontzas, E. Tylianakis and G. E. Froudakis, J. Am. Chem. Soc., 2009, 131, 13410–13414.
64. E. Klontzas, E. Tylianakis and G. E. Froudakis, J. Phys. Chem. C, 2009, 113, 21253–21257.
65. T. Stergiannakos, E. Tylianakis, E. Klontzas and G. E. Froudakis, J. Phys. Chem. C, 2010, 114, 16855–16858.
66. R. B. Getman, J. H. Miller, K. Wang and R. Q. Snurr, J. Phys. Chem. C, 2011, 115, 2066–2075.
67. E. Klontzas, E. Tylianakis and G. E. Froudakis, J. Phys. Chem. Lett., 2011, 2, 1824–1830.
68. E. Tsivion, S. P. Veccham and M. Head-Gordon, ChemPhysChem, 2017, 18, 184–188.
69. H. Tachikawa and T. Iyama, J. Phys. Chem. C, 2019, 123, 8709–8716.
70. R. Chakraborty, K. M. Carsch, D. E. Jaramillo, Y. Yabuuchi, H. Furukawa, J. R. Long and M. Head-Gordon, J. Phys. Chem. Lett., 2022, 13, 10471–10478.
71. Y. S. Bae and R. Q. Snurr, Microporous Mesoporous Mater., 2010, 132, 300–303.
72. E. Garrone, B. Bonelli and C. Otero Areán, Chem. Phys. Lett., 2008, 456, 68–70.
73. S. M. Chavan, O. Zavorotynska, C. Lamberti and S. Bordiga, Dalton Trans., 2013, 42, 12586–12595.
74. E. Tsivion and M. Head-Gordon, J. Phys. Chem. C, 2017, 121, 12091–12100.
75. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
76. N. Mardirossian and M. Head-Gordon, Mol. Phys., 2017, 115, 2315–2372.
77. J. Lee, A. Rettig, X. Feng, E. Epifanovsky and M. Head-Gordon, J. Chem. Theory Comput., 2022, 18, 7336–7349.
78. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184–32215.
79. N. Mardirossian and M. Head-Gordon, J. Chem. Phys., 2016, 144, 214110.
80. O. A. Vydrov and T. Van Voorhis, J. Chem. Phys., 2010, 133, 244103.
81. A. Najibi and L. Goerigk, J. Chem. Theory Comput., 2018, 14, 5725–5738.
82. G. Santra and J. M. Martin, Some observations on the performance of the most recent exchange-correlation functionals for the large and chemically diverse GMTKN55 benchmark, AIP Conf Proc., 2019, 2186, 030004.
83. B. Chan, P. M. W. Gill and M. Kimura, J. Chem. Theory Comput., 2019, 15, 3610–3622.
84. Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2008, 4, 1849–1868.
85. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
86. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
87. R. F. Ribeiro, A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2011, 115, 14556–14562.
88. E. Epifanovsky, A. T. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons and A. L. Dempwolff, et al. , J. Chem. Phys., 2021, 155, 084801.
89. R. Z. Khaliullin, E. A. Cobar, R. C. Lochan, A. T. Bell and M. Head-Gordon, J. Phys. Chem. A, 2007, 111, 8753–8765.
90. D. S. Levine, P. R. Horn, Y. Mao and M. Head-Gordon, J. Chem. Theory Comput., 2016, 12, 4812–4820.
91. D. S. Levine and M. Head-Gordon, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 12649–12656.
92. P. R. Horn, Y. Mao and M. Head-Gordon, J. Chem. Phys., 2016, 144, 114107.
93. P. R. Horn and M. Head-Gordon, J. Chem. Phys., 2015, 143, 114111.
94. N. L. Rosi, J. Eckert, M. Eddaoudi, D. T. Vodak, J. Kim, M. O'Keeffe and O. M. Yaghi, Science, 2003, 300, 1127–1129.
95. J. L. Rowsell and O. M. Yaghi, J. Am. Chem. Soc., 2006, 128, 1304–1315.
96. D. Dubbeldam, K. S. Walton, D. E. Ellis and R. Q. Snurr, Angew. Chem., Int. Ed., 2007, 46, 4496–4499.
97. K. Sumida, D. L. Rogow, J. A. Mason, T. M. McDonald, E. D. Bloch, Z. R. Herm, T.-H. Bae and J. R. Long, Chem. Rev., 2012, 112, 724–781.
98. J. L. Rowsell, E. C. Spencer, J. Eckert, J. A. Howard and O. M. Yaghi, Science, 2005, 309, 1350–1354.
99. K. Sillar, A. Hofmann and J. Sauer, J. Am. Chem. Soc., 2009, 131, 4143–4150.
100. F. M. Mulder, T. J. Dingemans, H. G. Schimmel, A. J. Ramirez-Cuesta and G. J. Kearley, Chem. Phys., 2008, 351, 72–76.
101. S. Bordiga, J. G. Vitillo, G. Ricchiardi, L. Regli, D. Cocina, A. Zecchina, B. Arstad, M. Bjørgen, J. Hafizovic and K. P. Lillerud, J. Phys. Chem. B, 2005, 109, 18237–18242.
102. H. Deng, S. Grunder, K. E. Cordova, C. Valente, H. Furukawa, M. Hmadeh, F. Gándara, A. C. Whalley, Z. Liu and S. Asahina, et al. , Science, 2012, 336, 1018–1023.
103. N. L. Rosi, J. Kim, M. Eddaoudi, B. Chen, M. O'Keeffe and O. M. Yaghi, J. Am. Chem. Soc., 2005, 127, 1504–1518.
104. P. D. Dietzel, Y. Morita, R. Blom and H. Fjellvåg, Angew. Chem., 2005, 117, 6512–6516.
105. T. Xiao and D. Liu, Microporous Mesoporous Mater., 2019, 283, 88–103.
106. S. R. Caskey, A. G. Wong-Foy and A. J. Matzger, J. Am. Chem. Soc., 2008, 130, 10870–10871.
107. A. F. Cozzolino, C. K. Brozek, R. D. Palmer, J. Yano, M. Li and M. DincaÌ, J. Am. Chem. Soc., 2014, 136, 3334–3337.
108. K. Lee, W. C. Isley, A. L. Dzubak, P. Verma, S. J. Stoneburner, L. C. Lin, J. D. Howe, E. D. Bloch, D. A. Reed, M. R. Hudson, C. M. Brown, J. R. Long, J. B. Neaton, B. Smit, C. J. Cramer, D. G. Truhlar and L. Gagliardi, J. Am. Chem. Soc., 2014, 136, 698–704.
109. E. Haldoupis, J. Borycz, H. Shi, K. D. Vogiatzis, P. Bai, W. L. Queen, L. Gagliardi and J. I. Siepmann, J. Phys. Chem. C, 2015, 119, 16058–16071.
110. K. Lee, J. D. Howe, L. C. Lin, B. Smit and J. B. Neaton, Chem. Mater., 2015, 27, 668–678.
111. A. L. Dzubak, L.-C. Lin, J. Kim, J. A. Swisher, R. Poloni, S. N. Maximoff, B. Smit and L. Gagliardi, Nat. Chem., 2012, 4, 810–816.
112. A. Kundu, K. Sillar and J. Sauer, Chem. Sci., 2020, 11, 643–655.
113. S. A. FitzGerald, D. Mukasa, K. H. Rigdon, N. Zhang and B. R. Barnett, J. Phys. Chem. C, 2019, 123, 30427–30433.
114. G. M. Su, H. Wang, B. R. Barnett, J. R. Long, D. Prendergast and W. S. Drisdell, Chem. Sci., 2021, 12, 2156–2164.
115. S. Biswas, M. Tonigold, M. Speldrich, P. Kögerler, M. Weil and D. Volkmer, Inorg. Chem., 2010, 49, 7424–7434.
116. D. Sengupta, P. Melix, S. Bose, J. Duncan, X. Wang, M. R. Mian, K. O. Kirlikovali, F. Joodaki, T. Islamoglu and T. Yildirim, et al. , J. Am. Chem. Soc., 2023, 145, 20492–20502.
117. D. E. Jaramillo, D. A. Reed, H. Z. Jiang, J. Oktawiec, M. W. Mara, A. C. Forse, D. J. Lussier, R. A. Murphy, M. Cunningham, V. Colombo, D. K. Shuh, J. A. Reimer and J. R. Long, Nat. Mater., 2020, 19, 517–521.
118. M. J. Sever and J. J. Wilker, Dalton Trans., 2004, 1061–1072.
119. H. Fei, J. Shin, Y. S. Meng, M. Adelhardt, J. Sutter, K. Meyer and S. M. Cohen, J. Am. Chem. Soc., 2014, 136, 4965–4973.
120. M. Hmadeh, Z. Lu, Z. Liu, F. Gándara, H. Furukawa, S. Wan, V. Augustyn, R. Chang, L. Liao and F. Zhou, et al. , Chem. Mater., 2012, 24, 3511–3513.
121. D. Rankine, T. D. Keene, C. J. Doonan and C. J. Sumby, Cryst. Growth Des., 2014, 14, 5710–5718.
122. S. J. Stoneburner, V. Livermore, M. E. McGreal, D. Yu, K. D. Vogiatzis, R. Q. Snurr and L. Gagliardi, J. Phys. Chem. C, 2017, 121, 10463–10469.
123. H. Chen and R. Q. Snurr, J. Phys. Chem. C, 2021, 125, 21701–21708.
124. S. J. Stoneburner and L. Gagliardi, J. Phys. Chem. C, 2018, 122, 22345–22351.
125. S. De, T. Devic and A. Fateeva, Dalton Trans., 2021, 50, 1166–1188.
126. D. W. Smithenry, S. R. Wilson and K. S. Suslick, Inorg. Chem., 2003, 42, 7719–7721.
127. K. S. Suslick, P. Bhyrappa, J.-H. Chou, M. E. Kosal, S. Nakagaki, D. W. Smithenry and S. R. Wilson, Acc. Chem. Res., 2005, 38, 283–291.
128. E.-Y. Choi, P. M. Barron, R. W. Novotney, C. Hu, Y.-U. Kwon and W. Choe, CrystEngComm, 2008, 10, 824–826.
129. L. Bonomo, M.-L. Lehaire, E. Solari, R. Scopelliti and C. Floriani, Angew. Chem., Int. Ed., 2001, 40, 771–774.
130. S. Muniappan, S. Lipstman, S. George and I. Goldberg, Inorg. Chem., 2007, 46, 5544–5554.
131. J. Wang, Y. Fan, Y. Tan, X. Zhao, Y. Zhang, C. Cheng and M. Yang, ACS Appl. Mater. Interfaces, 2018, 10, 36615–36621.
132. D. Feng, H.-L. Jiang, Y.-P. Chen, Z.-Y. Gu, Z. Wei and H.-C. Zhou, Inorg. Chem., 2013, 52, 12661–12667.
133. J. Ryou, G. Kim and S. Hong, J. Chem. Phys., 2011, 134, 234701.
134. G. J. Kubas, Chem. Rev., 2007, 107, 4152–4205.
135. R. C. Lochan, R. Z. Khaliullin and M. Head-Gordon, Inorg. Chem., 2008, 47, 4032–4044.
136. R. J. Comito, Z. Wu, G. Zhang, J. A. Lawrence, M. D. Korzyński, J. A. Kehl, J. T. Miller and M. Dincă, Angew. Chem., Int. Ed., 2018, 57, 8135–8139.
137. H. Shen, S. P. Veccham and M. Head-Gordon, J. Chem. Theory Comput., 2023, 19, 8624–8638.
138. P. L. Arnold, F. G. N. Cloke, P. B. Hitchcock and J. F. Nixon, J. Am. Chem. Soc., 1996, 118, 7630–7631.
139. F. Calderazzo, I. Ferri, G. Pampaloni, U. Englert and M. L. Green, Organometallics, 1997, 16, 3100–3101.
140. J. J. Talbot, R. Chakraborty, H. Shen and M. Head-Gordon, J. Phys. Chem. Lett., 2023, 14, 5432–5440.
141. K. A. Dill, S. Bromberg and D. Stigter, Molecular driving forces: statistical thermodynamics in biology, chemistry, physics, and nanoscience, Garland Science, 2010, p. 563.
142. H. Frost and R. Q. Snurr, J. Phys. Chem. C, 2007, 111, 18794–18803.
143. A. Lamaire, J. Wieme, S. M. Rogge, M. Waroquier and V. Van Speybroeck, J. Chem. Phys., 2019, 150, 094503.
144. V. Kapil, J. Wieme, S. Vandenbrande, A. Lamaire, V. Van Speybroeck and M. Ceriotti, J. Chem. Theory Comput., 2019, 15, 3237–3249.
145. D. Chandler and P. G. Wolynes, J. Chem. Phys., 1981, 74, 4078–4095.
146. C. P. Herrero and R. Ramirez, J. Phys.: Condens. Matter, 2014, 26, 233201.
147. V. Kapil, E. Engel, M. Rossi and M. Ceriotti, J. Chem. Theory Comput., 2019, 15, 5845–5857.
148. R. J. Comito, K. J. Fritzsching, B. J. Sundell, K. Schmidt-Rohr and M. Dincă, J. Am. Chem. Soc., 2016, 138, 10232–10237.

---

Footnote

† Electronic supplementary information (ESI) available: Coordinates for all converged geometries. See DOI: https://doi.org/10.1039/d3cp05540j

---