Hydrogen adsorption mechanism of MOF-74 metal–organic frameworks: an insight from first principles calculations

Abstract

The microscopic mechanism of the H₂ adsorption of two Mg-MOF-74 isoreticular frameworks, one with a benzenedicarboxylate (BDC) linker and the other with a dihydroxyfumarate (DHF) linker, were studied on the basis of density functional theory (DFT) method. Possible adsorption sites on the internal surface of the two MOFs were detected using ab initio molecular dynamics (AIMD) annealing simulations. The simulations were able to reproduce all adsorption sites which have been experimentally observed for the BDC-based M-MOF-74 frameworks with M = Ni and Zn. In descending order of binding strengths, they are the adsorption sites primarily induced by the open metal sites P₁, the oxygen atoms of the oxido groups P₂ and the aromatic rings P₃. The H₂–framework binding strengths were properly evaluated by taking into account the vibrational zero-point energy (ZPE) contribution. An additional type of adsorption sites induced by the oxygen atoms of the carboxyl groups P₄ is predicted for the Mg-MOF-74 framework. Two types of adsorption sites primarily induced by the open metal sites P₁ and oxygen atoms of the carboxyl groups P₂ were predicted for the DHF-based Mg-MOF-74 framework. Detailed analysis of the electron density showed that the electrostatic interaction of the H₂ molecule with the charge distribution of the local framework environment within a radius of ∼3.5 Å is a key factor to define adsorption positions and binding strength. The absence of the P₄ sites in the BDC-based Zn-MOF-74 framework is caused by the lower charge density at the oxygen atoms induced by less electro-positive metal. The substitution of the nonaromatic DHF linker for the aromatic BDC linker reduces the binding strength at the metal induced adsorption sites by 1.45 kJ mol⁻¹ due to the absence of the aromatic ring.

---

1. Introduction

In recent decades, metal–organic frameworks (MOFs) have emerged as an advanced generation of nanoporous materials.^1–3 They have crystalline structures that are composed of metal clusters linked to each other by organic molecules to form coordination networks. The large variation in chemical compositions and conformations of both metal-cluster nodes and organic linkers allows a great flexibility in designing and adjusting MOFs to achieve desired internal surface area, void volume as well as other physicochemical properties.² Consequently, designing MOFs for a wide range of applications such as catalysis,⁴ gas storage,⁵ gas sensors,⁶ and biomedicine⁷ etc. is a very active research area.

Numerous research studies have been devoted to designing MOFs for solid-state hydrogen storage applications.^8,9 One of the big challenges is the weak interaction between non-polar covalent hydrogen molecules and MOFs which typically can maintain high hydrogen uptake only at extreme conditions of low temperature and high pressure.⁸ Using Grand Canonical Monte Carlo (GCMC) simulations of hydrogen adsorption in the IRMOF series, Frost and Snurr suggested that a reasonable H₂ uptake at ambient condition could be obtained with an isosteric heat of adsorption of 15 kJ mol⁻¹ (0.155 eV).¹⁰ Particular interests have been paid for MOFs with open metal sites which can attract H₂ molecules via stronger electrostatic interactions. For example, the low loading heat of adsorption in Co-based SNU-15 MOF with open Co(II) sites is up to 15 kJ mol⁻¹.¹¹ However, its hydrogen uptake at T = 77 K and P = 1 bar is very low 0.74 wt% due to the low metal site density 0.0012 sites per ų.

Considerable attention has been paid for MOF-74 series since they have higher open metal site density 0.0045 sites per ų which corresponds to hydrogen uptake 2.5 wt% at T = 77 K and P = 1 bar.^12–15 The experimental values of low loading heat of H₂ adsorption of these materials are between 8 and 13 kJ mol⁻¹.¹³ GCMC simulations using classical force field demonstrated that the primary adsorption sites which are responsible for such loading heat of adsorption are induced by the electrostatic interaction between the positively charged open metal sites and the induced dipole moment of H₂ molecules.¹² This interaction was demonstrated to be strongly correlated with the ionic charge and the polarizability of the metal. The role of the linker was also revealed with the combination of experimental measurements and theoretical calculations. Kapelewski et al. were able to increase the low loading heat of adsorption by 1.5 kJ mol⁻¹ with a substitution of the original linker 2,5-dioxido-1,4-benzenedicarboxylate (BDC) using an isomer 4,6-dioxido-1,3-benzenedicarboxylate (mBDC).¹⁶ Their electronic structure calculations based on the density functional theory (DFT) suggest that the extra contribution of mBDC linker enhances the metal site–H₂ interaction. Moreover, the neutron powder diffraction (NPD) experiments evidence that the change from para (BDC linker) to meta (mBDC linker) configurations of the carboxylic acid functionalities and hydroxyl substituents of the linker lead to the change in secondary adsorption sites.^16,17 The hydrogen uptake of mBDC-Co/Ni-MOF-74 frameworks at P = 100 bar, T = 25 °C is up to 13 g L⁻¹ and 1.2 wt%, slightly higher than their BDC-based counterparts. These frameworks were believed to be among the best of hydrogen physisorptive adsorbents.¹⁷ However, they are still much lower than the targets of the US Department of Energy for practical fuel cell which require a volumetric uptake at least 50 g L⁻¹ in the temperature range −40 ÷ 60 °C at maximum delivery pressure 12 bar.¹⁸ Based on classical GCMC simulations, Witman and coworkers suggested that substituting the aromatic linker with a shorter aliphatic linker, 2,3-dihydroxyfumarate (DHF) can strongly enhance the volumetric uptake owing to the higher open-metal site density.¹⁹

In order to improve MOFs for practical hydrogen storage applications, it is crucial to understand the underlying mechanism of the H₂ adsorption. Various adsorption sites on the internal surface of MOF-74 frameworks have been reported by experimental works.^16,20,21 However, previous theoretical calculations primarily focused on the role of the open metal sites.^12,13 In this work, basing on ab initio molecular dynamics (AIMD) simulations, all possible adsorption sites including both metal-induced and linker-induced sites were predicted for the original BDC-based Mg-MOF-74 framework and its isoreticular DHF-based framework. Detail electronic structure analysis demonstrates that changing the electro-positivity of the metal affects the stability of linker-induced adsorption sites whereas the linker substitution gives significant influences on the binding strength of metal-induced adsorption sites.

2. Computational details

The crystal structures of MOF-74 family are based on rhombohedral lattice (Fig. 1a and b). Both BDC- and DHF-based frameworks belong to R3̄ space group. Their crystal structures were fully optimized on the base of static DFT calculations. There are 3 equivalent 1-dimensional (1D) void channels along c direction within a unit cell of each framework. Thus, adsorption sites in one channel are assumed to be the same as those in the others so that we only need to probe one channel.

Fig. 1 (a) and (b) represent the crystal structures of the BDC-Mg-MOF-74 and DHF-Mg-MOF-74 frameworks, respectively. Red, brown, cyan and white spheres represent oxygen, magnesium, carbon and hydrogen atoms, respectively. The insets are schematic representations of corresponding linkers, i.e. (–O)₂C₆H₂–1,4-(COO–)₂ linker derived from 2,5-dioxido-1,4-benzenedicarboxylic acid (BDC) in (a) and (–COO)₂C=C–2,3-(O–)₂ linker derived from 2,3-dihydroxyfumaric acid (DHF) in (b). Dash lines in the insets denote dangling bonds which can be attached to metal atoms. The carbon atoms of the aromatic ring are labeled C1 and the carbon atoms of the carboxyl groups are labeled C2. The oxygen atoms of the oxido groups are labeled O1 and the oxygen atoms of the carboxyl groups are labeled O2. (c) A schematic representation of the largest included sphere (LIS) and the largest free sphere (LFS). The red circle, green circle and black shaded area denote the LIS, LFS and the space filled by the framework atoms, respectively.

Molecular dynamics simulated annealing method has been adopted to locate minima of the potential energy surface of adsorption systems.^22–24 In this work, n hydrogen molecules were randomly put into one void channel with n = 12 for the DHF-based framework and 25 for the BDC-based framework. These systems are so-call framework–nH₂ systems. The total number of atoms in a unit cell is 132 in the case of the DHF-based framework–12H₂ system and 212 in the case of the BDC-based framework–25H₂ system. All framework–nH₂ systems were equilibrated at 300 K for 10 ps by AIMD simulations in NVT ensembles and then annealed from 300 K to 0 K for 10 ps. The H₂ molecules were expected to reach various adsorption sites on the potential energy surfaces induced by the frameworks at the end of the annealing processes. In order to evaluate the influence of initial conditions on the annealing processes, three initial configurations sampled from 300 K NVT ensembles were considered.

The frameworks with a single H₂ molecule at various adsorption sites, which are so-called single molecule adsorption models, were further optimized from the annealed structures with static DFT calculations. The electronic adsorption energies were estimated via the following formula:

ΔE⁰_i = E^fr+gas_i − E^gas − E^fr (1)

E denotes total energy at T = 0 K derived from static DFT calculations. The superscript 0 denotes the electronic contribution. The superscripts fr + gas, gas, and fr denote the single molecule adsorption system, one free gas molecule and bare framework respectively. Subscript “i” denotes type of adsorption site. The vibrational zero-point energy (ZPE) which purely originates from the quantum nuclear effect was taken into account for the total adsorption energy using the following formula:

ΔE_i = (E^fr+gas_i + ZPE^fr+gas_i) − (E^gas + ZPE^gas) − (E^fr + ZPE^fr) (2)

ZPE denotes the vibrational zero-point energy correction calculated using the following formula:²⁵

The sum is taken over all vibrational modes of the system, h is the Planck's constant and v_j is the frequency of the jth vibrational mode. The adsorption energies calculated using the single molecule adsorption model correspond to the zero loading adsorption energies.

First principles calculations were introduced within the framework of DFT implemented within the Vienna ab initio simulation package (VASP).^26,27 Exchange–correlation potentials were expressed within the generalized gradient approximation (GGA) using Perdew–Burke–Ernzerhof (PBE) functional.²⁸ The dispersion correction was added to original Kohn–Sham DFT equations in terms of atom-pairwise sum-over C₈R⁻⁸ potentials with cut-off distances and dispersion coefficients explicitly computed by S. Grimme et al. from first principles.²⁹ This method is commonly called DFT-D3. The valence electrons' wave functions, i.e. the wave functions of Mg 3s, C and O 2s, 2p and H 1s electrons were represented in plane wave basis set with cut-off energy of E_cut-off = 550 eV. Highly oscillating core parts were simulated with projector augmented wave (PAW) method.^30,31 Since both BDC-based and DHF-based Mg-MOF-74 frameworks have large unit cells of sizes ∼26.08 × 26.08 × 6.87 ų and ∼18.76 × 18.76 × 6.87 ų respectively,¹⁹ total energies were calculated on a Γ-center 1 × 1 × 3 Monkhorst–Pack mesh in reciprocal space. Both cut-off energy and k-mesh were chosen according to the convergence of the electronic adsorption energy of the strongest site within an error <1 meV per atom in comparison to the case using E_cut-off = 1000 eV and k-mesh of 4 × 4 × 12.

Self-consistent-field (SCF) convergence criterion for electronic structure optimization was 10⁻⁶ eV. For crystal structural optimization, all atomic positions were relaxed until all atomic forces smaller than 10⁻³ eV Å⁻¹. The vibrational frequencies required for calculating ZPE corrections were obtained from normal mode analysis within harmonic approximation. The elements of the Hessian matrix were calculated from the DFT forces using finite difference method with the atomic step size 0.015 Å. The frequencies were then calculated as the square roots of the eigenvalues of the Hessian matrix. For reducing computational costs, only the vibrations of hydrogen molecules were considered. Benchmarks on two adsorption sites of two different types in DHF-based framework showed that the ZPE contributions of the framework to the adsorption energies are less than 2 meV.

AIMD simulations were carried out within the framework of Born–Oppenheimer approximation.³² The temperature was controlled using Nosé–Hoover thermostat while the optimized volume of the unit cell of the frameworks were kept constant.³³ During the AIMD simulations, all displacements of the frameworks and gas molecules were allowed. The time step is 0.2 fs – 1/40 of the period of H–H stretching mode.

3. Results and discussion

3.1. Crystal structure and porous geometry

The porosity characteristics including accessible surface area per gram (ASA), accessible volume per gram (AV), largest included sphere (LIS) radius (R_I) and largest free sphere (LFS) radius (R_F) were examined using Voronoi-mesh method with Zeo++ package.³⁴ A schematic representation of the LIS and the LFS is shown in Fig. 1c. In Table 1, calculated structure parameters and porosity characteristics of the two isoreticular Mg-MOF-74 frameworks are summarized and compared with those derived from powder neutron diffraction (PND) data³⁵ and previous PBE calculations using the same basis set and pseudopotential method.²³ According to this table, the dispersion correction is important to reproduce experimental lattice constants of the BDC-based framework. It should be noted that the NPD data show a contraction of the lattice constant a and extension of the lattice constant c when the temperature increases from 20 K to ∼470 K.³⁵

Table 1 A summary of lattice parameters and porosity characteristics (accessible surface area (ASA), accessible volume (AV), largest included sphere (LIS) radius (R_I) and largest free sphere (LFS) radius (R_F)) of the BDC-based and DHF-based Mg-MOF-74 frameworks

	BDC-based Mg-MOF-74 framework			DHF-based Mg-MOF-74 framework
	PBE-D3/plane wave + PAW (this work)	PBE/plane wave + PAW^23	PND^35	PBE-D3/plane wave + PAW (this work)
a = b (Å)	26.06	26.22	25.92 (T = 20 K), 25.88 (T = 298 K)	18.80
c (Å)	6.86	6.96	6.86 (T = 20 K), 6.88 (T = 298 K)	6.84
ASA (m^2 g^−1)	1763	—	—	1042
AV (cm^3 g^−1)	0.340	—	—	0.088
RI (Å)	11.7	—	—	7.6
RF (Å)	10.9	—	—	6.2

---

The resulted unit cell volume change is a tiny oscillation ∼0.09% around the average value ∼3393 ų. The lattice constant changes are also very small within this temperature range, i.e. <0.3%. Hence, the optimized lattice unit cell can be used for the NVT AIMD simulation within the temperature range from 300 K down to 0 K.

There is no experimental data for the porosity characteristics except from the ASA which can be compared with Langmuir surface area calculated from the experimental N₂ adsorption isotherms. However, only experimental N₂ adsorption isotherms for the BDC-based Mn, Fe, Co and Ni-MOF-74 frameworks were reported with the deduced Langmuir surface areas 1797, 1535, 1432 and 1574 m² g⁻¹ respectively.¹⁶ The theoretical ASA value for the BDC-based Mg-MOF-74 framework 1763 m² g⁻¹ is well within the experimental range. The shorter linker of the DHF-based framework lead to smaller lattice constant a (b) and narrower void channel with smaller LIS and LFS radii but does not significantly change the lattice constant c. The narrower void channel, in turn, lead to the smaller ASA and AV. The ASA of the BDC-based framework is 1.69 times as large as that one of the DHF-based framework. Whereas the AV is 3.86 times larger accordingly.

3.2. Annealing processes

Fig. 2 shows the annealed structures of framework–nH₂ systems starting from three different initial configurations sampling within 300 K NVT ensembles. A general agreement between these annealing simulations is that there are five types of H₂ adsorption sites in the BDC-based framework and three types of adsorption sites in the DHF-based framework. In the Fig. 2, adsorption positions of the same type are represented by the same color of hydrogen atoms. Before going to analyse in detail of the adsorption positions, we first discuss the occupation temperature range within which H₂ molecules gradually occupy adsorption sites of each type. Table 2 is a summary of the occupation temperature ranges of various types with three different initial configurations. According to this, different initial configurations lead to a slight fluctuation ∼5 ÷ 15 K of the occupation temperature range of each type. This fluctuation can be considered as a result of the thermal fluctuation of the system during the annealing process. The initial configuration does not affect the occupation order between various types.

Fig. 2 Annealed structures of the BDC-based Mg-MOF-74–25H₂ system (on the left) and DHF-based Mg-MOF-74–12H₂ systems (on the right) starting from three different initial configurations sampling within the 300 K NVT ensemble. Red, brown, cyan and white spheres represent oxygen, magnesium, carbon and hydrogen atoms of the frameworks. Hydrogen atoms of the molecules adsorbed at P₁, P₂, P₃, P₄ and P₅ sites of the BDC-based framework are represented by blue, yellow, green, pink and white spheres respectively. Hydrogen atoms of the molecules adsorbed at P₁, P₂ and P₃ sites of the DHF-based framework are represented by blue, pink and white spheres respectively.

Table 2 A summary of adsorption temperature ranges in which various types of adsorption sites are populated during simulated annealing processes starting from 3 different initial configurations sampling from the 300 K NVT ensemble

	P1	P2	P3	P4	P5
BDC based	300–210 K	205–55 K	70–30 K	45–20 K	25–15 K
	300–215 K	220–55 K	85–30 K	55–20 K	30–20 K
	300–175 K	170–70 K	80–45 K	55–10 K	40–15 K
DHF based	300–170 K	260–50 K	50 K
	300–210 K	155–65 K	—
	300–210 K	150–20 K	20 K

---

For the BDC-based MOF–25H₂ system, during the simulated annealing processes, the blue sites in Fig. 2, start to be occupied first, and then yellow sites, green sites, pink sites and white sites which are called P₁, P₂, P₃, P₄ and P₅ sites respectively. Due to the symmetry of honeycomb-like lattice, there are 6 sites of each type, and hence totally 30 adsorption sites, in one void channel of a R3̄ unit cell. The P₁, P₂ and P₃ sites have been experimentally found for the BDC-based Zn- and Ni-MOF-74 frameworks via NPD data.^20,21 Interestingly, the NPD experiments showed that these sites can be fully occupied in the BDC-based Zn-MOF-74 framework at T down to 30 K (ref. 20) which is the same as the lowest temperature for the P₃ sites to be fully occupied in our simulated annealing processes. The P₅ sites were coincidentally observed via the NPD experiments for the BDC-based Zn-MOF-74 framework at T = 15 K. It should be noted that these sites were not experimentally observed for the Ni-based framework even at T = 4 K because, in this case, the NPD data were collected at low loadings (0.5 and 1.5H₂/Ni).²¹ The P₄ sites were not experimentally observed for the BDC-based Zn-MOF-74 framework at loadings up to 4.2H₂/Zn. The absence of P₄ sites in this framework will be discussed with the electronic structure analysis later. It should be noted that the number of H₂ molecules in the simulations, 25 molecules, is smaller than the total number of all adsorption sites, 30 sites. The first three types of adsorption sites have higher probability to be occupied so that they are fully occupied at the end of the annealing processes. The P₄ and P₅ sites are not fully occupied so the final temperatures of the corresponding occupation temperature ranges are the temperatures at which H₂ molecules populate P₄ and P₅ sites the last time.

For the DHF-based MOF–12H₂ system, the occupation order is the blue sites first, then yellow sites and white sites in the Fig. 2, mentioned as P₁, P₂ and P₃ sites respectively. The occupation temperature ranges of P₁ and P₂ sites shown in the Table 2 are similar to the P₁ and P₂ sites in the BDC-based framework, respectively. However, the dependence of the occupation temperature range of the P₂ sites on the initial configuration is stronger, denoting by a stronger thermal fluctuation ∼100 K. There are 6 P₁ and 6 P₂ sites but only one P₃ site, and hence 13 adsorption sites in one void channel of a unit cell while there are 12H₂ molecules. The P₁ sites are fully occupied first at temperature down to 170 K. In two of the simulated annealing processes, only 5 of 6 P₂ sites are occupied, the last H₂ molecule go to P₃ site instead of populating the last P₂ site. The P₃ site is occupied at the same time with the final occupation of P₂ sites at temperature down to 20 K.

In general, P₁, P₂, P₃ and P₄ sites of the BDC-based Mg-MOF-74 framework and P₁ and P₂ sites of the DHF-based framework are direct adsorption sites since they are directly attached to the frameworks. The P₅ sites of the BDC-based framework and P₃ sites of the DHF-based framework are indirect adsorption sites. In the following, direct adsorption sites will be discussed in detail to clarify the mechanism of the interaction between the frameworks and the H₂ molecules.

3.3. Adsorption positions and binding mechanisms

3.3.1. Adsorption positions. Fig. 3 shows the optimized positions of a single H₂ molecule at various adsorption sites in the two Mg-MOF-74 isoreticular frameworks. A P₁ site of the DOBDC-based and DHF-based frameworks shown in Fig. 3a and e is the nearest to an open metal site with a distance 2.41 Å and 2.39 Å respectively. This distance is in good agreement with the experimental value 2.45 Å derived from the NPD data collected at 77 K for the DBC-based Mg-MOF-74 framework.³⁶ At this site, the H₂ molecule is likely in side-on position with respect to the Mg²⁺ ion. The side-on position was shown to be common for a H₂ molecule when it is attached to an isolated positive ion.³⁷ Since the H₂ molecule has a quadrupole moment 2.21 × 10⁻⁴⁰ C m²,³⁸ the side-on position with respect to the positive ion is to optimize the positive charge–quadrupole interaction. However, with the ligand environment of the MOF-74 frameworks, the hydrogen molecule is not exactly in ideal side-on position with equal Mg–H distances for two hydrogen atoms of the molecule. It slightly tilts, giving rise to the difference in Mg–H distances 0.01 Å in the case of the BDC-based framework and 0.02 Å in the case of the DHF-based framework.

Fig. 3 Foreground illustrations of (a) P₁, (b) P₂, (c) P₃ and (d) P₄ adsorption positions in the BDC-Mg-MOF-74 framework. The P₁ and P₂ sites in the DHF-Mg-MOF-74 are illustrated in (e) and (f) respectively. Red, brown, cyan and white spheres represent oxygen, magnesium, carbon and hydrogen atoms of the framework respectively. P₁, P₂, P₃ and P₄ hydrogen molecules in the BDC-based framework are represented by blue spheres in (a), yellow spheres in (b), green spheres in (c) and pink spheres in (d) respectively. P₁ and P₂ hydrogen molecules in the DHF-based framework are represented by blue spheres in (e) and pink spheres in (f) respectively. Dash lines denote the connection between the adsorbed H₂ molecules and the nearest framework atoms. Dot lines in (b) denote the oxygen triangle close to the P₂ site of the BDC-based Mg-MOF-74 framework.

A P₂ site of the BDC-based Mg-MOF-74 framework shown in Fig. 3b is located above a triangle of three oxygen atoms which connect to the same metal ion, one from an oxido group called O1 and two from two carboxyl groups called O2 (see Fig. 1b). This site locates the nearest to O1 atoms with a distance 3.21 Å. The experimental values for this distance are 3.5 Å from the NPD data collected at 77 K.³⁶ At this site, the H₂ molecule is likely in head-on position with respect to the O1 atom. This is also the common position for a H₂ molecule when it is attached to an isolated negative ion in order to optimize the negative charge–quadrupole interaction.³⁷ The BDC ligand environment of the framework also gives rise to a slight tilt by an angle 6° from the ideal head-on position. A P₂ site of the DHF-based framework shown in Fig. 3f is the nearest to the oxygen atom of a carboxyl group O2 at a distance 2.47 Å. The H₂ molecule populate this site nearly at head-on position to the O2 atom with a significant tilting angle 22.5° induced by the DHF ligand environment.

A P₃ site of the BDC-based Mg-MOF-74 framework shown in Fig. 3c is above an aromatic ring of a linker at a distance 3.05 Å from the ring plane. This distance is in good agreement with the NPD value 3.08 Å for Co-based framework. It was shown that when a hydrogen molecule interacts with a molecule of a benzene dicarboxylic acid derivative, it prefers to bind head-on above the center of the aromatic ring at a distance ∼3 Å with its axes perpendicular to the ring plane.³⁷ However, in the BDC-based framework, the H₂ molecule is significantly inclined by an angle of 25° with the normal of the ring plane denoting significant influence of metal cluster nodes on this linker-induced adsorption site. A P₄ site of the BDC-based Mg-MOF-74 framework shown in Fig. 3d is the nearest to the oxygen atoms of a carboxyl group O2 at a distance 3.01 Å. Similar to the P₂ site, the H₂ molecule is inclined by an angle of 6° from the ideal head-on position to the O2 atom because of the BDC ligand environment.

3.3.2. Microscopic origin of the H₂ attractions. The deformation of electron density upon adsorption at various sites at 0 K was calculated using the following formula:ρ denotes the electron density, Δρ denotes the electron density deformation. The subscripts framework + H₂, framework and H₂ denote the values of framework with one adsorption site occupied, bared framework and of an isolated hydrogen molecule respectively. Fig. 4 shows the iso-surfaces of the electron density deformation at isovalue 0.0001|e| Å⁻³.

Fig. 4 Isosurfaces of the electron density deformations upon hydrogen adsorption onto (a) P₁, (b) P₂, (c) P₃ and (d) P₄ sites in the BDC-based Mg-MOF-74 framework, (e) P₁ and (f) P₂ sites in the DHF-based Mg-MOF-74 framework. The isovalues is 0.0001|e| Å⁻³. Yellow and green denote electron donated and loss surfaces respectively. Red, brown, cyan and white spheres represent oxygen, magnesium, carbon and hydrogen atoms of the frameworks respectively. P₁, P₂, P₃ and P₄ hydrogen molecules in the BDC-based framework are represented by blue, yellow, green and pink spheres respectively. P₁ and P₂ hydrogen molecules in the DHF-based framework are represented by blue and pink spheres respectively.

For a P₁ site in both frameworks, the electron density is increased in an extended area on the Mg side of the H₂ molecule and decreased on the other side. This electron density deformation clearly evidences that the positive charge of the Mg²⁺ cation polarizes the hydrogen molecule perpendicularly to its axis. Thus, the interaction of the Mg²⁺ charge with the hydrogen induce-dipole is the dominating force to hold H₂ molecules at a P₁ site in both frameworks. Additionally, there are electron loss in the area between the oxygen and magnesium atoms, but the electron enhancements in the p-shaped area around the oxygen atoms and on the aromatic ring carbon atoms called C1 (see Fig. 1a). These atoms are located within a radius 3.5 Å around the P₁ site. This indicates remarkable contribution of the surrounding ligand environment within the radius 3.5 Å to the hydrogen attraction of the P₁ site.

The electron density deformations upon hydrogen adsorption to P₂ and P₄ sites in the BDC-based framework and a P₂ site in the DHF-based framework demonstrate the dominating role of the interaction between the H₂ induced-dipole and oxygen ions. The electron donated and loss regions on the two heads of the hydrogen molecule indicate a polarization of the hydrogen molecule along its axis. Obvious electron density deformations can be seen around the nearest framework atoms, i.e. O1 for the P₂ site, O2 for the P₄ site in the BDC-based framework and O2 for the P₂ site in the DHF-based framework. There are also other noticeable electron density deformation regions. For the P₂ site of the BDC-based framework, they are the regions around the two other oxygen atoms O2 of the oxygen triangle and the carbon atoms of the related carboxyl group C2 at distance. For the P₄ site of the BDC-based framework, they are the regions around the H atom of the aromatic ring. For the P₂ site of the DHF-based framework, they are the region between two carbon atoms within the radius 3.5 Å. These deformations indicate the contributions of the surrounding ligand environment within the radius 3.5 Å. The electron deformation upon the H₂ adsorption at a P₃ site in the BDC-based framework shows an electron donated region above the aromatic ring and a polarization of the H₂ molecule along its axis. This demonstrates the interaction between the H₂ induced dipole and the electron cloud of the ring at this site.

3.3.3. The binding strength and the role of charge distribution. Table 3 summarizes the electronic parts of the adsorption energies or the non-corrected adsorption energies ΔE⁰, ZPE corrections and ZPE corrected adsorption energies ΔE of the direct adsorption sites. Negative values of adsorption energies mean that the bound states are more stable than unbound state. The contributions of vibrational ZPEs are positive, denoting a destabilization of the host–guest binding states. According to the adsorption energies, the descending order of binding strengths is P₁, P₂, P₃ and P₄ in the BDC-based framework and P₁ and P₂ in the DHF-based framework. It consistently follows the order of occupation during the simulated annealing processes.

Table 3 A summary of noncorrected adsorption energies ΔE⁰, vibrational zero point energy corrections ZPE, ZPE corrected adsorption energies ΔE at zero loading and noncorrected adsorption energies ΔE⁰_f, vibrational zero point energy corrections ZPE_f, ZPE corrected adsorption energies ΔE_f at full loading of the P₁ sites of direct adsorption sites in the BDC-based and DHF-based Mg-MOF-74 frameworks

	BDC-based Mg-MOF-74 framework				DHF-based Mg-MOF-74 framework
	P1	P2	P3	P4	P1	P2
ΔE^0 (eV)	−0.192	−0.119	−0.105	−0.099	−0.183	−0.128
ΔE (eV)	−0.121	−0.088	−0.066	−0.062	−0.106	−0.080
ZPE (eV)	0.070	0.031	0.038	0.037	0.077	0.048
ΔE^0f (eV)	−0.171	−0.110	−0.090	−0.074	—	—
ΔEf (eV)	−0.100	−0.057	−0.042	−0.027	—	—
ZPEf (eV)	0.071	0.053	0.048	0.047	—	—

---

As discussed above, the electronic origin of the H₂ attraction is the electrostatic interaction between the H₂ molecule and the surrounding charged fragments of the MOF-74 frameworks within the radius ∼3.5 Å. Thus, charge distribution should be quantitatively considered to explain for the binding strength at various adsorption sites. The atomic charges were calculated from the DFT-derived charge density using Bader charge analysis method which defined atomic charge as total charge within zero charge flux surface.^39–42 Table 4 show the calculated Bader atomic charges of the MOF-74 frameworks.

Table 4 Atomic charges of metal and oxygen atoms in the MOF-74 frameworks derived from the DFT electron density by Bader charge analysis

	BDC-based Mg-MOF-74 framework	DHF-based Mg-MOF-74 framework	BDC-based Zn-MOF-74 framework
Metal	+2|e|	+2|e|	+1.36|e|
O1	−1.93 ± 0.02|e|	−1.93 ± 0.02|e|	−1.19 ± 0.02|e|
O2	−1.88 ± 0.01|e|	−1.88 ± 0.01|e|	−1.16 ± 0.01|e|

---

Both P₁ sites in the BDC- and DHF-based framework are primarily based on the interaction of the positive Mg ion and its surrounding oxygens with the H₂ molecule. As shown in Table 4, the Bader atomic charges of these ions are the same for both frameworks. However, the binding strength of a P₁ site in the DHF-based framework is smaller than in the BDC-based framework by 0.015 eV – 1.45 kJ mol⁻¹. This difference can be assigned to the contribution of the aromatic ring in the BDC linker as denoted by the electron deformation on the ring, which is absent from the nonaromatic DHF linker. Similar linker-induced binding strength variation was also experimentally observed when the configuration of the BDC linker was changed from para to meta (mBDC linker). In this case, the extra contribution of the mBDC linker lead to an enhancement of the binding strength at the metal-induced site by 1.5 kJ mol⁻¹.¹⁶

In the BDC-based Mg-MOF-74 framework, the Bader charge of oxygen ion in the oxido group O1 is slightly more negative than that of the oxygen ion in the carboxyl group O2 by 0.05|e| as seen in Table 4. However, the binding strength of the O1 induced adsorption site P₂ is significantly larger than the O2 induced adsorption site P₄, 0.023 eV (2.22 kJ mol⁻¹). This is due to the fact that, within the radius 3.5 Å, there are more framework atoms surrounding the P₂ site (one O1, two O2 and one C2 atom) than the P₄ site (only one O2 and one H atom). The Bader charge of the O2 ion in the DHF-based Mg-MOF-74 framework is the same as that of the O2 atom and less negative than that of the O1 atom in the BDC-based Mg-MOF-74 framework. However, the binding strength of the O2 induced site P₂ in the DHF-based framework is larger than the O2 induced site P₄ and very close to the O1 induced site P₂ in the BDC-based framework. The reason is that this site is also surrounded by more framework atoms (O2 and carbon atoms) than the O2 induced site P₄ in the BDC-based framework.

It is worth noting from experiments of H₂ adsorption in the BDC-based Mg-MOF-74 framework that the heat of adsorption depend on the loading rate.¹³ The heat of adsorption at zero loading is a bit higher than at nonzero loading rate below 1H₂/metal. In this loading rate range, the P₁ sites provide the main contribution to the hydrogen uptake.^20,21 The slight reduction of the heat of adsorption can be assigned to the site–site interaction. In the present work, the adsorption energy of each type of adsorption sites was calculated at zero loading and full loading rate of P₁ sites (1H₂/metal). The absolute values of the ZPE corrected adsorption energy of a P₁ site at the zero loading and 1H₂/metal loading in the BDC-based Mg-MOF-74 framework are 0.121 eV and 0.100 eV, respectively. These values are in good agreement with the experimental data of zero loading and low loading enthalpies of adsorption 0.119 eV and 0.109 eV.¹⁵ At loading rate larger than 1H₂/metal, the main contributions to the experimental enthalpies come from the P₂ and P₃ sites. The absolute values of the ZPE corrected adsorption energies of a P₂ and P₃ site at 1H₂/metal loading are 0.057 eV and 0.042 eV respectively, which are around the experimental enthalpies of adsorption at loading rate >1H₂/metal ∼0.052 eV.

Concerning the absence of P₄ sites in the BDC-based Zn-MOF-74 framework suggested by the NPD experiment,²⁰ it is worth noting that the main origin of these sites is the O2 negative charge. The more negative the O2 ion is, the stronger the H₂ attraction is. On the other hand, Mg is more electro-positive than Zn so that the oxygen atoms surrounding Mg are expected to be more negatively charged than those surrounding Zn. The reducing in charge density at the O2 site is expected to reduce the binding strength at the P₄ site or even destabilize this site. In this work, the charge distribution of the BDC-based Zn-MOF-74 was examined using DFT derived electron density and Bader charge analysis. The results are also shown in Table 4. It is clear that the charge of Mg ion in the Mg-MOF-74 framework is +2|e|, larger that of Zn ion in the Zn-MOF-74 framework +1.36|e|. The ionic charges of neighbouring O2 atoms are −1.9 and −1.2|e| in the Mg- and Zn-based framework respectively. The ZPE corrected adsorption energy of a P₄ site in the BDC-based Zn-MOF-74 was also calculated at zero loading and 1H₂/metal loading. The zero loading value is −0.042 eV, 0.02 eV – 1.93 kJ mol⁻¹ weaker than the case of the Mg-MOF-74. At loading rate 1H₂/metal, the adsorption energy of a P₄ site is −0.027 eV and −0.007 eV for Mg- and Zn-MOF-74 framework, respectively. The predicted smaller values of the adsorption energy of the P₄ site agree with the fact that it has been absent in the Zn-based framework. However, due to the higher adsorption energy it would be interesting to check experimentally whether the H₂ molecules could be occupied at the P₄ sites in Mg-based framework.

The electronic adsorption energy at the ring induced site P₃ in the BDC-based MOF-74 framework is 0.105 eV (10.13 kJ mol⁻¹), higher than those of the BDC derivatives ∼4 kJ mol⁻¹.³⁷ The ZPE corrected adsorption energy of this site is −0.066 eV (6.37 kJ mol⁻¹), stronger than that of the ring site in the MOF-5 framework 5.1 kJ mol⁻¹.²² The enhancement of the binding strength between the ring and the H₂ molecule in the MOF-74 can be explained basing on the fact that adding electron-donated substituents such as –OH, –CH₃, –NH₂ can increase the binding strength.³⁷ In the case of MOF-74, the enhancement of H₂ attraction can come from the fact that there are two oxido substituents on the ring which directly connect to strong electron-donating metal atoms. Whereas, there is no such substituent on the BDC linker in the MOF-5.

4. Conclusions

In conclusion, the microscopic origin of H₂ adsorption onto the internal surface of the MOF-74 frameworks were explored on the base of AIMD simulations and DFT electronic structure calculations. The guest molecules show significant electrostatic interactions with the local framework environment within a radius 3.5 Å. Four types of direct adsorption sites on the internal surface of the BDC-based Mg-MOF-74 framework were specified. With shorter nonaromatic linker, the DHF-based MOF-74 framework provides only two types of direct adsorption sites. The primary adsorption sites P₁ in both frameworks are due to the interaction of the H₂ molecules with the positively charged open metal sites and their surrounding oxygen sites. In the case of the aromatic BDC linker, the contribution of the high electron density aromatic rings within the distance 3.5 Å from these sites enhances the binding strength by 1.45 kJ mol⁻¹ in comparison to the case of nonaromatic DHF linker. The secondary adsorption sites originate from the interactions of the H₂ molecule with the negatively charged oxygen sites and the high density electron clouds of the aromatic rings. The reduction of the electro-positivity of the metal can induce reductions of the negative charges of the oxygen sites which in turn destabilize the oxygen induced adsorption sites.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors would like acknowledge the financial supports of the World Bank and the Ministry of Science and Technology of Vietnam, grant number 13/FIRST/1.a/VNU1; and Vietnam National University Ho Chi Minh city, grant number TX-2020-50-01.

References

1. S. T. Meek, J. A. Greathouse and M. D. Allendorf, Metal-organic frameworks: a rapidly growing class of versatile nanoporous materials, Adv. Mater., 2011, 23, 249–267.
2. H. Furukawa, K. E. Cordova, M. O'Keeffe and O. M. Yaghi, The chemistry and applications of metal-organic frameworks, Science, 2013, 341(6149), 1230444.
3. S. Ahmad, A. Ansari, W. A. Siddiqi and M. K. Akram, Nanoporous Metal-Organic Framework, Mater. Res. Found., 2019, 58, 107–139.
4. A. Dhakshinamoorthy, Z. Li and H. Garcia, Catalysis and photocatalysis by metal organic frameworks, Chem. Soc. Rev., 2018, 47, 8134–8172.
5. R. B. Getman, Y.-S. Bae, C. E. Wilmer and R. Q. Snurr, Review and Analysis of Molecular Simulations of Methane, Hydrogen, and Acetylene Storage in Metal–Organic Frameworks, Chem. Rev., 2012, 112, 703–723.
6. H. Lee, T.-B. Nguyen, D.-K. Nguyen, J.-H. Kim, J.-Y. Kim, B. T. Phan and S. S. Kim, Gas Sensing Properties of Mg-Incorporated Metal–Organic Frameworks, Sensors, 2019, 19(15), 3323.
7. P. Horcajada, R. Gref, T. Baati, P. K. Allan, G. Maurin, P. Couvreur, G. Ferey, R. E. Morris and C. Serre, Metal–Organic Frameworks in Biomedicine, Chem. Rev., 2012, 112, 1232–1268.
8. M. P. Suh, H. J. Park, T. K. Prasad and D.-W. Lim, Hydrogen Storage in Metal–Organic Frameworks, Chem. Rev., 2012, 112, 782–835.
9. A. Ahmed, S. Seth, J. Purewal, A. G. Wong-Foy, M. Veenstra, A. J. Matzger and D. J. Siegel, Exceptional hydrogen storage achieved by screening nearly half a million MOFs, Nat. Commun., 2019, 10, 1568.
10. H. Frost and R. Q. Snurr, Design Requirements for Metal-Organic Frameworks as Hydrogen Storage Materials, J. Phys. Chem. C, 2007, 111, 18794.
11. Y. E. Cheon and M. P. Suh, Selective gas adsorption in a microporous metal–organic framework constructed of Co₄^II clusters, Chem. Commun., 2009, 45, 2296.
12. T. Pham, K. A. Forrest, R. Banerjee, G. Orcajo, J. Eckert and B. Space, Understanding the H₂ Sorption Trends in the M-MOF-74 Series (M = Mg, Ni, Co, Zn), J. Phys. Chem. C, 2015, 119, 1078–1090.
13. W. Zhou, H. Wu and T. Yildirim, Enhanced H₂ Adsorption in Isostructural Metal–Organic Frameworks with Open Metal Sites: Strong Dependence of the Binding Strength on Metal Ions, J. Am. Chem. Soc., 2008, 130, 15268–15269.
14. J. L. C. Rowsell and O. M. Yaghi, Effects of Functionalization, Catenation, and Variation of the Metal Oxide and Organic Linking Units on the Low-Pressure Hydrogen Adsorption Properties of Metal–Organic Frameworks, J. Am. Chem. Soc., 2006, 128, 1304–1315.
15. P. D. Dietzel, P. A. Georgiev, J. Eckert, R. Blom, T. Strassle and T. Unruh, Interaction of hydrogen with accessible metal sites in the metal–organic frameworks M2(dhtp) (CPO-27-M; M = Ni, Co, Mg), Chem. Commun., 2010, 46, 4962–4964.
16. 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. H. Gordon, C. M. Brown and J. R. Long, M2(m-dobdc) (M = Mg, Mn, Fe, Co, Ni) Metal–Organic Frameworks Exhibiting Increased Charge Density and Enhanced H₂ Binding at the Open Metal Sites, J. Am. Chem. Soc., 2014, 136, 12119–12129.
17. M. T. Kapelewski, T. Runčevski, J. D. Tarver, H. Z. H. Jiang, K. E. Hurst, P. A. Parilla, A. Ayala, T. Gennett, S. A. FitzGerald, C. M. Brown and J. R. Long, Record High Hydrogen Storage Capacity in the Metal–Organic Framework Ni2(m dobdc) at Near-Ambient Temperatures, Chem. Mater., 2018, 30(22), 8179–8189.
18. https://www.energy.gov/eere/fuelcells/hydrogen-storage.
19. M. Witman, S. Ling, A. Gladysiak, K. C. Stylianou, B. Smit, B. Slater and M. Haranczyk, Rational Design of a Low-Cost, High-Performance Metal–Organic Framework for Hydrogen Storage and Carbon Capture, J. Phys. Chem. C, 2017, 121, 1171–1181.
20. Y. Liu, H. Kabbour, C. M. Brown, D. A. Neumann and C. C. Ahn, Increasing the density of adsorbed hydrogen with coordinatively unsaturated metal centers in metal-organic frameworks, Langmuir, 2008, 24, 4772–4777.
21. C. M. Brown, A. J. R. Cuesta, J.-H. Her, P. S. Wheatley and R. E. Morris, Structure and spectroscopy of hydrogen adsorbed in a nickel metal-organic framework, Chem. Phys., 2013, 427, 3–8.
22. K. Sillar, A. Hofmann and J. Sauer, Ab Initio Study of Hydrogen Adsorption in MOF-5, J. Am. Chem. Soc., 2009, 131, 4143–4150.
23. P. Suksaengrat, V. Amornkitbamrung, P. Srepusharawoot and R. Ahuja, Density Functional Theory Study of Hydrogen Adsorption in a Ti-Decorated Mg-Based Metal–Organic Framework-74, ChemPhysChem, 2016, 17, 879–884.
24. M. Dixit, D. T. Major and S. Pal, Hydrogen adsorption in ZIF-7: a DFT and ab initio molecular dynamics study, Chem. Phys. Lett., 2016, 651, 178–182.
25. D. A. Mc Quarrie and J. D. Simon, Molecular Thermodynamics, University Science Books, Sausalito, CA, 1999.
26. G. Kresse and J. Furthmuller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
27. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev., 1999, 59, 1758.
28. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865.
29. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132, 154104.
30. P. E. Blochl, Projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
31. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
32. D. Marx and J. Hutter, Ab initio Molecular Dynamics: Basic Theory and Advanced Methods, Cambridge UP, 2012.
33. D. Frenkel and B. Smit, Understanding Molecular Simulation, Academic Press, London, 1996.
34. T. F. Willems, C. H. Rycroft, M. Kazi, J. C. Meza and M. Haranczyk, Algorithms and tools for high-throughput geometry-based analysis of crystalline porous materials, Microporous Mesoporous Mater., 2012, 149, 134–141.
35. W. L. Queen, C. M. Brown, D. K. Britt, P. Zajdel, M. R. Hudson and O. M. Yaghi, Site-Specific CO2 Adsorption and Zero Thermal Expansion in an Anisotropic Pore Network, J. Phys. Chem. C, 2011, 115, 24915–24919.
36. K. Sumida, C. M. Brown, Z. R. Herm, S. Chavan, S. Bordigad and J. R. Long, Hydrogen storage properties and neutron scattering studies of Mg2(dobdc)—a metal–organic framework with open Mg²⁺ adsorption sites, Chem. Commun., 2011, 47, 1157–1159.
37. R. C. Lochan and M. Head-Gordon, Computational studies of molecular hydrogen binding affinities: the role of dispersion forces, electrostatics, and orbital interactions, Phys. Chem. Chem. Phys., 2006, 8, 1357–1370.
38. D. E. Stogryn and A. P. Stogryn, Molecular multipole moments, Mol. Phys., 1966, 11, 371.
39. W. Tang, E. Sanville and G. Henkelman, A grid-based Bader analysis algorithm without lattice bias, J. Phys.: Condens. Matter, 2009, 21, 084204.
40. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, An improved grid-based algorithm for Bader charge allocation, J. Comput. Chem., 2007, 28, 899–908.
41. G. Henkelman, A. Arnaldsson and H. Jónsson, A fast and robust algorithm for Bader decomposition of charge density, Comput. Mater. Sci., 2006, 36, 354–360.
42. M. Yu and D. R. Trinkle, Accurate and efficient algorithm for Bader charge integration, J. Chem. Phys., 2011, 134, 064111.

---