Bo
Yang
*a,
Joshua I.
Wheeler
a,
Brett
Sorensen
a,
Robert
Steagall
a,
Taylor
Nielson
a,
Jianhua
Yao
b,
Jose
Mendez-Arroyo
*b and
Daniel H.
Ess
*a
aDepartment of Chemistry and Biochemistry, Brigham Young University, Provo, Utah 84602, USA. E-mail: dhe@chem.byu.edu; bo.yang@chem.byu.edu
bPhillips 66, Bartlesville, Oklahoma 74003, USA. E-mail: jose.e.mendez-arroyo@p66.com
First published on 27th April 2021
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 N2 and Ar adsorption isotherms. Our simulated average N2 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 1025 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.
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.9 For example, density functional calculations in combination with experiments revealed the thermodynamically most stable proton-topology of MOF NU-1000.10
MOF-808 is a Zr-based mesoporous material that has relatively high thermal and mechanical stability,5 which prompted structure modification and multiple studies in the application of adsorption, gas separation, and catalysis. MOF-808 is constructed by linking Zr6O4(OH)4 nodes (i.e., Zr6 node) with benzene-1,3,5-tricarboxylate (BTC) linkers, as shown in Fig. 1. Each Zr6 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 H1-NMR, IR, and XRD data11–13 have revealed that ligands such as formate and water can coordinate to a Zr6 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 Zr6O4(OH)4 nodes (shown as cyan polyhedron) and benzene-1,3,5-tricarboxylate linkers. Formates are coordinated to the node in the μ2-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 Zr6 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.13 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.13
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 N2 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. N2 and Ar adsorption isotherms were computed for selected models using grand canonical Monte Carlo (GCMC) simulations, and we obtained N2 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.
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 Zr6 node. The positioning of acetate capping groups around the Zr6 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.21 We used the M06-L22 density functional with SDD23 pseudopotentials and its corresponding basis set for Zr atoms and def2-TZVP24,25 basis set for all the other atoms. The M06-L functional was chosen because of its general accuracy in transition metal chemistry.26 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 N2 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.27 N2 and Ar molecules were described by the TraPPE force field.28 MOF-808 models were described by the GenericMOFs force field27 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.29 Partial atomic charges were computed for MOF atoms using the extended charge equilibration method.30 A total of 110000 GCMC cycles were performed to converge the adsorption isotherms with the first 10000 cycles for equilibrium and the latter 100000 cycles for ensemble averages.
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 κ1 or κ2. In models where a formate coordinates with two Zr atoms (μ2-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.02331 to 35.232 Å.32 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 Zr6 nodes that each bind with six formate ligands. In model I, all formates have a μ2-bridging configuration and bind with two adjacent Zr atoms. From model I to model II, three of the six μ2-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 μ2-bridging formates dissociate one leg from one of the two Zr atoms and extend into the pore, forming formates with κ1 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 PBEsol15,16 density functional with D317 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 N2 and Ar adsorption isotherms using models I–IV. A survey of the literature revealed several measurements of N2 adsorption at 77 K for MOF-80813,37–40 and one measurement of Ar adsorption at 87 K.41 Surprisingly, there is a relatively large deviation in experimentally measured maximum uptakes among different studies. Specifically, the N2 uptake ranges from ∼350 cm3 g−1 to ∼600 cm3 g−1 at p/p0 = 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 N2 adsorption at 77 K using model I with one missing linker per primitive unit cell has slightly higher (8 cm3 g−1) N2 uptake at p/p0 = 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/p0 region to obtain key information, such as the maximum uptake.
Good agreement was achieved between simulated and experimentally measured adsorption isotherms. For the N2 adsorption, the experimentally measured value varies between 350 and 600 cm3 (STP) g−1. Our calculated uptakes of 519, 502, 487, and 536 cm3(STP) g−1 at p/p0 = 0.3 were achieved for models I–IV, respectively. Interestingly, we noticed a decrease in N2 uptake when the formate ligands changed from binding with the Zr6 node in μ2-bridging configuration to the partially dissociated κ1 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 2832119, 2741040, 2654747, and 2862110 cm3 (STP) mol−1 achieved for models I–IV, respectively. The experimentally measured Ar uptake for MOF-808 is about 560 cm3 (STP) g−1 at p/p0 = 0.3.41 For the Ar adsorption using models I–IV, we obtained uptakes of 653, 632, 619 and 675 cm3 (STP) g−1, respectively, at p/p0 = 0.3. The same trend was observed for Ar adsorption using different coordination structures as for N2 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/p0 = 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. |
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 Zr6 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.13 As shown in Fig. 7, this model contains two κ2 chelating sulfate ligands, two hydroxo ligands, and six aqua ligands per Zr6 node. This model accounts the experimentally measured average molecular formula.12,13
To determine the acidity, we calculated the proton affinities (PAs) for all hydrogens of the Zr6 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.49 In our procedure, we calculated proton affinities using cluster models truncated from fully optimized periodic S-MOF-808 models (see Computational methods section).
PA = EMOF− + EH+ − EMOF–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.50 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−1 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−1 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−1 range can be roughly equated to 1025 order magnitude in acidity range. The PAs for the twelve aqua protons range from 299 to 336 kcal mol−1. The PAs for the four μ-hydroxo-protons range from 311 to 326 kcal mol−1. The average value of all the PAs for model A is 316 kcal mol−1.
Fig. 8 Distribution of proton affinities (PA, kcal mol−1) 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−1 and at most 4 kcal mol−1. 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.13
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−1. 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.
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 zeolites51 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 (R2 = 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 Zr6 node (Zr6–O(H)H) results in Zr-node stabilized hydroxide ligand ([Zr6–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 Zr6 node with one less aqua ligands than the starting Zr6 node as well as hydroxide (eqn (2)). EZr6, EOH−, and EZr6OH− are the energies of the dehydrated Zr6 node, OH−, and the conjugate base ([Zr6–OH]−). The dehydrated Zr6 node and the conjugate base are not optimized for the calculation of D.
D = EZr6 + EOH− − EZr6OH− | (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−1, 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−1) and heterolytic Zr–O bond dissociation energies (kcal mol−1) 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 Zr6 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 (H2O–Zr–OH2) and the other group includes substructures where a Zr atom coordinates with one aqua and one hydroxo ligand (H2O–Zr–OH). The PAs associated with aqua-protons in H2O–Zr–OH2 range from 299 to 318 kcal mol−1 and are always lower than those with the H2O–Zr–OH motif that range from 319 to 336 kcal mol−1. Indeed, these different ranges are not surprising given that a coordination structure motif of H2O–Zr–OH2 contains a more electron-deficient Zr center capable of better stabilizing the resulting hydroxide after proton loss (i.e. [H2O–Zr–OH]−) and this is in contrast to the H2O–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−1 computed using periodic DFT methods were reported.50 Superacids FSO3SbF5H (264), HSbF6 (265), HC5(CN)5 (262), HAlCl4 (266), HAlBr4 (265) have lower PAs (given in parentheses in kcal mol−1) than PAs of zeolites.45 Model A, as discussed earlier, has an averaged PA of 316 kcal mol−1 and lowest PA at 299 kcal mol−1. 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 Zr6 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.11 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−1, 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−1 (see the ESI† for comparison of all models), and therefore the highest acidity; Model B has the highest averaged PA of 327 kcal mol−1 and the lowest acidity. Interestingly, model K, which contains two bisulfate (SO4H) groups, has averaged PA of 309 kcal mol−1. The lack of acidity in model K is caused by the increased amount of hydroxo groups on the Zr6 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.
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 1025 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.
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 |
This journal is © The Royal Society of Chemistry 2021 |