Improving the Mechanical Stability of Metal-Organic Frameworks Using Chemical Caryatids.

Metal-organic frameworks (MOFs) have emerged as versatile materials for applications ranging from gas separation and storage, catalysis, and sensing. The attractive feature of MOFs is that, by changing the ligand and/or metal, they can be chemically tuned to perform optimally for a given application. In most, if not all, of these applications one also needs a material that has a sufficient mechanical stability, but our understanding of how changes in the chemical structure influence mechanical stability is limited. In this work, we rationalize how the mechanical properties of MOFs are related to framework bonding topology and ligand structure. We illustrate that the functional groups on the organic ligands can either enhance the mechanical stability through formation of a secondary network of nonbonded interactions or soften the material by destabilizing the bonded network of a MOF. In addition, we show that synergistic effect of the bonding network of the material and the secondary network is required to achieve optimal mechanical stability of a MOF. The developed molecular insights in this work can be used for systematic improvement of the mechanical stability of the materials by careful selection of the functional groups.

Introduction

Like any other material, metal–organic frameworks (MOFs), as an important class of porous materials with large diversity of pore shapes and sizes, and rich chemical functionalities must pass the stability criteria to be used in most practical applications.[1−3] Despite having superior performance for many applications, MOFs are vulnerable with respect to stability compared to the competing materials. For instance, due to the relatively weak metal–ligand coordination bonds, many MOFs are chemically unstable and have low endurance in different types of chemicals environments, e.g., acidic or basic environment.[3] Significant progress has been made in developing MOFs that are chemically stable, e.g., zirconium-based MOFs.[4] Since applications of MOFs often involve repetitive, cyclic temperature and pressure variations and capillary forces exerted by guest molecules, sufficient mechanical stability is of equal importance.[5,6] The mechanical stability for porous materials measures the stiffness of a material to withstand its pore size and structure under mechanical load. Clearly, deformations due to external pressure will disrupt pore shape and size, resulting in significantly reduced performance. In this study, we focus on strategies to improve the mechanical stability of a particular MOF.

The mechanical properties of materials vary by several orders of magnitude with changing atomic composition and/or crystal structure.[7−9] As the mechanical stiffness, i.e., modulus of elasticity, typically scales quadratically with the density,[10] mechanical stability is of particular importance for applications of low-density materials, such as MOFs.[6,11,12] For these materials special strategies are often required to improve their mechanical stability. Often these strategies are inspired by nature (e.g., wood and bones[13,14]) and involve a fractal and hierarchical design to make highly connected materials over multiple length scales.[15−17] Indeed, improving the mechanical stability of MOFs by tuning the chemistry has become an important focus of attention.[3,18−20] In analogy to the concept of high connectivity of the hierarchical design of materials, it has been shown that the MOFs with high degrees of framework interconnectivity, i.e., high coordination number of metal nodes, have improved the mechanical stability.[18,21] However, a particular MOF can not be easily replaced for all applications, and therefore, Kapustin et al. developed a strategy to retrofit a particular MOF by adding additional ligands to the framework.[20] This strategy is robust but limited to the MOFs that permit ligand installation.[22,23] In both cases, the mechanical stability is improved by increasing the connectivity of the bonding topology.

In this work, we explore the option of decorating the organic ligands of a MOF with functional groups. The significant progress in computational material science in in silico generation of MOFs[24,25] and reliable prediction of their mechanical properties[11,26] permits studying a large and diverse set of materials to extract structure–property relationships to design materials with enhanced mechanical stability. We show that the nonbonded interactions play an important role in the stiffness of the materials, and therefore, strategically placed functional groups can introduce extra framework connectivity via nonbonded interactions. This secondary network of nonbonded interactions can enhance the mechanical stability of the framework considerably. We use the term “chemical caryatids” for those functional groups that are contributing in carrying the mechanical load applied to the material. In addition, we show that the optimum mechanical stability of a MOF framework is obtained by the cooperative effect of the primary network, determined by the bonding topology, and the secondary network, which is governed by the nonbonded interactions.

Results and Discussion

In this work, we focus on zeolitic imidazolate frameworks (ZIFs), which are a special class of MOFs composed of four coordinated metals, typically zinc, with imidazolate (IM) derivative ligands. ZIFs are an ideal case study for our work because they all have the same coordination environment, but diverse bonding topologies and functional groups.[27,28] This allows us to focus on the effects of bonding topology and functional groups on the mechanical properties, while keeping coordination environment fixed, i.e., keeping the same metal node. In addition, because of the pioneering work of Cheetham and co-workers, ZIFs are among the very few MOFs for which systematic research has been done on their mechanical stability.[5,6,12] To characterize the mechanical properties of ZIFs, Cheetham and co-workers used nanoindentation to measure the Young’s modulus, i.e., the resistance of materials to the tensile stress.[6] These and related studies concluded that for these materials the mechanical properties can be described with the low density–stiffness correlation.[6,29−32] As these experiments require sufficiently large single crystals, the number of studied structures is relatively small compared to the total number of possible ZIFs. In this work, we expand the studied materials to, in addition to the known ZIF structures, a large set of in silico constructed materials using 50 different zeolite topologies[33] with four types of ligands. Such a large set of materials allows us to cover a representative range of bonding topologies and functional groups. The ligands used for in silico construction of materials include the commonly[27,34] used derivatives of IM shown in Figure .

Figure 1

Four different ligands used to construct hypothetical materials. (a) IM = imidazolate, (b) mIM = 2-methylimidazolate, (c) dcIM = dichloroimidazolate, and (d) nIM = 2-nitroimidazolate.

Four different ligands used to construct hypothetical materials. (a) IM = n class="Chemical">imidazolate, (b) mIM = 2-methylimidazolate, (c) dcIM = dichloroimidazolate, and (d) nIM = 2-nitroimidazolate.

Theoretically, mechanical properties of materials are described by their stiffness matrix.[35] Young’s and other moduli of elasticity, including bulk and shear modulus, which characterize material’s resistance to hydrostatic pressure and shear stress, respectively, can be extracted from the stiffness matrix. Since the mechanical properties of the materials in our study do not involve the breaking/formation of chemical bonds and other quantum effects, we used an approach based on a classical force field to compute the stiffness matrix for each material. The reliability of our force field is evaluated by comparison with the experimental and ab initio calculated values of Young’s modulus reported in the literature. Figure shows a comparable agreement between the ab initio and force field results with the experimental data, supporting the conclusion of our previous work that these classical force fields yield sufficiently reliable data on the mechanical properties of these materials.[36]

Figure 2

Young’s modulus versus density; for each material we plot the value along each of the three lattice principle axes. The filled markers with unique marker for each structure are used for those structures we can compare our force field (FF) with experimental (exp) or ab initio density functional theory (DFT) calculations, with the markers representing as follows: ◆, ZIF-8;[6,12] ●, ZIF-20;[6] ▼, ZIF-68;[6] ★, ZIF-4;[6,30,31] ■, ZIF-7;[6] and ▲, ZIF-zni.[6,30,32] The color coding is used to indicate the different ligands. If the density–stiffness correlation were perfectly obeyed, a principal component analysis would give a narrow cloud around the dashed line. The clouds derived from principal component analysis demonstrate the deviations for the different ligands. The complete set of data can be found in the SI.

If we focus on those materials in Figure for which experimental data are available, we observe the same low density–stiffness correlation as found experimentally.[6] However, if we include all our data, the picture becomes quite different. By expanding the chemistry and topology of ZIF structures, Figure shows large deviations from the density–stiffness correlation. Changing the underlying network topology and/or ligand can lead to larger variations in mechanical stability than changes in density, and in some cases, even reverse the trend. For instance, many ZIF structures with dcIM ligand have similar or lower stiffness in comparison to the structures in mIM and nIM ligand families, although they have higher density. A molecular-level explanation of these deviations is provided below, the understanding of which will allow us to exploit the chemical and topological features of a material to improve its mechanical stability.

The structures in Figure differ in their bonding topology and/or functional group of ligand. We introduce a computational approach to disentangle the effects of changes of the topology from changes of the ligand. To distinguish the role of the bonding topology on the mechanical properties, we first look at the mechanical stability of a simplified network of atoms composed of atomic bonding, and we refer to this network as the primary network. Several approaches have been used to define such a primary network.[21,37] Here, we define the primary network as the ZIF structure in the absence of nonbonded interactions. Since the ligands in our study only differ in their functional groups, the primary network of the structures with the same underlying network topology but different ligands are nearly identical. Hence, we expect similar mechanical properties for the structures with the same underlying network topology. Indeed, Figure a shows that all ZIFs with the same topology have similar bulk and shear modulus, and hence, superpose on each other.

Figure 3

Differentiating the contributions from bonding topology (primary network) and nonbonded interactions (secondary network) in mechanical properties. Considerable contribution from the secondary network is observed in some of the materials with functional groups. (a, b) Bulk modulus with respect to shear modulus of the materials computed without and with nonbonded interactions, respectively. (c) Bulk modulus of the materials versus the bulk modulus of them without nonbonded interactions. Dashed line represents identical properties computed with and without nonbonded interactions, i.e., no contribution from the secondary network. In all parts, each marker (open, half-filled, and filled) represents a unique underlying network topology while the colors represent the ligand.

Figure b,c shows the effects of switching on the nonbonded interactions where a large effect of functionalization on mechanical properties is observed. As there is no functional group on the IM ligand, it can be seen as bare backbone, and we see that the mechanical properties for this ligand are indeed dominated by the primary network. However, for the other ligands, functionalization can have a large effect on some topologies while on others surprisingly little. Moreover, although mIM, dcIM, and nIM exhibit observable contributions to the stiffness of ZIF structures in comparison to IM, depending on the topology one functional group might show greater enhancement. For instance, for LTA topology, mIM gives higher stiffness, while for GIS topology (ZIF-6[34,38]) dcIM has higher stiffness. Similar changes in the mechanical properties were observed experimentally for ZIFs with the same underlying net but different functional groups and were associated with the ligand–ligand interactions.[6] It is instructive to try to explain these deviations with a simple extension of the density–stiffness model. This model assumes a solid which has only nonbonded interactions, for example, a primitive cubic lattice with only nearest neighbor, (Lennard-Jones type) pairwise interactions. In this simple model, the only variable is the density-dependent nearest neighbor distance. The bulk modulus is given by the second derivative of the potential energy of the crystal with respect to isotropic deformations. The second derivative of the Lennard-Jones potential changes sign from positive to negative at ∼1.2 σ, where σ is the van der Waals radius (see the SI). As the second derivative for each pairwise interaction can be positive or negative depending on the nearest neighbor distance, the bulk modulus of this simple solid consists of a sum of positive or negative contributions, giving the well-known density–stiffness correlation. In a ZIF structure, however, there is a distribution of interatomic distances; some have a positive contribution (i.e., stiffening interactions), and some have a negative contribution (i.e., softening interactions) to bulk modulus. One can argue that this distribution depends on the topology and functional group. If we now assume that the contributions of the nonbonded interactions are independent of the contribution of the primary network, we can obtain a simple correction to the density–stiffness correlation by adding the sum of the contribution of the nonbonded interactions to the bulk modulus resulting from the primary network.

In Figure a, we plot the distribution of stiffening/softening contributions for the ZIFs with BCT topology for three different ligands. The BCT zeolite topology includes some known ZIFs, e.g., ZIF-1.[34,38] As expected, for IM, which has no functional group, this contribution is small. For the mIM and dcIM ligands, Figure a shows higher peaks in the stiffening regime which is consistent with the observed increase in mechanical stability due to functionalization. Figure b shows an example of two materials in which the distributions of the stiffening and softening contributions are nearly identical. For the BCT topology structure we observe the expected stiffening compared to the primary network. However, for AEI topology we observe only a small effect of the nonbonded interactions on the bulk modulus. This is where our simple correction to the density–stiffness correlation breaks down. This example illustrates that the contributions of nonbonded interactions and the primary network to the stiffness can be highly nonadditive. The reason for this nonadditive behavior becomes clear by introducing the concept of a secondary network.

Figure 4

Distribution of stiffening/softening nonbonded contributions for (a) structures with BCT topology and IM, mIM, and dcIM ligands, and (b) structures with mIM ligand and BCT and AEI topologies. The vertical axes represent the sum of second derivative of van der Waals (vdW) energy () plotted with respect to the interatomic distances normalized with vdW radii (R̅ = r/σ). The bulk moduli are 4.7, 9.4, and 7.1 for the BCT structures with IM, mIM, and dcIM ligands, and 6.3 and 7.0 for the AEI structures with IM and mIM ligands, respectively (values are in GPa). (c, d) Atomic representation and primary and secondary networks for the mIM ligand structures with BCT and AEI topologies, respectively. Details of ligands and metals were omitted in visualization of the primary and secondary networks for clarity. The primary net is demonstrated with red tubes and secondary net with cyan tubes; white, black, blue, and gray spheres represent H, C, N, and Zn atoms, respectively. The corresponding structures with IM ligand have the same primary net and no secondary net.

We define the secondary network by connecting pairs of atoms with nonbonded interactions that have a stiffening contribution to the bulk modulus. Figure c,d shows the primary (red tubes) and secondary (cyan tubes) networks for the BCT and AEI topologies, respectively. Both materials have a 3D percolating primary network, but the pronounced difference is in the secondary networks. For BCT the secondary network is percolating in all three dimensions, while for AEI topology it percolates only in one dimension, and there are no contributions in the other two dimensions. Inspection of the primary network of AEI topology shows that the weak spots are on the ligands while the backbone is relatively stiff. Figure d shows that the corresponding secondary network reinforces this stiff backbone, but not the links between the backbones. Hence, the secondary network is only supporting AEI topology in a direction in which the primary network is already strong. As the mechanical properties are dominated by the weakest link, we now understand why we see such a small effect of the secondary network on the mechanical properties. To have an effect, we need to add a functional group that would form a secondary network orthogonal to the current network which would significantly increase the bulk modulus. This type of synergy between the primary and secondary networks explains why some topologies show a large effect of functionalization, while for others this effect can be small.

It is interesting to apply our concept of primary and secondary networks to MOF-520-BPDC. Kapustin et al.[20] retrofitted the mechanically unstable MOF-520 by adding an additional linker to allow for its use at high pressures. This retrofitting procedure changes the underlying network topology from a fon net to a more connected skl net.[39] This improved mechanical stability can be explained in terms of changes of the primary network (see the SI). This form of topological tunability is very robust. However, it does rely on the ability to add extra linkers to support the weak spots of the primary network, which can be challenging from a chemical point of view for most materials.

Alternatively, the mechanical properties of MOFs can be tuned by creating a secondary network via ligand functionalization. The presence of such a secondary network can shed some light on the experimental observation on the amorphization of ZIFs.[40] Amorphization is directly related to the mechanical stability of these materials.[41] Cheetham and co-workers showed that ZIFs with the bare IM ligand amorphize relatively easily under pressure and heating, while the corresponding ZIFs with functionalization ligands required extreme conditions; specifically, they observed thermal amorphization only in ZIFs with the bare IM ligand.[40,42] These results are consistent with our molecular dynamics simulations (see the SI). Our analysis of the mechanical stability shows that “switching on” the secondary network in ZIF-3 and ZIF-4 improves the mechanical stability by as much as ∼80% in shear modulus of both structures, and 300% and 150% in their bulk modulus, respectively. Figure shows that for both ZIFs the functionalized structures have a secondary network that spans the entire unit cell in all three directions. Such increased mechanical stability explains why these materials are stable at conditions where the unsubstituted IM structure amorphizes.

Figure 5

(a, b) Atomic representation and the primary and secondary networks of ZIF-3 and ZIF-4 structures with mIM ligands, respectively. The corresponding structures with IM ligand have the same primary net and no secondary net. The bulk and shear moduli for ZIF-3 (ZIF-4) are 2.0 and 0.53 (3.1 and 0.80) for IM and 7.8 and 0.96 (7.6 and 1.49) for mIM structures, respectively (values are in GPa). The functional groups of the ligands form a secondary network which enhance the mechanical stability considerably. The primary net is demonstrated with red tubes and secondary net with cyan tubes; white, black, blue, and gray spheres represent H, C, N, and Zn atoms, respectively.

Conclusions

Our study shows that there are two strategies to improve the mechanical stability of a nanoporous material: modifying the primary and/or secondary network. Changing the primary network can be challenging as it requires the addition of extra linkers. In this respect the work of Kapustin et al.[20] is a remarkable, but exceptional, achievement. Functionalization of ligands to create or modify the secondary network, much like the caryatids holding up the porch of the Erechtheion on the Acropolis, might be a more generally applicable route. Our study shows that such a network, however, is only effective if it supports the weak points of the primary network.

It is interesting to envision how these results can be used from an experimental perspective. Suppose we have a particular MOF for a given application, but the mechanical stability needs to be improved. As the tools developed in this work are applicable to any MOF, we can determine the primary and secondary network of this material. If this analysis shows weak spots, a simple screening of different functional groups should give a clear prediction whether the mechanical properties of the material can be improved. As these functional groups may change the details of the pores, other computational tools should be used to ensure that these changes do not influence the performance of the modified material.

Methods

To compute the mechanical properties of a material we start with the crystal structure either from experimental or from an in silico predicted structure. The procedure of computing the mechanical properties relies on the assumption that the structure corresponds to the minimum energy configuration that is consistent with the force field used to describe the potential energy surface of the material. We developed a structural minimization procedure to efficiently obtain this minimum energy configuration for all materials. All calculations were carried out within the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular simulation package.[43] The VMD–Visualize Molecular Dynamics package was used for the structural figures and visualization of the primary and secondary networks.[44] No unexpected or unusually high safety hazards were encountered. Below we summarize the computational procedures that we have used. A more detail description can be found in the Supporting Information.

Hypothetical Material Generation

Each material was assembled with the ToBasCCo algorithm,[45] using a representative set of 50 zeolite topologies. Inputs into the program were the underlying networks, as obtained from the International Zeolite Association web site,[33] and two geometric building blocks; a 4-connected tetrahedral (Zn2+) and 2-connected imidazole type ligands. This procedure yielded 200 materials, i.e., 50 structures for each of the four types of ligands, IM, nIM, mIM, and dcIM. All the structures are available through the materials cloud Web site and the Supporting Information.

Structure Minimization

Simulated annealing algorithm was used to minimize lattice parameters and atomic sites using DREIDING force field[46] as implemented[36] in LAMMPS for all the structures. To avoid getting trapped in local minima we combined temperature annealing with expansion/relaxation cycles. The details of the algorithm and its efficiency are discussed in the Supporting Information.

Calculation of the Mechanical Properties

The moduli of elasticity were extracted from the stiffness tensor based on Voigt–Reuss–Hill averages. The stiffness tensor was calculated for the minimized structures based on the curvature of the potential energy surface with respect to lattice deformations estimated by second order polynomials (see the SI for details).