Theoretical prediction of the mechanical properties of zeolitic imidazolate frameworks (ZIFs) †

A good resistance against mechanical stress is essential for the utilization of metal – organic frameworks (MOFs) in practical applications such as gas sorption, separation, catalysis or energy conversion. Here, we report on the successful modi ﬁ cation of the mechanical properties of zeolitic imidazolate frameworks (ZIFs) achieved through a substitution of the terminal group. The mechanical modulus of SALEM-2 was found to signi ﬁ cantly improve when the – H groups at position 2 of the imidazole linkers were replaced with electron withdrawing groups ( – CHO, – Cl, or – Br). The charge distribution and electron density were analyzed to reveal the mechanism behind the observed variation of the elastic sti ﬀ ness. Furthermore, ZIF-I with a – I group at position 2 of the imidazole linkers was predicted to exhibit an excellent mechanical strength in our study and then prepared experimentally. The results indicate that an inconspicuous change of the structure of ZIFs, i.e. , additional groups strengthening the ZnN 4 tetrahedron, will lead to a sti ﬀ er framework.


Introduction
3][4][5][6][7] However, the mechanical instability of some MOFs limits their commercial and industrial application potential.For instance, a loss of porosity may occur during the shaping or post-synthetic processing of MOF powders through mechanical pelletization, extrusion or sintering (i.e., the formation of a solid mass through heating). 8,9ue to their exceptionally high chemical and thermal stability, zeolitic imidazolate frameworks (ZIFs) are among the most frequently investigated MOFs. 10 The earliest study on the mechanical stability of ZIFs was conducted by Tan et al., who performed single-crystal nanoindentation tests. 11They found the mechanical properties of ZIFs to be superior to the corresponding properties of other MOFs and to be mainly determined through the rigidity and bulkiness of the substituted imidazolate linkers.Further experimental and theoretical investigations revealed that ZIF-8 exhibits an exceptionally low shear modulus below 1 GPa, [12][13][14][15] which is much lower than values previously reported for single-crystalline extended solids.
In particular, the shear elastic constant (C 44 ) of the cubic ZIF-8 was found to decrease with the applied external pressure, an effect which is called shear mode soening. 16The weak resistance against shear deformation was later veried by the observation that ZIF-8 might undergo a rapid and catastrophic porosity collapse during ball-milling. 17,18he addition of modulator ligands is a common approach to improve the mechanical robustness of MOFs. 11,17,19However, this may induce a structural instability through guest-host interaction, 20 and/or reduce the framework's porosity.Overcoming the limitations imposed by this trade-off relationship between the elastic stiffness and the framework's porosity therefore remains the main challenge when aiming to improve the mechanical properties of MOFs.
2][23] However, the low number of reports suggests that it is much more difficult to experimentally determine the corresponding shear modulus. 12herefore, the theoretical computation of the mechanical modulus of ZIFs has become an important alternative, 12,15,16 although it remains challenging to obtain accurate values through computational methods. 24Yet a comparison of the mechanical modulus values obtained for different ZIFs might yield important insights.
In this work, we employed density functional theory (DFT) calculations to investigate the mechanical properties of four ZIFs of the same geometry in order to isolate the effect of linker functionalization on the strength of MOFs.In particular, we demonstrate that small changes in the linker molecules can have a drastic impact on the mechanical properties of the studied ZIFs.These results are expected to be helpful for designing future ZIF materials with excellent mechanical properties.

Theoretical methods
The topological unit of the ZIF structures studied in this work consists of two Zn tetrahedrons linked via an imidazole ring.Substituting the terminal group of imidazole has been demonstrated to be a good strategy to tune the properties of ZIFs without having to change its symmetry.Our calculations were performed for four types of ZIFs (Fig. 1), which have been experimentally investigated before and share the same sodalite (SOD) topology but contain different functional groups: SALEM-2, 25 in which the linker is an unfunctionalized imidazolate; ZIF-Cl, 26 with a chlorofunctionalized linker (cim ¼ 2-chloroimidazolate); ZIF-Br, 26 with a bromo-functionalized linker (bim ¼ 2-bromoimidazole); and ZIF-90, 27 where the linker (ica ¼ imidazolate-2-carboxaldehyde) contains a aldehyde group.Other terminal functional groups, such as stronger electron donating groups (-NH 2 ) or stronger electron withdrawing groups (-NO 2 ), were not included in this study because we decided to focus on real ZIF structures and exclude purely hypothetical frameworks.Although 2-nitroimidazole can be used to synthesize an SOD framework, the metal center is Co instead of Zn. 28 In this study, the calculations were performed utilizing the CASTEP code 29 based on density functional theory (DFT). 30,31ifferent exchange-correlation functions with or without dispersion interaction corrections were tested (Tables S1 and  S2 †).We found that utilizing a generalized gradient approximation functional with Perdew-Burke-Ernzerhof parametrization (GGA-PBE) 32 using Vanderbilt's ultraso pseudopotential 33 and including dispersion corrections improved the reproduction of the experimental lattice parameters for most of the ZIFs studied in this work.The geometry optimization was performed employing the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, 34 with the cutoff energy selected to be 500 eV. 3 Â 3 Â 3 Monkhorst-Pack k-point meshes 35 were chosen for the Brillouin zone sampling to achieve a total energy convergence lower than 1 Â 10 À5 eV per atom, a maximum Hellmann-Feynman force of 0.03 eV ÅÀ1 , a maximum ionic displacement of 0.001 Å and a maximum stress of 0.05 GPa.The neness of the cutoff energy and the k-mesh used in our DFT calculations were tested and the results are presented in Fig. S1 and S2.† The nite basis set correction method was employed to optimize the unit cell parameters.
The elastic constants c ij were calculated through c ij ¼ s i /3 j , where s and 3 are the elastic stress and strain, respectively, and the subscripts i and j denote the Cartesian coordinates of the considered structures. 36A 1% strain was chosen to guarantee a purely elastic response.For the optimization of the internal atomic freedoms, the criteria for convergence were selected in such a way that the difference in total energy was lower than 2 Â 10 À6 eV per atom, the maximum force was 0.006 eV ÅÀ1 and the maximum displacement was 2 Â 10 À4 Å.Although the computed mechanical modulus of the ZIFs depends on the selected exchange-correlation functions and the dispersion correction, the variation trend of the mechanical properties obtained for the different ZIFs proved to be consistent (Table S2 †).

Results and discussion
In this work, the mechanical properties were obtained through a tensorial analysis of the single-crystal elastic constants, C ij .For the cubic phase, three non-zero independent elastic constants, i.e., C 11 , C 12 , and C 44 (Table 1), were used to compute the mechanical properties.The coefficient C 11 (¼C 22 ¼ C 33 ) indicates the stiffness along the crystal axis a (b or c), the coefficient C 12 (¼C 21 ) corresponds to the tensile-tensile coupling between the crystal axes a and b, and the shear coef-cient C 44 (¼C 55 ¼ C 66 ) corresponds to the stiffness against angular distortion when the cell is subjected to shear strain.The elastic constants determined for the studied ZIFs, which are listed in Table 1, satisfy the fundamental elastic stability criteria, which imposes the following restrictions on a cubic crystal: C 11 > rC 12 r, C 11 + 2 C 12 > 0, and C 44 > 0. 37 The elastic stability criterion was based on the convexity requirement of the equilibrium free energy of a stress-free crystal under small strain uctuation, according to the theory of lattice dynamics. 38n order to satisfy above demand, the single-crystal elastic constant matrix must be positive denite and then the restrictions on any types of crystal can be obtained.
According to the results presented in Table 1, the elastic constants C 11 , C 12 and C 44 do not depend on the size of the terminal group.ZIF-90 exhibits a large group size (-CHO) but low C 11 and C 44 values.This indicates that it is the electronegativity of  the terminal group which mainly contributes to the deformation resistance of the ZIFs and not the geometry.Here, we therefore only focus on the contribution of the electronegativity of the terminal groups and neglect the effect of the geometrical shape.
The shear modulus G may serve as a measure of the rigidity of the framework against structural distortion due to external shear forces.In this work, the G values were computed using eqn (1). 39Young's modulus E is dened as the reciprocal of the elastic compliance constants (eqn (2)).The bulk modulus B represents the resistance of the structure against volumetric strain under hydrostatic pressure and can be computed through eqn (3). (1) where S ij is the compliance matrix, which is the inverse of the elasticity matrix (tensor).The elastic compliance constants (S ij , in GPa À1 ) of the isostructural ZIF-8 are listed in Table 2.The Voigt-Reuss-Hill (VRH) averaging method 39 was used to compute the shear moduli (G H ) in this work (Table 1).The Voigt value (G V ) assumes a uniform strain, the Reuss value (G R ) corresponds to a uniform stress, and the Voigt-Reuss-Hill value (G H ) is the average of the two.According to eqn (1), the shear modulus strongly depends on the shear coefficient C 44 .The larger the C 44 value, the higher the shear modulus (G).Table 1 shows that ZIF-Br exhibits the largest shear modulus, followed by ZIF-Cl, ZIF-90 and SALEM-2.
Similarly, the highest Young's modulus (E) was also obtained for ZIF-Br due to its higher C 11 value.Table 1 also reveals the same order, i.e., ZIF-Br, ZIF-Cl, ZIF-90 and SALEM-2 (from high to low), for the bulk modulus (B).Excellent overall mechanical properties (G, E and B) were obtained for ZIF-Br, followed by ZIF-Cl, ZIF-90 and SALEM-2 (Table 1).
In order to better understand the variation of the mechanical properties of ZIFs, we have analyzed their charge distribution.
The calculated Mulliken atomic charges (q) of the four ZIFs investigated in this study are summarized in Table 2.The results show that the charge distribution throughout all of the ZIFs hardly changes, except for the Zn and C1 atoms and the terminal groups.The ZIFs can be divided into two groups according to the electronegativity of the terminal groups.The rst group only consists of SALEM-2, which features the electron donating group -H (q ¼ 0.240).A conical shape of the electron density distribution was obtained for the C-H bond (the le circular shade in Fig. 2) in SALEM-2, indicating an electron transfer from the -H group to the 5-member ring.The other three ZIFs studied here form the second group as they contain an electron withdrawing group (EWG).The electron density distribution in the terminal group features a dumbbell shape.The -Cl group show the highest electron-withdrawing ability, followed by -Br and -CHO.The analysis of the results presented in Tables 1 and 2 revealed two interesting phenomena: rst, the mechanical strength can benet from the presence of an EWG located at position 2 of the imidazole ring.The weakest mechanical strength was obtained for SALEM-2 due to its electron donating group (EDG).Second, the mechanical strength of a ZIF becomes weaker if a strong EWG is present, such as the -CHO group in ZIF-90.
The charge transfer model can be used to explain the above phenomena.Table 2 shows that the Zn atoms are positively charged with q ranging from 1.52 for SALEM-2 to 1.35 for ZIF-Br, which indicates a reduction of the electron donation by the Zn atoms.While the charge of the N atoms remains almost constant, the electrons donated by the Zn atoms will eventually contribute to the charge of the C1 atoms (Table 2).The electron loss at the Zn-N bond lowers the local electron density (Fig. 2), thereby reducing the strength of the ZnN 4 tetrahedral, which ultimately determines the mechanical properties of the ZIFs. 40ble 2 Comparison of the atomic charges (q, e) of the studied ZIFs.The definition of the notation is provided in Fig. 1 Type Atomic charge In this regard, terminal groups with a weak electronwithdrawing ability can contribute to the electron accumulation in ZnN 4 , thus strengthening the ZIFs.The above relation between the electron-withdrawing ability of the terminal group and the framework strength could be further validated through investigating the mechanical properties of ZIF-I (Fig. 3), where the linker (iim ¼ 2-iodoimidazole) features a -I group with a comparatively weak electron-withdrawing ability (Table 2).As expected, ZIF-I exhibits a larger mechanical modulus (Fig. 3) than the other ZIFs (Table 1).We then successfully synthesized ZIF-I (ESI †) to prove that it is not just a hypothetical structure.
In previous studies, 40 we have argued that the axial stretching and shear deformation of ZIFs depend on the Zn-N bond length and the N-Zn-N angle, respectively.Here, we used the cluster model to quantify the enhancing effect of the terminal group on the bond length and angle in ZnN 4 (Table 3).The model and computational details can be found in the ESI.† The results presented in Table 3 clearly show that the force constants (K b and K a ) obtained for the Zn-N bond and the N-Zn-N angle stretching increase as the electron-withdrawing ability of the terminal groups decreases.The results calculated using the cluster model are consistent with the information given in the periodic table, again conrming that the strength of the ZnN 4 tetrahedral and thereby the mechanical properties of ZIFs can be adjusted by replacing the EWGs.

Conclusions
In summary, we have computed and assessed the complete elastic properties of the experimentally accessible ZIFs with SOD geometry by means of ab initio density functional theory (DFT).Our calculation results indicate that ZIF-Br with its weak electron withdrawing group (EWG) may exhibit excellent overall mechanical properties (G, E and B).In comparison, strong EWGs, such as -Cl in ZIF-Cl and -CHO in ZIF-90, can attract more electrons from the Zn atoms.These electrons then aggregate at the C atoms directly bonded to the EWGs.The resulting lower electron density in the ZnN 4 tetrahedrons weakens the bond and angle strength and thereby the mechanical strength of the ZIFs.Considering these results, a new synthesized ZIF structure (ZIF-I) with a terminal group that exhibits a low electron-withdrawing ability was modeled and predicted to possess even better mechanical properties than ZIF-Br.In this work the cluster model was used to compute the force parameters of the Zn-N bond and the N-Zn-N angle stretching, conrming the enhancing effect of the terminal groups.The substitution of the terminal groups therefore seems a promising strategy to obtain strong MOFs which retain their high porosity and are capable of resisting the pressure or stresses involved in practical applications.

Conflicts of interest
There are no conicts to declare.

Fig. 1
Fig. 1 Representation of (a) the unit cell of ZIFs with an SOD topology and (b) the basic unit of ZIFs with the different functional groups studied in this work.

Fig. 3
Fig. 3 Basic unit of ZIF-I (left) and its predicted mechanical properties (right).

Table 1
Comparison of the single-crystal elastic constants C ij (GPa), the elastic compliance constants S ij (GPa À1 ), and the values calculated for Young's modulus E (GPa), the shear modulus G (GPa) and the bulk modulus B (GPa) of the studied ZIFs Electron density distribution computed for the selected ZIFs.