Enabling Computational Design of ZIFs Using ... - ACS Publications

Sep 28, 2018 - Department of Materials Science and Metallurgy, University of Cambridge, 27 Charles Babbage Road, Cambridge CB3 0FS, U.K...
0 downloads 0 Views 2MB Size
Subscriber access provided by UNIV OF LOUISIANA

B: Fluid Interfaces, Colloids, Polymers, Soft Matter, Surfactants, and Glassy Materials

Enabling Computational Design of ZIFs Using ReaxFF Yongjian Yang, Yun Kyung Shin, Shichun Li, Thomas D. Bennett, Adri C.T. van Duin, and John C. Mauro J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b08094 • Publication Date (Web): 28 Sep 2018 Downloaded from http://pubs.acs.org on September 29, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Enabling Computational Design of ZIFs Using ReaxFF Yongjian Yang1, Yun Kyung Shin2, Shichun Li3,4, Thomas D. Bennett4, Adri C. T. van Duin1,2* and John C. Mauro1* 1 Department of Materials Science and Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, USA 2 Department of Mechanical and Nuclear Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, USA 3 Institute of Chemical Materials, China Academy of Engineering Physics, Mianyang 621900, P. R. China 4 Department of Materials Science and Metallurgy, University of Cambridge, 27 Charles Babbage Road, CB3 0FS, UK Corresponding Author: * [email protected] and [email protected] ABSTRACT: Classical force fields have been broadly used in studies of metal-organic framework crystals. However, processes involving bonds breaking or forming are prohibited due to the non-reactive nature of the potentials. With emerging trends in the study of ZIFs that include glass formation, defect engineering, and chemical stability, enhanced computational methods are needed for efficient computational screening of ZIF materials. Here, we present simulations of three ZIF compounds using a ReaxFF reactive force field. By simulating the melt-quench process of ZIF-4, ReaxFF can reproduce the atomic structure, density, thermal properties, and pore morphology of the glass formed (agZIF-4), showing remarkable agreement with experimental and first-principle molecular dynamics results. The predictive capability of ReaxFF is further exemplified in the melting of ZIF-62, where the balancing of electronic and steric effects of benzimidazolate yields a lower Tm. Based on electron withdrawing effect of –NO2 group, ReaxFF simulations predict that ZIF-77 has an even lower Tm in terms of Zn-N interaction but its low chemical stability makes it unsuitable as a glass former. Because of its low computational cost and transferability, ReaxFF will enable the computational design of ZIF materials by accounting for properties associated with disorder/defects.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 25

INTRODUCTION Metal−organic frameworks (MOFs) are reticular solids wherein inorganic nodes (e.g., metal ions: Co, Zn, Cu, etc., or cluster: Zn4O1, Al(OH)2O42, etc.) are linked via organic ligands in an infinite 3D array. Because the placement of these secondary building units (SBUs) within the framework is highly regular, the materials possess long-range order. Due to the nature of the chemical bonds between SBUs, MOFs can adopt a variety of crystalline structures. Among them, the zeolitic imidazolate frameworks (ZIFs) –– a subfamily of MOFs –– has recently attracted significant interest3–5. In ZIFs, tetrahedral metal ions, e.g., Zn and Co, are connected by imidazolate (Im) linkers with an average M-Im-M angle of 145°6, similar to the Si-O-Si angle preferred in zeolites; as such, these materials possess zeolitic structures. The industrial applications of ZIFs (and MOFs in general) include gas absorption and separation, catalysis, chemical sensing, and drug delivery. In the past decade, the main focus of research into metal-organic frameworks has been on the synthesis of new structures4,5,7. Recently, more effort is being devoted to understanding the chemical and physical flexibility of the structures, which7 are well documented in databases such as the Cambridge Structural Database (CSD)8 and the Computation-Ready, Experimental (CoRE) MOF database9. Two new promising trends emerging from these efforts are: (1) vitrification of the initially crystalline framework into a glass with a desired shape via a melt-quench process10,11, and (2) tailoring of MOF properties through nanoscale defect engineering7,12,13. One common issue shared by both processes is that they both involve reconfiguration of coordinative bonds in the material. Recently, the computational study of MOF structures, especially through molecular dynamics (MD) simulations, has become a valuable tool for both understanding of the framework and screening/optimization of the desired properties10,14–16. The set of force fields used for classical MD simulations of MOFs (including ZIFs) include UFF17, DREIDING18, BTW-FF19, UFF4MOF20,21, and QuickFF22,23. They have facilitated studies of various MOF features, including elasticity (e.g., Young’s modulus, bulk modulus, and Poisson’s ratio), phonon properties (e.g., Grüneisen parameters and thermal expansion)24,25, and sieving behavior26. However, due to the non-reactive nature of these potentials, no 2 ACS Paragon Plus Environment

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

bond breaking or reformation is allowed, thus precluding the study of MOFs at high temperature or under inelastic deformation, and also inhibiting the study of defect-containing or non-crystalline MOFs. While quantum mechanical techniques such as density functional theory (DFT) and first-principles molecular dynamics simulation (FPMD) can be applied to study bond breaking and reformation10,27, their high computational cost severely limits their potential usage for, e.g. high-throughput screening. In contrast, reactive force fields such as ReaxFF are more computationally economical and have already found applications in the study of water stability of MOFs28–30. Inspired by the emerging need for studies of ZIFs with defects or disordered structures, we demonstrate the capability of ReaxFF for accurately modeling the melt-quench process of three ZIF compounds.

Driven by recent discoveries of a selection of MOF glasses11,31, we first simulate ZIF-4, Zn(Im)2 (Im – C3H3N2-), a prototypical MOF glass former, whose crystalline/liquid/glass/amorphous solid structures and melting mechanism have been revealed very recently in experiments and molecular simulations10,11,15,32–36. Through comparison with this existing abundant data, we show that ReaxFF can reproduce accurately the most critical structural properties of ZIF-4. We then apply ReaxFF to the ZIF-62 framework, Zn(Im)1.75(bIm)0.25 (bIm – C7H5N2-). The framework shares similar crystal structure with ZIF-4 but has a much lower melting temperature, and an empirically observed ultrahigh glass-forming ability,32,37 due to a tradeoff between electronic and steric effects of the larger organic ligand. Lastly, we extrapolate from ReaxFF melting simulations that electron withdrawing groups such as –NO2 may lower the melting temperature of ZIFs. However, we also show that there is a trade-off with decomposition temperature, as ZIF-77, Zn(nIm)2 (nIm – C3H2N3O2), is shown to oxidize or decompose due to its lower chemical stability. New opportunities of ReaxFF in efficient computational design of ZIFs are discussed in terms of its promising applications in chemical stability and defect engineering of ZIFs.

METHODS 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 25

ReaxFF molecular dynamics simulations The ReaxFF force fields from Ref.38 has been used to simulate ZIF-4, ZIF-62 and ZIF-77 using LAMMPS39. A timestep of 0.25 fs is used for time integration. The initial configurations of ZIF-4, ZIF62, and ZIF-77 crystals are extracted from CSD8 (The Cambridge Structural Database) and are equilibrated at 10K using a Berendsen thermostat. ZIF-4 and ZIF-62 are then heated to 300 K in 2.5 ps with the temperature and the pressure controlled using a Nose-Hoover thermostat/barostat in the isothermalisobaric ensemble (NPT), followed by a further heating to 1500 K within 12.5 ps for melting. ZIF-77 is directly heated from 10 K to 900 K using the same heating rate. ZIF-4 sample is then quenched to 300 K within 12.5ps using an NPT ensemble, to form a glass. For the ReaxFF, a tolerance of 10-6 for the charge equilibrium40–42is used. The average pressure of the system is maintained at 1 atm. The sample contains 2176 atoms, 2368 atoms, and 3360 atoms for the ZIF-4, ZIF-62, and ZIF-77 simulations respectively. In terms of computational cost, the melt-quench simulation (25 ps) for 2176 atoms can be completed in about 1.5 hours using 8 Intel® Xeon® processors (2.2GHz) available in common workstations. Because of the low CPU-time cost, ReaxFF is suitable for simulating many ZIFs materials that contain C, H, O, N, and Zn, and can be extended to cover other elements if needed. We note that in Ref.10, because of the high computational cost of first-principles MD, the ZIF-glass model from reverse Monte Carlo modeling – whose unit cell contains thousands of atoms – cannot be directly simulated.

Energy minimization of ZIF-4 crystal structure The ZIF-4 structure from CSD8 with 272 atoms and without solvent was energy-minimized using the ReaxFF. The equation of state is obtained by varying the unit cell volume of the CSD ZIF-4 structure between -15% and +15% from which the bulk modulus can be calculated by fitting the ReaxFF energy/volume relation using a Birch-Murnaghan equation43. The N–Zn–N bond angle distribution is sampled from the optimized ZIF-4 structure.

4 ACS Paragon Plus Environment

Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Heat capacity and coordination analysis The isobaric heat capacity Cp along the MD trajectory is calculated from the change of enthalpy by finite difference10. Since the PV term has a much larger fluctuation than the internal energy U(T) but its change is less than 0.5% of that of U(T), we directly use the change of U(T)): ‫ܥ‬௣ ቀ

்భ ା்మ ଶ

ቁ=

ு(்మ )ିு(்భ ) ்మ ି்భ



௎(்మ )ି௎(்భ ) ்మ ି்భ

(1)

where T1 and T2 are two consecutive temperatures in the heating cycle. The raw data is block averaged every 0.25 ps to reduce noise. The Zn coordination in all of the ZIF samples is computed using a cut-off distance of 2.5 Å for the Zn–N bond which roughly corresponds to the first minimum in Zn–N pair distribution function.

X-ray and neutron scattering The total scattering structure factor F(Q) and the differential correlation function D(r) are calculated following the procedure in Ref.44. For neutron scattering, the scattering length for deuterium is used because a deuterated ZIF-4 sample was used in the experiments10. Due to the limitation of sample size, D(r) is calculated up to half of the smallest sample dimension. The experimental scattering data are extracted from Ref.10,35. The Zn–N pair distribution function is also computed using the atomic configuration along the MD trajectories upon heating and cooling.

Pore structure analysis The porosity of ZIF-4 crystal, liquid and glass (agZIF-4) was analyzed using the software Zeo++45–47, which uses Voronoi decomposition to compute the accessible and non-accessible volume/surface to a given sphere probe. The data from MD trajectories is averaged over measurement from three successive configurations. A probe radius of 1.2 Å is used for computing the pore characteristics listed in Table 1. For pore size distribution, a smaller radius (0.8 Å) is used for ZIF-4 glass because the connectivity be5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 25

tween pore pockets is poor in the fast-quenched sample. The pore structure is visualized using the OVITO software48.

Results and Discussion Structure and Stability of ZIF-4 In ZIF-4 (Figure 1a), each Zn2+ atom forms a tetrahedron by linking to four N atoms of the organic groups. The N-Zn-N angle distributions of the original CSD structure and the optimized structure are both clustered around the ideal tetrahedral angle 109.5°, ranging from 90° to 135° (Figure 1b)49. According to the potential energy landscape of the ZIF family50–52, ZIF-4 is a highly metastable polymorph of Zn(Im)2 chemistry, which is captured in the equation of state (EOS) of ZIF-4 crystal (Figure 1c). The equilibrium volume, ca. 4570 Å3, from the EOS is very close to the experimental unit cell volume 4383 Å3 34 and 4574/4430 Å3 from density functional calculation53,54. Two fittings to Birch-Murnaghan were conducted by including/excluding data points from higher volumetric strain beyond ~8%. Because the equation of state has a larger curvature near the bottom of the curves than the curvature on both sides, two bulk moduli, 1.99 GPa and 4.97 GPa are obtained from the fitting. While the modulus considering large volumetric strain is close to the experimental measurement, ca. 2.6 GPa by Bennett et al.34, both are within the same order of magnitude of experimentally measured bulk modulus on ZIF-4. The EOS of the highest density Zn(Im)2 polymorph was also calculated as shown in Figure S1 (Supporting Information). The ReaxFF derived bulk modulus (volume) for ZIF-zni from EOS is 14.9 GPa (7100 Å3), very close to the experimental values of 14.0 GPa (6883 Å3)55 and 15.6 GPa (7056 Å3) from density functional theory calculations53.

Melting and Quenching of ZIF-4

6 ACS Paragon Plus Environment

Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

According to earlier experimental work, heating ZIF-4 under argon, leads to the melting of the crystal structure at 865 K11. Cooling then leads to vitrification, with the product termed agZIF-4, in accordance with prior literature56. In order to examine the ZIF-4 liquid in the relatively short time scale accessible to FPMD, Gaillac et al.10 melted the structure at a much higher temperature, with an estimated melting point between 1000 K and 1500 K. Using a similar melting and quenching scheme (Figure 2a), we found the ZIF-4 crystalline structure simulated using ReaxFF, can melt upon heating to 1000 K.

Crystalline ZIF-4 (at 300 K), liquid ZIF-4 (at 1,500 K) and agZIF-4 (at 300 K) have distinctive atomic configurations (Figure 2b). For example, all Zn atoms in the ZIF-4 configuration are 4-fold coordinated. The liquid phase however has a broad distribution of Zn coordinations between 1-fold and 4-fold as a result of Zn-N bond breaking at high temperature10. Interestingly, the Zn coordination in the glass is dominated by 3- and 4-fold configurations, meaning that the 1-fold and 2-fold coordinations are not energetically favorable at low temperature and that the part of the framework with 3-fold Zn is less constrained than in the crystal phase. The glass phase thus comprises both tetrahedral and trigonal Zn centers, and hence the model for agZIF-4 in the current work is distinct from the glass model built using reverse Monte Carlo modelling10, where soft constraints mean that tetrahedral Zn centers are retained.

At room temperature, the density of the ZIF-4 configuration is 1.28 g/cm3 (Figure 3a), which is very close to the crystallographic density of 1.22 g/cm3, and density of 1.25 g/cm3 obtained from FPMD10. Upon heating, the density shows a small decrease to 1.2 g/cm3 at 1,500 K, consistent with the FPMD calculation under constant-pressure condition by Gaillac et al.10. It has been shown in previous work that ZIF-4 undergoes amorphization and recrystallization to the high density ZIF-zni polymorph upon heating and prior to melting33. However, due to the limited time scale in MD and the high heating/cooling rate used in MD, it can be difficult to observe such a phase transition. When cooling down to room temperature, the ZIF-4 liquid vitrifies to become a glass with a higher density of 1.59 g/cm3. 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 25

According to Bennett et al.32, agZIF-4 has a higher density of 1.62 g/cm3 compared to its crystalline counterpart. When the quench rate is 10 times slower, i.e., 9.6 K/ps, the obtained agZIF-4 has a slightly higher density of 1.62 g/cm3, as seen in Figure S2 (Supporting Information), where the glass transition temperature is around 800 K. To the best of our knowledge, the cooling rate effect on the density of agZIF-4 has not been reported in experiments. However, it has been shown that agZIF-62 (i.e. the melt quenched glass formed from ZIF-62) formed from a higher temperature ZIF-62 melt possesses a higher density, indicating a better structural relaxation at higher temperature32.

The heat capacity during the heating process (Figure 3b) resembles the results from FPMD10. It seems there are two peaks associated with the melting. The first one is located around 1000 K and the second around 1300 K. Although the corresponding two peaks in FPMD locate at higher temperature (at 1350 K and 1880 K), the value of the heat capacity from our simulation at each peak is close to those measurements in FPMD. Because the first peak is closer to the melting point of ZIF-4 (863 K), it is more likely to be correlated with the melting process of the sample.

Zn Coordination and Zn-N pair radial distribution function As shown previously, the distribution of Zn coordination environments is different in crystal, liquid, and glass phases. The change of Zn coordination upon heating and cooling is shown in Figure 3c and 3d, where it can be seen that the Zn in the sample remains 4-fold coordinated up to 600 K during heating. At higher temperature, the fraction of 3-fold coordinated Zn increases. The proportions of 4-fold Zn and 3fold Zn flip at 1000 K, which may account for the heat capacity peak at around the same temperature (Figure 3b). At temperatures above 1300 K, there is an apparent increase of 2-fold coordinated Zn as a result of Zn-N bond cleavage in the 3-fold Zn, whose concentration begins to drop. The second Zn-N bond-breaking in the Zn tetrahedra may account for the sharp peak seen in the heat capacity curve above 1200 K, where the energy cost is higher than to break the first Zn-N bond. At 1500K, the liquid sample 8 ACS Paragon Plus Environment

Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

is mostly comprised of 2- and 3-fold Zn coordinations. The change of Zn coordination, shown in Figure 3c, resembles those observed in FPMD simulations10 (dashed lines) and demonstrates the capability of ReaxFF for simulating this process. As the sample cools, the proportions of Zn with different coordinations reverse their changes upon heating, as shown in Figure 3d. There is an increase of 3-fold coordination accompanying a decrease of 2fold coordination, since 2-fold Zn is thermodynamically unfavorable. The 4-fold Zn only increases when temperature is lower than ~1100 K and rises up to 55% at room temperature. In the final glass sample, there are mainly 3- and 4-fold coordinated Zn configurations with a trivial amount of 2-fold Zn, consistent with the atomic configuration in Figure 2b. The change of Zn coordination upon cooling demonstrates the thermal history dependence of the glass topology, which is crucial in glass forming and property controlling. Both two57, and three58 coordinated zinc species have been reported, though in non-controlled atmospheres they are likely to be chemically unstable. Throughout the melting and quenching process, no 5-fold or higher coordination Zn2+ species were observed.

The melting of ZIF-4 has been shown to be driven by exchange of two neighboring Im groups connected with a Zn atom10,15. This is reflected in the Zn-N partial radial distribution (Figure 3e) upon heating. There is general thermal peak broadening when the temperature increases from 300 K to 1500 K, and no medium to long range order is observed beyond 900 K. The positions of the first two peaks are the same as the results from FPMD10. In contrast, there is a clear sharpening of the first two peaks in the liquid upon cooling (Figure 3f). A small peak locating at r = 3 Å appears in the glass, which may come from the under-coordination of Zn2+ as shown in Figure 3d. The medium to long range order in the crystal is not recovered in the agZIF-4. The non-crystallization of the liquid upon cooling is ascribed to the strong steric hindrance experienced by the imidazolate ligand, which prohibits crystallization.

Neutron and X-ray scattering of liquid and glass phases. 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 25

Because the atomic structures of ZIF-4 in the crystal, liquid, and glass phases have been characterized in experiments using neutron and X-ray scattering10,33, we compare the structures from ReaxFF MD simulations with the experimental results, in Figure 4. The structure factor from neutron scattering (Figure 4a) shows that both the ZIF-4 glass and liquid obtained from ReaxFF can accurately reproduce almost all the features that can be seen in the experimental neutron scattering, evidenced by the peak position (dashed lines) and their relative intensity up to a Q value of 40 Å-1. The peak positions in the X-ray scattering of the simulated ZIF-4 glass and liquid seem to shift towards smaller Q (Figure 4b); however, all the peaks existing in the experimental scattering data are reproduced in the simulation structure. Comparing the simulated differential correlation function DN(r) with the experimental neutron scattering data (Figure 4c), the simulated ZIF-4 glass and liquid have the same short range order (SRO) as evidenced by the dashed lines which reflect the intratetrahedral correlation between the imidazolate atoms and the Zn atom31,37. There exists some medium range order (MRO) in the neutron scattering of agZIF-4 beyond r = 7 Å, which is also weak in the ReaxFF-generated glass. In the pair distribution function DX(r) from X-ray scattering (Figure 4d), the SRO peaks of the glass and the liquid phases are generally reproduced in the corresponding simulations, with the peak at 1.33 Å from C-C bonds in the same position while the other major peaks shift a little towards higher r. In all comparisons, the simulated neutron scattering from the ReaxFF-generated glass/liquid structure reproduced most peaks observed in the experiments, which is clear evidence for some capability of ReaxFF in modeling the ZIF-4 glass/liquid system. According to Keen and Bennett35, the neutron data is more sensitive to the arrangement of the imidazole ions and the X-ray is more sensitive to the ZnN4 polyhedra, indicating that the ReaxFF can reproduce not only a reasonable structure from the ZnN4 polyhedra, but also capture the structural features of the imidazole units.

Pore configuration and pore size distribution

10 ACS Paragon Plus Environment

Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The porous architecture of ZIF materials is of high importance for applications in gas sorption/separation5. The pore characteristics of ZIF-4 from ReaxFF simulations are computed and listed in Table 1. At low temperature, the ZIF-4 crystal has a small amount of inaccessible volume, possibly due to torsion of the framework. The size and shape of the accessible porosity is well defined, owing to the crystallographic periodicity. The pores in the agZIF-4 are, in contrast, irregularly distributed throughout the framework and vary in size (Figure 5a). Unlike the crystalline phase, where periodic channels between pores make internal porosity accessible to external gas molecules, the majority in agZIF-4 are inaccessible. The cavity size distribution among the three ZIF-4 phases studied is shown in Figure 5b. The pore structure in the crystal phase has a bimodal distribution, consistent with study by Thornton et al.56 using both positron annihilation lifetime spectroscopy (PALS) and pore analysis on a structural model using the Polymatic code59. The pore size in the ZIF-4 liquid is broadly distributed between 3 and 7 Å and shifts toward smaller cavity size in the glass phase, in agreement with the PALS measurement of ZIF-4 glass56.

ZIF-62 and ZIF-77 To further demonstrate the predictive capability of the ReaxFF method, we have melted ZIF-62 [Zn(Im)1.75(bIm)0.25] (Fig. 6a). This structure shares the same network topology as ZIF-4, but with 1/8 of the imidazolate ligands substituted by benzimidazolate. It has been confirmed in experiments that ZIF-62 has a melting point around 708 K, ca 150 K lower than that of ZIF-432,37. Figure 6b shows that the first heat capacity jump occurs around 900 K, i.e. ca 100 K lower than that of ZIF-4. The lower insilico melting point has been further confirmed in the coordination distribution in Figure 6c where the crossover point between 3-fold coordination and 4-fold coordination locates also at 900K. The decrease in melting point temperature from ZIF-4 to ZIF-62 in our calculations is close to the relative decrease from experiments, which may prove useful for the purposes of pre-screening ZIFs to elucidate their glass-forming abilities. 11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 25

The effects of framework chemistry (linker or node identity), and topology on melting point are of interest in designing new ZIF materials which undergo melting. For example, the benzimidazolate ligand in ZIF-62 is expected to have a stronger affinity for Zn2+. The strengthening of the Zn–N bond would then be expected therefore to increase the melting temperature. Yet the opposite is true, and the reason may lie in that the strong bIm-Zn bonds could possibly weaken the Im-Zn bonds which instead facilitates melting. Another example, ZIF-77, is considered to demonstrate the predictability of ReaxFF in terms of its chemical stability at high temperature. Given the electronic-withdrawing effect of –NO2 group and its limited steric effect32, we expect that the addition of this group to Im, resulting in ZIF-77 (Figure 7a), would weaken the Zn–N bonds and decrease the melting point. The heat capacity result (Figure 7b) shows that the melting of ZIF-77 takes less energy than ZIF-4 and ZIF-62. Zn coordination with N in Figure 7c further shows that there are ~10% 3-fold coordinated Zn at 300 K atoms, indicative of possible spontaneous defects in ZIF-77. However, further analysis of radial distribution function and visualization of O-Zn pairs as seen in Figure 7d and 7e show that the O in the –NO2 group can form bonds with Zn at a temperature as low as 425 K. The amount of Zn oxidation increases with the temperature and it has been found in the melting simulation at temperature higher than 1140 K that part of the C-NO2 bonds also break leading to decomposition of ZIF-77. Based on our calculation, ZIF-77 is not suitable as a glass former due to its low chemical stability at high temperature.

As discussed previously, reactive force fields such as ReaxFF enable the study of processes that involve bond breaking and formation, which is not possible using classical nonreactive potentials. Because of the low computational cost of ReaxFF (see the Methods section) compared to quantum-level simulations, new opportunities are opened for the study of ZIF materials through efficient computational screening. For example, Zheng et al.60 used ReaxFF to show a reduced yield strength of MOF-5 resulted from dislocation. Although the elastic properties of many existing ZIFs have been well documented in 12 ACS Paragon Plus Environment

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the literature61–63, there is a lack of studies on their plastic response in tension/compression, which are crucial in industrial applications. Another application of ReaxFF is to assess the chemical stability of different ZIF compounds in atmospheres including water, CO2, H+, etc. Because degradation of ZIFs in these environments may affect their functionality27, a better understanding of the mechanism is needed before applications of ZIFs. With the promising application of ReaxFF in these aspects, together with the existing database of ZIFs8,9, and its low computational cost and expandability to include other elements, we expect the ReaxFF will promote a better and more efficient design of novel ZIF materials.

CONCLUSIONS In summary, the melt-quench behavior of ZIF-4 was studied by molecular dynamics simulations using ReaxFF. The ReaxFF simulations accurately reproduce the atomic structures of ZIF-4 crystal, liquid, and glass observed in experiments or first-principles MD models. The heat capacity and the change of Zn coordination during heating in ReaxFF simulation reproduces the majority of results from the firstprinciples MD as performed by Ref.10. The retention of undercoordinated Zn during quenching of the glass is interesting and deserves further study in experiments. Porosity analysis of the ZIF-4 crystal and ZIF-4 glass have shown remarkable agreement with measurements from positron annihilation lifetime spectroscopy and from simulation.

The predictive capability of ReaxFF has been further exemplified by simulating the melting process of ZIF-62. Specifically, the electronic effect of benzimidazolate group is demonstrated by heating ZIF-62, where a lower melting point is observed, consistent with recent experiments. With the same force field, ZIF-77 is simulated for which although a lower melting temperature is predicted, according to the electron withdrawing effect of -NO2, the low chemical stability of ZIF-77 at high temperature would make it unsuitable as a glass former. Our ReaxFF simulations of ZIF-4, ZIF-62, and ZIF-77 are the first to use

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 25

ReaxFF to study the glass-forming process of ZIFs. ReaxFF can be used to study other ZIFs that contain C, H, O, N, and Zn.

Emerging research areas such as liquid/glass ZIFs, inelastic deformation, and defect structures require use of reactive force fields that allow for bond breaking and formation. While nonreactive potentials are useful in the study of elastic properties, they are not suitable for application to these emerging areas. Given the low computational cost of ReaxFF compared with other techniques such as first-principles MD and density functional theory, and its expandability to cover more complex systems, ReaxFF is therefore perhaps a suitable framework to enable the more efficient computational design of new ZIF materials.

ASSOCIATED CONTENT Supporting Information The supporting information includes equation of state of ZIF-zni, and density of ZIF-4 at slower cooling rate.

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected] and [email protected] Notes The authors declare no competing financial or non-financial interests. ACKNOWLEDGMENT

14 ACS Paragon Plus Environment

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Y.Y. and J.C.M. acknowledge the Institute for CyberScience Advanced CyberInfrastructure (ICS-ACI) at the Pennsylvania State University for providing computing resources. S.L. acknowledges funding from China Scholarship Council (CSC). T.D.B acknowledges the Royal Society for a University Research Fellowship (UF150021).

REFERENCES (1) Li, H.; Eddaoudi, M.; O’Keeffe, M.; Yaghi, O. M. Design and Synthesis of an Exceptionally Stable and Highly Porous Metal-Organic Framework. nature 1999, 402 (6759), 276. (2) Férey, G.; Serre, C. Large Breathing Effects in Three-Dimensional Porous Hybrid Matter: Facts, Analyses, Rules and Consequences. Chem. Soc. Rev. 2009, 38 (5), 1380–1399. (3) Hayashi, H.; Cote, A. P.; Furukawa, H.; O’Keeffe, M.; Yaghi, O. M. Zeolite A Imidazolate Frameworks. Nat. Mater. 2007, 6 (7), 501. (4) 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. 2006, 103 (27), 10186–10191. (5) Banerjee, R.; Phan, A.; Wang, B.; Knobler, C.; Furukawa, H.; O’keeffe, M.; Yaghi, O. M. HighThroughput Synthesis of Zeolitic Imidazolate Frameworks and Application to CO2 Capture. Science 2008, 319 (5865), 939–943. (6) Phan, A.; Doonan, C. J.; Uribe-Romo, F. J.; Knobler, C. B.; O’keeffe, M.; Yaghi, O. M. Synthesis, Structure, and Carbon Dioxide Capture Properties of Zeolitic Imidazolate Frameworks. Acc Chem Res 2010, 43 (1), 58–67. (7) Bennett, T. D.; Cheetham, A. K.; Fuchs, A. H.; Coudert, F.-X. Interplay between Defects, Disorder and Flexibility in Metal-Organic Frameworks. Nat. Chem. 2017, 9 (1), 11. (8) Allen, F. H. The Cambridge Structural Database: A Quarter of a Million Crystal Structures and Rising. Acta Crystallogr. B 2002, 58 (3), 380–388. (9) Chung, Y. G.; Camp, J.; Haranczyk, M.; Sikora, B. J.; Bury, W.; Krungleviciute, V.; Yildirim, T.; Farha, O. K.; Sholl, D. S.; Snurr, R. Q. Computation-Ready, Experimental Metal–Organic Frameworks: A Tool to Enable High-Throughput Screening of Nanoporous Crystals. Chem. Mater. 2014, 26 (21), 6185–6192. (10) Gaillac, R.; Pullumbi, P.; Beyer, K. A.; Chapman, K. W.; Keen, D. A.; Bennett, T. D.; Coudert, F.X. Liquid Metal–Organic Frameworks. Nat. Mater. 2017, 16 (11), 1149. (11) Bennett, T. D.; Tan, J.-C.; Yue, Y.; Baxter, E.; Ducati, C.; Terrill, N. J.; Yeung, H. H.-M.; Zhou, Z.; Chen, W.; Henke, S.; et al. Hybrid Glasses from Strong and Fragile Metal-Organic Framework Liquids. Nat. Commun. 2015, 6. (12) Dissegna, S.; Epp, K.; Heinz, W. R.; Kieslich, G.; Fischer, R. A. Defective Metal-Organic Frameworks. Adv. Mater. 2018, 1704501. (13) Fang, Z.; Bueken, B.; De Vos, D. E.; Fischer, R. A. Defect-Engineered Metal–Organic Frameworks. Angew. Chem. Int. Ed. 2015, 54 (25), 7234–7254. (14) Boyd, P. G.; Lee, Y.; Smit, B. Computational Development of the Nanoporous Materials Genome. Nat. Rev. Mater. 2017, 2 (8), 17037. (15) 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. 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 25

(16) 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 (12), 124110. (17) Rappé, A. K.; Casewit, C. J.; Colwell, K.; Goddard Iii, W.; Skiff, W. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114 (25), 10024–10035. (18) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94 (26), 8897–8909. (19) Bristow, J. K.; Tiana, D.; Walsh, A. Transferable Force Field for Metal–Organic Frameworks from First-Principles: BTW-FF. J. Chem. Theory Comput. 2014, 10 (10), 4644–4652. (20) Coupry, D. E.; Addicoat, M. A.; Heine, T. Extension of the Universal Force Field for Metal– Organic Frameworks. J. Chem. Theory Comput. 2016, 12 (10), 5215–5225. (21) Addicoat, M. A.; Vankova, N.; Akter, I. F.; Heine, T. Extension of the Universal Force Field to Metal–Organic Frameworks. J. Chem. Theory Comput. 2014, 10 (2), 880–891. (22) Vanduyfhuys, L.; Vandenbrande, S.; Verstraelen, T.; Schmid, R.; Waroquier, M.; Van Speybroeck, V. QuickFF: A Program for a Quick and Easy Derivation of Force Fields for Metal-Organic Frameworks from Ab Initio Input. J. Comput. Chem. 2015, 36 (13), 1015–1027. (23) Vanduyfhuys, L.; Vandenbrande, S.; Wieme, J.; Waroquier, M.; Verstraelen, T.; Van Speybroeck, V. Extension of the QuickFF Force Field Protocol for an Improved Accuracy of Structural, Vibrational, Mechanical and Thermal Properties of Metal–Organic Frameworks. J. Comput. Chem. 2018. (24) Bristow, J. K.; Skelton, J. M.; Svane, K. L.; Walsh, A.; Gale, J. D. A General Forcefield for Accurate Phonon Properties of Metal–Organic Frameworks. Phys. Chem. Chem. Phys. 2016, 18 (42), 29316–29329. (25) Heinen, J.; Dubbeldam, D. On Flexible Force Fields for Metal–Organic Frameworks: Recent Developments and Future Prospects. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2018, e1363. (26) Liu, B.; Smit, B. Molecular Simulation Studies of Separation of CO2/N2, CO2/CH4, and CH4/N2 by ZIFs. J. Phys. Chem. C 2010, 114 (18), 8515–8522. (27) Han, C.; Zhang, C.; Tymińska, N.; Schmidt, J. R.; Sholl, D. S. Insights into the Stability of Zeolitic Imidazolate Frameworks in Humid Acidic Environments from First-Principles Calculations. J. Phys. Chem. C 2018, 122 (8), 4339–4348. (28) Han, S. S.; Choi, S.-H.; van Duin, A. C. Molecular Dynamics Simulations of Stability of Metal– Organic Frameworks against H 2 O Using the ReaxFF Reactive Force Field. Chem. Commun. 2010, 46 (31), 5713–5715. (29) Liu, X. Y.; Pai, S. J.; Han, S. S. ReaxFF Molecular Dynamics Simulations of Water Stability of Interpenetrated Metal–Organic Frameworks. J. Phys. Chem. C 2017, 121 (13), 7312–7318. (30) Greathouse, J. A.; Allendorf, M. D. The Interaction of Water with MOF-5 Simulated by Molecular Dynamics. J. Am. Chem. Soc. 2006, 128 (33), 10678–10679. (31) Bennett, T. D.; Cheetham, A. K. Amorphous Metal–Organic Frameworks. Acc. Chem. Res. 2014, 47 (5), 1555–1562. (32) Bennett, T. D.; Yue, Y.; Li, P.; Qiao, A.; Tao, H.; Greaves, N. G.; Richards, T.; Lampronti, G. I.; Redfern, S. A.; Blanc, F.; et al. Melt-Quenched Glasses of Metal–Organic Frameworks. J. Am. Chem. Soc. 2016, 138 (10), 3484–3492. (33) 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 MetalOrganic Framework. Phys. Rev. Lett. 2010, 104 (11), 115503. (34) 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 (28), 7983–7985. (35) Keen, D. A.; Bennett, T. D. Structural Investigations of Amorphous Metal-Organic Frameworks Formed Via Different Routes. Phys. Chem. Chem. Phys. 2018. 16 ACS Paragon Plus Environment

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(36) 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 (13), 3067– 3071. (37) Qiao, A.; Bennett, T. D.; Tao, H.; Krajnc, A.; Mali, G.; Doherty, C. M.; Thornton, A. W.; Mauro, J. C.; Greaves, G. N.; Yue, Y. A Metal-Organic Framework with Ultrahigh Glass-Forming Ability. Sci. Adv. 2018, 4 (3). (38) Shin, Y. K.; Ashraf, C. M.; Van Duin, A. C. Development and Applications of the ReaxFF Reactive Force Field for Biological Systems. In Computational Materials, Chemistry, and Biochemistry: From Bold Initiatives to the Last Mile; in press; Springer, 2017. (39) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117 (1), 1–19. (40) Rappe, A. K.; Goddard III, W. A. Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95 (8), 3358–3363. (41) Nakano, A. Parallel Multilevel Preconditioned Conjugate-Gradient Approach to Variable-Charge Molecular Dynamics. Comput. Phys. Commun. 1997, 104 (1–3), 59–69. (42) Aktulga, H. M.; Fogarty, J. C.; Pandit, S. A.; Grama, A. Y. Parallel Reactive Molecular Dynamics: Numerical Methods and Algorithmic Techniques. Parallel Comput. 2012, 38 (4–5), 245–259. (43) Birch, F. Finite Elastic Strain of Cubic Crystals. Phys. Rev. 1947, 71 (11), 809. (44) Keen, D. A. A Comparison of Various Commonly Used Correlation Functions for Describing Total Scattering. J. Appl. Crystallogr. 2001, 34 (2), 172–177. (45) 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. Microporous Mesoporous Mater. 2012, 149 (1), 134–141. (46) Martin, R. L.; Smit, B.; Haranczyk, M. Addressing Challenges of Identifying Geometrically Diverse Sets of Crystalline Porous Materials. J. Chem. Inf. Model. 2011, 52 (2), 308–318. (47) 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. Mol. Graph. Model. 2013, 44, 208–219. (48) Stukowski, A. Visualization and Analysis of Atomistic Simulation Data with OVITO–the Open Visualization Tool. Model. Simul. Mater. Sci. Eng. 2010, 18 (1), 015012. (49) Beake, E.; Dove, M.; Phillips, A.; Keen, D.; Tucker, M.; Goodwin, A.; Bennett, T.; Cheetham, A. Flexibility of Zeolitic Imidazolate Framework Structures Studied by Neutron Total Scattering and the Reverse Monte Carlo Method. J. Phys. Condens. Matter 2013, 25 (39), 395403. (50) Baburin, I. A.; Leoni, S. The Energy Landscapes of Zeolitic Imidazolate Frameworks (ZIFs): Towards Quantifying the Presence of Substituents on the Imidazole Ring. J. Mater. Chem. 2012, 22 (20), 10152–10154. (51) Galvelis, R.; Slater, B.; Cheetham, A. K.; Mellot-Draznieks, C. Comparison of the Relative Stability of Zinc and Lithium-Boron Zeolitic Imidazolate Frameworks. CrystEngComm 2012, 14 (2), 374–378. (52) Lewis, D. W.; Ruiz-Salvador, A. R.; Gómez, A.; Rodriguez-Albelo, L. M.; Coudert, F.-X.; Slater, B.; Cheetham, A. K.; Mellot-Draznieks, C. Zeolitic Imidazole Frameworks: Structural and Energetics Trends Compared with Their Zeolite Analogues. CrystEngComm 2009, 11 (11), 2272– 2276. (53) 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 (2), 375–382. (54) Adhikari, P.; Xiong, M.; Li, N.; Zhao, X.; Rulis, P.; Ching, W.-Y. Structure and Electronic Properties of a Continuous Random Network Model of an Amorphous Zeolitic Imidazolate Framework (a-ZIF). J. Phys. Chem. C 2016, 120 (28), 15362–15368. (55) Spencer, E. C.; Angel, R. J.; Ross, N. L.; Hanson, B. E.; Howard, J. A. Pressure-Induced Coop17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(56) (57)

(58) (59) (60)

(61)

(62)

(63)

Page 18 of 25

erative Bond Rearrangement in a Zinc Imidazolate Framework: A High-Pressure Single-Crystal X-Ray Diffraction Study. J. Am. Chem. Soc. 2009, 131 (11), 4022–4026. Thornton, A.; Jelfs, K.; Konstas, K.; Doherty, C.; Hill, A.; Cheetham, A.; Bennett, T. Porosity in Metal–Organic Framework Glasses. Chem. Commun. 2016, 52 (19), 3750–3753. Ellison, J. J.; Ruhlandt-Senge, K.; Hope, H. H.; Power, P. P. Synthesis and Characterization of the New Selenolate Ligand-SeC6H3-2, 6-Mes2 (Mes= C6H2-2, 4, 6-Me3) and Its Two-Coordinate Zinc and Manganese Derivatives: Factors Affecting Bending in Two-Coordinate Metal Complexes with Aryl-Substituted Ligands. Inorg. Chem. 1995, 34 (1), 49–54. Darensbourg, D. J.; Niezgoda, S. A.; Draper, J. D.; Reibenspies, J. H. Trigonal-Planar Zinc (II) and Cadmium (II) Tris (Phenoxide) Complexes. Inorg. Chem. 1999, 38 (6), 1356–1359. Abbott, L. J.; Hart, K. E.; Colina, C. M. Polymatic: A Generalized Simulated Polymerization Algorithm for Amorphous Polymers. Theor. Chem. Acc. 2013, 132 (3), 1334. Zheng, B.; Fu, F.; Wang, L. L.; Wang, J.; Du, L.; Du, H. Effect of Defects on the Mechanical Deformation Mechanisms of Metal–Organic Framework-5: A Molecular Dynamics Investigation. J. Phys. Chem. C 2018, 122 (8), 4300–4306. Tan, J. C.; Cheetham, A. K. Mechanical Properties of Hybrid Inorganic–Organic Framework Materials: Establishing Fundamental Structure–Property Relationships. Chem. Soc. Rev. 2011, 40 (2), 1059–1080. Tan, J. C.; Bennett, T. D.; Cheetham, A. K. Chemical Structure, Network Topology, and Porosity Effects on the Mechanical Properties of Zeolitic Imidazolate Frameworks. Proc. Natl. Acad. Sci. 2010, 107 (22), 9938–9943. Tan, J.-C.; Civalleri, B.; Lin, C.-C.; Valenzano, L.; Galvelis, R.; Chen, P.-F.; Bennett, T. D.; Mellot-Draznieks, C.; Zicovich-Wilson, C. M.; Cheetham, A. K. Exceptionally Low Shear Modulus in a Prototypical Imidazole-Based Metal-Organic Framework. Phys. Rev. Lett. 2012, 108 (9), 095502.

1 2 3 4 5 6

18 ACS Paragon Plus Environment

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figures

Figure 1. (a) Unit cell viewed along the c axis of ZIF-4 crystal. Big gray atoms are Zn, blue atoms are N, small gray atoms are C and white atoms are H. (b) N-Zn-N angle distribution. (c) Equation of state from minimization including bulk modulus and equilibrium volume from the fitting of the curves.

Figure 2. (a) Melt-quenching process of ZIF-4. Atomic configurations of the system (2176 atoms) at three moments (1,2 and 3) are shown on the bottom. Whereas all Zn atoms in ZIF-4 are 4-fold coordinated, there is a distribution of coordination associated with Zn atoms in both the liquid and the meltquenched glass, as denoted by the red arrows. For clarity, only a thin layer of the liquid and the glass samples are shown in 2 and 3. Large gray balls are Zn atoms, small gray balls are C atoms and blue balls are N atoms.

ACS Paragon Plus Environment

19

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 25

Figure 3. (a) Variation of density with temperature upon heating and cooling of the ZIF-4 sample. (b) Evolution of the heat capacity with temperature. The experimental melting point of ZIF-4 at 865 K is indicated by the blue vertical line. (c) and (d) Distribution of zinc coordination numbers as a function of temperature upon heating (c) and cooling (d), from 0-fold coordinated to 5-fold. 5-fold coordination is close to the y = 0 axis. The dashed lines in (c) are FPMD results10. (e) Zn-N pair distribution function upon heating. (f) Zn-N pair distribution function upon cooling. The FPMD results in (b) and (c) are reproduced with permission from Ref.10.

ACS Paragon Plus Environment

20

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4. Structure factors (a and b) and differential correlation function (c and d) from neutron and Xray scattering. Black and green curves are scattering data extracted from existing experimental data. The data, i.e., the spurious ripples, below 1 Å in the black curve are not reliable. Blue curves are ZIF-4 glass at 300K. Red curves are ZIF-4 liquid at 1500K. Green curves correspond to ZIF-4 glass quenched from 1500K to 300K, and blue curves correspond to crystalline ZIF-4 glass at 300K. The results from simulation have been rescaled for comparison with the experiments. The top three curves are shifted for clarity. The experimental results are reproduced with permission from Ref.10,35.

ACS Paragon Plus Environment

21

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 25

Figure 5. (a) Pore structure and atomic configuration of ZIF-4 crystal and glass. Free volume is represented by the golden areas. (b) Pore size distribution of ZIF-4 crystal, liquid, and glass using probe of different sizes (0.01 Å for glass and 1.2 Å for the other two configurations).

Figure 6. (a) Atomic configuration of crystalline ZIF-62 at 300 K. (b) Evolution of the heat capacity with temperature. The vertical blue line stands for the melting point. (c) Distribution of Zn coordination number as a function of temperature upon heating.

Figure 7. (a) Atomic configuration of ZIF-77 crystal at 300 K using the current ReaxFF. (b) Evolution of the heat capacity with temperature. (c) Distribution of Zn coordination number as a function of temperature upon heating. (d) Radial distribution function (RDF) of O-Zn in ZIF-77. (e) O-Zn bonds in ZIF-77 at ~680 K. A bond cutoff of 2.1 Å is used for O-Zn bond length. The atom color scheme in (a) and (e): large gray balls are Zn atoms, small gray balls are C atoms, white balls are H atoms, blue balls are N atoms, and red balls are O atoms. ACS Paragon Plus Environment

22

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ACS Paragon Plus Environment

23

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 25

Table 1. Porosity analysis for the ReaxFF generated ZIF-4 crystal, liquid and glass simulated with Zeo++45–47. A probe radius of 1.2 Å is used. ASA (m2/g)

ISA (m2/g)

AV (cm3/g)

IV (cm3/g)

Di (Å)

Df (Å)

Dif (Å)

ZIF-4 crystal from CSD8

1353



0.05



5.0

2.6

5.0

ZIF-4 crystal (~300 K)

1114

162

0.04

0.005

5.4

2.8

5.1

ZIF-4 liquid (~1500 K)

1384

78

0.06

0.002

6.5

2.7

6.1

ZIF-4 glass (~300 K)



295



0.007

5.4

1.8

4.9

ASA: accessible surface area. ISA: inaccessible surface area. AV: accessible volume. IV: inaccessible volume. Di: the largest included sphere. Df: the largest free sphere. Dif: the largest included sphere along the free sphere path (Dif).

ACS Paragon Plus Environment

24

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table of Contents (TOC Graphic)

ACS Paragon Plus Environment

25