Deciphering the controlling factors for phase transitions in zeolitic imidazolate frameworks

ABSTRACT Zeolitic imidazolate frameworks (ZIFs) feature complex phase transitions, including polymorphism, melting, vitrification, and polyamorphism. Experimentally probing their structural evolution during transitions involving amorphous phases is a significant challenge, especially at the medium-range length scale. To overcome this challenge, here we first train a deep learning-based force field to identify the structural characteristics of both crystalline and non-crystalline ZIF phases. This allows us to reproduce the structural evolution trend during the melting of crystals and formation of ZIF glasses at various length scales with an accuracy comparable to that of ab initio molecular dynamics, yet at a much lower computational cost. Based on this approach, we propose a new structural descriptor, namely, the ring orientation index, to capture the propensity for crystallization of ZIF-4 (Zn(Im)2, Im = C3H3N2−) glasses, as well as for the formation of ZIF-zni (Zn(Im)2) out of the high-density amorphous phase. This crystal formation process is a result of the reorientation of imidazole rings by sacrificing the order of the structure around the zinc-centered tetrahedra. The outcomes of this work are useful for studying phase transitions in other metal-organic frameworks (MOFs) and may thus guide the development of MOF glasses.


Supplementary Text
Ab initio MD simulations.The data set for training the deep learning potential were first generated using ab initio MD simulations following the procedure employed in Ref. [1].The input files for ab initio MD simulations are adapted from the Supplementary file in Ref. [1].A series of systems containing C, H, N, and Zn elements were selected as the training sets, including both the crystalline and disordered structures of ZIF-4, SALEM-2, ZIF-8, ZIF-62, elementary substances, and compounds.Note that, two different compositions of ZIF-62 are selected, i.e., ZIF-62-0.5 and ZIF-62-0.25,corresponding to Zn(Im)1.5(bIm)0.5 and Zn(Im)1.75(bIm)0.25,respectively.An overview of the systems used for generating the data set is provided in Supplementary Table S1.
All the ab initio calculations were carried out at the DFT level [2] using the Quickstep module of the CP2K package [3] with the hybrid Gaussian and plane wave method (GPW) [4].As suggested in Ref. [5], the basis functions are mapped onto a multi-grid system with the default number of four different grids with a plane-wave cutoff for the electronic density to be 600 Ry, and a relative cutoff of 40 Ry to ensure the computational accuracy.The AIMD trajectories at elevated temperatures were obtained in the NVT ensemble with a timestep of 0.5 fs.The temperature was controlled using the Nosé-Hoover thermostat [6].The exchange-correlation energy was calculated using the Perdew-Burke-Ernzerhof (PBE) approximation [7], and the dispersion interactions were handled by utilizing the empirical dispersion correction (D3) from Grimme [8].The pseudo potential GTH-PBE combined with the corresponding basis sets DZVP-MOLOPT-SR-GTH were employed to describe the valence electrons [9].
Deep learning force field construction.The deep learning force field was constructed using the DeePMD-kit package.The training dataset consists of two parts: (1) the initial training data from the trajectory of AIMD simulations and the single energy calculation of the crystalline structures (i.e., ZIF-4, SALEM-2, ZIF-8, ZIF-62-0.5)under deformations; and (2) the expanded dataset realized by single energy calculation using the active machine learning method implemented in the DP-GEN package [10].The detailed constituents and configurations of training data are summarized in Supplementary Fig. S1 and S2 and Supplementary Tables S1 and S2.The test set is adopted from the trajectory of AIMD simulations of ZIF-62-0.25.We employed a commonly used network structure [11,12] for deep learning as follows: (1) the embedding network containing 3 layers with the number of neurons per layer equal to 25, 50, and 100, respectively; (2) the fitting network containing 3 layers with the number of neurons per layer all equal to 240.The descriptor "se_e2_a", which contains both angular and radial information of atomic configurations, was used as the model of force field.The training process included two processes, i.e., the initial training and final training.Both processes included 6,000,000 steps of iterations.The loss function for the training process included the energy, force, and virial terms.For the initial training, the learning rate dynamically changed from 1e-3 to 1e-9.The atomic interaction beyond 0.6 Å were treated using the ZBL repulsive interactions to avoid the overlap of atoms at high temperatures, and the prefactors of the energy, force, and virial terms dynamically changed from 0.2 to 1, 500 to 1, and 0.02 to 0.2, respectively.The final training process removed the ZBL interaction and the prefactors of all terms were set to 1 as the learning rate changed from 1e-5 to 1e-8.
We also tried the descriptor "se_e2_r", which only included the radial information of atomic configurations to train the model.However, the resulting force field cannot reproduce the glass structure as well as the descriptor "se_e2_a" (Supplementary Fig. S21a).We also tried different cutoffs for interatomic interaction (i.e., 5.0, 6.0, 7.0 Å), but the simulation result is not significantly influenced by the cutoffs (Supplementary Fig. S21b).Following the literature results, a cutoff for interatomic interaction was set to be 6 Å to achieve a balance between accuracy and computational efficiency.The comparisons between the results of energy, force, and stress obtained with the optimized DP model and the AIMD calculations are provided in Supplementary Figs.S22-S24.
More details of MD simulations.The initial configuration of the ZIF-4 crystal was obtained from the Cambridge Structural Database [13], which was further replicated into a 2×2×2 supercell of 2,176 atoms.This system size is commonly used in the MD simulations of ZIF glasses using ReaxFF when studying the structural features upon glass formation [14,15].To exclude the possibility of significant size effects in this work, we also attempted the use of a much larger system size with a 4×4×4 supercell of 17,408 atoms.However, the resulting glass structures did not exhibit a significant dependence on the system size (Supplementary Fig. S25) and we therefore use systems with 2,176 atoms in the remainder of this work.Other crystalline phases of Zn(Im)2 (i.e., SALEM-2, ZIF-zni, ZIF-6, and ZIF-10) are dynamically equilibrated under 300 K and zero pressure for 12.5 ps before structural analysis.The transferability of this DLFF is verified by simulating the newly discovered ZIF glasses with the cyano-functionalized imidazolate linkers [16].The structural properties and comparison with experiments confirm the accuracy and generalizability of the DLFF (Supplementary Fig. S26).
Note that the thermodynamic and structural properties of the structure at each state were analyzed based on the trajectory equilibrated at the specified temperature for 50 ps in the NVT ensemble with a dump interval of 0.25 ps to ensure proper statistical averaging.During the simulation, periodical boundary conditions (PBC) were applied in all directions.A timestep of 0.25 fs was used for the motion integration.The configurations obtained from the simulations were visualized using the OVITO package [17].
Local entropy of Zn atoms.Since the two-body excess entropy based on the pair correlation function accounts for about 90% of the configurational entropy, the local entropy of Zn atoms was calculated using the formula from Ref. [18] as follows: where  is the system density, kB is the Boltzmann constant, rm is the upper limit for integration, g(r) is the radial function centered at the ith atom, rij is the distance between ith and jth atoms, and σ is the broadening parameter.Here, rm was selected to be 5 Å, and σ to be 0.2.To increase the statistics of the analysis, the entropy calculations were based on the trajectory of equilibration.
Ring orientation index of the imidazole rings.The ring orientation index (ROI) of the imidazole rings was calculated based on the stereographic projection method [19], which enables the mapping of the three-dimensional vectors at the surface of a sphere onto a plane.Specifically, we used the following steps to calculate ROI as shown schematically shown in Fig. 4: (1) After identifying all the imidazole rings, the positions of the ring members were used to calculate the normal vector of each ring; (2) A unit sphere was constructed to include all the normal vectors; (3) For each projected plane (i.e., yz-, xy-, and xz-planes), the unit sphere was divided into two equal hemispheres, while all the normal vectors were transformed at the surface of one hemisphere since each normal had two opposite vectors; (4) Each normal vector was connected with the pole of the opposite hemisphere, and the intersection point at the projected plane was defined as the stereographic projection of the corresponding normal vector; and (5) The pole figure of each projected plane was obtained by calculating the local density of the points, while the distribution 6 histogram of each pole figure was calculated by binning the number of the points as a function of the central angle θ.Here, the bin size to calculate the histogram was selected as 2°.We note that the selection of bin size influences the value of ROI without changing the relative correlation (Supplementary Fig. S27).This analysis was also carried out on the trajectory of equilibration to increase statistics.The ROI value of each system was then finally calculated as, where h(x) is the distribution histogram.The overall ROI value was obtained by averaging over the three projected planes.We note that the ROI is found not to be influenced by the selected system size as the pole figure is irrespective of the simulated system size (see Supplementary Fig. S17).

Figure S2 .
Figure S2.Selected configurations of atomic structures containing C, H, N, and Zn elements for the ab initio simulations.Zn, C, H, and N are colored green, grey, white, and blue, respectively.The four crystalline structures on the top are also used as the initial structures for the ab initio simulations after perturbation and active learning simulations using DPGEN.

Figure S3 .
Figure S3.(a) Total and (b) Zn-N partial pair distribution functions of ZIF-4 glass from the melt quench simulation using the DLFF trained from different datasets.DLFF ver. 1 (the present datasets), and DLFF ver. 2 (extended datasets).Specifically, the extended datasets include the present datasets and additional datasets of the SALEM-2 crystal with methane molecules.

Figure S4 .
Figure S4.Comparison of experimental (X-ray) and MD-simulated differential correlation functions D(r) of ZIF-4 glass using ReaxFF and DLFF.The experimental data is from Ref.[20].

Figure S6 .
Figure S6.Zn-N partial pair distribution functions of ZIF-4 in crystalline, melt, and glassy states.Inset: Partial magnification to highlight the overlapped peaks in the molten state.

Figure S7 .
Figure S7.Zn-N partial pair distribution functions of ZIF-4 during the melting process of 400 K, 1000 K, and 1500 K at varying pressures (from 0 to 1.0 Gpa).Inset: Partial magnification to highlight the overlapped peaks.

Figure S8 .
Figure S8.Potentials of mean force along the Zn−N coordinate for the structure at 1500 K and varying pressures in the unit of kT.

Figure S9 .
Figure S9.Distribution of N-Zn-N bond angle of (a) ZIF-4 at 400 K and (b) ZIF-4 at 1000 K under varying pressures.

Figure S10 .
Figure S10.Distribution of entropy of Zn atoms in (a) ZIF-4 in crystalline, melt, and glassy states processed under zero pressure, and (b) different crystalline structures of Zn(Im)2.

Figure S11 .
Figure S11.Distribution of entropy of Zn atoms in ZIF-4 glasses processed under varying pressures.For comparison, the entropy distribution of ZIF-4 crystal equilibrated at 300 K is also provided as the green dashed line.

Figure S12 .
Figure S12.Comparison of structural properties of ZIF-zni crystal and ZIF-4 glass prepared under zero pressure.(a) Total pair distribution function.(b) Distribution of entropy of Zn atoms

Figure S15 .
Figure S15.Orientation distribution histograms of imidazole ring in SALEM-2 crystal projected in all three directions.

Figure S16 .
Figure S16.Orientation distribution histograms of imidazole ring in ZIF-4 at (b) 400 K and zero pressure and (c) 1500 K and zero pressure projected in all three directions.

Figure S18 .
Figure S18.The density evolution as a function of temperature during (a) melting and (b) quenching process under varying pressures.

Figure S20 .
Figure S20.(a) Total and (b) Zn-N partial pair distribution functions, and (c) N-Zn-N bond angle distribution function of ZIF-4 glass quenched under different cooling rates.

Figure S21 .
Figure S21.Pair distribution functions of the ZIF-4 glasses simulated using the DLFF with different models.(a) The "se_e2_a" and "se_e2_r" descriptors with a cutoff of 6 Å.(b) The "se_e2_a" descriptor with different cutoffs.

Figure S22 .
Figure S22.Comparison of the energy of different ZIF structures calculated by DFT and predicted by the present ML force field for both training sets of (a) SALEM-2, (b) ZIF-62-0.25,(c) ZIF-8, and (d) ZIF-4, and (e) test sets of ZIF-62-0.25.

Figure S23 .
Figure S23.Comparison of the atomic force components of the unit cell of ZIF-62 calculated by DFT and predicted by the present ML force field for both training and test sets at (a) x-direction, (b) y-direction, and (c) z-direction.

Figure S24 .
Figure S24.Comparison of the stress components of the unit cell of ZIF-62 calculated by DFT and predicted by the present ML force field for both training and test sets.(a) xx, (b) yy, (c) zz, (d) x€(e) xz, and (f) yz components.

Figure S27 .
Figure S27.Calculated (a) ring orientation distribution histogram and (b) ring orientation index as a function of the bin size of the ring orientation distribution histogram.

Table S2 .
Details of the explored data set obtained from the single energy calculation using CP2K.