Evaluation of ZIF-8 flexible force fields for structural and mechanical properties

Metal–Organic Frameworks (MOFs) offer considerable potential for applications in adsorption due to their large pore volumes and surface areas. Studies on mechanical stability of MOFs are scarce. Seminal experimental work has shed a new light on the role that elastic constants play in establishing the structural stability of the prototypical ZIF-8 MOF

In a previously published computational study, Ortiz et al. [26] monitored the mechanical properties of ZIF-8 as a function of hydrostatic pressure and found that the elastic constant corresponding to the shear modulus ( 44 ) rapidly drops as the pressure approaches ∼0. 4 GPa, resulting in the amorphization of the framework.The weak resistance against shear deformation (shear-mode softening) was later experimentally verified during ball-milling [27,28].Interestingly, ZIF-8 is one of the few MOFs for which the elastic constants have been experimentally measured, via Brillouin scattering [29].It was found that the shear modulus of the crystal is remarkably low under atmospheric https://doi.org/10.1016/j.micromeso.2022  conditions (≤1 GPa), which is the lowest yet reported for a singlecrystalline extended solid.Moreover, its elastic deformation mechanism was linked to the pliant ZnN 4 tetrahedra.
To properly exploit ZIF-8 for its proposed technical applications, it is essential to gain an atomistic understanding of the underlying mechanisms that control its flexible behaviour.For this, molecular simulations are a valuable tool.For computing structural, gas adsorption and gas diffusion properties, classical flexible force fields (FFs) have been developed to provide sufficiently good predictions [30].However, the accurate description of the interactions between the organic and inorganic parts of the framework remains a challenge.Over the past decade a number of flexible force fields have been proposed for ZIF-8, the earliest of which [31][32][33][34][35] make use the functional baseline provided by the generic GAFF [36] force field and combine UFF [37] and AMBER [36] parameters, modified in an ad-hoc manner, to reproduce experimental observables.In terms of results, these FFs mostly focus on the prediction of gas adsorption and diffusion properties and were initially validated by comparing structural properties, such as the lattice parameter of the unit cell at atmospheric conditions, with crystallographic data.In recent years, more sophisticated and complex ab-initio derived force fields have been proposed.In 2019, R. Schmid and colleagues [38] extended the MOF-FF methodology [39], which focuses on accuracy instead of transferability, and proposed a force field capable of reproducing the subtle swing-effect behaviour, as well as certain structural and dynamical properties of ZIF-8.In contrast, the one proposed by Weng and Schmidt [40] trades accuracy for transferability and is able to describe several of the ZIF-8 polymorphs.Systematic comparisons and evaluations of these force fields for the prediction of mechanical properties are lacking.
In a recent study, Maul et al. [41] investigated the mechanical response of ZIF-8 using quantum-mechanical calculations over a range of pressures and found that at  > 0.2 GPa an anisotropic response along ⟨111⟩ and ⟨100⟩ is observed, characterized by nonlinear behaviour of the  11 and  12 elastic constants.These results were contrasted to classical simulation data obtained by Ortiz et al. [26] using the Zhang [33] force field, which was unable to describe said behaviour.In general, we found that systematic comparisons and evaluations of classical force fields for the prediction of mechanical properties are lacking.
In this work, five published flexible force fields for ZIF-8 were implemented and validated in both the RASPA [42] and LAMMPS [43] simulation packages.Molecular dynamics simulations in the isothermalisobaric ensemble were conducted to determine the temperature and pressure dependence of the framework lattice parameters, and the 0 K elastic constants at different values of pressure were calculated using the eigenmode-following technique [44,45], which are compared with existing experimental measurements and quantum-mechanical calculations.The obtained results provide insight into the relationship between structural and elastic properties and the chosen parametrization, and suggest further adoption of ''top-down'' philosophies to the construction of future classical force fields would be beneficial to their accuracy in this respect.Our two-code approach led us to identify that the calculation of 0 K elastic constants is highly dependent on finding true (global) energy minimum configurations, and therefore the chosen minimization procedure is essential.We found that eigenmode-following fared best in this respect when compared to methodologies that rely on other minimization procedures, such conjugate gradient and steepest descent.

Methodology
The term ''force field'' refers to the functional form and parameter sets of the energy equation.A force field consists of two parts: 1.An analytical expression for the inter-atomic potential energy  () as a function of the atomic coordinates . 2. Definitions and parameters for ''atom types'', classified based on bonding and environment.Parameters are assigned based on the atom types involved.
The Generalized AMBER Force Field (GAFF) energy equation is defined as [46]: The three main categories of classical simulation methodologies are Molecular Dynamics (MD), and Monte Carlo (MC) and Molecular Mechanics (MM).In a Molecular Dynamics (MD) simulation, all atoms move according to Newton's equations of motion and the forces computed on the atoms.The forces are given by the gradients of the potential energy with respect to the internal coordinates of the molecule.The potential energy surface  of a (periodic) system can be Taylor-expanded around a configuration  of the system where  = is usually referred to as the Hessian matrix.The superscript  denotes the transpose of a vector or matrix.The generalized Hessian, having dimensions (3 + 6) × (3 + 6) with  being the number of atoms, is given as with the force constant matrix   (3 × 3) and the Born term   (6 × 6) being the second order derivative of the internal energy with respect to position and strain , respectively.The strain (and also stress tensor) is symmetric and can be simplified to a 6-dimensional vector using Voigt notation: The Born term accounts for distortions of the lattice,   and   are cross-terms.The 0 K elastic tensor is defined as the derivative of the stress with respect to the strain [47,48] at zero gradient  = 0 and can be expressed in terms of the generalized Hessian [49] The relaxation term arises when more than one particle is present in the unit cell [49].When the system is strained the atoms need to relax relative to one another, because before and after the strain applied the system must be in a state of zero net force.The state at which elastic constant are computed must a true (global) energy minimum: all forces (also on the cell) must be zero and all eigenvalues of the generalized Hessian matrix must be positive.The systematic study of the lattice stability was done by Born and Huang [50], who formulated the stability criteria in terms of the elastic constants   , by expanding the internal crystal energy in a power series in the strain and imposing the convexity of the energy.This criterion expresses the fact that any mechanical strain must increase the total mechanical energy of a system at equilibrium (resulting from the requirement that the eigenvalues of the elastic stiffness matrix  must be positively defined).The eigenvectors are the deformation modes.The Born elastic stability criteria restrict what values various moduli may have if the crystal is to be stable.In the case of external pressure the relevant free energy is  and the relevant moduli are   .In the cubic crystal system the stability criteria are (Voigt notation) [ The flexible force fields considered in this work are the following: • FF5 : T. Weng and J. R. Schmidt (2019) [40] They were implemented in both the LAMMPS [43] and RASPA [42] classical simulation packages and subsequently validated.This was done via comparison of the average lattice parameters at atmospheric conditions, which are sensitive to the parameters and the details of the implementations.The sets of functional forms, corresponding parameters and additional implementation details for all force fields are included in the supporting information document.Other existing force fields such as the ones proposed by Krokidas et al. [35] and Schmid et al. [38] were left out of the scope of this study, the former because of multiple ambiguities in the definitions of the potential energy terms in the original paper and the latter because of the inability to properly implement the specified non-bonded potentials in RASPA.
The average lattice parameters of the ZIF-8 crystal structure at different fixed values of temperature and pressure were calculated using Molecular Dynamics (MD) in the isothermal-isobaric ensemble (NPT) with LAMMPS.The simulations were performed in a supercell of ZIF-8 consisting of 2 × 2 × 2 unit cells.Periodic boundary conditions were applied in all three dimensions.Long-range electrostatic interactions were treated via the Ewald method.The simulation timestep used for the integrator was 1 fs.For each framework atom, a ''scaled 1-4'' policy was taken into account; that is, both non-bonded interactions (VdW and electrostatic) between couples of bonded atoms (1)(2) or between atoms bonded to a common atom (1-3) were excluded, while the interaction between atoms separated by two others (1)(2)(3)(4) were scaled according to the information provided by each force field.
In order to employ the two-code approach for the analysis and comparison of structural and mechanical properties, it is necessary to first verify that the FF implementations provide congruent results between the codes.For this, the term-by term contributions to the energy were calculated with both codes, for the same atomic configuration.The results for LAMMPS and RASPA for a chosen initial configuration per energy contribution were calculated and are included in the supporting information document (Table S9).For all implemented force fields the energy terms of each contribution match up to at least the fourth decimal in the case of bonded and VdW interactions and up to the second Fig. 3. Lattice parameter of the unit cell predicted by the implemented force fields.MD data obtained at a temperature of 300 K. Available experimental data at  < 0.4 GPa (before framework amorphization) represented by crosses [24,56].decimal in the case of electrostatic interactions.The differences in the latter might be due to either different implementations of the Ewald summation in the codes or the finite precision of the specified atom positions.We found that often the literature force field descriptions are defined unclear and/or ambiguously.Significant reverse-engineering and trial and error was required to resolve all ambiguities.

Temperature dependence
To examine the effect of temperature on the size and shape of the ZIF-8 unit cell, a series of MD simulations were performed at a pressure of 1 atm between the ranges of 100-900 K for all implemented force fields.As can be seen in Fig. 2 for FF1, and in Figure S3 in the supporting information for all evaluated force fields, it was found that the lattice parameter depends only weakly on temperature, confirming their high thermal stability as experimentally observed.Also in agreement with experiments, is that the unit cell becomes larger with increasing temperature (positive thermal expansion) [53].FF1 qualitatively matches the data of Zhou et al. [53] quite well, albeit lower in value.As can be seen in Figure S3, all other evaluated force fields do better in this respect, providing closer estimates.
In Fig. 2 we also compared the Conjugate Gradient [54] (CG) method, with a combination of Steepest Descent [55] (SD) and CG methods, and with the more advanced (but also more expensive) eigenmode following [44,45] method.A crucial requirement is that the MD results for low temperatures must converge to the 0 K result obtained with energy minimization.Eigenmode following minimization consistently displays this behaviour (Figure S3), while the SD and CG methods gave erroneous results.We confirmed that these structures do have zero forces, but do not pass the requirement of having all positive eigenvalues of the generalized Hessian matrix.That is, the CG and SD/CG methods were not able to find true minimum energy structures.This has important implications for the reliability and accuracy of elastic constants.

Unit cell size as a function of pressure
Next, we examine the effect of pressure on the ZIF-8 unit cell lattice parameter.A series of MD simulations were performed at a temperature of 300 K (Fig. 3 and Figure S4 in the supporting information) over a wide range of pressures using all implemented force fields.It was found that force fields 2 through 5 were able to accurately reproduce the experimental geometric properties of the ZIF-8 unit cell at low values of pressure ( < 0.4 GPa).These force fields qualitatively reproduce the experimental amorphization of the system that occurs at higher values of pressure.FF1 both underestimates the lattice parameter at lower values of pressure and fails to reproduce the pressure-induced phase change.

Elastic constants as a function of pressure
Evaluating the elastic constant dependency on pressure allows us to shed more light on the amorphization mechanism.Fig. 4 shows the computed elastic constants as a function of pressure at 0 K for all implemented force fields.When compared to each other, the force fields provide not only different quantitative results, but more importantly different qualitative results as well.
At zero pressure, the agreement between M06-2X DFT and experimental values from ambient temperature and pressure measurements is good, except for a slight overestimation of  12 .Force fields 3 and 4 compare reasonably well with experimental and DFT values, while FF1, FF2 and FF5 differ significantly in one or more elastic constant value.We can observe that while all the studied force fields present a negative linear dependency of  44 throughout the pressure range; with FF5 providing exceptional quantitative estimations, none of them are able to properly describe the non-linear behaviour of  11 and  12 at  > 0.2 predicted by DFT calculations.
FF3 and FF4, while underestimating the amorphization pressure, provide qualitatively similar results and come closest to predict the EC trends described by the DFT results.FF1 provides a consistent overestimation  44 over the pressure range, leading to a predicted amorphization pressure that is too high (i.e.FF1 is stable up to and over 1 GPa), which is congruent with results shown in Fig. 3.These three force fields correctly suggest that the negative pressure dependence of  44 is the clear direct route to instability, that is, via the breaking of the third Born criteria (Eq.( 14)) or shear mode softening, in agreement with the experimental finding.In contrast, both FF2 and FF5 qualitatively show another amorphization mechanism, namely Eq. ( 13), where  11 becomes larger than  12 + 2, which for the former FF occurs at pressures lower than 0.1 GPa.Owing to the sharp increase of the  12 elastic constant at  > 0.2, FF5 presents ambiguity regarding which stability criteria is violated first and therefore the amorphization mechanism related to the structural transformation of the framework.
In Table 1 we compare the computed 0 K elastic constants and elastic moduli obtained with RASPA to the reported simulated and experimental values.Quantitatively, at zero pressure, FF2 and FF3 Fig. 4. 0 K elastic constants obtained using the mode-following technique over a range of pressures.Black data points indicate the C11 elastic constants, red the C12 elastic constants and blue the C44 elastic constants.The simulated values are compared to the Brillouin scattering measurements (filled symbols with green borders) [29] conducted at 295 K as well as B3LYP DFT (crosses) [29] and M06-2X DFT (open symbols) [41] results.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)provide similar 0 K results and come the closest to the reported DFT and experimental values.FF1 overestimates both  11 and  12 while FF2 greatly underestimates  11 .FF5 provides an exceptional estimation of  44 and significantly underestimates  11 and  12 .Of note is the fact that the reported 0 K values for FF5 are in heavy disagreement with the ones obtained in this work.This discrepancy might be related to the difference in methodologies employed to calculate them.

Elastic constants using SD/CG vs. mode-following minimization
To determine how the calculated elastic constants depend on the energy minimization algorithm, the SD/CG methods used in LAMMPS and eigenmode-following as implemented in RASPA were compared.For each force field, the output from the 0 K conjugate gradient LAMMPS routine was used as the input for RASPA's minimization.It Fig. 5. Eigenmode-following minimization of ZIF-8 using FF1, starting from a configuration previously minimized using the conjugate gradient algorithm.Initially the structure presents 170 negative eigenvalues.Notice that the number of negative eigenvalues is reduced to zero, then increases (between steps 38 and 51) before decreasing again.The increase corresponds to a structural change in the framework which the SD and CG methods are unable to induce.
was found that the former is not able to find true minimum energy configurations (Fig. 2) due to the identification of negative eigenvalues (corresponding to first-order saddle-points).Furthermore we noted that in most cases, to find relaxed configurations resulting in zero negative eigenvalues, the CG energy minimization routine requires either a preceding NPT MD simulation at very-low temperatures or the implementation of additional iterative minimization runs of different types (such as steepest descent).This last methodology was found not to be self-consistent across force fields.The problem is that these methods have no inherent mechanisms to systematically remove negative eigenvalues.Fig. 5 shows how the eigenmode-following technique quickly drives the system to zero eigenvalues.Note that the number can increase again during the minimization which often signals large-scale structural changes that need to be overcome.Fig. 6 shows the final energies obtained from the CG-based methods compared to the iteration steps of the eigenmode following.What becomes clear is that methods based on energies or forces are not able to find (a) the lowest energy, (b) the correct structure at zero K.We note that, even though the difference in lowest energy is small, the obtained unit cells differ in size and in atomic structure.When comparing the resulting energy values and unit-cell dimensions obtained using the CG-based routines with the ones obtained with eigenmode-following, it becomes clear that former are unable to efficiently relax both the framework atom's positions and the simulation box at the same time, which might produce problems when used to minimize the system after the box deformations are applied (in the elastic constants calculation routine).
In order to compute elastic constants with LAMMPS we employed modified versions of the manual-suggested routine.These rely on userdefined deformation parameters and make use of numerical differentiation to obtain derivatives.In contrast, in RASPA, the elastic constants are computed via Eq.( 5) using analytical derivatives of the functional forms of the potentials.Since RASPA gives zero Kelvin elastic constants basically at machine precision, we can test how sensitive the numerical procedure is to the deformation parameters.The elastic constants were re-computed using the LAMMPS CG-based routine using different box deformation parameters, with magnitudes ranging from 5e−5 GPa (corresponding to <0.01%change in simulation box side dimensions) to 0.1 GPa (about 10%).We found that the predicted elastic constant values severely depend on the chosen parameter (Fig. 7).Furthermore, no systematic dependency of the elastic constants estimate with respect to deformation parameter was identified in Fig. 7(b).

Conclusions
We evaluated a set of literature classical force fields to test whether they can accurately reproduce elastic constants and thermal expansion for ZIF-8.The force fields not only gave different quantitative results, but more importantly gave different qualitative results, both when compared to each other and, in the case of elastic constants, to reported DFT-calculated values.The reason is that the majority of these force fields have been developed ''bottom-up'', starting from atomic bonding parameters, parameters for bending and torsion, and next long-range interactions like Van der Waals and electrostatics.It is therefore not surprising that any reproduction of system properties is fortuitous.We find that this applies in particular when examining these properties as a function of pressure, i.e. the force fields differ significantly in the prediction of the pressure where amorphization or collapse starts to occur.Future work on force fields should therefore include this type of information in the parametrization.Although the force fields provide mostly acceptable quantitative estimations of elastic constants at very low pressures ( < 0.1 GPa), they are unable to reproduce the nonlinear behaviour of  11 and  12 at  > 0.2 predicted by DFT calculations and, in the case of FF2 and FF5, unequivocally point to the correct amorphization mechanism.In this study, we used both the RASPA and LAMMPS simulation packages and compared the implementation of the literature force fields for all energy terms separately.We found that often the literature force field descriptions are defined unclear and/or ambiguously in the supporting information.Significant reverseengineering and trial and error was required.RASPA and LAMMPS gave identical results for the unit cell volumes as a function of temperature and pressure.However, we found that zero Kelvin unit cells structures Fig. 6.Eigenmode-following minimization of ZIF-8 using FF1, starting from a configuration previously minimized using the conjugate gradient algorithm.The convergence of system energy is presented in red.The minimum energy obtained using CG and SG iterative routines are pictured with teal and blue lines, while the one corresponding to CG minimization after a low-temperature PT simulation is presented in purple.(a) total range (b) zoomed-in portion where the corresponding lattice parameter values are 16.7122Å (purple), 16.725 Å (blue) and 16.711 Å (red).The eigenmode-following technique relaxes both the cell shape and the atoms at the same time and produces the configuration with the lowest energy.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 7. (a) Box plot of 0 K elastic constant values obtained using the CG routine for FF1.The blue dots correspond to the elastic constant values obtained using different boxdeformation parameters; those outside of the whiskers are considered outliers.The red crosses represent the elastic constants obtained using the eigenmode-following technique, a self-consistent method that does not require user-defined parameters.(b) Bar plot of the  11 value obtained for each box-deformation parameter using the CG routine.The red line represents the  11 value obtained using eigenmode-following. and elastic constants depend on the procedure.The conjugate gradient minimization routine implemented in LAMMPS is often unable to find the true energy minimum and zero K structure and therefore, when employed in elastic constant calculations, provides erroneous values.The eigenmode-following minimization in RASPA is guaranteed to find the correct zero K structure, and therefore a comparison could be made.This methodology is important in MOFs, and we showed that unit cell obtained by RASPA is the same as extrapolated from finite temperature MD data to zero Kelvin, while the minimized unit cells from conjugate gradient are often erroneous.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.
.112406 Received 30 October 2022; Received in revised form 7 December 2022; Accepted 11 December 2022 Size of the unit cell of ZIF-8 predicted by the model proposed by Zheng et al.

Table 1
Comparison of Elastic constants (  ), Shear modulus (), Young's modulus () and Bulk modulus () found in the literature and those reported in this work.All values in GPa.