The Atomic Structure and Mechanical Properties of ZIF-4 under High Pressure: Ab Initio Calculations

The effects of pressure on the structural and electronic properties and the ionic configuration of ZIF-4 were investigated through the first-principles method based on the density functional theory. The elastic properties, including the isotropic bulk modulus K, shear modulus G, Young’s modulus E, and Poisson’s ratio ν of the orthorhombic-type structure ZIF-4 were determined using the Voigt–Reuss–Hill averaging scheme. The results show that the ZIF-4 phase is ductile according to the analysis of K/G and Cauchy pressure. The Debye temperatures obtained from the elastic stiffness constants increase with increasing pressure. Finally, the pressure-dependent behaviors of the density of states and ionic configuration are successfully calculated and discussed.


Introduction
Zeolitic imidazolate frameworks (ZIFs) are a unique class of porous hybrid materials known as a subset of metal-organic frameworks (MOFs), and have been studied extensively over the past decade due to their potential applications in gas sorption and separation, catalysis, drug delivery, and sensing applications [1][2][3][4][5]. Especially in the field of environmental protection, ZIFs and their composite structures have also shown considerable potential in research into the capture of greenhouse gases recently [6][7][8][9]. The ZIF structure consists of tetrahedral metal ions, typically Zn 2+ or Co 2+ , linked by organic bridging ligands, which are derived from imidazolate anions (Im = C 3 N 2 H 3 − ) [10]. They have significantly softer mechanical properties compared with their inorganic cousins and may structurally collapse upon heating [11][12][13][14], pressurization [15][16][17], ball-milling [18][19][20], or under an electric field [21] to form amorphous frameworks [22], which possess the same short-range connectivity as their crystalline counterparts [23,24]. Among the various ZIFs, ZIF-4 (Zn(C 3 N 2 H 3 − ) 2 ) has attracted lots of interest because of its unique behaviors, such as remarkable structural flexibility [25], and has great potential in the field of small-molecule adsorption and separation [26,27]. In particular, the complex phase transition behavior of ZIF-4 under high pressure [28] and its vitrification have attracted a great deal of attention from researchers [29]. Bennett et al. showed that amorphization occurs at various hydrostatic pressure values depending on the presence or absence of solvent molecules in the pores of the material, as well as the nature of the hydrostatic medium. The evacuated ZIF-4 demonstrated amorphization at very low pressure (0.35 GPa). At the same time, the presence of DMF (dimethylformamide, C 3 H 7 NO) molecules in its pores could shift amorphization to higher pressure and lead to an intermediate monoclinic crystalline phase (ZIF-4-I) [15]. Recent in situ XRD (X-ray diffraction) studies have shown that ZIF-4 undergoes multiple crystal-crystal transitions before undergoing an irreversible amorphous transition [28]. The elastic behavior of ZIF-4 and its high-pressure phase ZIF-zni an irreversible amorphous transition [28]. The elastic behavior of ZIF-4 and its highpressure phase ZIF-zni is well explained by the simulation research of Tan et al. [30]. Furthermore, the op-cp phase transitions (i.e., breathing transitions) of ZIF-4 driven by mechanical pressure were identified and taken into account. The change in pore size suggested potential applications of the gas storage of this functional material [31]. Lately, Vervoorts et al. have applied a cutting-edge HPXRD (high-pressure X-Ray diffraction) setup to probe the response of ZIF-4 at low hydrostatic pressures. An initial correlation between the mechanical properties of ZIF-4 and mid-range structural parameters, such as free volume, is presented [32].
The structure-properties relationship of ZIF-4 is plagued by several obstacles, such as small pressure needs and small structure characterization. Computation simulation provides a powerful tool to probe the microscopic details of structure evaluation under pressure [33,34]. For complex structures such as ZIFs, ab initio methods without empirical parameters can reveal deeper correlations between electronic structural and mechanical properties [35][36][37][38]. To our knowledge, there is no systematic study of the mechanical properties-structure relationship of ZIF-4 with varying pressure up to now. Figure 1 shows the schematic structure of ZIF-4 at p = 0, 0.5, and 1 GPa. Additionally, the insert diagrams show the average Zn-N bond lengths at varied pressures. The variation in the local environment of Zn is very limited, as can be seen by the change in Zn-N bond length, and Zn(Im)4 still has a tetrahedral structure at different pressures. Changes in mechanical properties under pressure can often be explained by structural features in the mid-range. To further investigate the response to different pressures, the present work is conducted to clarify the phase stability, elastic constants, and electronic properties of ZIF-4 under high pressure of up to 1 GPa.

Results
The model of the ZIF-4 structure is firstly fully optimized under 0 GPa. The calculated lattice parameters are very close to the available experimental data [15] and other DFT calculation works [39]. From Table 1, the structure at zero pressure has a difference of about 3% on the b-axis. The main reason is the overestimation of the lattice parameters by the GGA-PBE method [40,41]. Such systematic errors would not have a practical impact on the reliability of this work, considering that all calculations were performed under the same parameters and criteria. The pressure dependences of lattice parameter ratios and volume ratios are illustrated in Figure 2a. The corresponding results are listed in Table 1.

Results
The model of the ZIF-4 structure is firstly fully optimized under 0 GPa. The calculated lattice parameters are very close to the available experimental data [15] and other DFT calculation works [39]. From Table 1, the structure at zero pressure has a difference of about 3% on the b-axis. The main reason is the overestimation of the lattice parameters by the GGA-PBE method [40,41]. Such systematic errors would not have a practical impact on the reliability of this work, considering that all calculations were performed under the same parameters and criteria. The pressure dependences of lattice parameter ratios and volume ratios are illustrated in Figure 2a. The corresponding results are listed in Table 1. The a-axis slightly increases as the pressure increases from 0.3 to 0.5 GPa, indicating that ZIF-4 has negative linear compressibility (NLC) in this pressure range. A recent study has shown that ZIF-75 exhibits similar negative linear compressibility under 0.1 GPa [42]. The NLC effect has also been investigated in more MOF-related studies, which are thought to be caused by the contraction of metal-ligand bonds and the tilting of ligands [43]. Drastic volume change occurs in structures at 0.3 to 0.4 GPa, where the cell volume is reduced from 4224.366 Å 3 to 3629.236 Å 3 . This mutation corresponds to a change in the structural characteristics of the middle range. It has also been shown experimentally that at such low pressures, the major cause of ZIF-4 compression is the contraction of free volume in the mid-range structure [32]. The a-axis slightly increases as the pressure increases from 0.3 to 0.5 GPa, indicating that ZIF-4 has negative linear compressibility (NLC) in this pressure range. A recent study has shown that ZIF-75 exhibits similar negative linear compressibility under 0.1 GPa [42]. The NLC effect has also been investigated in more MOF-related studies, which are thought to be caused by the contraction of metal-ligand bonds and the tilting of ligands [43]. Drastic volume change occurs in structures at 0.3 to 0.4 GPa, where the cell volume is reduced from 4224.366 Å 3 to 3629.236 Å 3 . This mutation corresponds to a change in the structural characteristics of the middle range. It has also been shown experimentally that at such low pressures, the major cause of ZIF-4 compression is the contraction of free volume in the mid-range structure [32].   The calculated values of transverse sound velocity Vt, longitudinal sound velocity Vl, average wave velocity Vm, and Debye temperature θD from 0 GPa to 1 GPa are listed in Table S1 (see Supplementary Materials) and Figure 2b. The Debye temperature increases with increasing pressure. It is difficult to compare calculated and experimental data since The calculated values of transverse sound velocity V t , longitudinal sound velocity V l , average wave velocity V m , and Debye temperature θ D from 0 GPa to 1 GPa are listed in Table S1 (see Supplementary Materials) and Figure 2b. The Debye temperature increases with increasing pressure. It is difficult to compare calculated and experimental data since there are no experimental data regarding ZIF-4 under the high pressure of the Debye temperature. The present results can be served as a prediction for future experiments. In Debye's theory, the θ D is the temperature of a crystal's highest normal mode of vibration, i.e., the highest temperature that can be achieved due to a single normal vibration. The longitudinal sound velocity V l is insensitive to pressure until 0.4 GPa, and thereafter it displays a positive pressure dependence. Similarly, the transverse sound velocity V t and average wave velocity V m also undergo a transition at 0.4-0.5 GPa. Combined with the dramatic decrease in structural parameters occurring at 0.4 GPa, we suggest that this may correspond to the reversible phase transition process of ZIF-4 found in the experiments [15]. Combined with the experimentally measured phenomenon of increasing sound velocity during the vitrification of ZIF-4, we suggest that this change correlates with the degree of ZIF-4 amorphization. The discussion of the structural transformation and mechanical property changes in this pressure range will be continued below. The elastic stiffness constants of solids are essential parameters because they can provide a link between the mechanical properties and dynamic information concerning the forces operating in solids, especially for the stability and stiffness of materials. Elastic constants reflect the resistance of a crystal to externally applied stress. Study of the elastic properties of materials is essential to understand the chemical bonds and the cohesion of materials. Thus, it is significant to investigate elastic constants to understand the mechanical properties of the ZIF-4 phase at different pressures. The elastic constants of a single crystal were obtained through geometry optimization. The ZIF-4 phase belongs to the orthorhombic crystal, which has nine independent elastic constants, i.e., C 11 , C 22 , C 33 , C 44 , C 55 , C 66 , C 12 , C 13 , and C 23 . In this work, the calculated elastic constants of ZIF-4 at different pressures are listed in Table S2. According to Mouhat et al., the necessary and sufficient Born criteria of mechanical stability for orthorhombic structures are as follows [44]: (1) In our results, all structures are mechanically stable up to 1 GPa according to these stability criteria.
Among the nine independent elastic constants, C 11 , C 22 , and C 33 characterize the incompressibility along the x, y, and z directions, respectively. For the independent elastic constants obtained with GGA-PBE potential, the ratio of values calculated under high pressure to those under zero pressure, C ij /C ij (0 GPa), are depicted in Figure 3a. It can be seen from Table S2 and Figure 3a that before 0.4 GPa, C 22 is the largest among C 11 , C 22, and C 33 . However, C 11 becomes larger than the other two constants when the pressure increases above 0.4 GPa. This indicates that below 0.4 GPa, ZIF-4 is more incompressible along the b-axis (y direction). The pressure leads to changes in the structure so that it is more difficult to be compressed along the a-axis (y direction) above 0.4 GPa. C 11 , C 12, and C 13 are more sensitive to the change in pressure than the other elastic constants. Additionally, the gradient of the constants increases between 0.4 GPa and 0.5 GPa, especially for C 11 , which indicates a fundamental change in the structure during this pressure range.
Similarly, the K/K 0 , G/G 0, and E/E 0 ratios are also calculated to ascertain the influence of pressure on elastic and mechanical properties, as illustrated in Figure 3b. The detailed values of K, G, and E are listed in Table S3. For the mechanical moduli of ZIF-4, the strengthening effect can be queued as E > K > G with increasing pressure. At different hydrostatic pressures, the shear modulus G is always smaller than bulk modulus K, which indicates the shear instability of orthorhombic ZIF-4. This is consistent with previous works that have claimed that the shear mode collapse is the main reason for pressure-induced amorphization [23,[45][46][47].
The hardness and brittleness of the compounds are represented as the ratio (K/G) between bulk modulus and shear modulus, according to Pugh's criterion [48]. Materials with good ductility usually have a high K/G value. The compound with a larger K/G ratio (>1.75) is usually ductile; otherwise, a K/G ratio < 1.75 is considered brittle. Cauchy pressure is another criterion used to determine structural brittleness/ductility. According to Pettifor's rule [49], materials with large positive Cauchy pressures have more metallic bonds and thus become more ductile. If Cauchy pressures of the materials are strongly negative, they possess more angular bonds and thus exhibit more brittleness. The results of the K/G ratio and the Cauchy pressure of ZIF-4 are listed in Table S3 and shown in Figure 3c. It is shown that ZIF-4 is brittle at zero pressure, but it becomes ductile with external pressure going up from 0.1 to 0.4 GPa. It is followed by a decrease in K/G ratio that occurs in the pressure range of 0.3-0.5 GPa. The K/G ratio drops below 1.75 again at 0.5 GPa, indicating that ZIF-4 exhibits brittleness at this pressure. As the external pressure increases from 0.5 to 1 GPa, the K/G ratio again increases with external pressure. At pressures greater Molecules 2023, 28, 22 5 of 11 than 0.6 GPa, ZIF-4 undergoes a second brittle-ductile transition and the K/G ratio increases above the critical value of 1.75.
Consistent with the trend of the change in the K/G ratio, the Cauchy pressure for ZIF-4 increases with pressure in both the 0-0.3 GPa and 0.4-1 GPa hydrostatic pressure change ranges. The trends of Cauchy pressure suggest an increase in the metallicity of bonding in the structure and ductility. However, the corresponding Cauchy pressure at 0.3 GPa is higher than at 0.4-0.8 GPa, suggesting that the change in bonding in the ZIF-4 structure with pressure is discontinuous. In conjunction with the structural changes discussed previously, the drastic volume change in the structure that occurs at 0.3-0.4 GPa essentially corresponds to a partial chemical bonding transition. Similarly, the K/K0, G/G0, and E/E0 ratios are also calculated to ascertain the influenc of pressure on elastic and mechanical properties, as illustrated in Figure 3b. The detaile values of K, G, and E are listed in Table S3. For the mechanical moduli of ZIF-4, th strengthening effect can be queued as E > K > G with increasing pressure. At differen hydrostatic pressures, the shear modulus G is always smaller than bulk modulus K, whic indicates the shear instability of orthorhombic ZIF-4. This is consistent with previou works that have claimed that the shear mode collapse is the main reason for pressur induced amorphization [23,[45][46][47].
The hardness and brittleness of the compounds are represented as the ratio (K/G between bulk modulus and shear modulus, according to Pugh's criterion [48]. Materia with good ductility usually have a high K/G value. The compound with a larger K/G rati (>1.75) is usually ductile; otherwise, a K/G ratio < 1.75 is considered brittle. Cauch pressure is another criterion used to determine structural brittleness/ductility. Accordin to Pettifor's rule [49], materials with large positive Cauchy pressures have more metall bonds and thus become more ductile. If Cauchy pressures of the materials are strongl negative, they possess more angular bonds and thus exhibit more brittleness. The resul of the K/G ratio and the Cauchy pressure of ZIF-4 are listed in Table S3 and shown i Figure 3c. It is shown that ZIF-4 is brittle at zero pressure, but it becomes ductile wit external pressure going up from 0.1 to 0.4 GPa. It is followed by a decrease in K/G rat that occurs in the pressure range of 0.3-0.5 GPa. The K/G ratio drops below 1.75 again 0.5 GPa, indicating that ZIF-4 exhibits brittleness at this pressure. As the external pressur In addition to the K/G ratio and Cauchy pressure, the Poisson's ratio is also commonly used to evaluate the materials' mechanical response under external stress. Typically, a larger Poisson's ratio usually corresponds to good ductility. The value of ν for covalent materials is small (ν = 0.1), and for ionic materials, a typical value is 0.25. For central force solids, the value of ν is between 0.25 and 0.5 [50]. The calculated results of Poisson's ratio for ZIF-4 are listed in Table S3  The lattice parameters and mechanical properties discussed above indicate that for ZIF-4, 0.4 GPa is a "watershed" in hydrostatic pressure rises from 0 to 1 GPa. To reveal changes in the lattice parameters and mechanical properties near 0.4-0.5 GPa, the structures of ZIF-4 at 0.1, 0.4, 0.5, and 1 GPa were further characterized. The radial distribution function (RDF) in Figure 4 shows a mid-range structural shift at 0.4-0.5 GPa. When the pressure is smaller than 0.4 GPa, the compression of the structure corresponds to the contraction of the pores, and the difference in the distribution of the Zn-Zn bond lengths is slight. When the pressure increases to 0.4-0.5 GPa, the change in the form of the midrange connections is reflected in the further reduction in Zn-Zn bonds. The Zn-Zn bond lengths at 0.9 GPa and above do not change further than at 0.5 GPa. In the short-range structure presented by the Zn-N pair RDF in Figure 4a, there is little change under different pressures. From this, we deduce that the anomaly of decreasing ductility and plasticity of the ZIF-4 structure with increasing pressure is due to the transformation in the way the [ZnN 4 ] tetrahedra are connected. The change in the mid-range order explains the sudden change in the mechanical properties of silicates with similar [SiO4] tetrahedral network structures [33,51].

Discussion
In this work, we have investigated the structural evolution and mechanical properties of ZIF-4 at low hydrostatic pressures (<1 GPa). For this purpose, the atomic structure and elastic constants of ZIF-4 have been investigated concerning external pressure ranging from 0 GPa to 1 GPa using the DFT calculation method. The phase stability, mechanical moduli, and Debye temperature are discussed. The elastic constants, mechanical moduli, and Debye temperature showed increases with increasing pressure. Furthermore, we analyzed the geometrical characteristics of the structure and found that the form of the connection of the [ZnN4] tetrahedra in the structure changes at 0.4 GPa. The pressure leads to a twisting of the Im ligands. This mid-range structural shift leads to the creation of intra-molecular forces between ligands, in addition to a reduction in free volume. Affected by the new mid-range structures, anomalies in the ZIF-4 elastic constants are produced at around 0.4 GPa. These include sudden increases in elastic moduli, negative linear compressibility, and the transition from ductile to brittle as reflected in K/G. This change in mid-range structure due to the twisting of polyhedral units has also been investigated in other simulations of the mechanical properties of tetrahedral network structures. An example is the effect of the [SiO4] tetrahedral connection on the anomalous mechanical properties of SiO2 glasses at extreme high pressures [33]. This suggests that it is instructive to study changes in the mid-range local structure of amorphous structures.
It is worth noting that errors between theoretical and experimental values are

Discussion
In this work, we have investigated the structural evolution and mechanical properties of ZIF-4 at low hydrostatic pressures (<1 GPa). For this purpose, the atomic structure and elastic constants of ZIF-4 have been investigated concerning external pressure ranging from 0 GPa to 1 GPa using the DFT calculation method. The phase stability, mechanical moduli, and Debye temperature are discussed. The elastic constants, mechanical moduli, and Debye temperature showed increases with increasing pressure. Furthermore, we analyzed the geometrical characteristics of the structure and found that the form of the connection of the [ZnN 4 ] tetrahedra in the structure changes at 0.4 GPa. The pressure leads to a twisting of the Im ligands. This mid-range structural shift leads to the creation of intra-molecular forces between ligands, in addition to a reduction in free volume. Affected by the new midrange structures, anomalies in the ZIF-4 elastic constants are produced at around 0.4 GPa. These include sudden increases in elastic moduli, negative linear compressibility, and the transition from ductile to brittle as reflected in K/G. This change in mid-range structure due to the twisting of polyhedral units has also been investigated in other simulations of the mechanical properties of tetrahedral network structures. An example is the effect of the [SiO 4 ] tetrahedral connection on the anomalous mechanical properties of SiO 2 glasses at extreme high pressures [33]. This suggests that it is instructive to study changes in the mid-range local structure of amorphous structures.
It is worth noting that errors between theoretical and experimental values are inevitable due to limitations in the accuracy of the theoretical simulation method and the size of the model used. Our analysis is based on the central idea that varied mechanical properties correlate with mid-range structural evolution. For amorphous or other structures that are not at the lowest point on the energy landscape map, the greatest difference in static structure from their crystalline state is reflected in the mid-range order. This difference causes the response of the structure under stress or other external fields to diverge from that of the crystalline state. For the crystal-amorphous transition in MOF and other complex systems, larger-scale ab initio simulations based on more advanced computers or algorithms are the focus of further research. Another possible research prospect is the establishment of correlations across multi-scale simulations, as has been completed in catalytic simulations [52,53]. By using structural features of smaller spatial and temporal scales as input, the establishment of correlations between simulation tools at different scales will have a positive effect on the experimental-simulation combination. Further exploration of the effects of different polyhedral unit connections on structural properties at different length scales will help to establish a clearer link between atomic-scale simulations and macroscopic-scale performance.

Materials and Methods
The first-principles method, based on density functional theory (DFT), was used in the geometry optimization and total energy calculation of the structures of ZIF-4. The original ZIF-4 structure was obtained from The Cambridge Crystallographic Data Centre (CCDC) [54] with the number VEJYUF01 [55], and the solvent molecules were removed. The projector augmented wave (PAW) method was implemented in the Vienna ab initio Simulation Package (VASP) code [56]. The electronic exchange and correlation potential were expressed using the Perdewe-Burkee-Ernzerh (PBE) version of the generalized gradient approximation (GGA) [57]. The cut-off energy for the plane wave basis set was 600 eV for all structures. Structural relaxations were performed to a tolerance of 1 × 10 −7 eV/atom in the total energy, yielding average forces of 1 × 10 −3 eV/Å. The stress level of the final equilibrium structure was less than 0.1 GPa. The relaxation of all the structures imposed no restrictions on the volume and lattice vectors. Since the unit cells were large, containing nearly 400 hundred atoms, only one K-point at Γ (0, 0, 0) was used.
The effects of external pressure were enforced by the stress tensor (PSTRESS keyword) corresponding to a selected pressure value in the range between 0 GPa and 1 GPa. Energies and enthalpies were calculated at multiple pressure points starting from the experimental structure, but both the atomic positions and the unit cell parameters relaxed during energy optimization under the corresponding pressures.
The single-crystal elastic constants, C ij , of the elasticity matrix (tensor) were then computed by performing six finite distortions of the lattice and deriving the elastic constants from the strain-stress relationship with the code in VASP [58,59].

of 11
G R = 15[4(S 11 + S 22 + S 33 ) + 3(S 44 + S 55 + S 66 ) − 4(S 12 + S 13 + S 23 )] −1 (7) where S ij represents the elastic compliance matrices. The Hill bulk modulus (K H ) and shear modulus (G H ) can be given as: The Young's modulus (E) and Poisson's ratio (ν) of the polycrystalline aggregate can be calculated via K H and G H , as follows: In order to comprehensively estimate the elastic anisotropy of these materials, the universal elastic anisotropy index A U is adopted, accounting for both the shear and bulk contributions [63,64]: The Debye temperature (θ D ) of materials plays an important role in many physical properties of solids, such as specific heat, elastic stiffness constants, and melting temperature. The experimental value of a solid can usually be calculated from the sound velocity. The Debye temperature of materials can be deduced from elastic constants. The Debye temperature of ZIF-4 was calculated using the average sound velocity (V m ) with the following equation [65]: where h is Planck's constant, k is Boltzmann's constant, N A is Avogadro's number, ρ (=M/V) is the density, M is the molecular weight, and n is the number of atoms per formula unit. The average wave velocity (V m ) in the polycrystalline material can be calculated with the following equation: where V t and V l are the transverse and longitudinal sound velocity, respectively, which are obtained from the values of Hill's bulk modulus K and shear modulus G from Navier's equation [65]: