# Computation of Finite Temperature Mechanical Properties of Zeolitic Imidazolate Framework Glasses by Molecular Dynamics

Nicolas Castel, François-Xavier Coudert

► To cite this version:

Nicolas Castel, François-Xavier Coudert. Computation of Finite Temperature Mechanical Properties of Zeolitic Imidazolate Framework Glasses by Molecular Dynamics. Chemistry of Materials, 2023, 35 (10), pp.4038-4047. 10.1021/acs.chemmater.3c00392. hal-04094425

HAL Id: hal-04094425

https://hal.science/hal-04094425

Submitted on 11 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Computation of Finite Temperature Mechanical
Properties of Zeolitic Imidazolate Framework
Glasses by Molecular Dynamics

Nicolas Castel†,‡ and François-Xavier Coudert*,†

†Chimie ParisTech, PSL Research University, CNRS, Institut de Recherche de Chimie
Paris, 75005 Paris, France

‡École des Ponts, 77420 Marne-la-Vallée, France

E-mail: fx.coudert@chimieparistech.psl.eu
Abstract

Mechanical properties of amorphous phases of metal-organic frameworks (MOF), such as MOF glasses, are difficult to determine experimentally. Moreover, computational characterization is limited by the level of theory chosen for the description of interatomic interactions and is often computationally expensive. In this work, we have extensively investigated the computation of finite temperature mechanical properties of ZIF-4 in the crystal and glass phases. We critically assessed computational methodologies including *ab initio* molecular dynamics, reactive force fields, and classical force fields, based on a variety of glass models. We find that ZIF-4 glasses have a larger bulk modulus than the crystal and confirm previous studies that the density is larger for the glass phases. Moreover, we confirm in the case of zeolitic imidazolate framework (ZIF) glasses the relationship between density and bulk modulus, showing that obtaining models of correct density is key to the prediction of physical properties for these systems.

Introduction

Because of their wide diversity and tunability, both structurally and chemically, metal-organic frameworks (MOFs) have been proposed over the past two decades for a large variety of applications such as gas storage, fluid separation, and heterogeneous catalysis.¹ Out of this vast family of nano- and meso-porous materials, composed of inorganic nodes connected by organic linkers to form three-dimensional architectures, several classes demonstrated high chemical and thermal stability,² such as zeolitic imidazolate frameworks (ZIFs).³ However, while their mechanical stability is also essential for MOFs to fully achieve their potential in industrial-scale processes – in particular for processing (extrusion and pellet formation) –,
the study of how these materials respond to mechanical stress is comparatively still emerg-
ing.⁴,⁵ Particularly promising in this regard, is the growing family of non-crystalline MOFs,
coming as a large variety of MOF glasses and MOF gels (forming aerogels or monoliths
upon drying) and which can be formed by various production routes.⁶ While conserving
many of the intrinsic advantages of their crystalline counterparts, these states yield the po-
tential of increased mechanical robustness and would allow for greater ease of processing,
notably by circumventing the performance drop due to the necessary densification of the
MOF powders.⁷,⁸

While a series of experimental techniques have successfully been employed to determine
the mechanical properties of crystalline MOFs, there is a lack of studies on the amorphous
phases.⁹ Some methods such as those relying on high-pressure X-ray diffraction¹⁰ are not
straightforwardly applicable to disordered materials, and many others require large bulk
glass samples which are challenging to prepare.⁹,¹¹ Fortunately, the same mechanical prop-
erties can be studied by computational simulations,¹² which have been used extensively
as a complement to experiments on MOF crystals and can lead to systematic studies of
the structure–property relationships.¹³,¹⁴ Despite the challenges in the determination of the
framework structure of amorphous states at the microscopic scale, several methodologies
are available to generate the atomistic models which are a prerequisite of any computa-
tional determination of the mechanical properties.¹⁵ One family of approaches to computing
finite temperature mechanical properties consists in using Molecular Dynamics (MD), ei-
ther to mimic *in silico* the high-pressure experiments¹⁶⁻¹⁸ or by using the strain-fluctuation
method.¹⁹,²⁰ Unlike most computational works which led to the determination of mechanical
properties at 0 K, often using the stress-strain approach,²¹ the temperature dependence of
these methods allows for a better prediction and understanding of the mechanical strength
and stability of the studied materials, as the mechanical properties at operational tempera-
tures may deviate significantly from those computed at 0 K.²² It is key for disordered MOFs
which can display complex pressure and temperature behaviors at both low and high tem-
perature, with the occurrence of multiple phase transitions.²³,²⁴

In this work, we applied a range of methods and MD schemes – *ab initio*, reactive force
fields, classical force fields – to extract the finite temperature mechanical properties of ZIF
crystals and glasses. We first report how *ab initio* MD can be used to obtain reference
values, which have never been determined for a MOF glass in the literature before. We then
investigate two MD schemes with lower computational cost: classical and reactive force fields
(FF). We conclude that reactive force fields are inapplicable to low-density models, while
classical force fields reproduce with good agreement *ab initio* results. We find that ZIF-4
glasses have a larger bulk modulus than the crystal, and confirm previous studies that the
density is larger for the glass phases. In addition, we confirm the persistent relationship
between density and bulk modulus, visible across every model, MD scheme and method.

**Systems & Computational methods**

Systems

Our systems throughout this study belong to two different phases of ZIF-4: the crystal and
the melt-quenched glasses.²⁵ The first amorphous MOF discovered and studied in detail,²⁶
ZIF-4 has since been the subject of numerous works and can be seen as a prototypical
amorphous ZIF system.¹⁵ It is built up from Zn²⁺ metal nodes and imidazolate (Im) organic
linkers, which are organized in the crystalline state as Zn(Im)₄ tetrahedra linked by Zn–N
coordinative bonds as illustrated on Figure 1. The amorphous phases can be formed from the
parent crystal by a variety of experimental methods, including ball milling, melt-quenching,
or pressure-induced amorphization.⁶,⁹

Out of the many finite temperature mechanical properties that characterize anisotropic
materials, in this study we will mostly focus on the bulk modulus (K) which characterizes
the variation of the volume (V) of a solid under uniform hydrostatic pressure (P) and is
defined as K = -V (∂P/∂V)ₜ. This property is the most readily available experimentally,
<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of the assembly of ZIF-4 as a three-dimensional network of Zn(Im)₄ tetrahedra. Reproduced from Ref. 15. Copyright 2022 American Chemical Society.</figcaption></figure>

and provides key information relating to applications, including about the MOF stability
during the shaping process<sup>27</sup> or pressure swing adsorption (PSA) cycles.<sup>4</sup>

To this date, multiple experimental and computational values for the ZIF-4 crystal have been determined and are listed in Table 1. *K* is comprised in the 1.4–2.7 GPa range for systems at a density (*ρ*₀) around the crystallographic density³ of 1.22 g cm⁻³, corresponding to a phase called open-pore (op) in the two latest experimental works that have focused on low-pressure behavior.²⁴,²⁸ Comparison between experiments and computational works is not straightforward as the measurements are complicated by subtle differences between MOFs from different batches (e.g. unintended defects incorporation during synthesis and different degrees of activation) and are dependent on the experimental setup (e.g. pressure transmitting medium).²⁸,²⁹ Additionally, computational values reported so far are either computed at 0 K³⁰⁻³² or performed at a low level of theory,²⁰,²³ with no *ab initio* value at finite temperature reported in the literature for the ZIF-4 crystal.

Though some mechanical properties of amorphous ZIF-4 such as the Young's modulus³⁵ or the fracture toughness³⁶ have been determined,³⁷ studies are significantly scarcer than for crystals. There are only two reported values of a bulk modulus, which stems from computational studies with the stress-strain approach performed on the same Continuous Random
<table><caption>Table 1: Bulk moduli ($K$) and densities at zero pressure ($ρ_0$) of crystalline ZIF-4, as reported in previous experimental (top panel) and computational (bottom panel) works. With the exception of values computed by stress-strain (at 0 K), they are measured or computed at room temperature. The two most recent experimental works identified two low-pressure phases, open-pore (op) and closed-pore (cp), which are both reported.</caption><thead><tr><th></th><th>Method</th><th>Ref.</th><th>K (GPa)</th><th>ρ<sub>0</sub> (g cm<sup>-3</sup>)</th><th>Phase</th></tr></thead><tbody><tr><td rowspan="3">Exp.</td><td>High-Pressure</td><td>33</td><td>2.6 ± 0.1</td><td>1.21</td><td></td></tr><tr><td rowspan="2">Crystallography</td><td rowspan="2">28</td><td>2.01 ± 0.05</td><td>1.22</td><td>(op)</td></tr><tr><td>4.39 ± 0.20</td><td>1.56</td><td>(cp)</td></tr><tr><td rowspan="2">Mercury intrusion</td><td rowspan="2">24</td><td>1.42</td><td>1.22</td><td>(op)</td></tr><tr><td>4.88</td><td>1.53</td><td>(cp)</td></tr><tr><td rowspan="4">Comp.</td><td rowspan="2">Strain-fluctuation with<br>classical FF<sup>34</sup></td><td>20</td><td>2.58</td><td>1.31</td><td></td></tr><tr><td>23</td><td>2.69</td><td>1.31</td><td></td></tr><tr><td rowspan="3">Stress-strain from<br>first principles<sup>21</sup></td><td>30</td><td>2.41</td><td>1.16</td><td></td></tr><tr><td>31</td><td>1.54</td><td>1.03</td><td></td></tr><tr><td>32</td><td>1.76</td><td>1.18</td><td></td></tr></tbody></table>

Network (CRN) model. A first work³⁸ found a bulk modulus of 8.88 GPa for a density of 1.07 g cm⁻³, while a second one³⁹ reported $K = 4.47$ GPa and $ρ_0 = 0.99$ g cm⁻³. However, the physical realism of the CRN model has not been thoroughly demonstrated in the first place,¹⁵ notably as it features a surprisingly low density.⁴⁰ Additionally, the mismatch in reported $K$ values was left unexplained in the latest study, leaving us with no reliable value.

## Molecular dynamics simulations

We describe here summarily the different methods for all types of molecular dynamics (MD) simulations used in this work. Additionally, representative input files for the MD simulations are available online in our data repository at https://github.com/fxcoudert/citable-data
Ab initio molecular dynamics

DFT-based MD simulations were performed using the Quickstep module⁴¹ of the CP2K software package,⁴² using parameters already fine-tuned for ZIF-4 in previous works.²³,²⁵,⁴³ We used the hybrid Gaussian and plane wave method GPW,⁴¹ the exchange-correlation energy was evaluated in the PBE approximation,⁴⁴ and the dispersion interactions were treated at the DFT-D3 level.⁴⁵ The multigrid system was set up with four different grids, a plane-wave cutoff for the electronic density of 600 Ry, and a relative cutoff of 40 Ry. Valence electrons were described by double-ζ valence polarized basis sets and norm-conserving Goedecker-Teter-Hutter⁴⁶ pseudopotentials, all adapted for PBE (DZVP-GTH-PBE) for H, C, and N or optimized for solids (DZVP-MOLOPT-SR-GTH) in the case of Zn.

The simulations were performed in the constant-volume (N, V, T) ensemble with a fixed size and shape of the unit cell. A timestep of 0.5 fs was used in the MD runs; the temperature was controlled by velocity rescaling⁴⁷ with a time constant of 1000 fs. Volumetric deformation consisted in an instantaneous and isotropic volume change by 2% from the previous volume, followed by a sufficient equilibration time of around ~100–200 ps.

Reactive force field

We employed a parametrization of ReaxFF, a flavor of reactive force fields,⁴⁸ developed for ZIF materials,⁴⁹ and used the same parameters as in our previous work.⁴⁰

ReaxFF simulations were performed using LAMMPS.⁵⁰,⁵¹ A timestep $\delta t$ of 0.25 fs was used in the MD runs and, the temperature (and pressure when applicable) were controlled using a Nosé-Hoover thermostat (and barostat). Temperature and pressure damping parameters $T_{damp}$ and $P_{damp}$ were fixed at 100 fs and 1000 fs respectively. Constant-pressure (N, P, T) simulations were performed using either an isotropic or a flexible cell, respectively, with the LAMMPS keywords `iso` and `tri`. Volumetric deformation in the (N, V, T) ensemble was a continuous process lasting $\tau_{\Delta V}$ = 250 ps, during which each dimension of the box changed linearly with time from its initial to final value to achieve a volume reduction by
1.5%. To reduce finite size effects, a (2 × 2 × 2) supercell of ZIF-4 with 2,176 atoms was simulated.

## Classical force field

We employed MOF-FF for ZIFs$^{52,53}$ and molsys, the python package behind MOF-FF for the parametrization of systems.$^{54}$ No adaptation was required for the ZIF-4 crystal, and two FF adaptations of MOF-FF for ZIFs were employed for the glasses: Zn3tetra and Zn3trig. In Zn3tetra, every three-coordinated Zn atom (denoted zn3_n3 in molsys atypes format) was treated as a four-coordinated Zn (denoted zn4_n4). In Zn3trig, zn3_n3 were treated as zn4_n4 in all their interactions with other atoms, except for the N–Zn–N angle which was changed from the tetrahedral angle (arccos(-1/3) ≈ 109.47°) to 120°, using the same spring constant.

As the identification of the building units (Zn atoms and imidazolates) of the glass models could not be performed by molsys, due to large fluctuations in the bond angles and lengths as well as the presence of defects, it was carried out with our python library aMOF,$^{40}$ interfaced with molsys. aMOF is available online on https://github.com/coudertlab/amof, where the version used for this work is v1.0.0.

Classical MD simulations were performed using LAMMPS,$^{50}$ using identical procedures, system size and parameters than with ReaxFF excepted for the four following parameters: $\delta t = 1$ fs, $T_{damp} = 1$ ps, $P_{damp} = 10$ ps and $\tau_{\Delta V} = 500$ ps.

## Extraction of mechanical properties

Two approaches were used in this work to obtain finite temperature mechanical properties, focusing on the example of the bulk modulus $K$, defined as:$^{55}$

$$K = -V \left( \frac{\partial P}{\partial V} \right)_T \qquad (1)$$
Values of K reported in this work correspond to the equilibrium volume.

## Finite difference methods

The finite difference approach consists in first establishing the P-V relationship, and then fitting it with an equation of state (EoS). P-V data is generated here by running multiple MD simulations, enforcing either various values of pressure (stress) or volume (strain), and measuring the other variable. We will talk about the finite strain difference method if the MD simulations are performed in the constant-volume (N, V, T) ensemble by enforcing the volume V, with the resulting pressure P(V) of the equilibrated system being measured. Alternatively, if they are performed in the constant-pressure (N, P, T) ensemble, with the volume V(P) being measured for a given P, we talk about the finite stress difference method.

In both cases, the well-behaved region of the P-V data was fitted with the second order Birch–Murnaghan EoS$^{55,56}$

$$P = 3K f_E (1 + 2f_E)^{5/2} \qquad (2)$$

where $V_0$ is the volume at zero pressure and $f_E$ the Eulerian strain defined as

$$f_E = \frac{1}{2} \left[ \left( \frac{V_0}{V} \right)^{2/3} - 1 \right] \qquad (3)$$

In addition to providing K, this EoS fit returns $V_0$ and thus the density at zero pressure $\rho_0$.

We systematically tried to use the third order Birch–Murnaghan EoS, and found that it leads to overfitting in the large majority (> 80%) of cases due to the small number of P-V data points, therefore it is not reported in this work.

## Strain-fluctuation method

Mechanical properties at a given temperature T can also be evaluated from the fluctuations of a system at equilibrium, simulated under a constant stress in the (N, σ, T) ensemble.$^{19}$
The bulk modulus can be obtained directly from the fluctuations of the volume:¹⁹

$$ \langle (\Delta V)^2 \rangle = \frac{k_B T \langle V \rangle}{K} \qquad (4) $$

where $k_B$ is the Boltzmann constant, $\Delta V = V - \langle V \rangle$, and $\langle X \rangle$ denotes the time average of any quantity $X$.

With this method, it is possible to obtain more than the bulk modulus, and access anisotropic mechanical properties. The elastic stiffness tensor $C$ can be obtained from the fluctuations of the unit cell matrix $\mathbf{h}$ through the following relation:¹⁹,²⁰

$$ \left( \frac{k_B T}{V} \right) C^{-1}_{ijkl} = \langle \epsilon_{ij} \epsilon_{kl} \rangle - \langle \epsilon_{ij} \rangle \langle \epsilon_{kl} \rangle \qquad (5) $$

In this equation, $C_{ijkl}$ are the components of the fourth-order tensor $C$, called elastic constants, and $\epsilon$ is the unit cell strain defined by:⁵⁷

$$ \epsilon = \frac{1}{2} \left( (\mathbf{h}_0^T)^{-1} \mathbf{h}^T \mathbf{h} \mathbf{h}_0^{-1} - 1 \right) \qquad (6) $$

with $\mathbf{h}_0$ the reference unit cell, corresponding to the first frame of the equilibrated trajectory.

From this elastic stiffness tensor, multiple mechanical properties — bulk modulus $K$, Young's modulus $E$, shear modulus $G$, and Poisson's ratio $\nu$ — were computed using the ELATE code.⁵⁸ Values reported in this work were obtained with the Hill averaging scheme.⁵⁹

The computation of $C$ was performed using the Python library `aMOF`, in which we adapted and implemented a code developed for previous works.²⁰,²³

## Structural properties

Structural analyses were performed using the parameters used and further detailed in a previous study,⁴⁰ using the Python library `aMOF`.

Unless explicitly stated in the caption, properties are averaged over an $(N, V, T)$ trajec-
tory of 1 ns for classical MD simulations and 60 ps for *ab initio* simulations. For classical MD, frames are taken every 1 ps for the radial distribution function (RDF) and bond angle distribution, and every 50 ps for pore statistics. For *ab initio* MD, those intervals are respectively of 0.5 fs and 0.5 ps.

The potential of mean force (PMF)⁶⁰ is computed from the RDF $g(r)$ through the relation

$F(r) = -k_B T (\ln g(r) - \ln g_{max})$, where $g_{max}$ is the maximum of $g(r)$ over $r$.

N–Zn–N angle distribution is computed by taking a cutoff radius of 2.5 Å for Zn–N dis-
tances, a value determined based on the Zn–N PMF, and validated in previous *ab initio* ²⁵,⁴³,⁶¹
and ReaxFF⁴⁰ studies.

The total pore volume is computed on individual frames using the Zeo++ software, ⁶²⁻⁶⁴
as the sum of accessible and non-accessible volume with a helium probe of radius 1.2 Å.

**Results & Discussion**

***Ab initio*** *molecular dynamics*

*Ab initio* methods, and in particular density functional theory (DFT)-based methods, allow for a full description of the electronic state of physical systems at the quantum chemical level. They are commonly used to compute mechanical properties at 0 K using the stress–strain approach,²¹ and have been applied to many MOFs, including for the ZIF-4 crystal³⁰⁻³² (see Table 1). Combining the precision of *ab initio* methods with molecular dynamics allows the determination of finite temperature mechanical properties and can account for potential phase transitions that may occur under the application of both temperature and pressure.²³,²⁴,²⁸ Already applied to molten salts¹⁸,⁶⁵ and metals,¹⁷ *ab initio* molecular dynamics (AIMD) has never — to our knowledge — been used to this aim for any MOF system. Despite its significant computational cost, AIMD can provide much needed reference values for ZIF glasses, where there are no experimental results currently available.

In prior work, AIMD has successfully been used to model the melting of the ZIF-4 crys-
tal,²⁵ and to produce atomistic configurations of melt-quenched glasses.⁴³ All these models have been generated and equilibrated in the (N, V, T) ensemble, and thus are all at the same density (taken to be the crystallographic density 1.22 g cm⁻³). Although this difficulty in capturing changes in density has been identified as an important limitation of the *ab initio* methodology,²⁵,⁴⁰ such glass models still represent the most chemically accurate atomistic description of ZIF-4 glasses published to date in the literature.¹⁵

Here, we have computed the bulk modulus *K* for four of these systems: the crystal and three melt-quenched glasses, chosen to preserve the diversity found in the coordination environments. As only three different discrete values were found for the Zn–N coordination numbers of the 10 original glass models, we chose one model per value. This restriction was guided by the significant computational cost of *ab initio* MD.

Finite strain difference method

The finite difference approach to compute *K* consists in producing *P* − *V* data, and fitting the results with an equation of state (EoS). Widely performed experimentally,²⁸,³³ it is also within the reach of computational methods by running multiple MD simulations at various values of pressure or volume. This computational approach has been used on multiple materials such as molten salts,¹⁸,⁶⁵ ceramics⁶⁶ or silicates,¹⁶ and can be used with MD schemes ranging from *ab initio*¹⁷,¹⁸,⁶⁵ to classical¹⁶,⁶⁶ and reactive⁶⁷ force fields.

Due to the difficulty of performing (N, P, T) *ab initio* simulations of soft porous crys-
tals<sup>40,68</sup> and to the reduced computational cost of (N, V, T) MD runs, we employed the finite
strain difference method. All *ab initio* simulations were performed in the (N, V, T) ensem-
ble to enforce a volume *V* and compute the resulting pressure *P*(*V*) of the equilibrated
system. Successive simulations at different volumes were performed, consisting of isotropic
volume change by 2% from the previous volume, followed by an equilibration time of around
~ 100–200 ps. Equilibration was monitored by examining the pressure as a function of time
with a sufficiently large moving average of 5 ps, as illustrated on Figure 2. Convergence was
determined by the relative absence of drift in the pressure averaged over a larger period of
time of 50 ps, which led to longer simulations for larger deformations. The final value of the
pressure for a given volume was taken as the average over the last 50 ps of the simulation.
To estimate the uncertainty due to the fluctuations and limited simulation time, the last
50 ps were divided into blocks of 5 ps. The average pressure of each block was computed,
and a standard deviation was calculated on these 10 values. Both the bulk modulus *K* and
the extrapolated density at zero pressure *ρ*₀ were obtained by fitting the *P* − *V* data with
the second order Birch-Murnaghan EoS.⁵⁵,⁵⁶ Results are shown on Table 2. The *P* − *V* plot
with the standard deviations as well as the fitted EoS is shown on Figure 3 in the case of
the crystal. We note that despite the large pressure fluctuations of Figure 2, inherent to the
computation of stresses in molecular dynamics simulations, we obtain a reasonable averaged
pressure for each volume.

<figure><img src="image_3.png" /><figcaption>Figure 2: Pressure as a function of time for different deformations of the crystal. A moving average over 5 ps is used and deformation is defined as relative volume change to the reference volume of the crystal.</figcaption></figure>

The same procedure was then used for the glass models, although the complex energy
landscapes of ZIF systems — prone to polyamorphism under the application of temperature⁶⁹
or pressure⁷⁰,⁷¹ — meant that it was not straightforward to observe a well-behaved *P*(*V*)
regime for a large range of volumes. As illustrated on Figure S1, a structural rearrangement
<figure><img src="image_4.png" /><figcaption>Figure 3: Pressure as a function of volume for the crystal, fitted with the second order Birch-Murnaghan Equation of State (red). Each volume corresponds to a (N, V, T) simulation with fixed deformation. Pressure is computed as the average over the last 50 ps of simulations. Error bars are the standard deviations of the pressure averaged with a rolling window of 5 ps.</figcaption></figure>

<table><caption>Table 2: Bulk moduli *K* and densities at zero pressure *ρ*₀ for the ZIF-4 crystal and glasses obtained by AIMD⁴³ with the finite strain difference method. Glass values are averaged over 3 models.</caption><thead><tr><th></th><th>Crystal</th><th>Glass</th></tr></thead><tbody><tr><td><i>K</i> (GPa)</td><td>1.39</td><td>2.43 ± 0.09</td></tr><tr><td><i>ρ</i><sub>0</sub> (g cm<sup>-3</sup>)</td><td>1.29</td><td>1.34 ± 0.03</td></tr></tbody></table>

of the system caused by a large strain can lead to a change in the final pressure, although this structural change is not visible in local structural properties.

## Reference values

As shown on Table 2, we obtain a bulk modulus of $K_{\text{crystal}} = 1.39$ GPa, in reasonable agreement with previous experimental$^{24,28,33}$ and computational$^{30-32}$ works. It is in the lower range of values (see Table 1), but further interpretation is limited, as the precise value determined depends on the computational methodology and – for experiments – also on subtle differences between MOFs from different batches or on the experimental setup.$^{28,29}$ However, the comparison between two different systems within the same methodology is meaningful, and
we find that the ZIF-4 glasses are less soft with $K_{glass} = 2.43 \pm 0.09$ GPa. This new result complements a previous nanoindentation study³⁵ that found that ZIF-4 melt-quenched glasses showed larger Young’s modulus ($E_{glass} = 8.2$ GPa) than in their parent crystalline phase ($E_{crystal} = 4.6$ GPa).⁷²

Interestingly, this finite strain difference approach also yields the density at zero pressure $\rho_0$ which is of $1.29 \text{ g cm}^{-3}$ for the crystal and of $1.34 \pm 0.03 \text{ g cm}^{-3}$ for the glass. We note that the glass density is higher than for the crystal, as expected from pycnometric measurements³⁵ ($\rho_{glass} = 1.63 \text{ g cm}^{-3}$ and $\rho_{crystal} = 1.50 \text{ g cm}^{-3}$) and from a recent CO₂ physisorption study⁷³ ($\rho_{glass} = 1.38 \text{ g cm}^{-3}$ which allows direct comparison to the crystallographic density³ $\rho_{crystal} = 1.22 \text{ g cm}^{-3}$).⁴⁰ However, the difference is not as large as what could be expected from these experiments and may come from the *ab initio* methodology, which led to glass models of the same density as the crystal before any deformation step.

In spite of the limitations discussed above, we find that *ab initio* simulations provide a reference value for the bulk modulus at room temperature of ZIF-4 glasses, which may be thought of as a lower bound due to the low density of the particular glass models used. We note that such a quantity had never been determined before at the quantum chemical level. By comparing two states of similar densities, it further demonstrates the influence of the system topology on the mechanical properties of ZIF-4 systems. It complements previous results on the high-density states which found the closed-pore ZIF-4 phase bulk modulus ($\sim 4.4$ to $4.9$ GPa)²⁴,²⁸ to be much lower than its polymorph ZIF-zni ($\sim 14$ GPa),⁷⁴ while they all show similar densities ($1.53$ to $1.56 \text{ g cm}^{-3}$).

## Reactive force fields

With a reference value now at hand, it is possible to develop and validate more computation-
ally efficient alternatives for the determination of finite temperature mechanical properties.
Although classical MD simulations are routinely used for crystalline MOFs,²⁰ they are un-
able by design to simulate bond breaking or reformation, a key mechanism in the formation
of ZIF-4 glasses by melt-quenching.²⁵ To develop a method able to screen the mechanical properties of a large number of ZIF glasses, it would be convenient and consistent to use the same MD scheme for the glass model generation and the subsequent determination of its properties. An option is the use of reactive force fields, which have been proposed to generate models of ZIF glasses by melt-quenching with the development of a ReaxFF parametrization for ZIFs.⁴⁹ They have previously been used to compute mechanical properties³⁶,⁷⁵,⁷⁶ and are a natural candidate for the task at hand. Able to explore larger spatial and time scales, they could in principle be used to obtain the bulk modulus with several methods while reducing the finite size effects inherent to AIMD. However, we showed in a previous work the atypical structural properties of the glasses obtained with this approach,⁴⁰ and there is still no in-depth validation of the mechanical properties obtained with this force field in the existing literature.¹⁵

In order to offer a meaningful comparison to other MD schemes, we consider several glass models: a ReaxFF glass obtained in our previous work,⁴⁰ the same three *ab initio* glass models studied in the previous section⁴³ and a glass obtained by Reverse Monte Carlo (RMC) modeling.²⁵,⁷⁷ These models, as well as the crystalline system, went through a careful preparation process⁴⁰ and are thoroughly presented in the supporting information.

## Finite strain and finite stress difference methods

To allow direct comparison to our *ab initio* results, we first computed *K* with the finite strain difference method. The reduced computational cost allowed for smoother volumetric deformations (1.5% change over 250 ps) and a longer equilibration (250 ps, shown on Figure S3) than was possible with AIMD. This approach led to well-behaved *P* − *V* plots and to the values presented in Table 3.

Although of reasonable density, the ReaxFF crystal shows a larger bulk modulus by comparison with AIMD and experimental results. The bulk moduli and densities reported for the *ab initio* glasses do not contrast with the values computed for the crystal as clearly
than with AIMD. The glass models which display higher *K*, represent states of considerably higher densities than the *ab initio* glass model and than the experimentally determined crystallographic density of ZIF-4 glasses,⁷³ and are therefore not directly comparable.

<table><caption>**Table 3:** Bulk moduli *K* and densities at zero pressure *ρ*₀ for multiple ZIF-4 models obtained with the finite strain difference method with ReaxFF.</caption><thead><tr><th></th><th>Crystal</th><th>ReaxFF glass</th><th>RMC glass</th><th>Ab initio glass</th></tr></thead><tbody><tr><td><i>K</i> (GPa)</td><td>3.15</td><td>10.36</td><td>6.29</td><td>3.57 ± 0.75</td></tr><tr><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td>1.26</td><td>1.67</td><td>1.55</td><td>1.31 ± 0.03</td></tr></tbody></table>

An alternative option is the use of the finite stress difference method, in which the MD simulations are performed in the constant-pressure (*N*, *P*, *T*) ensemble at a given *P*, before measuring the volume *V*(*P*). It is illustrated on Figure S4. This approach is more commonly used in conjunction with classical and reactive force fields.¹⁵,⁶⁷ We evaluate the method using two possible variants, with an isotropic or a flexible cell (respective LAMMPS keywords *iso* and *tri*), as a previous work⁴⁰ showed that this choice can impact the description of the glass.

While this method yields results for every studied system, as reported on Table 4, the significant mismatch between values computed with different cell types, which should not exist for the crystal, severely questions the results and illustrate once again the tendency of ReaxFF to densify systems.⁴⁰ Additionally, the only two low-density states, the crystal and the *ab initio* glasses with an isotropic cell, have lower values than what was obtained with AIMD, the opposite of what yielded the finite strain difference method with ReaxFF. This seriously discourages using the finite stress difference method with ReaxFF.

## Strain-fluctuation method

Another approach to obtain mechanical properties at a given temperature *T* is to evaluate them from the fluctuations of the unit cell of an equilibrated system simulated under a constant-stress (*N*, *σ*, *T*) ensemble.¹⁹ Unlike the previous finite difference methods which
<table><caption>Table 4: Bulk modulus *K* and density at zero pressure *ρ*₀ for multiple ZIF-4 models obtained with the finite stress difference method with ReaxFF, either using an isotropic (iso) or a flexible (tri) cell.</caption><thead><tr><th></th><th>cell</th><th>Crystal</th><th>ReaxFF glass</th><th>RMC glass</th><th>Ab initio glass</th></tr></thead><tbody><tr><td rowspan="2">K (GPa)</td><td>iso</td><td>0.91</td><td>7.64</td><td>1.83</td><td>1.87 ± 0.75</td></tr><tr><td>tri</td><td>6.11</td><td>10.60</td><td>7.89</td><td>4.07 ± 1.52</td></tr><tr><td rowspan="2">ρ₀ (g cm<sup>-3</sup>)</td><td>iso</td><td>1.31</td><td>1.74</td><td>1.50</td><td>1.36 ± 0.06</td></tr><tr><td>tri</td><td>1.61</td><td>1.68</td><td>1.57</td><td>1.46 ± 0.08</td></tr></tbody></table>

only yield one property at a time, it leads to the estimation of the entire tensor of second-order elastic constants $C_{ijkl}$ which in turns is linked to mechanical properties such as Young's modulus, shear modulus, and Poisson's ratio.⁵⁸ This approach, already used for the ZIF-4 crystal with a classical force field,²⁰,²³ requires long equilibration times (~5 – 10 ns). While out of reach of AIMD simulations, it is tractable with ReaxFF. The bulk modulus can also be directly computed from the fluctuations of the volume,¹⁹ a second method we used as a validation.

We performed (*N*, *P*, *T*) simulations with a flexible cell (LAMMPS keyword *tri*), corresponding to the (*N*, *σ*, *T*) ensemble required by the method, until convergence of the volume (see Figure 4) and elastic constants (see Figure S5), both computed over the last 2.5 to 3 ns. It required a total equilibration time comprised between 5 ns and 22 ns, as detailed in the supporting information.

From the resulting bulk moduli, shown on Table 5, we see that both strain-fluctuation methods are consistent. While such a long ReaxFF equilibration only yields highly densified states, with large *K* as a consequence, the results are consistent with the high-density states explored with finite difference methods.

This additional method confirms the very direct relationship between large *K* and *ρ*₀, apparent on Figure 5, which is found almost irrespective of the model — crystal or glass — and of the method used to compute *K*. It is therefore difficult to separate the influence of the topology (e.g. between different glasses) from the density, or to study low-density glass
<figure><img src="image_5.png" /><figcaption>Figure 4: Volume as a function of time during a long 5 ns equilibration in the (N, P, T) ensemble with a flexible cell for the crystal. Elastic constants were computed on the red part of the figure. The orange line shows a moving average over 100 ps.</figcaption></figure>

<table><caption>Table 5: Bulk moduli *K* and densities at zero pressure *ρ*₀ for multiple ZIF-4 models obtained with two strain-fluctuation methods with ReaxFF.</caption><thead><tr><th></th><th>Crystal</th><th>ReaxFF glass</th><th>RMC glass</th><th><i>Ab initio</i> glass</th></tr></thead><tbody><tr><td><i>K</i> - volume fluctuations (GPa)</td><td>5.65</td><td>12.73</td><td>7.47</td><td>4.30 ± 0.64</td></tr><tr><td><i>K</i> - elastic constants (GPa)</td><td>7.58</td><td>12.85</td><td>7.62</td><td>4.77 ± 0.88</td></tr><tr><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td>1.61</td><td>1.68</td><td>1.56</td><td>1.46 ± 0.07</td></tr></tbody></table>

models. As a consequence, we cannot recommend the use of the ReaxFF parametrization for ZIFs⁴⁹ to compute mechanical properties, especially at finite temperature, and suggest that investigations which intend to employ it should first go through a detailed validation.

## Classical Force Fields

Considering the strong limitations of the ReaxFF parametrization for ZIF glasses, no computationally efficient MD scheme published to this date in the literature can both generate ZIF melt-quenched glasses and provide their finite temperature mechanical properties. In this section, we look into a possible strategy to separate these two aspects and investigate the use of classical force fields (FF) to analyze the crystal and glass models generated by AIMD²⁵,⁴³ studied in the first section. Owing to the computational efficiency of classical MD,
<figure><img src="image_6.png" /><figcaption>Figure 5: Bulk modulus *K* as function of the density at zero pressure *ρ*₀ for multiple ZIF-4 models and methods with ReaxFF.</figcaption></figure>

all 10 configurations of the melt-quenched glasses were simulated. The preparation of these
systems is detailed in the supporting information.

Adaptation of the Force Fields

Even if a few studies used classical MD to study the amorphization of ZIFs,²⁰,⁷⁸ every classical FF for ZIFs published to this date has been developed for crystals¹⁵ and should thus first be adapted and validated before studying ZIF glasses. In our study, we chose to use the recently developed MOF-FF for ZIFs by Dürholt et al.⁵³ This force field showed good transferability across a series of ZIF polymorphs, indicating that it can handle various network topologies, and led to mechanical properties in line with earlier FFs.⁷⁹ Developing a new FF for glasses from scratch is beyond the scope of this study, but the MOF-FF methodology for an automatic parametrization of FFs would make such a development achievable in future work.⁵²

While ZIF-4 crystals are made only of Zn(Im)₄ tetrahedra (see Figure 1), a number of three-coordinated Zn ions can be found in ZIF-4 glasses, which are not readily parametrized in MOF-FF. On the N–Zn–N bond angle distribution of the *ab initio* glass shown on Figure 6, we see that these three-coordinated Zn deviate from the 109° angle of the tetrahedra. To reproduce the glass properties and test the sensitivity of our results to our adaptation of
MOF-FF, we introduced and evaluated two adaptations of the FF, which we name **Zn3tetra** and **Zn3trig** which respectively shift the N–Zn–N angle of the three-coordinated Zn towards 109° or 120° as shown on Figure 6. As it features no under-coordinated Zn atoms, the ZIF-4 crystal is simulated with the original FF.

<figure><img src="image_7.png" /><figcaption>Figure 6: Distribution of the N–Zn–N angle for (a) every Zn atom and (b) only three-coordinated Zn atom, for the MOF-FF crystal (blue) and glasses (red and orange) compared to the *ab initio* crystal (green) and glass (grey).</figcaption></figure>

In order to evaluate the validity of these FF adaptations, we compared the structural characteristics of ZIF-4 systems simulated with MOF-FF to the original *ab initio* models, using the same metrics as in our previous work on ReaxFF.⁴⁰ An investigation of the local order with the partial radial distribution functions (RDF) and potentials of mean force (PMF) (see Figure S6), evidences an excellent reproduction of interatomic distances and similar energy landscapes close to the energy minima. The region between the two minima in PMF
is qualitatively different, as MOF-FF is unable by design to simulate non-bonded but nearby Zn-N atoms. Despite differences in the angle distribution of the three-coordinated Zn atoms, Figure 6 shows that the total distribution is similar for both FF adaptations, which reproduce the wider distribution of angles for the glasses than for the crystal. Although affected by the small difference in geometry imposed by the FF adaptations, both FF adaptations lead to glass models of higher porosity than the crystal, as shown on Figure S7. Finally, as classical MD cannot simulate bond breaking, the coordination and topology (e.g. ring statistics) are automatically preserved. Despite not being developed for ZIF glasses, both FF adaptations reasonably reproduce the structural properties of the *ab initio* systems, at least compared to the only force field available in the literature for amorphous ZIFs, namely ReaxFF.⁴⁰

## Validation with AIMD results

We investigate the applicability of these classical FF adaptations by computing finite temperature bulk moduli, with the same methods we used on ReaxFF, starting with the finite strain difference method (illustrated on Figure S8) to allow direct comparison with AIMD results. The parameters used for each method are detailed in the supporting information.

<table><caption>**Table 6: Bulk modulus *K* and density at zero pressure *ρ*₀ for ZIF-4 crystal and *ab initio* glasses obtained with the finite strain difference method with two different adaptations of MOF-FF.**</caption><thead><tr><th></th><th>Crystal</th><th>Glass w/ Zn3trig</th><th>Glass w/ Zn3tetra</th></tr></thead><tbody><tr><td><i>K</i> (GPa)</td><td>1.48</td><td>3.01 ± 1.15</td><td>2.46 ± 0.65</td></tr><tr><td>ρ<sub>0</sub> (g cm<sup>-3</sup>)</td><td>1.22</td><td>1.29 ± 0.06</td><td>1.31 ± 0.05</td></tr></tbody></table>

From the results, summarized in Table 6, we see that both FF adaptations lead to larger bulk moduli for the glasses than the crystal. As shown on Figure 7, *K* values are in reasonable agreement with AIMD, where ReaxFF failed to differentiate the crystal from the glasses. The densities are also consistent with both AIMD and experiments.

We note, however, that there is an important variance in the bulk moduli of the glasses,
<figure><img src="image_8.png" /><figcaption>Figure 7: Bulk modulus *K* for ZIF-4 crystal and *ab initio* glasses with the finite strain difference method and different MD schemes. Average and standard deviation for classical FF are computed over every glass and FF adaptation.</figcaption></figure>

with potentially divergent values with respect to the FF adaptation for the same glass model,
as shown on Figure S9. It highlights the importance of using a large enough number of glass
configurations to obtain average mechanical properties for models with a small system size.

Exploring other methods

We then investigated the use of the (*N*, *P*, *T*) ensemble with our two FF adaptations, again
with both an isotropic and a flexible cell as done for ReaxFF (see Figure S10). We see on
the results presented in Table 7 that no system reaches a state of unreasonably high density
(unlike with a reactive force field). In particular, MOF-FF handles better the additional
degree of freedom brought by the free shape of the cell than does ReaxFF, although a flexible
cell still leads to a larger densification, particularly for the crystal. We note that both FF
adaptations give similar average results, with still a significant variance over the different
glass models. Finally, *K* values are consistent with the finite strain difference method, and
the previously observed trends of *K*<sub>glass</sub> > *K*<sub>crystal</sub> and increasing *K* with respect to *ρ*<sub>0</sub> are
found again.

Finally, to obtain additional mechanical properties, we applied the strain-fluctuation but
found that it fails to highlight the difference in *K* between the crystal and the glasses as
shown on Figure 8, which may be explained by the unreasonably large density of the crystal
<table><caption>Table 7: Bulk modulus *K* and density at zero pressure *ρ*₀ for ZIF-4 crystal and *ab initio* glasses obtained with the finite stress difference method with two different adaptations of MOF-FF, either using an isotropic (iso) or a flexible (tri) cell.</caption><thead><tr><th></th><th>cell</th><th>Crystal</th><th>Glass w/ Zn3trig</th><th>Glass w/ Zn3tetra</th></tr></thead><tbody><tr><td rowspan="2">K (GPa)</td><td>iso</td><td>1.00</td><td>2.16 ± 0.89</td><td>2.24 ± 0.94</td></tr><tr><td>tri</td><td>1.80</td><td>2.77 ± 0.98</td><td>2.59 ± 1.04</td></tr><tr><td rowspan="2">ρ₀ (g cm<sup>-3</sup>)</td><td>iso</td><td>1.21</td><td>1.33 ± 0.08</td><td>1.31 ± 0.06</td></tr><tr><td>tri</td><td>1.39</td><td>1.35 ± 0.08</td><td>1.35 ± 0.06</td></tr></tbody></table>

equilibrated in the (*N*,*P*,*T*) ensemble with a flexible cell. The strain-fluctuation method with MOF-FF is detailed in the supporting information.

<figure><img src="image_9.png" /><figcaption>Figure 8: Bulk modulus *K* for ZIF-4 crystal and glasses with different methods and two different adaptations of MOF-FF.</figcaption></figure>

Despite this challenge, all finite difference methods are overall consistent, as illustrated on Figure S9, albeit with a large variability across glass configurations. Both FF adaptations yield similar results, suggesting they are not ill adapted and that this procedure is not very sensitive to such a small perturbation in the MOF-FF parametrization. We thus showed that the finite difference methods with classical FFs are ready to be deployed, and that the development of a specific MOF-FF parametrization for amorphous ZIFs could prove useful.
Conclusions and perspectives

In this work, we have extensively investigated the computation of finite temperature mechan-
ical properties of ZIF-4 in the crystal and glass phases. We have compared different molecular
dynamics schemes and computational methods on a series of glass models. This is important
because such mechanical properties are difficult to access experimentally, especially for the
glasses.

Demonstrating for the first time the applicability and accuracy of the finite strain differ-
ence method with *ab initio* molecular dynamics to MOF systems, we provide the first reliable
value for the bulk modulus of an amorphous ZIF in the literature. We find that ZIF-4 glasses
have a larger bulk modulus than the crystal, and confirm previous results which found the
density to be larger for the glass phase. Not only does this study provide intrinsically valuable
data, it yields reference values for the computation of finite temperature mechanical proper-
ties, making it possible to validate alternative methods with less computationally expensive
MD schemes.

We have investigated the use of ReaxFF, the only alternative MD scheme present in the
literature that could be used to both generate melt-quenched glass models and study their
mechanical properties. We have reported and analyzed a tendency of ReaxFF simulations
to densify the systems, particularly if they are performed in the constant-pressure (N, P, T)
ensemble, making it all but inapplicable to low-density models. Additionally, we have shown
that the bulk moduli reported are primarily a function of the density and not of its phase or
topology. All these observations suggest that mechanical properties obtained from the use
of the ReaxFF force field for ZIFs should be interpreted with caution, and makes a strong
case for the use of alternative methodologies, or for the further optimization of the ReaxFF
force field.

Finally, we examined the use of classical force fields for the computation of mechanical
properties of models created through other methodologies, demonstrating that all finite dif-
ference methods are consistent, and that the values are in good agreement with AIMD. We
have thus shown that this method is ready to be deployed, and that the development of a
specific MOF-FF parametrization for amorphous ZIFs could prove useful.

However, this last approach requires glass models to be first generated by other modeling
strategies and by itself cannot lead to a systematic screening of the mechanical properties of
numerous ZIF glasses. A promising integrated strategy is the development of machine learnt
(ML) potentials,⁸⁰,⁸¹ that could constitute a new generation of specific and accurate reactive
potentials which could generate glasses by melt-quenching. The methods presented in this
work to compute their finite temperature mechanical properties could then be carried out
with the same ML potentials, using for validation the AIMD reference value determined in
this study.

**Supporting Information Available**

System description and methodological details for ReaxFF and classical MD, strain-fluctuations
method with MOF-FF, volume, pressure and elastic constant convergence, *P* − *V* data, ad-
ditional mechanical properties, radial distribution functions, potentials of mean force, angle
distributions, porous volumes.

Acknowledgement

We thank Kim Jelfs and Tom Bennett for discussions and ongoing collaboration on this topic.
We acknowledge access to high-performance computing platforms provided by GENCI grant
A0130807069.

**References**

(1) Furukawa, H.; Cordova, K. E.; O'Keeffe, M.; Yaghi, O. M. The Chemistry and Applications of Metal-Organic Frameworks. *Science* **2013**, *341*, 673.
(2) Healy, C.; Patil, K. M.; Wilson, B. H.; Hermanspahn, L.; Harvey-Reid, N. C.; Howard, B. I.; Kleinjan, C.; Kolien, J.; Payet, F.; Telfer, S. G. et al. The thermal stability of metal-organic frameworks. *Coord. Chem. Rev.* **2020**, *419*, 213388.

(3) Park, K. S.; Ni, Z.; Côté, A. P.; Choi, J. Y.; Huang, R.; Uribe-Romo, F. J.; Chae, H. K.; O'Keeffe, M.; Yaghi, O. M. Exceptional chemical and thermal stability of zeolitic imidazolate frameworks. *Proc. Natl. Acad. Sci. U.S.A.* **2006**, *103*, 10186–10191.

(4) Yang, K.; Zhou, G.; Xu, Q. The elasticity of MOFs under mechanical pressure. *RSC Adv.* **2016**, *6*, 37506–37514.

(5) Redfern, L. R.; Farha, O. K. Mechanical properties of metal-organic frameworks. *Chem. Sci.* **2019**, *10*, 10666–10679.

(6) Bennett, T. D.; Horike, S. Liquid, glass and amorphous solid states of coordination polymers and metal-organic frameworks. *Nature Rev. Mater.* **2018**, *3*, 431–440.

(7) Bennett, T. D.; Cheetham, A. K. Amorphous Metal–Organic Frameworks. *Acc. Chem. Res.* **2014**, *47*, 1555–1562.

(8) Nandasiri, M. I.; Jambovane, S. R.; McGrail, B. P.; Schaef, H. T.; Nune, S. K. Adsorption, separation, and catalytic properties of densified metal-organic frameworks. *Coord. Chem. Rev.* **2016**, *311*, 38–52.

(9) Fonseca, J.; Gong, T.; Jiao, L.; Jiang, H.-L. Metal–organic frameworks (MOFs) beyond crystallinity: amorphous MOFs, MOF liquids and MOF glasses. *J. Mater. Chem. A* **2021**, *9*, 10562–10611.

(10) Katrusiak, A. High-pressure crystallography. *Acta Cryst.* **A** **2008**, *64*, 135–148.

(11) Varshneya, A. K.; Mauro, J. C. In *Fundamentals of Inorganic Glasses (Third Edition)*, third edition ed.; Varshneya, A. K., Mauro, J. C., Eds.; Elsevier, 2019; pp 187–214.
(12) Coudert, F.-X.; Fuchs, A. H. Computational characterization and prediction of metal-organic framework properties. *Coord. Chem. Rev.* **2016**, *307*, 211–236.

(13) de Jong, M.; Chen, W.; Angsten, T.; Jain, A.; Notestine, R.; Gamst, A.; Sluiter, M.; Krishna Ande, C.; van der Zwaag, S.; Plata, J. J. et al. Charting the complete elastic properties of inorganic crystalline compounds. *Sci Data* **2015**, *2*, 345.

(14) Chibani, S.; Coudert, F.-X. Systematic exploration of the mechanical properties of 13
621 inorganic compounds. *Chem. Sci.* **2019**, *10*, 8589–8599.

(15) Castel, N.; Coudert, F.-X. Atomistic Models of Amorphous Metal–Organic Frameworks. *J. Phys. Chem. C* **2022**, *126*, 6905–6914.

(16) Matsui, M. Molecular dynamics simulation of structures, bulk moduli, and volume thermal expansivities of silicate liquids in the system CaO-MgO-Al₂O₃-SiO₂. *Geophys. Res. Lett.* **1996**, *23*, 395–398.

(17) Ono, S.; Brodholt, J. P.; David Price, G. Elastic, thermal and structural properties of platinum. *J. Phys. Chem. Sol.* **2011**, *72*, 169–175.

(18) Andersson, D.; Beeler, B. Ab initio molecular dynamics (AIMD) simulations of NaCl, UCl₃ and NaCl-UCl₃ molten salts. *J. Nucl. Mater.* **2022**, *568*, 153836.

(19) Parrinello, M.; Rahman, A. Strain fluctuations and elastic constants. *J. Chem. Phys.* **1982**, *76*, 2662–2666.

(20) Bouëssel du Bourg, L.; Ortiz, A. U.; Boutin, A.; Coudert, F.-X. Thermal and mechanical stability of zeolitic imidazolate frameworks polymorphs. *APL Mater.* **2014**, *2*, 124110.

(21) Perger, W.; Criswell, J.; Civalleri, B.; Dovesi, R. Ab-initio calculation of elastic constants of crystalline systems with the CRYSTAL code. *Comput. Phys. Commun.* **2009**, *180*, 1753–1759.
(22) Wang, Y.; Wang, J. J.; Zhang, H.; Manga, V. R.; Shang, S. L.; Chen, L.-Q.; Liu, Z.-K. A first-principles approach to finite temperature elastic constants. *J. Phys.: Condens. Matter* **2010**, *22*, 225404.

(23) Widmer, R. N.; Lampronti, G. I.; Anzellini, S.; Gaillac, R.; Farsang, S.; Zhou, C.; Belenguer, A. M.; Wilson, C. W.; Palmer, H.; Kleppe, A. K. et al. Pressure promoted low-temperature melting of metal-organic frameworks. *Nature Mater.* **2019**, *18*, 370–376.

(24) Henke, S.; Wharmby, M. T.; Kieslich, G.; Hante, I.; Schneemann, A.; Wu, Y.; Daisenberger, D.; Cheetham, A. K. Pore closure in zeolitic imidazolate frameworks under mechanical pressure. *Chem. Sci.* **2018**, *9*, 1654–1660.

(25) Gaillac, R.; Pullumbi, P.; Beyer, K. A.; Chapman, K. W.; Keen, D. A.; Bennett, T. D.; Coudert, F.-X. Liquid metal-organic frameworks. *Nature Mater.* **2017**, *16*, 1149–1154.

(26) Bennett, T. D.; Goodwin, A. L.; Dove, M. T.; Keen, D. A.; Tucker, M. G.; Barney, E. R.; Soper, A. K.; Bithell, E. G.; Tan, J.-C.; Cheetham, A. K. Structure and Properties of an Amorphous Metal-Organic Framework. *Phys. Rev. Lett.* **2010**, *104*, 2272.

(27) Vervoorts, P.; Stebani, J.; Méndez, A. S. J.; Kieslich, G. Structural Chemistry of Metal–Organic Frameworks under Hydrostatic Pressures. *ACS Mater. Lett.* **2021**, *3*, 1635–1651.

(28) Vervoorts, P.; Hobday, C. L.; Ehrenreich, M. G.; Daisenberger, D.; Kieslich, G. The Zeolitic Imidazolate Framework ZIF-4 under Low Hydrostatic Pressures. *Z. Anorg. Allg. Chem.* **2019**, *645*, 970–974.

(29) Collings, I. E.; Goodwin, A. L. Metal–organic frameworks under pressure. *J. Appl. Phys.* **2019**, *126*, 181101.
(30) Tan, J.-C.; Civalleri, B.; Erba, A.; Albanese, E. Quantum mechanical predictions to elucidate the anisotropic elastic properties of zeolitic imidazolate frameworks: ZIF-4 vs. ZIF-zni. *CrystEngComm* **2015**, *17*, 375–382.

(31) Ryder, M. R.; Tan, J.-C. Explaining the mechanical mechanisms of zeolitic metal-organic frameworks: revealing auxeticity and anomalous elasticity. *Dalton Trans.* **2016**, *45*, 4154–4161.

(32) Shi, Z.; Weng, K.; Li, N. The Atomic Structure and Mechanical Properties of ZIF-4 under High Pressure: Ab Initio Calculations. *Molecules* **2023**, *28*, 22.

(33) Bennett, T. D.; Simoncic, P.; Moggach, S. A.; Gozzo, F.; Macchi, P.; Keen, D. A.; Tan, J.-C.; Cheetham, A. K. Reversible pressure-induced amorphization of a zeolitic imidazolate framework (ZIF-4). *Chem. Commun.* **2011**, *47*, 7983.

(34) Zhang, L.; Hu, Z.; Jiang, J. Sorption-Induced Structural Transition of Zeolitic Imidazolate Framework-8: A Hybrid Molecular Simulation Study. *J. Am. Chem. Soc.* **2013**, *135*, 3722–3728.

(35) Bennett, T. D.; Yue, Y.; Li, P.; Qiao, A.; Tao, H.; Greaves, N. G.; Richards, T.; Lampronti, G. I.; Redfern, S. A. T.; Blanc, F. et al. Melt-Quenched Glasses of Metal–Organic Frameworks. *J. Am. Chem. Soc.* **2016**, *138*, 3484–3492.

(36) To, T.; Sørensen, S. S.; Yue, Y.; Smedskjaer, M. M. Bond switching is responsible for nanoductility in zeolitic imidazolate framework glasses. *Dalton Trans.* **2021**, *50*, 6126–6132.

(37) Horike, S.; Ma, N.; Fan, Z.; Kosasang, S.; Smedskjaer, M. M. Mechanics, Ionics, and Optics of Metal–Organic Framework and Coordination Polymer Glasses. *Nano Lett.* **2021**, *21*, 6382–6390.
(38) Adhikari, P.; Li, N.; Rulis, P.; Ching, W.-Y. Deformation behavior of an amorphous zeolitic imidazolate framework – from a supersoft material to a complex organometallic alloy. *Phys. Chem. Chem. Phys.* **2018**, *20*, 29001–29011.

(39) Shi, Z.; Hartati, S.; Arramel, A.; Li, N. Unraveling the bond structure, porosity, and mechanical properties amorphous ZIF-4 and its topological equivalents: Large scale ab initio calculations. *APL Mater.* **2023**, *11*, 021103.

(40) Castel, N.; Coudert, F.-X. Challenges in Molecular Dynamics of Amorphous ZIFs Using Reactive Force Fields. *J. Phys. Chem. C* **2022**, *126*, 19532–19541.

(41) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. *Comput. Phys. Commun.* **2005**, *167*, 103–128.

(42) Kühne, T. D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V. V.; Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R. Z.; Schütt, O.; Schiffmann, F. et al. CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations. *J. Chem. Phys.* **2020**, *152*, 194103.

(43) Gaillac, R.; Pullumbi, P.; Bennett, T. D.; Coudert, F.-X. Structure of Metal–Organic Framework Glasses by *Ab Initio* Molecular Dynamics. *Chem. Mater.* **2020**, *32*, 8004–8011.

(44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. *Phys. Rev. Lett.* **1996**, *77*, 3865–3868.

(45) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate *ab initio* parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. *J. Chem. Phys.* **2010**, *132*, 154104.
(46) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials.
*Phys. Rev. B* **1996**, *54*, 1703–1710.

(47) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling.
*J. Chem. Phys.* **2007**, *126*, 014101.

(48) Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M. et al. The ReaxFF reactive force-field: development, applications and future directions. *npj Comput. Mater.* **2016**, *2*, 15011.

(49) Yang, Y.; Shin, Y. K.; Li, S.; Bennett, T. D.; van Duin, A. C. T.; Mauro, J. C. Enabling Computational Design of ZIFs Using ReaxFF. *J. Phys. Chem. B* **2018**, *122*, 9616–9624.

(50) Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in 't Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D. et al. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. *Comput. Phys. Commun.* **2022**, *271*, 108171.

(51) Aktulga, H.; Fogarty, J.; Pandit, S.; Grama, A. Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques. *Parallel Computing* **2012**, *38*, 245–259.

(52) Bureekaew, S.; Amirjalayer, S.; Tafipolsky, M.; Spickermann, C.; Roy, T. K.; Schmid, R. MOF-FF - A flexible first-principles derived force field for metal-organic frameworks.
*Phys. Status Solidi B* **2013**, *250*, 1128–1141.

(53) Dürholt, J. P.; Fraux, G.; Coudert, F.-X.; Schmid, R. Ab Initio Derived Force Fields for Zeolitic Imidazolate Frameworks: MOF-FF for ZIFs. *J. Chem. Theo. Comput.* **2019**, *15*, 2420–2432.

(54) Schmid, R.; Schmitz, G.; Dürholt, J. P.; Keupp, J.; Amabile, R. cmc-tools. https://github.com/MOFplus/cmc-tools, (accessed 2023-04-28).
(55) Angel, R. J. Equations of State. *Rev. Mineral. Geochem.* **2000**, *41*, 35-59.

(56) Birch, F. Finite Elastic Strain of Cubic Crystals. *Phys. Rev.* **1947**, *71*, 809-824.

(57) Ray, J. R. Elastic constants and statistical ensembles in molecular dynamics. *Comput. Phys. Rep.* **1988**, *8*, 109-151.

(58) Gaillac, R.; Pullumbi, P.; Coudert, F.-X. ELATE: an open-source online application for analysis and visualization of elastic tensors. *J. Phys.: Condens. Matter* **2016**, *28*, 275201.

(59) Hill, R. The Elastic Behaviour of a Crystalline Aggregate. *Proc. Phys. Soc. A* **1952**, *65*, 349-354.

(60) Chandler, D. *Introduction to Modern Statistical Mechanics*; Oxford University Press, 1987.

(61) Gaillac, R.; Pullumbi, P.; Coudert, F.-X. Melting of Zeolitic Imidazolate Frameworks with Different Topologies: Insight from First-Principles Molecular Dynamics. *J. Phys. Chem. C* **2018**, *122*, 6730-6736.

(62) Willems, T. F.; Rycroft, C. H.; Kazi, M.; Meza, J. C.; Haranczyk, M. Algorithms and tools for high-throughput geometry-based analysis of crystalline porous materials. *Micro. Meso. Mater.* **2012**, *149*, 134-141.

(63) Pinheiro, M.; Martin, R. L.; Rycroft, C. H.; Jones, A.; Iglesia, E.; Haranczyk, M. Characterization and comparison of pore landscapes in crystalline porous materials. *J. Molec. Graph. Model.* **2013**, *44*, 208-219.

(64) Pinheiro, M.; Martin, R. L.; Rycroft, C. H.; Haranczyk, M. High accuracy geometric analysis of crystalline porous materials. *CrystEngComm* **2013**, *15*, 7531.
(65) Nam, H.; Bengtson, A.; Vörtler, K.; Saha, S.; Sakidja, R.; Morgan, D. First-principles molecular dynamics modeling of the molten fluoride salt with Cr solute. *J. Nucl. Mater.* **2014**, *449*, 148–157.

(66) Ortiz, A.; Sánchez-González, J.; González-Méndez, L.; Cumbrera, F. Determination of the thermal stability and isothermal bulk modulus of the ZrO₂ polymorphs at room temperature by molecular dynamics with a semi-empirical quantum-chemical model. *Ceram. Int.* **2007**, *33*, 705–709.

(67) Padilla Espinosa, I. M.; Barua, N.; Mohan, R. V. Hydrostatic compression and pressure phase transition of major Portland cement constituents – Insights via molecular dynamics modeling. *Cement* **2022**, *7*, 100017.

(68) Haigis, V.; Belkhodja, Y.; Coudert, F.-X.; Vuilleumier, R.; Boutin, A. Challenges in first-principles NPT molecular dynamics of soft porous crystals: A case study on MIL-53(Ga). *J. Chem. Phys.* **2014**, *141*, 064703.

(69) Schröder, C. A.; Baburin, I. A.; van Wüllen, L.; Wiebcke, M.; Leoni, S. Subtle polymorphism of zinc imidazolate frameworks: temperature-dependent ground states in the energy landscape revealed by experiment and theory. *CrystEngComm* **2013**, *15*, 4036–4040.

(70) Widmer, R. N.; Lampronti, G. I.; Chibani, S.; Wilson, C. W.; Anzellini, S.; Farsang, S.; Kleppe, A. K.; Casati, N. P. M.; MacLeod, S. G.; Redfern, S. A. T. et al. Rich Polymorphism of a Metal–Organic Framework in Pressure–Temperature Space. *J. Am. Chem. Soc.* **2019**, *141*, 9330–9337.

(71) Katsenis, A. D.; Puškarić, A.; Štrukil, V.; Mottillo, C.; Julien, P. A.; Užarević, K.; Pham, M.-H.; Do, T.-O.; Kimber, S. A. J.; Lazić, P. et al. In situ X-ray diffraction monitoring of a mechanochemical reaction reveals a unique topology metal-organic framework. *Nat Commun* **2015**, *6*, 413.
(72) Bennett, T. D.; Keen, D. A.; Tan, J.-C.; Barney, E. R.; Goodwin, A. L.; Cheetham, A. K. Thermal Amorphization of Zeolitic Imidazolate Frameworks. *Angew. Chem. Int. Ed.* **2011**, *50*, 3067–3071.

(73) Frentzel-Beyme, L.; Kolodzeiski, P.; Weiß, J.-B.; Schneemann, A.; Henke, S. Quantification of gas-accessible microporosity in metal-organic framework glasses. *Nature Commun.* **2022**, *13*, 705.

(74) Spencer, E. C.; Angel, R. J.; Ross, N. L.; Hanson, B. E.; Howard, J. A. K. Pressure-Induced Cooperative Bond Rearrangement in a Zinc Imidazolate Framework: A High-Pressure Single-Crystal X-Ray Diffraction Study. *J. Am. Chem. Soc.* **2009**, *131*, 4022–4026.

(75) To, T.; Sørensen, S. S.; Stepniewska, M.; Qiao, A.; Jensen, L. R.; Bauchy, M.; Yue, Y.; Smedskjaer, M. M. Fracture toughness of a metal–organic framework glass. *Nature Commun.* **2020**, *11*, 357.

(76) Du, T.; Sørensen, S. S.; Zhou, Q.; Bauchy, M.; Smedskjaer, M. M. Accessing a Forbidden Disordered State of a Zeolitic Imidazolate Framework with Higher Stiffness and Toughness through Irradiation. *Chem. Mater.* **2022**, *34*, 8749–8759.

(77) McGreevy, R. L. Reverse Monte Carlo modelling. *J. Phys. Cond. Matter* **2001**, *13*, R877–R913.

(78) Weng, T.; Schmidt, J. R. Flexible and Transferable ab Initio Force Field for Zeolitic Imidazolate Frameworks: ZIF-FF. *J. Phys. Chem. A* **2019**, *123*, 3000–3012.

(79) Acuna-Yeomans, E.; Gutierrez-Sevillano, J.; Calero, S.; Dubbeldam, D. Evaluation of ZIF-8 flexible force fields for structural and mechanical properties. *Micro. Meso. Mater.* **2023**, *348*, 112406.
(80) Behler, J. First Principles Neural Network Potentials for Reactive Simulations of Large Molecular and Condensed Systems. *Angew. Chem. Int. Ed.* **2017**, *56*, 12828–12840.

(81) Eckhoff, M.; Behler, J. From Molecular Fragments to the Bulk: Development of a Neural Network Potential for MOF-5. *J. Chem. Theory Comput.* **2019**, *15*, 3793–3809.
TOC Graphic

<figure><img src="image_10.png" /></figure>