Computational determination of coordination structure impact on adsorption and acidity of pristine and sulfated MOF-808

Abstract

Metal–organic frameworks (MOFs) are composed of metal nodes connected by organic linkers. With the massive amount of possible metal and linker combinations, it is critical to develop structure–property relationships. However, a major impediment to developing these relationships is that MOFs often have multiple metal coordination structures in a single crystal. Here, we report periodic and cluster density functional theory (DFT) calculations analyzing the coordination structures of MOF-808 and sulfated MOF-808 and their physical and chemical properties. For MOF-808, we determined coordination structures by comparing computationally determined lattice constants with experimental values and then used these coordination structures to simulate N₂ and Ar adsorption isotherms. Our simulated average N₂ and Ar uptakes agree very well with the experimental values. For the sulfated MOF-808, which has been proposed to be a superacidic material, we determined the impact of coordination structure on acidity. Surprisingly, our results based on calculated proton affinities suggest a 10²⁵ range in acid site strength depending on the coordination structure, with only a few sites having high acidity. The use of vibrational frequencies and other property-based methods for determining relative acidity fail. Our analysis revealed that the acidity of sulfated MOF-808 is unlikely controlled by a sulfate–water hydrogen bond that was previously proposed. Instead, we show that there is a strong correlation between conjugate base stability and proton affinity to rationalize acidity of protons in a single coordination structure model as well as across all coordination structures.

---

1. Introduction

Metal–organic frameworks (MOFs) are crystals that consist of metal nodes connected by organic ligand linkers. The large number of metal nodes and organic linkers provides an almost overwhelming possible diversity of MOFs with varying sizes, shapes, topologies, and physical and chemical properties.^1,2 In addition to this diversity, MOFs have ordered and tunable structures that provide opportunities for post-synthesis modifications.³ With unique physical and chemical properties, MOFs are beginning to be targeted for many industrial applications, including gas storage, water purification, separation, toxic chemical decomposition, and catalysis.^4–6

While experimental and computational studies have begun to understand MOF structure–property relationships,^5,7,8 because of the large chemical diversity of MOFs and varying physical and chemical properties, there remains a large gap in knowledge about how the coordination structure of MOFs control particular properties, especially gas adsorption and acidity related to catalysis. The term coordination structure is used here to describe different combinations and arrangements of organic and inorganic ligands that coordinate with the metal centers in the nodes. Complicating the development of coordination-structure-to-property relationships is that MOFs are often non-uniform and have multiple periodic coordination structures in a single crystal. Therefore, the experimentally observed properties of a MOF are potentially the results of co-existent coordination structures. Surprisingly, while thousands of MOFs have been experimentally synthesized, there are only a few computational studies that have systematically examined MOF coordination structures and resulting properties.⁹ For example, density functional calculations in combination with experiments revealed the thermodynamically most stable proton-topology of MOF NU-1000.¹⁰

MOF-808 is a Zr-based mesoporous material that has relatively high thermal and mechanical stability,⁵ which prompted structure modification and multiple studies in the application of adsorption, gas separation, and catalysis. MOF-808 is constructed by linking Zr₆O₄(OH)₄ nodes (i.e., Zr₆ node) with benzene-1,3,5-tricarboxylate (BTC) linkers, as shown in Fig. 1. Each Zr₆ node cluster coordinates with six BTC linkers, and there are twelve open coordination sites on each node that can be utilized for binding other ligands. Due to the open metal sites in MOF-808, experimental H¹-NMR, IR, and XRD data^11–13 have revealed that ligands such as formate and water can coordinate to a Zr₆ node in different configurations. This results in many different coordination structures throughout the crystal of MOF-808. To date, no systematic computational analysis of MOF-808 coordination structures has been reported. Therefore, it remains unknown how local structure controls local and global MOF-808 properties.

Fig. 1 MOF-808 with Zr₆O₄(OH)₄ nodes (shown as cyan polyhedron) and benzene-1,3,5-tricarboxylate linkers. Formates are coordinated to the node in the μ²-bridging configuration.

A sulfated MOF-808 can be obtained through post-synthetic modification by exposing the crystals to aqueous sulfuric acid.^11,13 Sulfate ligands are incorporated into the Zr₆ nodes through the occupation of the aforementioned open metal sites of MOF-808 through Zr ligand substitution. Because there are a variety of possible Zr–sulfate coordination modes, similar to MOF-808, sulfated MOF-808 is a composite of multiple different coordination structures. Sulfated MOF-808 has been reported to be super acidic and acts as a potent Brønsted acid catalyst. For example, Yaghi demonstrated that the acidity of sulfated MOF-808 induces the dimerization of isobutene (2-methyl-1-propene) to isooctene.¹³ The proposed superacidity of sulfated MOF-808 was briefly analyzed using a single coordination structure and acidity proposed to arise from a sulfate ligand interacting with a neighboring aqua ligand through hydrogen bonding.¹³

Here we report periodic and cluster density functional theory (DFT) modeling of MOF-808 and sulfated MOF-808 to directly determine the impact of different coordination structures on N₂ and Ar adsorption and acidity properties as well as build general structure–property relationships. Importantly, this work demonstrates the importance and nontrivial analysis of many coordination structures to analyze and understand the origin of physical and chemical properties for MOF-808. For MOF-808, the most probable coordination structures were determined by comparing the lattice constants of models with experimentally measured lattice constants of MOF-808. N₂ and Ar adsorption isotherms were computed for selected models using grand canonical Monte Carlo (GCMC) simulations, and we obtained N₂ and Ar uptake results that agree with the experimental measurements. For sulfated MOF-808, we find that the originally proposed hydrogen bond interaction in a single coordination structure likely does not significantly control acidity. Instead, the number of aqua groups and their relative positions affect the stability of conjugate base, which results in a strong correlation between conjugate base stability and proton affinity in a single coordination structure model as well as across all coordination structures. The typical use of vibrational frequencies does not adequately reproduce proton affinities.

2. Computational methods

All MOF-808 and sulfated MOF-808 periodic models were optimized using DFT calculations as implemented in the CP2K program.¹⁴ We used the PBEsol^15,16 functional with the addition of the D3 damped dispersion terms of Grimme,¹⁷ the DZVP-MOLOPT basis set for Zr atoms, and TZVP-MOLOPT basis set for C, H, O, and S atoms. A plane wave cut-off energy of 200 Ry, and pseudopotentials for core electrons (as formulated by Geodecker et al.¹⁸) were used for all atoms. The PBEsol functional has been widely applied for calculating geometric properties of bulk solids and crystals, and generally gives accurate results of physical properties, such as lattice constants.^19,20 All periodic models were optimized using triclinic primitive unit cells to facilitate fast calculations. Optimized primitive unit cells were transformed to cubic unit cells to acquire the average lattice constants, which were then compared with experimental values.

Cluster models were constructed from optimized MOF-808 and sulfated MOF-808 periodic structures. The construction involved the truncation of BTC linkers and replacing them with acetate capping groups around one Zr₆ node. The positioning of acetate capping groups around the Zr₆ node was done so that the spatial constraint imposed on the node is the same as in the periodic structure. The position of the methyl-C atom in each acetate is not optimized during geometry optimizations so that the imposed constraints mimic those in a periodic system. All cluster models were optimized using the Gaussian 16 program.²¹ We used the M06-L²² density functional with SDD²³ pseudopotentials and its corresponding basis set for Zr atoms and def2-TZVP^24,25 basis set for all the other atoms. The M06-L functional was chosen because of its general accuracy in transition metal chemistry.²⁶ All cluster models were optimized to their respective ground states which were confirmed by frequency analysis.

In addition to DFT calculations, we have also performed GCMC simulations on the DFT-optimized periodic geometries to obtain adsorption properties. Specifically, isotherms for N₂ and Ar adsorptions at 77 and 87 K, respectively, were simulated using selected MOF-808 models. All GCMC simulations were performed using the RASPA 2.0 program.²⁷ N₂ and Ar molecules were described by the TraPPE force field.²⁸ MOF-808 models were described by the GenericMOFs force field²⁷ and were assumed rigid. The Lorentz–Berthelot mixing rule was used to calculate the Lennard-Jones parameters between different types of atoms. The Lennard-Jones interactions in the simulation were calculated using a spherical cut-off of 12 Å with long-range correction. The long-range electrostatic interactions were evaluated using Ewald summation.²⁹ Partial atomic charges were computed for MOF atoms using the extended charge equilibration method.³⁰ A total of 110 000 GCMC cycles were performed to converge the adsorption isotherms with the first 10 000 cycles for equilibrium and the latter 100 000 cycles for ensemble averages.

3. Results and discussions

3.1 MOF-808 and gas adsorption

We developed 18 different periodic boundary models of pristine (i.e., no defects) MOF-808, each with a different coordination structure. All the MOF-808 models have the same backbone framework structure with Zr₆ nodes interconnected by BTC linkers with each node binding with six linkers (Fig. 2a). All the Zr atoms in a Zr₆ node have a formal oxidation state of +4 and can form eight bonds in total with surrounding ligands and linkers in a square anti-prism geometry. Surrounding each individual Zr atom, in addition to the BTC linkers and μ-oxo, μ-hydroxo groups within the Zr₆ node, a maximum of two other ligands can be introduced at positions α and β shown in Fig. 2b. The types of ligands used in each model were based on the reported synthetic procedures and experimental XRD crystal data.¹¹ Specifically, we consider three ligands including aqua (–OH₂), hydroxo (–OH), and formate (–OOCH).

Fig. 2 (a) Visualization of the Zr₆ node with six BTC coordinating linkers. (b) A top-down view of the Zr₆ node where the BTC linkers were removed for clarity. Each Zr atom can bind two ligands at locations α and β. C atoms are shown in grey, H atoms in white, O atoms in red, and Zr atoms in cyan.

The different coordination structures of each model were designed by changing the number and coordination configuration of ligands. A formate ligand can bind with either one or two Zr atoms. In models where a formate coordinates to a single Zr, the formate can chelate through either κ¹ or κ². In models where a formate coordinates with two Zr atoms (μ²-formate), there is a Zr–formate–Zr bridging structure motif. A list of the models and their coordination geometries is given in Table S1 of the ESI.†

The experimentally reported MOF-808 crystal data suggest cubic unit cells with averaged lattice constants ranging from 35.023³¹ to 35.232 Å.³² The variation in lattice constants results from MOF-808 being a mixture of coordination structures that each give a slightly different set of lattice constants. Therefore, we examined 18 coordination structure models to determine if a single model fits this lattice constants range or if a combination of coordination structures provides a better fit. Of the 18 models, the four models I–IV shown in Fig. 3 provide average lattice constants of 35.282, 35.131, 35.033, and 35.270 Å, respectively, that are closest to the experimental lattice values.

Fig. 3 Top-down views of MOF-808 models I–IV where the BTC linkers were removed for clarity. C atoms are shown in grey, H in white, O in red, and Zr in cyan.

The structural differences between the selected four models can be considered through a series of steps of bond breaking and new bond forming involving the formate ligands. Model I, II, and III all contain only Zr₆ nodes that each bind with six formate ligands. In model I, all formates have a μ²-bridging configuration and bind with two adjacent Zr atoms. From model I to model II, three of the six μ²-bridging formates dissociate one leg from one of the two Zr atoms and bind, instead, with μ-hydroxo groups on the node through hydrogen bonding. From model II to model III, two of the remaining three μ²-bridging formates dissociate one leg from one of the two Zr atoms and extend into the pore, forming formates with κ¹ configuration. From model I to model IV, one formate is replaced with a hydroxo ligand.

MOF-808 models I–IV have calculated lattice constants ranging from 35.033 to 35.282 Å. To verify our choice of models, we performed a higher level calculation for model I which has averaged lattice constants of 35.282 Å. The higher level calculation uses PBEsol^15,16 density functional with D3¹⁷ dispersion correction and up to 500 eV plane-wave basis set in the Vienna Ab initio Simulation Package (VASP).^33–36 Due to the very large computational cost, we did not perform the higher level calculation for all models. This high-level calculation decreased the average lattice constant of model I from 35.282 to 35.182 Å, which is lower than the largest experimentally measured lattice constant, suggesting that model I is a suitable candidate for representing the MOF-808 structure. We also note that the experimentally measured lattice constants might be affected by residual water within the MOF-808 framework.^11,31 Therefore, we examined the effect of water on physical properties of MOF-808. Several water-containing models were constructed by adding water molecules to the tetrahedral cage of model I and then re-optimized using DFT methods in CP2K. Structural properties, including lattice constants and BET surface areas were calculated (see Fig. 4). Importantly, both properties decrease as the number of water molecules increase within the tetrahedral cage of MOF-808.

Fig. 4 Changes in lattice constants and BET surface area as the number of water molecule increases (x-axis) within the tetrahedral pore of the MOF-808 model I.

To determine the impact of MOF model coordination structure on physical properties, we simulated both N₂ and Ar adsorption isotherms using models I–IV. A survey of the literature revealed several measurements of N₂ adsorption at 77 K for MOF-808^13,37–40 and one measurement of Ar adsorption at 87 K.⁴¹ Surprisingly, there is a relatively large deviation in experimentally measured maximum uptakes among different studies. Specifically, the N₂ uptake ranges from ∼350 cm³ g⁻¹ to ∼600 cm³ g⁻¹ at p/p₀ = 0.3. It is possible that the experimentally observed deviation in adsorption capacities is due to synthetic imperfections. These imperfections are potentially the result of missing linkers/nodes, a mixture of different coordination structures, or other structural deformations.^42–44 Our simulated N₂ adsorption at 77 K using model I with one missing linker per primitive unit cell has slightly higher (8 cm³ g⁻¹) N₂ uptake at p/p₀ = 0.3 than that simulated using pristine model I.

Fig. 5 shows the results of simulated isotherms for the MOF-808 models I–IV. Data points are simulated mainly at the low p/p₀ region to obtain key information, such as the maximum uptake.

Fig. 5 Calculated N₂ (top panel) and Ar (bottom panel) adsorption isotherms of the MOF-808 model I (black), model II (blue), model III (yellow), and model IV (red). The black curves are fitted to the model I uptakes.

Good agreement was achieved between simulated and experimentally measured adsorption isotherms. For the N₂ adsorption, the experimentally measured value varies between 350 and 600 cm³ (STP) g⁻¹. Our calculated uptakes of 519, 502, 487, and 536 cm³(STP) g⁻¹ at p/p₀ = 0.3 were achieved for models I–IV, respectively. Interestingly, we noticed a decrease in N₂ uptake when the formate ligands changed from binding with the Zr₆ node in μ²-bridging configuration to the partially dissociated κ¹ configuration. We also noticed an increase in the uptake when a hydroxo ligand replaces a formate. Note, the observed trend also holds in per mole basis with uptakes of 2 832 119, 2 741 040, 2 654 747, and 2 862 110 cm³ (STP) mol⁻¹ achieved for models I–IV, respectively. The experimentally measured Ar uptake for MOF-808 is about 560 cm³ (STP) g⁻¹ at p/p₀ = 0.3.⁴¹ For the Ar adsorption using models I–IV, we obtained uptakes of 653, 632, 619 and 675 cm³ (STP) g⁻¹, respectively, at p/p₀ = 0.3. The same trend was observed for Ar adsorption using different coordination structures as for N₂ adsorption.

We examined the effect of water content within MOF-808 on adsorption. The water-containing models were constructed by adding water molecules to the tetrahedral cage of model I. The calculated Ar uptakes at p/p₀ = 0.3 are shown in Fig. 6 and suggest that the amount of adsorbed Ar decreases with the increase of water content within MOF-808 as the water content reduces the surface area of the model (Fig. 4).

Fig. 6 Change in calculated Ar uptake as the number of water molecules increases within the tetrahedral pore of the MOF-808 model I.

3.2 Sulfated MOF-808 coordination structure and acidity

The sulfated MOF-808 can be synthesized by treating pristine MOF-808 with aqueous sulfuric acid.^11,13 The sulfated MOF-808 is then “activated” by heating it under dynamic vacuum.^11,13 Here, we refer to the sulfated and activated MOF as “S-MOF-808”. Similar to pristine MOF-808, the S-MOF-808 crystal contains Zr₆ nodes with different coordination structures. These coordination structures may differ in terms of the number and the relative position of attached ligands, in particular coordinating sulfate dianions (SO₄). For S-MOF-808, Yaghi¹¹ reported average molecular formula of Zr₆O₄(OH)₄(BTC)₂(SO₄)_2.3(OH)_1.4(OH₂)_3.1(DMF)_0.4. By considering the average molecular formula and the synthesis procedures,¹¹ we constructed nine periodic S-MOF-808 models that differ by (a) the number of sulfate ligands per node (from 1 to 3), (b) the bonding coordination configuration of the sulfate ligand to include μ²-bridging, κ²-chelating, and κ¹ bonding accompanied by hydrogen bonding, and (c) the relative positions of sulfate and hydroxo ligands.

S-MOF-808 models were constructed with two approaches. First, alterations were made directly to the optimized MOF-808 models according to the experimentally measured average molecular formula.^11,13 We altered the pristine MOF-808 model using various combinations of ligands to capture potential coordination structures that could exist around Zr₆ node. The second approach to construct S-MOF-808 models involved using experimental XRD structures. The advantage of this latter approach is that it retains possible structural nuances that can be potentially be lost as artifacts in the cell optimization procedure. A detailed list of studied models and their optimized coordinates is provided in the ESI.† Because of the activated form of S-MOF-808 there are no direct adsorption measurements to compare coordination structure models to, and therefore, rather than analyzing a single coordination structure, we compared all calculated coordination structures and their acidity by calculating proton affinities.

To begin, we first analyzed the acidity of coordination structure model A that was proposed by Yaghi and Head-Gordon.¹³ As shown in Fig. 7, this model contains two κ² chelating sulfate ligands, two hydroxo ligands, and six aqua ligands per Zr₆ node. This model accounts the experimentally measured average molecular formula.^12,13

Fig. 7 Unit cell of model A and one Zr₆ node showing two sulfate, two hydroxo, and six aqua ligands. Aqua-protons are labeled H1 to H12 and μ-hydroxo-protons are labeled μH1 to μH4 (μH1 is hidden behind the cluster). BTC linkers have been removed for clarity. O atoms are shown in red, S in yellow, Zr in cyan, C in grey, and H in white.

To determine the acidity, we calculated the proton affinities (PAs) for all hydrogens of the Zr₆ node in model A. The PA measures the energy requirement to separate the O–H bond into a conjugate anion-proton pair (defined in eqn (1)). This type of analysis has been applied to both organic and inorganic acids,^45,46 zeolites,^47,48 and MOFs.⁴⁹ In our procedure, we calculated proton affinities using cluster models truncated from fully optimized periodic S-MOF-808 models (see Computational methods section).

PA = E_MOF⁻ + E_H⁺ − E_MOF–H (1)

We also explored ensemble-averaged proton affinities for hydrogens in model A. For a given hydrogen, the ensemble-averaged proton affinity (〈PA〉) takes into account all the possible binding locations on the conjugate anion that are in thermodynamic equilibrium.⁵⁰ The PA as defined in eqn (1) does not consider equilibrium proton binding locations; however, is computationally less demanding than 〈PA〉 because the PA only requires the optimization of one MOF–H complex while the 〈PA〉 requires additional optimization of all structures. For two tested protons in model A, the computed 〈PA〉 values are within 0.6 kcal mol⁻¹ to the PAs computed using eqn (1). This very small energy difference gave us confidence to generally use PA values rather than ensemble-averaged values. We have also compared 〈PA〉 values with PA values for MOF-808 model XIV and also found nearly identical values (see the ESI†).

For model A, we computed 16 proton affinities in total for the deprotonation of 12 protons from six aqua ligands and four protons from four μ-hydroxo groups. Fig. 8 plots the calculated PAs for model A and categorizes them into five kcal mol⁻¹ ranges. This revealed that there is an extremely wide range of acidities just for this single coordination structure model. This greater than 35 kcal mol⁻¹ range can be roughly equated to 10²⁵ order magnitude in acidity range. The PAs for the twelve aqua protons range from 299 to 336 kcal mol⁻¹. The PAs for the four μ-hydroxo-protons range from 311 to 326 kcal mol⁻¹. The average value of all the PAs for model A is 316 kcal mol⁻¹.

Fig. 8 Distribution of proton affinities (PA, kcal mol⁻¹) of model A. Protons are labelled following Fig. 7. In each column PA increases from bottom to top of bar graph.

The two protons from the same aqua ligand gave similar PAs. For the six aqua ligands in model A, the PA difference between two protons from the same aqua ligand is on average 1 kcal mol⁻¹ and at most 4 kcal mol⁻¹. This similarity in proton affinities for both hydrogens of an aqua ligand is perhaps surprising given that the 3D structural and bonding environment are different. But more importantly, this observation impacts the interpretation of the previous computational assessment of model A where the superacid label of this MOF was proposed to arise from a single proton of one aqua ligand and its hydrogen bond with a neighboring sulfate ligand.¹³

Fig. 9 shows an aqua ligand with two protons (H1 and H2, also see Fig. 7) that have a different chemical environment. H1 was previously assigned to be highly acidic due to the hydrogen bond interaction with the neighboring sulfate ligand. H2 is not involved in a hydrogen bonding interaction. The interaction between atom H1 and the sulfate ligand does induce an elongated O1–H1 bond. We computed the O1–H1 bond length of 1.006 Å and O1–H2 bond length of 0.968 Å. However, despite the difference in their surrounding chemical environment, and the difference in O–H bond lengths, proton H1 and H2 have identical PA values of 318 kcal mol⁻¹. Therefore, the water–sulfate hydrogen bond interaction highlighted in Fig. 9 is not the source of superacidity in S-MOF-808. Just as important, this aqua ligand does not possess the most acidic hydrogens in model A. Instead, the most acidic protons (H9 and H10, Fig. 7) come from an aqua ligand that does not interact with sulfate ligands and bind with a Zr atom that is coordinating with another aqua ligand.

Fig. 9 The model A cluster with a highlighted region showing the interaction between adjacent aqua and sulfate ligands. Bond lengths of O1–H1 and O1–H2 are 1.006 and 0.968 Å, respectively. C atoms are shown in grey, H in white, O in red, sulfate in yellow, and Zr in cyan.

With the observation that the most acidic hydrogen in model A is not apparent from intra-node hydrogen bonding, we wondered if conventional analysis of O–H bond stretching frequencies, ν(OH), would correlate with the calculated proton affinities and begin to provide an origin of acidity. O–H stretching frequencies correlate well with PA in zeolites⁵¹ and some other crystalline materials,^52,53 and is a commonly used acidity descriptor. However, ν(OH) values fail to capture the relative proton acidity for S-MOF-808 and hydrated MOF-808. We calculated ν(OH) values, O–H bond lengths, and PAs for three aqua ligands from model A and six aqua ligands from fully hydrated MOF-808 (details are given in ESI†). As shown in Fig. 10(a), there is a clear linear correlation (R² = 0.996) between ν(OH) and the O–H bond length. However, Fig. 10(b) shows that there is no correlation between ν(OH) and PA.

Fig. 10 Calculated O–H stretching frequencies plotted versus (a) O–H bond length and (b) proton affinity.

With the inability of vibrational frequencies to correlate with acidity, we desired to directly analyze S-MOF-808 conjugate base relative stabilities that control acidity. The deprotonation of a proton from an aqua ligand on the Zr₆ node (Zr₆–O(H)H) results in Zr-node stabilized hydroxide ligand ([Zr₆–OH]⁻). The stability of this Zr–hydroxide conjugate base can be expressed in terms of the heterolytic Zr–O bond energy D. The heterolytic dissociation results in the dehydrated Zr₆ node with one less aqua ligands than the starting Zr₆ node as well as hydroxide (eqn (2)). E_Zr6, E_OH⁻, and E_Zr6OH⁻ are the energies of the dehydrated Zr₆ node, OH⁻, and the conjugate base ([Zr₆–OH]⁻). The dehydrated Zr₆ node and the conjugate base are not optimized for the calculation of D.

D = E_Zr6 + E_OH⁻ − E_Zr6OH⁻ (2)

For model A, the linear correlation between each PA and the corresponding heterolytic dissociation energy is displayed in Fig. 11. For protons H1 and H2 shown in Fig. 9, the heterolytic Zr–O dissociation energies of their respective conjugate bases are 103 and 105 kcal mol⁻¹, which explains their essentially identical PA values. This inverse relationship between proton acidity and conjugate base stability has been further verified using other S-MOF-808 models which are summarized in ESI.†

Fig. 11 Linear correlation plot between proton affinities (kcal mol⁻¹) and heterolytic Zr–O bond dissociation energies (kcal mol⁻¹) for hydrogens in model A. H1 and H2 are the hydrogen atoms shown in Fig. 9.

As expected, the relative stability of the conjugate base ([Zr–OH]⁻) is significantly impacted by the coordination ligands around the Zr atom. As described earlier, for each Zr atom of the Zr₆ node, there are two open coordination sites, and these coordination sites can be occupied by aqua, hydroxo, or sulfate ligands. Here, we compare two groups of sub-structures in model A; one group includes substructures where a Zr atom coordinates with two aqua ligands (H₂O–Zr–OH₂) and the other group includes substructures where a Zr atom coordinates with one aqua and one hydroxo ligand (H₂O–Zr–OH). The PAs associated with aqua-protons in H₂O–Zr–OH₂ range from 299 to 318 kcal mol⁻¹ and are always lower than those with the H₂O–Zr–OH motif that range from 319 to 336 kcal mol⁻¹. Indeed, these different ranges are not surprising given that a coordination structure motif of H₂O–Zr–OH₂ contains a more electron-deficient Zr center capable of better stabilizing the resulting hydroxide after proton loss (i.e. [H₂O–Zr–OH]⁻) and this is in contrast to the H₂O–Zr–OH motif that would result in a less stable [HO–Zr–OH]⁻ conjugate base.

With the wide range of acidities found in model A, and the identification of hydrogens more acidic than previously known, we wanted to determine if S-MOF-808 should be considered a superacid. Therefore, we decided to compare the proton affinities for model A to zeolite acids and known superacids. For zeolites MFI, BEA, FER, MOR, CHA, and FAU, PAs of 287.0 ± 2.6 kcal mol⁻¹ computed using periodic DFT methods were reported.⁵⁰ Superacids FSO₃SbF₅H (264), HSbF₆ (265), HC₅(CN)₅ (262), HAlCl₄ (266), HAlBr₄ (265) have lower PAs (given in parentheses in kcal mol⁻¹) than PAs of zeolites.⁴⁵ Model A, as discussed earlier, has an averaged PA of 316 kcal mol⁻¹ and lowest PA at 299 kcal mol⁻¹. Because model A's PAs are significantly higher than that of both zeolites and traditional molecular superacids, it is unlikely that model A provides superacidic protons. Note that the PAs reported for zeolites and superacids were computed using different DFT functionals and basis sets than ours; however, despite the differences in methods, we believe these PAs are reasonable to make a qualitative comparison with S-MOF-808.

As we discussed earlier, a S-MOF-808 crystal contains different coordination structures, and therefore, proton affinity based on a single coordination structure is insufficient to give a full picture of the acidity of S-MOF-808. First, we compare models with different number of sulfate groups, namely models A, B, and C. Model A, as discussed earlier, has two sulfate ligands per Zr₆ node. Models B and C have coordination structures that are based on model A. Model B has one sulfate ligand per node and was constructed by replacing one sulfate ligand of every node in model A with an aqua and a hydroxo ligand. Model C has three sulfate ligands per node and was constructed by replacing an aqua and a hydroxo ligand per node in model A with a sulfate group. Both models B and C were optimized with periodic boundary conditions, and their PAs were evaluated using corresponding cluster models. The choice of number of sulfate groups is based on the experimentally measured molecular formula of S-MOF-808 which suggests an average of 2.3 sulfate ligands per node.¹¹ We observed an increase in average acidity of S-MOF-808 models with the increase of sulfate groups per node. Specifically, models B, A, C have averaged PAs of 327, 316, 309 kcal mol⁻¹, respectively, which indicate progressively increasing acidity. Among the six neutral models we have studied, model J, which contains three sulfate groups, has the lowest averaged PA of 306 kcal mol⁻¹ (see the ESI† for comparison of all models), and therefore the highest acidity; Model B has the highest averaged PA of 327 kcal mol⁻¹ and the lowest acidity. Interestingly, model K, which contains two bisulfate (SO₄H) groups, has averaged PA of 309 kcal mol⁻¹. The lack of acidity in model K is caused by the increased amount of hydroxo groups on the Zr₆ node due to charge balance, resulting in the loss of electron deficient Zr centers. We noticed that none of the neutral coordination structures give PA values that would suggest superacidity compared to classical superacids.

4. Conclusions

MOF-808 is a Zr-based MOF with relatively high chemical and mechanical stability and has often been for applications in chemical adsorption and gas separation. In this study, we applied periodic and cluster DFT calculations to determine the impact of MOF-808 and sulfated MOF-808 coordination structures on physical and chemical properties. For MOF-808, many calculated coordination structures were compared to experimental lattice constants. Our calculated lattice constants for multiple coordination structures, and N₂ and Ar adsorption isotherms, agree very well with experimental values.

For the sulfated MOF-808 that was previously proposed to be superacidic, we determined the impact of coordination structure on acidity. Surprisingly, our results based on proton affinities suggest a large distribution of acid site strength equivalent to 10²⁵ orders of magnitude depending on the coordination structure with only few sites having high acidity. Furthermore, the data suggests that vibrational frequencies and other property-based methods for determining relative acidity fail. Our analysis revealed that the acidity of sulfated MOF-808 is unlikely controlled by a sulfate–water hydrogen bond that was previously proposed. Instead, we showed that there is a strong, likely causal, correlation between conjugate base stability and proton affinity to rationalize acidity of protons in a single coordination structure model as well as across all coordination structures.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Brigham Young University and the Office of Research Computing for computational resources. We thank Phillips 66 Company for financial support.

Notes and references

1. Y. G. Chung, E. Haldoupis, B. J. Bucior, M. Haranczyk, S. Lee, H. Zhang, K. D. Vogiatzis, M. Milisavljevic, S. Ling, J. S. Camp, B. Slater, J. I. Siepmann, D. S. Sholl and R. Q. Snurr, J. Chem. Eng. Data, 2019, 64, 5985–5998.
2. S. Yuan, J.-S. Qin, C. T. Lollar and H.-C. Zhou, ACS Cent. Sci., 2018, 4, 440–450.
3. Z. Wang and S. M. Cohen, Chem. Soc. Rev., 2009, 38, 1315–1329.
4. Y. Bai, Y. Dou, L.-H. Xie, W. Rutledge, J.-R. Li and H.-C. Zhou, Chem. Soc. Rev., 2016, 45, 2327–2367.
5. A. Corma, H. García and F. X. L. i Xamena, Chem. Rev., 2020, 110, 4606–4655.
6. A. U. Czaja, N. Trukhan and U. Müller, Chem. Soc. Rev., 2009, 38, 1284–1293.
7. M. D. Allendorf and V. Stavila, CrystEngComm, 2015, 17, 229–246.
8. L. Chen, X. Zhang, X. Cheng, Z. Xie, Q. Kuang and L. Zheng, Nanoscale Adv., 2020, 2, 2628–2647.
9. J. L. Mancuso, A. M. Mroz, K. N. Le and C. H. Hendon, Chem. Rev., 2020, 120, 8641–8715.
10. N. Planas, J. E. Mondloch, S. Tussupbayev, J. Borycz, L. Gagliardi, J. T. Hupp, O. K. Farha and C. J. Cramer, J. Phys. Chem. Lett., 2014, 5, 3716–3723.
11. J. Jiang, F. Gándara, Y.-B. Zhang, K. Na, O. M. Yaghi and W. G. Klemperer, J. Am. Chem. Soc., 2014, 136, 12844–12847.
12. H. Furukawa, F. Gándara, Y.-B. Zhang, J. Jiang, W. L. Queen, M. R. Hudson and O. M. Yaghi, J. Am. Chem. Soc., 2014, 136, 4369–4381.
13. C. A. Trickett, T. M. Osborn Popp, J. Su, C. Yan, J. Weisberg, A. Huq, P. Urban, J. Jiang, M. J. Kalmutzki, Q. Liu, J. Baek, M. P. Head-Gordon, G. A. Somorjai, J. A. Reimer and O. M. Yaghi, Nat. Chem., 2019, 11, 170–176.
14. J. Hutter, M. Iannuzzi, F. Schifmann and J. VandeVondele, Rev. Comput. Mol. Sci., 2014, 4, 15–25.
15. L. A. Constantin, J. P. Perdew and J. M. Pitarke, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 075126.
16. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
17. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
18. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
19. M. Ropo, K. Kokko and L. Vitos, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 195445.
20. G. I. Csonka, J. P. Perdew, A. Ruzsinszky, P. H. T. Philipsen, S. Lebègue, J. Paier, O. A. Vydrov and J. G. Ángyán, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 155107.
21. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 (Revision B.01), Gaussian, Inc., Wallingford, CT, 2016.
22. Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101.
23. D. Andrae, U. Häußermann, M. Dolg, H. Stoll and H. Preuß, Theor. Chim. Acta, 1990, 77, 123–141.
24. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
25. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
26. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2009, 11, 10757–10816.
27. D. Dubbeldam, S. Calero, D. E. Ellis and R. Q. Snurr, Mol. Simul., 2016, 42, 81–101.
28. J. J. Potoff and I. J. Siepmann, AIChE J., 2001, 47, 1676–1682.
29. P. P. Ewald, Ann. Phys., 1921, 369, 253–287.
30. C. E. Wilmer, K. C. Kim and R. Q. Snurr, J. Phys. Chem. Lett., 2012, 3, 2506–2511.
31. W. Zhang, A. Bu, Q. Ji, L. Min, S. Zhao, Y. Wang and J. Chen, ACS Appl. Mater. Interfaces, 2019, 11, 33931–33940.
32. S. Lin, Y. Zhao, J. K. Bediako, C.-W. Cho, A. K. Sarkar, C.-R. Lim and Y.-S. Yun, Chem. Eng. J., 2019, 362, 280–286.
33. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
34. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
35. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–610.
36. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
37. R. C. Klet, Y. Liu, T. C. Wang, J. T. Hupp and O. K. Farha, J. Mater. Chem. A, 2016, 4, 1479–1485.
38. J. Xu, J. Liu, Z. Li, X. Wang and Z. Wang, J. Mater. Sci., 2019, 54, 12911–12924.
39. K. M. Choi, H. M. Jeong, J. H. Park, Y.-B. Zhang, J. K. Kang and O. M. Yaghi, ACS Nano, 2014, 8, 7451–7457.
40. J. Zhang, S. B. Peh, J. Wang, Y. Du, S. Xi, J. Dong, A. Karmakar, Y. Ying, Y. Wang and D. Zhao, Chem. Commun., 2019, 55, 4727–4730.
41. K. Healey, W. Liang, P. D. Southon, T. L. Church and D. M. D'Alessandro, J. Mater. Chem. A, 2016, 4, 10816–10819.
42. G. C. Shearer, S. Chavan, J. Ethiraj, J. G. Vitillo, S. Svelle, U. Olsbye, C. Lamberti, S. Bordiga and K. P. Lillerud, Chem. Mater., 2014, 26, 4068–4071.
43. G. C. Shearer, S. Chavan, S. Bordiga, S. Svelle, U. Olsbye and K. P. Lillerud, Chem. Mater., 2016, 28, 3749–3761.
44. D. S. Sholl and R. P. Lively, J. Phys. Chem. Lett., 2015, 6, 3437–3444.
45. I. A. Koppel, P. Burk, I. Koppel, I. Leito, T. Sonoda and M. Mishima, J. Am. Chem. Soc., 2000, 122, 5114–5124.
46. M. M. Meyer, X.-B. Wang, C. A. Reed, L.-S. Wang and S. R. Kass, J. Am. Chem. Soc., 2009, 131, 18050–18051.
47. K.-P. Schröder, J. Sauer, M. Leslie, C. Richard, A. Catlow and J. M. Thomas, Chem. Phys. Lett., 1992, 188, 320–325.
48. J. Sauer, Faraday Discuss., 2016, 188, 227–234.
49. M. D. de Mello, G. Kumar, T. Tabassum, S. K. Jain, T.-H. Chen, S. Caratzoulas, X. Li, D. G. Vlachos, S.-I. Han, S. L. Scott, P. Dauenhauer and M. Tsapatsis, Angew. Chem., 2020, 132, 13362–13368.
50. A. J. Jones and E. Iglesia, ACS Catal., 2015, 5, 5741–5755.
51. S. Bordiga, C. Lamberti, F. Bonino, A. Travert and F. Thibault-Starzyk, Chem. Soc. Rev., 2015, 44, 7262–7341.
52. E. E. Platero, M. P. Mentruit, C. O. Arean and A. Zecchina, J. Catal., 1996, 162, 268–276.
53. M. Mihaylov, S. Andonova, K. Chakarova, A. Vimont, E. Ivanova, N. Drenchev and K. Hadjiivanov, Phys. Chem. Chem. Phys., 2015, 17, 24304–24314.

---

Footnote

† Electronic supplementary information (ESI) available: Periodic MOF-808 models and S-MOF-808 models (CIF); tabulated summary of MOF-808 and S-MOF-808 models, additional discussion on ensemble-averaged proton affinity, and additional correlation diagram between proton affinity and Zr–O bond dissociation energy (PDF). See DOI: 10.1039/d1ma00330e

---