Reliably Modeling the Mechanical Stability of Rigid and Flexible Metal-Organic Frameworks.

Over the past two decades, metal-organic frameworks (MOFs) have matured from interesting academic peculiarities toward a continuously expanding class of hybrid, nanoporous materials tuned for targeted technological applications such as gas storage and heterogeneous catalysis. These oft-times crystalline materials, composed of inorganic moieties interconnected by organic ligands, can be endowed with desired structural and chemical features by judiciously functionalizing or substituting these building blocks. As a result of this reticular synthesis, MOF research is situated at the intriguing intersection between chemistry and physics, and the building block approach could pave the way toward the construction of an almost infinite number of possible crystalline structures, provided that they exhibit stability under the desired operational conditions. However, this enormous potential is largely untapped to date, as MOFs have not yet found a major breakthrough in technological applications. One of the remaining challenges for this scale-up is the densification of MOF powders, which is generally achieved by subjecting the material to a pressurization step. However, application of an external pressure may substantially alter the chemical and physical properties of the material. A reliable theoretical guidance that can presynthetically identify the most stable materials could help overcome this technological challenge. In this Account, we describe the recent research the progress on computational characterization of the mechanical stability of MOFs. So far, three complementary approaches have been proposed, focusing on different aspects of mechanical stability: (i) the Born stability criteria, (ii) the anisotropy in mechanical moduli such as the Young and shear moduli, and (iii) the pressure-versus-volume equations of state. As these three methods are grounded in distinct computational approaches, it is expected that their accuracy and efficiency will vary. To date, however, it is unclear which set of properties are suited and reliable for a given application, as a comprehensive comparison for a broad variety of MOFs is absent, impeding the widespread use of these theoretical frameworks. Herein, we fill this gap by critically assessing the performance of the three computational models on a broad set of MOFs that are representative for current applications. These materials encompass the mechanically rigid UiO-66(Zr) and MOF-5(Zn) as well as the flexible MIL-47(V) and MIL-53(Al), which undergo pressure-induced phase transitions. It is observed that the Born stability criteria and pressure-versus-volume equations of state give complementary insight into the macroscopic and microscopic origins of instability, respectively. However, interpretation of the Born stability criteria becomes increasingly difficult when less symmetric materials are considered. Moreover, pressure fluctuations during the simulations hamper their accuracy for flexible materials. In contrast, the pressure-versus-volume equations of state are determined in a thermodynamic ensemble specifically targeted to mitigate the effects of these instantaneous fluctuations, yielding more accurate results. The critical Account presented here paves the way toward a solid computational framework for an extensive presynthetic screening of MOFs to select those that are mechanically stable and can be postsynthetically densified before their use in targeted applications.

Introduction

Metal–organic frameworks (MOFs), defined as hybrid materials that are composed of inorganic moieties connected through organic ligands to form highly ordered and often crystalline materials, have emerged as important contenders in a variety of applications, including gas storage[1] and heterogeneous catalysis.[2] However, to compete with existing materials such as zeolites, reliable scale-up production methods for MOFs need to be developed.[3] One of the key challenges for this scale-up is the shaping of MOF powders, as many MOFs cannot be synthesized on a large scale as single crystals and hence need to be densified. The most popular industrial densification procedure subjects the material to elevated pressures of up to a few hundred megapascals.[3] However, as indicated in Figure a, many MOFs undergo phase transitions or become amorphous when exposed to pressures of this magnitude. As a result, the densification procedure may degrade their physical and chemical properties.[4−7] Hence, the synthesis of MOFs with elevated mechanical stability becomes increasingly important for their further industrialization.

Figure 1

(a) Critical pressures to induce phase transition or amorphization for some selected flexible and rigid MOFs. (b–d) Schematic overviews of the three computational methods to determine the mechanical stability of MOFs.

A thorough mechanical characterization of MOFs is not only crucial for a scale-up toward technological applications but is also instructive to comprehend the microscopic parameters driving mechanical flexibility and amorphization in these materials.[8] Flexible MOFs exhibit structural deformations between metastable phases of the material under influence of, among others, temperature, pressure, or guest adsorption.[9,10] This dynamic behavior may endow the material with extraordinary performance, surpassing the performance of rigid materials. These intriguing properties open new possibilities for a myriad of applications, including efficient gas storage, sensing devices, and shock absorbers.[11,12] While pressure-induced phase transitions typically take place in the pressure range of 10–100 MPa (orange data points in Figure a), also more extreme cases exist for which the long-range crystalline order is lost upon application of elevated pressures in the gigapascal regime, resulting in amorphous materials (black data points in Figure a).[13] While this process may be irreversible and inevitably leads to a deterioration of the porosity, MOF amorphization can also be beneficial, for instance in the efficient trapping of radioactive species or drugs.[14,15]

Hence, whether or not the intended MOF application directly relies on the mechanical response of the material on pressure, the stability of these materials is a crucial challenge in the further development of MOFs for practical applications.[8] This is exemplified by some outstanding experimental reviews of the mechanical stability of MOFs[16−18] as well as the advent of mercury intrusion porosimetry as a reliable technique to characterize the response of MOFs upon pressure treatment.[12,19,20] In contrast, theoretical investigations on the mechanical stability of MOFs have emerged only recently. Moreover, they are often limited in scope to one material or one specific technique, varying from valuable yet less accurate graph-based methods[21] to the more accurate and expensive models employed in this work. However, a solid theoretical framework is indispensable for the mechanical characterization of MOFs. Such a framework not only may help to identify the most stable materials, but also would provide microscopic insight into the factors influencing this stability, facilitating the synthesis of MOFs with defined functions for specific applications.[10,11,22]

In principle, one could study the stability of MOFs by performing molecular dynamics (MD) simulations under experimental conditions of pressure and temperature.[23] A versatile set of MD techniques are available, mostly relying on a force field (FF) description of the potential energy surface, approximating the interparticle interactions with predefined analytical functions.[24] However, additional technical issues may arise during these straightforward simulations. As noted in a previous contribution, MD simulations suffer from large pressure fluctuations, which have a profound effect on the MOF’s properties.[25] For instance, for the flexible MIL-53(Al), we showed that these pressure fluctuations may lead to an underestimation of the transition pressure by several tens of megapascals. As this underestimation is on the same order of magnitude as the transition pressure of the material, this approach would lead to incorrect conclusions about the relative stability of the different MIL-53(Al) phases.

To overcome these issues, more advanced thermodynamic models have been developed. Recently, three such computational methods, schematically depicted in Figure b–d, have been proposed and applied successfully on MOFs by the Coudert group and ourselves. The first approach investigates the effect of small-amplitude deformations on the stability of a given crystalline phase, resulting in a fundamental set of necessary and sufficient criteria.[26] These so-called Born stability criteria have only recently been applied to study the pressure-induced amorphization of ZIF-8 and UiO-66,[27,28] yielding predictions in close agreement with experiment. A second method aims to relate the MOF’s flexibility to the variation of the mechanical moduli, notably the Young and shear moduli, along different directions. This technique distinguishes between flexible MOFs, for which the anisotropy is larger than 100, and rigid MOFs, for which the anisotropy is smaller.[29] Finally, a third technique is based on the construction of complete pressure-versus-volume equations of state, from which the relevant mechanical properties can be determined straightforwardly.[25] These equations of state have been applied to study the flexibility[25,30,31] and rigidity[25,28] of various MOFs, showcasing their potential to study both elastic and inelastic deformations.

Despite the development of these three promising theoretical methods, they have not yet been critically assessed with respect to their mutual consistency or their limitations when applied on a broad range of materials. As a consequence, it is unclear to which extent and for which materials one can rely on these concepts, hampering their widespread use. However, if found to be reliable, these methods could spark a revolution within application-driven MOF synthesis. Such a solid and reliable computational framework would enable a presynthetic theoretical screening of MOFs, identifying the most stable materials in a fraction of the time required to develop a reliable synthetic protocol. In this Account, we carry out a critical assessment of the three methods for four different MOFs and five different phases, covering the whole MOF spectrum stretching from the mechanically rigid MOFs UiO-66(Zr)[32] and MOF-5(Zn)[33] to the prototype flexible MOFs MIL-47(V)[34] and MIL-53(Al).[35] On the basis of this assessment, guidelines are provided to identify which methods can be applied to reliably predict the rigidity or flexibility crucial for current technological applications.

Theoretical Methods To Assess Mechanical Stability

The three computational techniques discussed in this Account are schematically illustrated in Figure b–d and discussed below. We refer the reader to the Supporting Information for additional technical background and the computational details.

Born Stability Criteria

First stated by Born in his 1940 paper,[26] the Born stability criteria form a set of necessary and sufficient criteria that determine whether a given unstressed material is stable. These stability criteria are obtained by requiring that the energy of the lattice increases for any infinitesimal strain ε imposed on the unit cell in its equilibrium volume V0. To second order, this energy is given byThese criteria therefore require the second-order elastic stiffness tensor C to be positive-definite in equilibrium. In Voigt notation, this tensor is defined asThese thermodynamic criteria were later generalized also to account for systems subject to an arbitrary external Cauchy stress σ = P1 + σa, where P is the isotropic pressure and σa is the deviatoric stress.[36] For a purely isotropic loading σ = P1, the unstressed second-order elastic stiffness tensor C is replaced by its stressed analogue B, given by[37]To confirm the stability of a material at a given temperature T and subject to a given pressure P, the positive-definiteness of the second-order stiffness tensor B needs to be confirmed, for instance by verifying that all of its eigenvalues are positive or that all of its leading principal minors are positive (Sylvester’s criterion). This procedure has to be repeated for a set of temperatures and/or pressures to determine the experimental conditions for which one of the Born stability criteria is first violated.

For a general crystal system, the second-order elastic stiffness tensor C consists of 21 independent elastic constants. For more symmetric systems, C becomes sparser, leading to a more comprehensive set of Born stability criteria (see Section S1 in the Supporting Information and the comprehensive work of Mouhat and Coudert[37]). For instance, for the orthorhombic crystal system, encountered in the large-pore phases of MIL-47(V) and MIL-53(Al), only nine independent elastic constants exist, such that the stability under tensile and shear deformations can be separated. The six Born stability criteria under a hydrostatic pressure P then readin whichare the three leading principal minors of B. As shown in Section S1, these six criteria correspond to stability under three independent tensile deformations (BSC > 0) and under three shear deformations (C – P > 0).

For the highly symmetric cubic system, encountered in UiO-66(Zr) and MOF-5(Zn) (see Figures a and 3a), only three independent elastic constants exist. The positive-definiteness of B under a hydrostatic pressure P then leads to the following three Born stability criteria:These three criteria encompass (i) an isotropic tensile deformation, (ii) a twofold-degenerate axial tensile deformation, and (iii) a threefold-degenerate shear deformation, which are schematically depicted in the insets of Figures d and 3d.

Figure 2

UiO-66(Zr): (a) unit cell with directions of highest and lowest Young modulus E; (b) directional Young modulus E at 300 K and 0 MPa; (c) P(V) profile at 300 K; (d) Born stability criteria as a function of the applied pressure at 300 K; (e) anisotropies of the mechanical moduli as a function of the applied pressure at 300 K, compared with DFT[38] and FF optimizations. Unstable regions are shaded in red.

Figure 3

MOF-5(Zn): (a) unit cell with directions of highest and lowest Young modulus E; (b) directional Young modulus E at 300 K and 0 MPa; (c) P(V) profile at 300 K; (d) Born stability criteria as a function of the applied pressure at 300 K; (e) anisotropies of the mechanical moduli as a function of the applied pressure at 300 K, compared with DFT[42] and FF optimizations. Unstable regions are shaded in red.

Anisotropy of Mechanical Properties

A second indicator of the mechanical stability of crystals was proposed in 2012 by the Coudert group.[29] In this method, the directional Young modulus E and shear modulus G are calculated from the second-order stiffness tensor C. These mechanical moduli are defined aswhere and are unit vectors.[29] Subsequently, the anisotropy of each of these properties is calculated as the ratio of its maximum value to its minimum value, which varies from 1 (for a perfectly isotropic material) to infinity. This anisotropy correlates with the flexibility of the material: for flexible MOFs, anisotropy factors are found to be larger than 100, while significantly lower values, near unity, are associated with rigid MOFs.[29] These anisotropy values are typically only calculated at one temperature and pressure, yielding no further quantitative information on the thermodynamic conditions under which the material becomes unstable.

Pressure-versus-Volume Equations of State

A last thermodynamic method to determine the mechanical stability of crystals was recently proposed by the authors of this Account.[25] Key in this method is the construction of pressure-versus-volume equations of state. To obtain this equation of state for a given material, MD simulations at a predefined temperature T are performed for a series of volume points within the range of interest (see Figure d). During these simulations, the volume V of the crystal is kept constant, whereas the cell shape can fluctuate—but in such a way that the ensemble-averaged deviatoric stress ⟨σa⟩ vanishes, leading to the dedicated (N, V, σa = 0, T) ensemble.[25] This ensemble was specifically designed so that the large fluctuations in the internal pressure, which are present in the (N, P, σa = 0, T) ensemble and generally lead to poor convergence and even premature phase transitions for flexible materials, no longer influence the volume of the material, which is kept fixed in the (N, V, σa = 0, T) ensemble. With the assumption of mechanical equilibrium, the well-defined ensemble-averaged internal pressure ⟨P⟩ exerted by the material on its environment is exactly the pressure the environment can exert on the material at this volume. Hence, the obtained ⟨P(V)⟩ profile corresponds with the macroscopic pressure-versus-volume equation of state.

From the pressure-versus-volume equation of state, the structure of the material at any given pressure can be predicted by determining the intersection(s) between the equation of state and a horizontal line at the required pressure (see Section S3). All intersections with a negative ∂P/∂V slope are thermodynamically stable and can hence be obtained experimentally by applying the given pressure. Transition pressures, i.e., those limiting pressures above or below which a metastable phase disappears, can directly be identified as the minima and maxima of the theoretical equation of state.

Loss of Crystallinity in Rigid MOFs: UiO-66(Zr) and MOF-5(Zn)

As a first test for the three theoretical methods, their performance on UiO-66(Zr) and MOF-5(Zn), schematically shown in Figures a and 3a, is discussed. As these materials are known for their mechanical rigidity, the models need to perform well even under pressures in the gigapascal range to accurately describe this rigidity. For UiO-66(Zr), its exceptional stability can be traced back to the high linker coordination number (12) and the strong Zr–O bonds. The material retains its crystalline structure under high pressures (>1 GPa),[20,38] high temperatures,[39] and harsh acidic conditions,[40] even when structural linker defects are present.[41] MOF-5(Zn) is less stable, for instance decomposing in humid air,[4] because of the weaker Zn–O bonds—zinc is less oxophilic than zirconium—and the lower coordination number (6).[33] Nevertheless, experimental and computational studies indicate that the material is still endowed with an appreciable mechanical stability[42,43] but yield no identification of the weakest modes in MOF-5(Zn) at elevated temperatures.

As a first computational validation of this stability, the directional Young moduli at 300 K and 0 MPa for UiO-66(Zr) and MOF-5(Zn) are depicted in Figures b and 3b, respectively. The maximum Young modulus is encountered along the body diagonal for UiO-66(Zr) and along the crystal axes for MOF-5(Zn), in both cases coinciding with the positions of the organic ligands and hence the directions of the highest atom density. The minimum Young modulus, in contrast, points toward the largest pores in both materials. As even this minimum Young modulus still amounts to about 10–30 GPa (see Figures S1 and S2), both MOFs are qualitatively expected to be rigid. While these Young moduli are among the highest encountered in MOFs, they are rather small compared with values encountered in conventional materials (see Figure 2 of ref (8).).

To quantify the pressures at which these materials become mechanically unstable, pressure-versus-volume equations of state were generated (see Figures c and 3c). UiO-66(Zr) exhibits an unstable region (∂P/∂V > 0) below about 8450 Å3, which can be reached by applying pressures higher than 1.83 GPa, while MOF-5(Zn) exhibits an unstable region below about 17 300 Å3, which can be reached by applying an appreciably lower pressure of 189 MPa. For UiO-66(Zr), this instability was revealed to coincide with a broadening of the peaks in the radial distribution function at a pressure that agrees well with the experimental loss of crystallinity when accounting for the defects present in experimental samples.[20,28] This same indicator of loss of short-range order is also observed in MOF-5(Zn) (see Figure S7), indicating that also for MOF-5(Zn) the instability is correlated with a loss of crystallinity.

This same loss-of-crystallinity pressure is also predicted by determining the Born stability criteria at several pressures. For the cubic UiO-66(Zr) and MOF-5(Zn) (space group Fm3̅m), only three independent Born stability criteria exist (see Figures d and 3d). At an extrapolated value of about 1.84 GPa, C11 – C12 – 2P is the first Born criterion to be violated in UiO-66(Zr), hence indicating that the material fails in an axial tensile deformation at a pressure very close to the one obtained from the P(V) profile. For MOF-5(Zn), the Born stability criterion corresponding to a shear deformation is the smallest at 0 MPa and decreases monotonically, which might give an indication that this is the weakest mode.[44] However, MD simulations at pressures close to instability show that the first Born criterion to fail corresponds to an isotropic compression at an extrapolated value of 190 MPa, whereas the criterion corresponding to a shear deformation remains positive, even for pressures larger than 190 MPa.

Finally, the anisotropies of the Young and shear moduli were calculated (Figures e and 3e). For UiO-66(Zr), FF geometry optimizations yield anisotropy factors in the range 1.2–1.3, similar to those values found by earlier reported DFT studies,[38] revealing its rigidity. To account for possible temperature effects, which may profoundly affect the structure of MOFs,[30,44,45] these anisotropy factors have also been determined at 300 K via FF MD simulations. For UiO-66(Zr), this temperature increase leads to a slight increase in anisotropy, but UiO-66(Zr) is still classified as a rigid MOF, even when the pressure is increased close to the point of instability. For MOF-5(Zn), a similar picture emerges: while the FF-optimized results are thrice as large as the DFT results of Bahr et al.,[42] FF MD simulations at 300 K predict anisotropy values of about 1–2, which increase slowly as the pressure converges to the loss-of-crystallinity pressure. These values are even smaller than those obtained for UiO-66(Zr) and clearly classify the material as rigid.

Phase Transformations in Flexible MOFs: MIL-47(V) and MIL-53(Al)

The grand challenge in the description of flexible MOFs exists in accurately determining those pressures that induce phase transitions, on the order of 10–100 MPa. MIL-47(VIV)[34] and the isoreticular MIL-53(Al),[35] which is obtained by substituting the vanadium–oxide chain with an aluminum–hydroxide chain, were selected as case studies. They are schematically shown in Figures a and 5a. While MIL-53 is the archetypical flexible MOF,[46−48] exhibiting transitions between a large-pore (LP) phase and a closed-pore (CP) phase at low pressures of 13–18 MPa,[12] transitions to the CP phase of MIL-47(V) are induced only by pressures of about 100 MPa.[19] While MIL-53(Al) remains in its CP phase even after the applied pressure reverts to atmospheric pressure, giving rise to potential applications as shock absorbers,[12] MIL-47(V) returns to its LP phase when the pressure is decreased and is hence endowed with the potential to act as an efficient shock dampener.[30] An accurate determination of these transition pressures is vital for the further development of flexible MOFs.

Figure 4

MIL-47(V): (a) unit cell with directions of highest and lowest Young modulus E; (b) directional Young modulus E at 300 K and 0 MPa; (c) P(V) profile at 300 K; (d) Born stability criteria as a function of the applied pressure at 300 K; (e) anisotropies of the mechanical moduli as a function of the applied pressure at 300 K, compared with DFT[29] and FF optimizations. Open symbols in (d) and (e) denote simulations with a shortened simulation time. Unstable regions are shaded in red.

Figure 5

LP phase of MIL-53(Al): (a) unit cell with directions of highest and lowest Young modulus E; (b) directional Young modulus E at 300 K and −100 MPa; (c) P(V) profile at 300 K with indication of the two metastable phases at 0 MPa; (d) Born stability criteria as a function of the applied pressure at 300 K; (e) anisotropies of the mechanical moduli as a function of the applied pressure at 300 K, compared with DFT[29] and FF optimizations. Open symbols in (d) and (e) denote simulations with a shortened simulation time. Unstable regions are shaded in red.

As these MOFs are isoreticular, the directional Young moduli of the LP phases of MIL-47(V) (Figure b) and MIL-53(Al) (Figure b) are very similar. Because of the weaker bonds in the inorganic chain of MIL-53(Al), however, its Young modulus in the y direction is slightly smaller, as evidenced by the different maximum E values in these figures. The MIL-53(Al) CP phase is denser, resulting in a stronger material with a larger Young modulus (Figure b). As for the rigid materials, the directions of highest Young modulus coincide with the positions of the linkers, whereas the lowest Young modulus points into the pores.

Figure 6

CP phase of MIL-53(Al): (a) unit cell with directions of highest and lowest Young modulus E; (b) directional Young modulus E at 300 K and 0 MPa; (c) P(V) profile at 300 K with indication of the two metastable phases at 0 MPa; (d) Born stability criteria as a function of the applied pressure at 300 K; (e) anisotropies of the mechanical moduli as a function of the applied pressure at 300 K, compared with DFT[49] and FF optimizations. Unstable regions are shaded in red.

The flexibility of MIL-47(V) is evidenced by the two metastable branches in the pressure-versus-volume equation of state (Figure c), in contrast to the single stable branch for the aforementioned rigid materials. The LP (V ≳ 1400 Å3) and CP (V ≲ 1100 Å3) branches are separated by an unstable region for intermediate volumes. Starting from the LP structure at 0 MPa, increasing the pressure results in mechanical instability when the pressure exceeds 123 MPa, in close agreement with the experimental LP-to-CP transition pressure. Moreover, this pressure profile also reveals the CP-to-LP transition pressure, amounting to 56 MPa at 300 K, correctly indicating hysteresis.[19] Similarly, the MIL-53(Al) pressure-versus-volume equation of state, shown in Figures c and 6c, exhibits a local maximum (PLP→CP) at 29.6 MPa and a local minimum (PCP→LP) at −182 MPa. While the LP-to-CP transition pressure is in rather good agreement with experiment,[12] the CP-to-LP transition pressure has not yet been measured experimentally, as it would require a negative pressure or pulling, for instance by embedding the MOF in a membrane. Finally, it should be noted that the profiles for the LP and CP phases (Figures c and 6c) are identical, indicating that information about both phases can be extracted from the same set of MD simulations.

Since the LP phases of MIL-47(V) and MIL-53(Al) belong to the orthorhombic crystal system (space group Pmma), the six Born stability criteria are no longer degenerate but can be divided into three tensile and three shear deformations (denoted by orange and blue diamonds, respectively, in Figures d and 5d). As observed from Figure d, the first Born stability criterion tending to fail in MIL-47(V) corresponds to a tensile deformation at an extrapolated transition pressure of 118 MPa. In contrast, the first Born stability criterion to fail for the LP phase of MIL-53(Al) is C66 – P. This corresponds to a shear deformation on the yz plane along the y direction, parallel to the relatively weak inorganic chain. Apart from this difference, the ordering and magnitudes of the other five Born stability criteria are shared between the LP phases of MIL-53(Al) and MIL-47(V). Compared with the cubic UiO-66(Zr) and MOF-5(Zn), the three tensile deformations can no longer be easily visualized, as they correspond to the determinants of the three leading minors of the stiffness tensor. Moreover, the stability analysis needs to be repeated for every phase, and large pressure fluctuations prohibit MD simulations close to instability.[25] Hence, the Born stability criteria need to be extrapolated to obtain the pressure at which the material becomes unstable, leading to larger inaccuracies.

To study the mechanical stability of the low-symmetric monoclinic MIL-53(Al) CP phase (space group P2/c) using the Born stability criteria, the positive-definiteness of the stiffness tensor was validated by determining its eigenvalues. As shown in Figure d, the smallest eigenvalue vanishes at a negative pressure of about −185 MPa, when the CP phase becomes unstable, in favor of the LP phase. It should be noted that no fixed eigenmode can be assigned along the whole range of pressures since the eigenmodes of the stiffness tensor are pressure-dependent, impeding visualization of the weakest eigenmode at the pressure of instability. As a result, the Born stability criteria are less suited for low-symmetric materials.

Finally, the anisotropies of the mechanical properties of the LP phases of MIL-47(V) and MIL-53(Al) at 300 K (Figures e and 5e) are one to two orders of magnitude larger than those of the rigid UiO-66(Zr) and MOF-5(Zn), in agreement with earlier DFT optimizations.[29] For the MIL-53(Al) CP phase, the anisotropy in the shear modulus is one order of magnitude smaller, while the anisotropy in the Young modulus is on the same order of magnitude in both MIL-53(Al) phases. All of these calculated anisotropy values indicate that these MOFs are flexible, and they increase substantially for pressures close to the transition pressure.

Closing Theoretical Guidelines and Perspectives

Based on the presented case studies, we can now formulate a comprehensive set of theoretical guidelines to computationally assess the flexibility or rigidity of MOFs. When one is solely interested in classifying the materials as either rigid or flexible, the anisotropy of the Young and shear moduli was found to be a reliable qualitative criterion. Whereas the exact values of these moduli and their anisotropies depend on the accuracy of the adopted level of theory, this technique correctly predicts the general flexibility behavior of the material under pressure. However, this criterion proved inadequate to reveal more quantitative information regarding the critical pressure at which the mechanical instability takes place and does not provide any microscopic or macroscopic insight into the origin of the instability.

To extract quantitative information on this instability, either the Born stability criteria or the pressure-versus-volume equations of state can be adopted, depending on the additional information one wishes to extract. If one wishes to shed light on the macroscopic weakest mode of deformation, the Born stability criteria can be employed. However, only for highly symmetric materials are these deformations easily interpretable. Moreover, large pressure fluctuations during these simulations may induce premature phase transitions,[25] impeding an accurate determination of these criteria near the mechanical instability, as one needs to rely on extrapolation. Especially for flexible materials, these two shortcomings prohibitively hamper the application of the Born stability criteria.

In contrast, these shortcomings are mitigated when pressure-versus-volume equations of state are used, resulting in more accurate results in a shorter time frame. Moreover, these profiles provide microscopic insight into the mode of deformation, for instance identifying the loss of crystallinity in rigid materials, and can be integrated to reveal the relative stability of the metastable phases. In contrast to the two aforementioned techniques, only a single equation of state needs to be constructed for a given material, even if multiple metastable phases are present, as demonstrated here for MIL-53(Al).

This critical assessment leads to the conclusion that only the Born stability criteria and the pressure-versus-volume equations of state provide reliable, quantitative information necessary for an in silico high-throughput screening to identify MOFs that are sufficiently stable to be shaped postsynthetically. Whereas both methods provide additional macroscopic or microscopic insight for rigid materials, only the latter procedure offers a facile way to characterize flexible materials. In the framework of a widely applicable procedure for the mechanical characterization of MOFs, this technique hence holds the highest potential to promote the further development of MOFs for technological applications.

While pressure-versus-volume equations of state can predict with a high reliability the mechanical rigidity or flexibility of a material and hence promote application-driven MOF synthesis, the accuracy of these computational methods is determined by the adopted level of theory. The level of theory should be able to reliably reproduce the potential energy surface, including the correct identification of all relevant metastable phases and their relative stability. First-principles-based methods yield a more accurate description of the potential energy surface, but they are often computationally too demanding for large screening studies. Moreover, these simulations need to be carried out at operating temperatures, as the stability in MOFs is often determined by a subtle interplay between entropic and dispersion interactions. Hence, at least for the near future, the further development of this computational protocol will be closely intertwined with the development of accurate force fields, with the common goal to bring the breakthrough of MOFs in industrial applications closer to reality.