Flexible and Transferable ab Initio Force Field for Zeolitic Imidazolate

Subscriber access provided by ECU Libraries

Article

A Flexible and Transferable Ab-initio Force Field for Zeolitic Imidazolate Frameworks: ZIF-FF Tingting Weng, and Jordan R. Schmidt J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b12311 • Publication Date (Web): 05 Mar 2019 Downloaded from http://pubs.acs.org on March 11, 2019

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 40 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

A Flexible and Transferable Ab-initio Force Field for Zeolitic Imidazolate Frameworks: ZIF-FF Tingting Weng and J.R. Schmidt∗ Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin–Madison, Madison, Wisconsin, 53706, United States E-mail: [email protected]

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

Abstract We have developed a transferable ab-initio intramolecular force field for zeolitic imidazolate frameworks (ZIFs), “ZIF-FF”, that is capable of quantitatively describing the structural properties and relative stabilities of ZIFs. In contrast to nearly all prior force fields, ZIF-FF properly describes the relative stability of ZIF polymorphs, a crucial element in ZIF nucleation and crystal growth. Beginning with a general Amber force field (GAFF), Zn-related force field parameters were optimized against dispersion-corrected DFT-calculated properties using a genetic algorithm. We validated the resulting force field by examining bond and angle distributions, phonon density of states, mechanical properties, diffusion properties and via modeling a ZIF amorphization process. Furthermore, we find that ZIF-FF is transferable, successfully describing relative stability of various ZIF surface structures, as well as the densities of ZIFs with diverse functionalized linkers.

Introduction Metal-organic frameworks (MOFs), a fascinating class of nanoporous material, have attracted considerable research interest due to their widespread potential applications 1 in gas storage, 2,3 separations, 4,5 catalysis, 6,7 drug delivery, 8,9 and many other fields. 10 In particular, zeolitic imidazolate frameworks (ZIFs), MOFs with tetrahedral metal cations (e.g. Zn, Co) bridged by imidazolate (Im) linkers, exhibit high thermal and chemical stability 11 and thus have gained a particular focus in terms of their synthesis, properties, and applications. 12 ZIFs are topologically isomorphic with zeolite structures due to the similarity in both their connectivity and of the bridging Zn-Im-Zn / Si-O-Si angles. As with zeolites, the large number of possible topologies, metal cations, and linker functionalizations have led to a vast variety of ZIF structures, opening up the possibility of rational design of ZIFs for specific applications. Yet targeted synthesis of ZIFs with a particular structure or polymorph is 2 ACS Paragon Plus Environment

Page 2 of 40

Page 3 of 40 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


challenging: subtle changes in the reaction parameters often leads to the formation of alternative crystalline or amorphous products. 13,14 Therefore it is of great importance to develop predictive computational models to simulate ZIF crystallization, which would shed light on the ZIF nucleation/crystallization mechanism(s) and guide targeted synthesis. Considering the significant polymorphism exhibited by ZIFs, an important prerequisite for any realistic simulation of ZIF nucleation and/or crystal growth is the development of accurate force field to describe intra-molecular flexibility of ZIFs and to energetically differentiate various ZIF structures. Yet existing ZIF force fields can be largely grouped into two categories: rigid, inter-molecular force fields, focusing on ZIF-adsorbate interactions (and often used to computed adsorption isotherms via Grand Canonical Monte Carlo); 15–19 or flexible ZIF force fields, which typically adapt inter-molecular terms from standard, generic force fields (e.g. UFF, 20 Dreiding, 21 Amber 22 ), and which are often employed where ZIF flexibility may be important (e.g. modeling diffusion within a ZIF framework). Examples in this latter category include the work of Hertäg et al., 23 who utilized parameters derived on the basis of electronic structure calculations on ZnN4-clusters (from metalloproteins), in conjunction with Amber force fields, to investigate CH4 and H2 diffusion in a ZIF-8 membrane. Their flexible force field reproduced experimentally-measured adsorption selectivity, however, could not reproduce the equilibrium bond lengths, angles, and torsions measured by via x-ray diffraction (XRD). A similar approach was utilized by Zheng et al. 24 to develop a flexible force field to investigate CO2 gas diffusion within ZIF-8. The resulting force field was able to reproduce the experimental equilibrium crystal structure determined via XRD. However, Gee and Sholl 25 found that this modified version of general Amber force field (GAFF) was incapable of capturing the relative energetics of various ZIF polymorphs, or of accurately describing the Zn-N bond and N-Zn-N angle distributions. In related work, Jiang and co-workers modified the Amber force field with Zn-related bonded parameters fit to the lattice constants of ZIF-8 (or cluster analogues), to reproduce 3 ACS Paragon Plus Environment


its structural, mechanical, and thermophysical properties, and to simulate gas adsorption and diffusion within the ZIF-8 framework. 26–28 Later, ab initio data were used by Krokidas et al. with a [Zn(MeIm)4]2– cluster model to fit bonded parameters in flexible force field. 29 Gabrieli et al. employed force matching of ZIF-8, using a fully periodic model system. 30 In each of these cases, these force fields were parameterized for specific ZIF systems, rather than a spectrum of ZIF topologies. As such, they should not be expected to reproduce either the structures nor the relative energetics of alternative ZIF topologies. In general, very little work has focused on developing transferable force fields for ZIFs or other MOFs that accurately reflect the relative stability (or “strain energy”) of various polymorphs. In one of the few such cases, Schmid and co-workers developed a flexible force field, MOF-FF, based on ab initio reference data (cluster calculations) for a variety of MOFs; 31,32 this methodology was subsequently extended and applied to covalent organic frameworks (COFs). 33,34 Later, D¨ urholt et al. derived a coarse-grained variant of MOFFF, MOF-FF-CGNB, which was used to compute the relative stability of tbo and pto topologies. 35 Yet to our best knowledge, there have not been any attempts to develop ab initio force fields based on strain energy analysis of ZIFs. Our present goal is thus to develop a transferable and flexible (non-rigid)) force field capable of quantitatively describing the structural and mechanical properties of a wide variety of ZIF structures, while simultaneously reproducing their relative stability. Periodic density functional theory (DFT) calculations on various ZIF structures are utilized to determine their relative stability, i.e. strain energy. These stabilities, along with the ZIF structures, provide input for optimizing our force field. Because intra-molecular bonded terms account for the majority of the strain associated with the distortions of bonds, angles, and dihedrals that characterize various ZIF topologies, our force field will mainly focus on fitting intramolecular bonded terms. We first choose our model system and perform electronic structure calculations to op4 ACS Paragon Plus Environment

Page 4 of 40

Page 5 of 40 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


timize structures and compute strain energy. Next we present our force field development methodology, optimizing Zn-related parameters using a genetic algorithm (GA). We validate the resulting force field, ZIF-FF, by examining bond and angle distributions among various polymorphs, phonon density of states, and mechanical properties. We then utilize ZIF-FF to model the pressure-induced amorphization of ZIF-8, making direct comparison against experiment. Finally, we demonstrate the transferability of the resulting force field by examining predicted ZIF surface stability and the bulk structure of various functionalized ZIFs.

Methods Model Systems We choose 25 polymorphs of the Zn(MeIm)2 (MeIm = 2-methylimidazolate) system 36 as models for our ab-initio force field development. Among these 25 polymorphs, at least the SOD and DIA topologies are directly synthetically accessible (for the MeIm linker), while many of the other topologies are accessible for alternative choices of metal cation or differing linker functionalization. In particular, the SOD topology is known as ZIF-8, and is a prototypical ZIF and an excellent candidate for various applications. 37,38 Examples of four different polymorphs are shown in Fig. 1. The relative stability of each polymorph can be characterized by a strain energy, defined as the relative energy per formula unit (Zn(MeIm)2), measured relative to the SOD topology:

∆Ei =

E(Zi ) E(ZSOD ) − Zi ZSOD

where Zi is the number of formula units in the system with polymorph i.


(1)


Figure 1: Examples of four ZIF polymorphs: SOD, QTZ, DFT, LTA.

Reference Calculations Electronic structure calculations were performed to provide benchmarks for the ab-initio force field development. The initial structures of the 25 Zn(MeIm)2 (MeIm = 2-methylimidazolate) polymorphs were all obtained from the Cambridge Structural Database (CSD), and were derived from experimental measurements (SOD 12,39 and DIA 40 ) or computationally (other topologies – yielding structures analogous to known ZIFs) by Baburin and Leoni. 36 The Vienna Atomistic Simulation Package (VASP) 41–44 was utilized to conduct periodic density functional theory (DFT) calculations of these structures. These calculations utilized a projector-augmented wave (PAW) 45,46 treatment of core electrons in conjunction with a 600 eV energy cutoff. All structural optimizations utilized the PBE generalized gradient approximation exchange-correlation functional 47 with both lattice constants and atom positions relaxed. Dispersion corrections were included in all calculations via an empirical Grimme D3 correction. 48 Note that inclusion of dispersion is essential to reproduce experimental ZIF densities (see Table S1 and Fig. S1 in the Supporting Information). For each polymorph structure, the k-point grid was converged separately with an energy tolerance of 10−6 eV, and forces are converged to a tolerance of 0.01 eV/˚ A. The optimized structures that provided the basis for force field development are included in the Supporting Information.


Page 6 of 40

Page 7 of 40 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


Force Field Development Several prior works have developed intra-molecular force field for ZIF-8 based on the Amber force field. 23,24,26–28,49–51 We choose an existing modified version of general Amber force field (GAFF) 24,50,52 as the starting point for further optimization and refinement. The GAFF energy expression is given by Eq. (2).

Etot = Ebond + Eangle + Eproper + Eimproper + Evdw + Eelec X X X = Kr (r − r0 )2 + Kθ (θ − θ0 )2 + Kφ [1 + cos(nφ − φ0 )] bonds

+

X

proper

angles

Kξ [1 + cos(nξ − ξ0 )] +

improper

X C6 r6

−

(2)

C12 X qi qj + r12 r2

Fig. 2 illustrates the secondary building unit (SBU) of the Zn(MeIm)2 systems, with each atom type labeled. GAFF has been shown to give good agreement with experiment for both the equilibrium structural parameters and diffusion of CO2 in ZIF-8. 24 However, it was shown by Gee and Sholl 25 that relative energies of different topologies, as computed using GAFF, are in conflict with results shown by 36 using DFT(PBE) calculations. We find a similar result (vide infra): that unoptimized GAFF is incapable of accurately describing the strain energies calculated by DFT(PBE-D3) among different polymorphs, and also yields poor agreement with a number of related structural properties. Nonetheless, we take GAFF as our starting point for further targeted optimized.

Bonded Parameters To effective capture the relative energetics of the ZIF polymorphs (i.e. ZIF strain energy), our force field development will mainly focus on the bonded terms. Based on the distribution of bond distances and angles found in the DFT reference structures among the 25 polymorphs (at 0K, see Fig. 3 and Fig. 4), we find that only Zn-related bonds and angles exhibit broad



Figure 2: SBU of Zn(MeIm)2 systems with labeled atom types. distributions, while all other values are extremely consistent among the various polymorphs. In other words, it is Zn-related parameters that distinguish among various polymorphs. Therefore, our target is to fit Zn-related bonded terms for force field (FF) on the basis of first principles reference data of periodic model systems. In contrast, linker (MeIm–) parameters have only a minor effect on strain energy analysis. However, to further reduce the error in strain energy analysis and to better reproduce phonon frequencies of linkers (discussed in phonon spectrum section), we substitute GAFF linker parameters with those obtained by Gabrieli et al. 30 via force matching. These optimized linker intra-molecular parameters (which utilize the related CHARMM functional form) are included in Table S3 in the Supporting Information.

Non-bonded Parameters We directly adopt non-bonded parameters from GAFF force field. 25 In GAFF, Amber parameters with LJ potential and Lorenz-Berthelot combination rules were utilized. Van der Waals interactions were calculated using a cutoff radius of 13 ˚ A and include long-range tail correction. Electrostatic interactions were treated via the particle-particle particle-mesh (PPPM) 53 method with same cutoff radius (13 ˚ A). Point charges were assigned to each atom using the Density Derived Electrostatic and Chemical (DDEC) method by Manz and


Page 8 of 40

Page 9 of 40 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


Figure 3: All 7 types of bond distributions of DFT reference structures with 25 polymorphs. The sum of each bin distribution has been normalized to 1 and bin width is set to be 0.01 ˚ A.

Figure 4: All 11 types of angle distributions of DFT reference structures with 25 polymorphs. The sum of each bin distribution has been normalized to 1 and bin width is set to be 1°.



Page 10 of 40

Sholl. 54,55 We found that the resulting framework charges of all 25 polymorphs were very similar. As a result, we can use only one set of point charges (DDEC charges of ZIF-8 with SOD polymorph) for all 25 polymorph structures. In addition, we parameterized an alternative version of the ZIF-FF force field that is compatible with the inter-molecular force field developed by McDaniel and Schmidt 18 using symmetry-adapted perturbation theory (SAPT). This SAPT-based FF was originally successful in reproducing gas adsorption in rigid ZIF frameworks. To extend to flexible framework systems, we replace original charge method based on cluster model with DDEC charges and refit charge associated Aelec parameters. The resulting set of bonded and nonbonded parameters parameters are shown in Table S5 in the Supporting Information. This modified force field parameters generate similar results in strain energy (see Fig. S3 in the Supporting Information).

Fitting Objectives We parameterized ZIF-FF to reproduce a wide variety of structural and energetic properties, using the DFT-calculated values for the 25 polymoprhs as our benchmark values. To capture relative energetic stability of the ZIF polymorphs, we utilize the strain energy per formula unit measured relative to the SOD polymorph. Based on this definition, we can calculate the first fitting objective χstrain , the root mean square deviation (RMSD) of strain energy, shown in Eq. (3). Here N equals total number of polymorphs, i.e. 25.

rP χstrain =

FF i (∆Ei

− ∆EiDF T )2 N

(3)

The second objective is to reproduce accurate structural information. For polymorph i, the structural deviation relative to the DFT reference structure is calculated as RMSD of


Page 11 of 40 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


atomic positions δi in Eq. (4), where Mi equals the total number of atoms in the structure T with polymorph i, rDF represents the atom position for atom j in polymorph i DFT referji

ence structure, and r0 Fji F denotes the atom position for atom j in polymorph i FF optimized structure; both positions are expressed in direct coordinates to eliminate any effect due to slight differences in the DFT vs. FF lattice constants. The second fitting objective, χposition is thus calculated as the average of RMSD δi for each polymorph i, v u Mi uX T 2 |r0 Fji F − rDF | /Mi δi = t ji

(4)

ji

P χposition =

i δi

N

(5)

The final fitting objective is related to RMSD of lattice constants (see Eq. (6)), where Li = (ai , bi , ci ) represents lattice constants for the structure with polymorph i, rP χlattice =

i

T 2 |LFi F − LDF | i N

(6)

Optimization Algorithm With multiple non-linear Zn-related parameters to optimize, we use a genetic algorithm (GA) to optimize and fit our force filed parameters. GAs are well suited for nonlinear multiple optima optimization problems with strongly correlated adjustable parameters. 56 In the present work, we utilize a Python implementation of GA in DEAP. 57 We have interfaced a serial version of DEAP with the LAMMPS package. 58 The evaluation of the fitness of an individual parameter set comprises the following basic steps: (1) FF geometry optimization



with both atom position and lattice constants relaxation of the reference structure based on the selected functional form and a trial set of Zn-related parameters; (2) data collection of strain energy, atom positions, and lattice constants for all 25 topologies; (3) calculation of RMSD of strain energies, atom positions and lattice constants (i.e. χstrain , χposition , χlattice between a trial solution and the reference data; (4) computation of the fitness function F (F = χstrain + χposition + χlattice ) for each set of trial parameters to optimize. Note that the numerical value χstrain (in units of kcal/mol) is approximately an order magnitude larger than χposition and χlattice (both in units of ˚ A). We choose to retain this weighting, implying that the χstrain has the largest weighting during the optimization process.

Results and Discussion Force Field Parameterization We optimized the Zn-related parameters in conjunction with a GAFF functional form, against the metrics described above. To further improve the accuracy of our force field, the parameters of linker MeIm– from GAFF were also replaced with optimized parameters derived by Gabrieli et al. 30 on the basis of force matching; all Zn-related parameters were reoptimized with GA. Note that the GA-optimized parameters are well reproducible based on the initial conditions we choose. The original ranges of Zn-related parameters are listed in Table S2 in SI. The resulting optimized bonded parameters are listed in Table 1, with the original GAFF 50 and force matched parameters also listed for comparison. Notably, both the force matched and our GA re-optimized (ZIF-FF) parameters show similar trends as compared to the original GAFF (note especially angle force constants Kθ ). For the nonbonded interactions, we utilize the DDEC charges with van der Waals parameters taken from GAFF. 50 The complete set of all parameters are summarized in Table S3 the Supporting Information. It is important to note that the fit is sensitive to the choice of the 1-4 scaling 12 ACS Paragon Plus Environment

Page 12 of 40

Page 13 of 40 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


factor for Coulomb interaction. Therefore, the 1-4 scaling factor for Coulomb interaction in ZIF-FF was also optimized, yielding a value of 0.6874; this value differs slightly from the conventional choice for GAFF (0.8333). Table 1: ZIF-FF Zn-related Parameters Bond Zn-N

Kr (kcal · mol−1 · ˚ A−2 ) 74.50 (78.50, 67.62)

a

r0 (˚ A) 2.024 (2.011, 2.017)

Angle Kθ (kcal · mol−1 · rad−2 ) θ(o ) Zn-N-C2 11.36 (32.48, 12.23) 126.95 (126.4, 127.028) Zn-N-C1 14.43 (48.68, 12.15) 126.85 (128.33, 125.721) N-Zn-N 11.09 (35.24, 5.24) 109.42 (109.48, 109.360) −1 Proper Kφ (kcal · mol ) n φ(o ) H2-C2-N-Zn 1.056 (2.325, 0.00) 2 180 C2-C2-N-Zn 1.416 (2.325, 0.82) 2 180 C2-N-Zn-N 0.021 (0.00, 0.10) 3 180 C3-C1-N-Zn 0.227 (5.00, 0.00) 2 180 Zn-N-C1-N 0.614 (5.00, 1.26) 2 180 C1-N-Zn-N 0.026 (0.00, 0.00) 3 0 Improper N-Zn-C1-C2

Kξ (kcal · mol−1 ) 0.056 (1.1, 0.00)

n 2

ξ(o ) 180

a. Values given in parenthesis are for comparison and are from from GAFF 50 and the work by Gabrieli et al., 30 respectively.

To assess the quality of our fit, we examine three properties (strain energy, crystal structures and density) that are directly associated with our three fitting objectives. Compared with the DFT reference data, the original, unoptimized GAFF cannot accurately capture the relative stability of the 25 ZIF polymorphs (see Fig. 5). Optimizing only the Zn-related bond parameters (while leaving the linker parameters at their original GA-FF value) yields an intermediate “GA-GAFF” force field that greatly reduces the error in the strain energy; the detailed Zn-related parameters of GA-GAFF is shown in Table S4 in the Supporting Information. To further improve the accuracy, we also replaced the linker 13 ACS Paragon Plus Environment


Figure 5: Comparison of the FF- and DFT-calculated in strain energy (∆E) for the 25 ZIF (Zn(MeIm)2) polymorphs. Data obtained using force filed parameters marked as GAFF (in green) is taken from ref. 50, those marked as GA-GAFF (in purple) utilize GA optimized Zn-related parameters and those marked as ZIF-FF (in orange) use parameters developed in this work.


Page 14 of 40

Page 15 of 40 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


parameters with those previously obtained via a force matching method, 30 yielding the final ZIF-FF force field and an additional slight reduction in the error, and excellent agreement with the DFT-calculated strain energies. (Note that the force-matched linker parameters alone cannot accurately capture the relative stability of the ZIFs; see Fig. S2 in the Supporting Information).

Figure 6: RMSD per atom of the 25 polymorph crystal structures as optimized by both GAFF and ZIF-FF from the DFT reference structures. We also compare the fidelity of ZIF-FF predicted crystal structures to the DFT benchmarks. The RMSD of GAFF / ZIF-FF vs. DFT-predicted atomic positions is shown in Fig. 6. We find that ZIF-FF yields far more accurate crystal structures as compared to the original GAFF. In particular, the average RMSD of the atomic positions for ZIF-FF is only 0.2 ˚ A, much lower than the corresponding value for GAFF (0.5 ˚ A). The optimized force field also dramatically reduces the maximum errors observed from more than 1.0 ˚ A(GAFF, e.g. unc, mmt, lta polymorphs) to less than 0.6 ˚ A(ZIF-FF). Finally, we examine the predicted densities of the various polymorphs as a proxy for the accuracy of the lattice constants. As shown in Fig. 7, the mean average (unsigned) density error for ZIFF-FF is 1.0 %, with little systematic error (mean signed error of only -0.4 %), and a maximum error of less than 5 % for all structures. This is a substantial decrease 15 ACS Paragon Plus Environment


Figure 7: Signed density error of the varous Zn(MeIm)2 polymorphs as compared to the DFT references. from the original GAFF values of 3.6% and -2.4 %, respectively. In addition, GAFF exhibits several outlier structures with errors approaching 15 %, while ZIF-FF yields a far more balanced performance.

Force Field Validation Structural Properties We further validated ZIF-FF by examining its performance via its prediction of a number of additional properties, none of which were directly included in the parameterization process. We first examined the various local environments around the Zn atoms found in the reference structures in terms of bond lengths (Zn-N) and bond angles (N-Zn-N, Zn-N-C1, Zn-N-C2) at 0K. Note that the 25 different ZIF crystal polymorphs represent different topologies, and thus each structure exhibits different equilibrium bond lengths and angles. We compared the distribution of these values against our DFT reference calculations. As shown in Fig. 8, ZIF-FF provides a distribution of Zn-N bond lengths much closer to DFT distribution as compared to GAFF. It is particularly notable that characteristic bond


Page 16 of 40

Page 17 of 40 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


Figure 8: Comparison of the distribution of bond (Zn-N) and angle (N-Zn-N, Zn-N-C1, Zn-N-C2) distributions as obtained by DFT (green), GAFF (blue), and ZIF-FF (red). The distributions are normalized.



length shifts to larger values by approximately 0.05 ˚ A, and amount which far exceeds the change in the equilibrium value of the associated bond stretching parameter (r0 ). This seems to imply that the original GAFF force field exhibited significant artificial bond “strain” that has been relieved in the optimized force field. In contrast, the N-Zn-N angle distributions of all three methods (DFT, GAFF and ZIF-FF) largely overlap, indicating that both ZIF-FF and GAFF can effective capture this crucial structural detail. Overall, these results suggest that ZIF-FF is capable of reproducing the ensemble of local structures around the Zn atoms found in the various polymorphs.

Phonon Spectrum The phonon spectrum plays an important role in the finite-temperature structural and mechanical properties of ZIFs. 59 We thus also compare the phonon spectra of ZIF-8 calculated by force field and DFT. For both force field and DFT, dynamical matrices at Γ-point were built from the supercell force constants following optimization of the structure, and then phonon vibrational frequencies were calculated as the square root of the eigenvalues of the mass weighted dynamical matrix. The partial density state (PDOS), which decomposes the DOS into various elemental contributions, was also calculated to provide further insight into the origins of various spectral contributions. The phonon DOS and PDOS for ZIF-8 (SOD topology) were calculated using Phonopy and plotted in Fig. 9. Note that all calculations utilized the harmonic approximation; although this approximation is certainly questionable, particularly at very low phonon frequency, its use consistently across all cases still allows for comparison across various force field / DFT methods. Compared to GAFF, good agreement is observed between DFT and ZIF-FF, suggesting that the phonon vibrational properties of ZIF-8 are well reproduced by ZIF-FF. In addition, the PDOS spectrum of Zn is notable improved, particularly at very low frequencies, confirming that ZIF-FF does enhance the Zn-related parameters. 18 ACS Paragon Plus Environment

Page 18 of 40

Page 19 of 40 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


Figure 9: Density of states (DOS) and Zn-projected partial density of states (PDOS) calculated by DFT, GAFF, and ZIF-FF. Mechanical Properties We also validated the predicted mechanism properties (elastic constants, shear modulus, bulk modulus) of ZIF-8 against literature values. For ZIF-8, a cubic crystal, three independent elastic constants (C11 , C12 , and C44 ) were computed and listed in Table 2. These elastic constants satisfy the fundamental elastic stability criteria for a cubic crystal: C11 > |C12 |, C11 + 2C12 > 0, and C44 > 0. 60 The shear modulus (G) represents the framework’s rigidity against structural distortion imposed by an external shear forces. In contrast, the bulk modulus (K), the inverse of compressibility, represents the resistance of the structure against volumetric strains under a hydrostatic pressure. In order to directly connect with corresponding measurements on polycrystalline samples, we utilize an isotropic assumption. In this case, the shear modulus and bulk modulus were computed via Voigt–Reuss–Hill (VRH) averages.(See Eq. (7) and Eq. (8), for cubic systems.) The Voigt scheme assumes a uniform strain, while the Reuss averaging correspond to a uniform stress; the Voigt–Reuss–Hill values are the average of these two value. 61 The VRH values are also listed in Table 2.



Page 20 of 40

Gmax = (C11 − C12 )/2, Gmin = C44 GV = (C11 − C12 + 3C44 )/5, GR = 5/[4(S11 − S12 + 3S44 ]

(7)

GV RH = (GV + GR )/2 where Sij is the compliance matrix, i.e. the inversion of the elasticity tensor Cij .

KV RH = KV = KR = (C11 + 2C12 )/3

(8)

From Table 2, we find that elastic constants, shear modulus and bulk modulus computed by ZIF-FF match the experimental data quite well. Consistent with conclusion of Tan et al., 62 ZIF-8 exhibits an exceptionally low shear modulus, implying a propensity for shearinduced amorphization. Alternatively, the bulk modulus, K can also be computed directly using its definition, through Eq. (9).

K = −V0

dP dV

(9)

where V and V0 are the crystal volumes at actual pressure P and reference pressure P0 , respectively. Using a third-order Birch-Murnaghan equation of state (EOS), the calculated K was 7.07 GPa, close to experimental value 6.52(35) GPa based on in-situ high pressure Xray diffraction experiment using a identical analysis method. 63 Both the experimental and simulated values calculated using this EOS method are smaller than the respective VRH values, which may be due to the higher pressure perturbation applied in EOS method and isotropic assumption used in VRH scheme.


Page 21 of 40 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


Table 2: Comparison of elastic constants(Cij ), shear modulus(G), and bulk modulus(K) of ZIF-8 by different techniques in units of GPa Technique Brillouin Scattering(295K) 62 Zhang2013 27 FF(300K) 64 ZIF-FF(0K) ZIF-FF(295K)

C11 C12 C44 Gmin 9.523 6.865 0.967 0.97 11.21 7.49 2.75 2.75 9.33 6.45 1.37 1.37 9.06 5.92 1.20 1.20

Gmax 1.33 – 1.44 1.57

GV RH 1.10 – 1.40 1.34

KV RH 7.75 8.00 7.41 6.96

Diffusion We also validated ZIF-FF by simulating methane diffusion in ZIF-8. We utilized dynamically corrected transition state theory (dcTST), 65,66 where the slow diffusion of adsorbates is modeled as a cage-to-cage hopping process, in conjunction with umbrella sampling (US) to estimate the associated free energy hopping barrier. The computational details are similar to the work by Verploegh et al., 65 but instead utilizing ZIF-FF for the framework force field. The free energy curve for methane diffusion at 308.15 K in both flexible and rigid ZIF-8 are shown in Fig. 10. From Fig. 10, we can see that ZIF-FF provides a free energy barrier of 22 kJ/mol, slightly lower than the 25 kJ/mol generated by Zhang FF. 27 For comparison, we also computed the diffusion free energy barrier for the ZIF-8 rigid framework (holding the framework atoms fixed at their equilibrium positions). The associated free energy barrier for the rigid case is much larger, 32 kJ/mol, demonstrating that the framework flexibility imbued by ZIF-FF can lead to important consequences for diffusion.

Amorphization We also validated our flexible ZIF-FF force field via application to the pressure-induced amorphization of ZIF-8. Amorphous MOFs typically share similar local structure with their crystalline counterparts, while lacking the long-range order characteristic of crystalline solids. 67 Amorphous MOFs have gained significant recent attention due to their potential 21 ACS Paragon Plus Environment


Figure 10: Free energy curve for methane diffusion at 308.15 K in both flexible and rigid ZIF-8 framework.


Page 22 of 40

Page 23 of 40 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


applications, including those that may exploit a crystalline-amorphous transition, such as reversible gas storage via pressure-induced amorphization, 68 irreversible trapping of harmful gas species via framework collapse, 69 or MOF-based drug delivery. 70 It is known that ZIF-8 can be amorphized via pressure 63 or ball-milling. 71 Those experiments seem to suggest that local structure is preserved during the amorphization process. We thus assume that no bond breaking takes place during amorphization, and that as such we can utilize the non-reactive ZIF-FF to simulate the pressure-induced amorphization process. In fact, transient bond breaking may occur, which might lower the amorphization threshold. Ortiz et al. have shown that the pressure-induced amorphization of ZIF-8 arises due to shear mode softening; 64 a similar low shear modulus is also observed using ZIF-FF. To simulate the amorphization process, constant stress and constant temperature (NsT) simulations were performed using the LAMMPS software package. The temperature was fixed at 300 K using Nose-Hoover thermostat, with a damping coefficient of 100 fs. The pressure was controlled via a Nose-Hoover barostat (damping time of 1000 fs), with the pressure steadily increased at a rate of 100 Pa/fs from ambient pressure to 1.5 GPa over the course of the simulation. Prior to the NsT simulation, a constant pressure and constant temperature (NPT) simluation at 300 K and ambient pressure was run for 1 ns to obtain the equilibrium structure. During the NsT simulation, the density trajectory, as well as representative snapshots of the crystalline or amorphous ZIF-8 structures, are shown in Fig. 11. During the pressure ramp, two sharp density increases at around 0.81 GPa and 1.1 GPa were observed, indicating two phase transitions, first from a crystalline to an intermediate amorphous structure, and later to a high-density amorphous structure. This assignment is supported by analysis of the associated simulation snapshots in Fig. 11b and via analysis of the structure factor (shown in Fig. 12). During the reverse pressure ramp, analogous transitions take place at lower pressures (around 0.26 GPa and 0.20 GPa respectively), showing 23 ACS Paragon Plus Environment


significant hysteresis and only transient presence of the intermediate amorphous phase. Using the two legs of the hysteresis curve, we estimate an amorphization pressure between 0.20-0.81 GPa, consistent with the experimental value 0.34 GPa observed during pressureinduced amorphization by Chapman et al.. 63 Cao et al. reported similar abrupt increase in density for the amorphization of ZIF-8 during ball-milling. 71 Note that the transition pressure observed during the increasing leg (0.81 GPa) is larger than that observed from experiment. This discrepancy may be due to potential transient bond breaking in experimental conditions, as well as kinetic factors induced by the rapid pressure ramp utilized in the simulation (and consistent with the observed hysteresis). Additionally, the presence or formation of defects in ZIF-8 would lead to a less strained framework, perhaps further lowering the amorphization pressure under experimental conditions. The structure factor, S(q), was computed to further characterize the amorphous structure. (See supporting information for complete computational details; the low q range of S(q) is computed using the method shown in the supporting information in the work. 72 ) Because the density of ZIF-8 continues to increases with pressure even after the amorphous transition, we compare the structure factor at the same density (1.52 g/cm3 ) as observed during ball-milling, 71 using the structures derived from the decreasing pressure leg of the simulation; the respective structure factors are compared in Fig. 12. The structure factor of intermediate (transient) amorphous ZIF-8 with a density of 1.40 g/cm3 is also shown (in red dash line) in Fig. 12, and is clearly distinct from both the high-density amorphous phase and the original crystalline structure. We find that the simulated S(q) for the high-density amorphous ZIF-8 matches well with experimental data, demonstrating that ZIF-FF provides a faithful reproduction of both the crystalline and amorphous phases.


Page 24 of 40

Page 25 of 40 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


3

4

2

1

5

(a)

ZIF-8

1

2

3

4

5

(b)

Figure 11: Density of ZIF-8 as a function of pressuring during pressure-induced amorphization, modeled using ZIF-FF. The top panel shows density change with pressure, while the bottom panel shows representative structures at the marked time points. 25 ACS Paragon Plus Environment


Figure 12: Structure factor of amorphous ZIF-8 with density of 1.52 g/cm3 and intermediate amorphous ZIF-8 with density of 1.40 g/cm3 obtained from constant stress simulation using ZIF-FF. Experimental data are reproduced from ref. 71 via ball-milling for 300 min. 71 The dashed vertical lines indicate the peak positions of structure factor from experimental data. The curves for simulated amorphous ZIF-8 and intermediate amorphous ZIF-8 are shifted in y-axis in order to differentiate from experimental one.


Page 26 of 40

Page 27 of 40 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


Transferability ZIF Surface Structures Our ZIF-FF parameters were fitted to bulk structures and strain energies. Here, we examine whether ZIF-FF yields transferable predictions for the structure and energetics of ZIF surfaces and interfaces. Comparison of the relative stability of surfaces/interfaces is generally made on the basis of surface energy (defined in Eq. (10)), which accounts both for the energy to break bonds to form the interface, as well as the energetics of the surface relaxation. Since ZIF-FF is a non-reactive force field, we focus exclusive on the latter by calculating the relative surface energies ∆γ of six different Zn(MeIm)2 slabs with same number of broken Zn-N bonds within the unit cell. The configuration of these six slabs (cag(010), crb(110), dft(001), irl(010), lon(010), pcl(010)) are shown as in Fig. S4 in the SI.

γ=

Eslab − N · Ebulk A

(10)

We choose lon(010) as our reference and calculate all surface energies, ∆γi , relative to this reference value; see Eq. (11).

∆γi = γi − γlon(010)

(11)

There is significant ambiguity regarding the stable termination of ZIF interfaces. Moh et al. performed step height analysis of AFM deflection images and indicated a surface termination motif of [Zn(MeIm)3]– under solvothermal growth condition. 73 Tian et al. utilized X-ray photoelectron spectroscopy (XPS) to identify the presence of termination with carbonates, water/hydroxides, uncoordinated amines and methylimidazole of ZIF-8 thin film. High-resolution TEM study by Zhu et al. revealed a ’armchair’-type rather than ’zigzag’27 ACS Paragon Plus Environment


type termination with 2-methylimidazole (HMeIm) capping ZIF-8 (110) surfaces. 75 Given the various possible terminations (likely reflecting differences in synthetic / environmental conditions), we examined two different surface terminations, [Zn(MeIm)3]– and [Zn(MeIm)4]2–, using ZIF-FF. In all cases, a 15 ˚ Avacuum gap was included and, with the bottom two layers of the ZIF slab fixed to the bulk structure. As shown in Fig. 13, ZIF-FF faithfully reproduces the relative surface strain energies of the six slabs in their varying terminations as compared to the DFT reference, while the GAFF results exhibit substantially poorer agreement. The success of ZIF-FF when applied to these very different (interfacial) conditions bodes well for the transferability of the model and its applications to modeling ZIF interfaces.

Functionalized ZIFs We also explore the transferability of the force field model to a variety of other ZIFs that differ in their linker functionalization (shown in Fig. 14), as well as in their topology. For SOD polymorphs, three different linkers were computed and compared. For simplicity, we employed standard GAFF parameters for the different linker functional groups, while using ZIF-FF parameters for the other terms. We then compared the predicted density of the ZIFs with DFT benchmark values. As shown in Fig. 14, the FF-computed density show good agreement with respect to both DFT and experiments, supporting the contention that ZIF-FF is not only transferable not only across topologies, but also across functionalization. Note that the slight discrepancy of densities between ZIF-FF and experiment in ZIF-72, which is slightly larger compared to other systems, is likely due to relatively inaccurate parameters of Cl from GAFF.


Page 28 of 40

Page 29 of 40 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


Figure 13: The comparison of ZIF-FF and DFT in the relative surface energy of six slabs, with both [Zn(MeIm)3]– and [Zn(MeIm)4]2– termination.



Figure 14: ZIF-FF predicted densities of various functionalized ZIFs, as compared to both DFT and experiment

Conclusions We have developed a simple and effective ab-initio intra-molecular force field for a variety of ZIFs, ZIF-FF, based upon periodic DFT reference data. Starting with standard GAFF parameters and functional forms, the Zn-containing bonded terms were optimized to reproduce DFT-calculated structure (atom positions, lattice constants) and energetic (strain energy) data. The resulting force fields not only reproduces the relative stability of wide range of ZIF polymorphs, but also exhibits excellent agreement with a wide array of other properties, including the experimentally-observed crystalline-amorphous transformation of ZIF-8 and ZIF surface structures. We anticipate that that ZIF-FF will facilitate future simulation studies of a wide array of ZIF bulk and interfacial processes.


Page 30 of 40

Page 31 of 40 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


Supporting Information Available Table S1 and Figure S1 comparing experimental and DFT densities of ZIFs with various topologies; Table S2 providing initial Zn-related parameter ranges for GA; Table S3 giving complete force field parameters of ZIF-FF; Table S4 giving Zn-related parameters used in GA-GAFF; Figure S2 comparing FF and DFT strain energy of 25 ZIF polymorphs using FF by Gabrieli et al. 30 and ZIF-FF; Table S5 providing non-bonded and Zn-related bonded force field parameters of SAPT-FF; Figure S3 comparing FF and DFT strain energy of 25 ZIF polymorphs using GAFF, ZIF-FF and SAPT-FF; details of structure factor calculation; Figure S4 showing six Zn(MeIm)2 slabs with different topologies; geometry files of DFToptimized ZIF reference systems with 25 polymorphs; example of LAMMPS input files using ZIF-FF in ZIF-8 (3 × 3 × 3) supercell. • Filename: supplementary.pdf • Filename: zif meth structures.zip • Filename: lammps input.zip This material is available free of charge via the Internet at http://pubs.acs.org/.



Acknowledgement This work was supported by Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy, under award DESC0014059. Computational resources were provided by the UW- Madison Center for High Throughput Computing (CHTC) in the Department of Computer Sciences. The CHTC is supported by UW-Madison, the Advanced Computing Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin Institutes for Discovery, and the National Science Foundation, and is an active member of the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy’s Office of Science. J.R.S is a Camille Dreyfus Teacher-Scholar. We thank Prof. Thomas Bennett for providing experimental raw data for structure factor of amorphous ZIF-8 by ball-milling for 300 min. We also thank Prof. McDaniel for providing code to calculate structure factor in low q region.


Page 32 of 40

Page 33 of 40 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


References (1) Furukawa, H.; Cordova, K. E.; O’Keeffe, M.; Yaghi, O. M. The Chemistry and Applications of Metal-Organic Frameworks. Science 2013, 341, 1230444. (2) Murray, L. J.; Dinc˘a, M.; Long, J. R. Hydrogen Storage in Metal–Organic Frameworks. Chem. Soc. Rev. 2009, 38, 1294–1314. (3) Chaemchuen, S.; Kabir, N. A.; Zhou, K.; Verpoort, F. Metal–Organic Frameworks for Upgrading Biogas via CO2 Adsorption to Biogas Green Energy. Chem. Soc. Rev. 2013, 42, 9304–9332. (4) Tanh Jeazet, H. B.; Staudt, C.; Janiak, C. Metal–Organic Frameworks in Mixed-Matrix Membranes for Gas Separation. Dalton Trans. 2012, 41, 14003–14027. (5) Qiu, S.; Xue, M.; Zhu, G. Metal–Organic Framework Membranes: from Synthesis to Separation Application. Chem. Soc. Rev. 2014, 43, 6116–6140. (6) Liu, J.; Chen, L.; Cui, H.; Zhang, J.; Zhang, L.; Su, C.-Y. Applications of Metal–Organic Frameworks in Heterogeneous Supramolecular Catalysis. Chem. Soc. Rev. 2014, 43, 6011–6061. (7) Farrusseng, D.; Aguado, S.; Pinel, C. Metal–Organic Frameworks: Opportunities for Catalysis. Angew. Chem. Int. Ed. 2009, 48, 7502–7513. (8) Rojas, S.; Wheatley, P. S.; Quartapelle-Procopio, E.; Gil, B.; Marszalek, B.; Morris, R. E.; Barea, E. Metal–Organic Frameworks as Potential Multi-Carriers of Drugs. CrystEngComm 2013, 15, 9364–9367. (9) Horcajada, P.; Gref, R.; Baati, T.; Allan, P. K.; Maurin, G.; Couvreur, P.; Férey, G.; Morris, R. E.; Serre, C. Metal–Organic Frameworks in Biomedicine. Chem. Rev. 2012, 112, 1232–1268. (10) Falcaro, P.; Ricco, R.; Doherty, C. M.; Liang, K.; Hill, A. J.; Styles, M. J. MOF Positioning Technology and Device Fabrication. Chem. Soc. Rev. 2014, 43, 5513–5560. (11) Park, K. S.; Ni, Z.; Côté, A. P.; Choi, J. Y.; Huang, R.; Uribe-Romo, F. J.; Chae, H. K.; O’Keeffe, M.; Yaghi, O. M. Exceptional Chemical and Thermal Stability of Zeolitic Imidazolate Frameworks. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 10186–91. (12) 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, 58–67. ˇ (13) Katsenis, A. D.; Puˇskarić, A.; Strukil, V.; Mottillo, C.; Julien, P. A.; Uˇzarević, K.; Pham, M.-H.; Do, T.-O.; Kimber, S. A. J.; Lazić, P. et al. In Situ X-ray Diffraction Monitoring of a Mechanochemical Reaction Reveals a Unique Topology Metal-Organic Framework. Nat. Commun. 2015, 6 . 33 ACS Paragon Plus Environment


(14) Van Vleet, M. J.; Weng, T.; Li, X.; Schmidt, J. R. In Situ, Time-Resolved, and Mechanistic Studies of Metal-Organic Framework Nucleation and Growth. Chem. Rev. 2018, 118, 3681–3721. (15) D¨ uren, T.; Bae, Y.-S.; Snurr, R. Q. Using Molecular Simulation to Characterise Metal–Organic Frameworks for Adsorption Applications. Chem. Soc. Rev. 2009, 38, 1237. (16) Pérez-Pellitero, J.; Amrouche, H.; Siperstein, F. R.; Pirngruber, G.; Nieto-Draghi, C.; Chaplais, G.; Simon-Masseron, A.; Bazer-Bachi, D.; Peralta, D.; Bats, N. Adsorption of CO2, CH4, and N2 on Zeolitic Imidazolate Frameworks: Experiments and Simulations. Chem. Eur. J. 2010, 16, 1560–1571. (17) McDaniel, J. G.; Yu, K.; Schmidt, J. R. Ab Initio, Physically Motivated Force Fields for CO2 Adsorption in Zeolitic Imidazolate Frameworks. J. Phys. Chem. C 2011, 116, 1892–1903. (18) McDaniel, J. G.; Schmidt, J. R. Robust, Transferable, and Physically-Motivated Force Fields for Gas Adsorption in Functionalized Zeolitic Imidazolate Frameworks. J. Phys. Chem. C 2012, 116, 14031–14039. (19) McDaniel, J. G.; Schmidt, J. R. Physically-Motivated Force Fields from SymmetryAdapted Perturbation Theory. J. Phys. Chem. A 2013, 117, 2053–2066. (20) Rappé, A. K.; Casewit, C. J.; Colwell, K. S.; 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, 10024–10035. (21) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: A Generic Force Field for Molecular Simulations. Journal of Physical Chemistry 1990, 94, 8897–8909. (22) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. An All Atom Force Field for Simulations of Proteins and Nucleic Acids. J. Comput. Chem. 1986, 7, 230–252. (23) Hertäg, L.; Bux, H.; Caro, J.; Chmelik, C.; Remsungnen, T.; Knauth, M.; Fritzsche, S. Diffusion of CH4 and H2 in ZIF-8. J. Membr. Sci. 2011, 377, 36–41. (24) Zheng, B.; Sant, M.; Demontis, P.; Suffritti, G. B.; Journal, T.; Chemistry, P.; Zheng, B.; Sant, M.; Demontis, P.; Suffritti, G. B. Force Field for Molecular Dynamics Computations in Flexible ZIF-8 Framework. J. Phys. Chem. C 2012, 116, 933–938. (25) Gee, J. A.; Sholl, D. S. Characterization of the Thermodynamic Stability of Solvated Metal-Organic Framework Polymorphs Using Molecular Simulations. J. Phys. Chem. C 2013, 117, 20636–20642. (26) Hu, Z.; Zhang, L.; Jiang, J. Development of a Force Field for Zeolitic Imidazolate Framework-8 with Structural Flexibility. J. Chem. Phys. 2012, 136 . 34 ACS Paragon Plus Environment

Page 34 of 40

Page 35 of 40 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


(27) Zhang, L.; Hu, Z.; Jiang, J. Sorption-Induced Structural Transition of Zeolitic Imidazolate Framework-8: A Hybrid Molecular Simulation Study. J. Am. Chem. Soc. 2013, 135, 3722–3728. (28) Zhang, L.; Wu, G.; Jiang, J. Adsorption and Diffusion of CO2 and CH4 in Zeolitic Imidazolate Framework-8: Effect of Structural Flexibility. J. Phys. Chem. C 2014, 8788–8794. (29) Krokidas, P.; Castier, M.; Moncho, S.; Brothers, E.; Economou, I. G. Molecular Simulation Studies of the Diffusion of Methane, Ethane, Propane, and Propylene in ZIF-8. J. Phys. Chem. C 2015, 119, 27028–27037. (30) Gabrieli, A.; Sant, M.; Demontis, P.; Suffritti, G. B. Fast and Efficient Optimization of Molecular Dynamics Force Fields for Microporous Materials: Bonded Interactions via Force Matching. Micropor. Mesopor. Mat. 2014, 197, 339–347. (31) Tafipolsky, M.; Schmid, R. Systematic First Principles Parameterization of Force Fields for Metal-Organic Frameworks using a Genetic Algorithm Approach. J. Phys. Chem. B 2009, 113, 1341–1352. (32) Bureekaew, S.; Amirjalayer, S.; Tafipolsky, M.; Spickermann, C.; Roy, T. K.; Schmid, R. MOF-FF - A Flexible First-Principles Derived Force Field for Metal-Organic Frameworks. Phys. Status Solidi B 2013, 250, 1128–1141. (33) Schmid, R.; Tafipolsky, M. An Accurate Force Field Model for the Strain Energy Analysis of the Covalent Organic Framework COF-102. J. Am. Chem. Soc. 2008, 130, 12600–12601. (34) Amirjalayer, S.; Snurr, R. Q.; Schmid, R. Prediction of Structure and Properties of Boron-Based Covalent Organic Frameworks by a First-Principles Derived Force Field. J. Phys. Chem. C 2012, 116, 4921–4929. (35) D¨ urholt, J. P.; Galvelis, R.; Schmid, R. Coarse Graining of Force Fields for Metal–Organic Frameworks. Dalton Trans. 2016, 45, 4370–4379. (36) 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, 10152–10154. (37) Lu, G.; Hupp, J. T. Metal Organic Franeworks as Sensors A ZIF-8 Based Fabry Perot Device as a Selective Sensor for Chemical Vapours and Gases. J. Am. Chem. Soc. 2010, 132, 7832–7833. (38) Tran, U. P. N.; Le, K. K. A.; Phan, N. T. S. Expanding Applications of Metal-Organic Frameworks: Zeolite Imidazolate Framework ZIF-8 as an Efficient Heterogeneous Catalyst for the Knoevenagel Reaction. ACS Catal. 2011, 1, 120–127. 35 ACS Paragon Plus Environment


(39) Huang, X.; Jiepeng, Z.; Xiaoming, C. [Zn(bim)2] · (H2O)1 · 67: A Metal-Organic OpenFramework with Sodalite Topology. Chin. Sci. Bull. 2003, 48, 1531–1534. (40) Shi, Q.; Chen, Z.; Song, Z.; Li, J.; Dong, J. Synthesis of ZIF-8 and ZIF-67 by SteamAssisted Conversion and an Investigation of Their Tribological Behaviors. Angew. Chem. Int. Ed. 2011, 50, 672–675. (41) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558–561. (42) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the LiquidMetal–Amorphous-Semiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251–14269. (43) Kresse, G.; Furthm¨ uller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186. (44) Kresse, G.; Furthm¨ uller, J. Efficiency of Ab-initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mat. Sci. 1996, 6, 15–50. (45) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979. (46) Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B 1999, 59, 1758–1775. (47) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (48) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (49) Zheng, B.; Pan, Y.; Lai, Z.; Huang, K. W. Molecular Dynamics Simulations on Gate Opening in ZIF-8: Identification of Factors for Ethane and Propane Separation. Langmuir 2013, 29, 8865–8872. (50) Gee, J. A.; Chung, J.; Nair, S.; Sholl, D. S. Adsorption and Diffusion of Small Alcohols in Zeolitic Imidazolate Frameworks ZIF-8 and ZIF-90. J. Phys. Chem. C 2013, 1–17. (51) Wu, X.; Huang, J.; Cai, W.; Jaroniec, M. Force Field for ZIF-8 Flexible Frameworks: Atomistic Simulation of Adsorption, Diffusion of Pure Gases as CH4, H2, CO2 and N2. RSC Adv. 2014, 4, 16503–16511. (52) Wang, J.; Wolf, R.; Caldwell, J. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 56531 . (53) Hockney, R.; Eastwood, J. Computer Simulation Using Particles; McGraw-Hill,New York, 1981. 36 ACS Paragon Plus Environment

Page 36 of 40

Page 37 of 40 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


(54) Manz, T. A.; Sholl, D. S. Improved Atoms-in-Molecule Charge Partitioning Functional for Simultaneously Reproducing the Electrostatic Potential and Chemical States in Periodic and Nonperiodic Materials. J. Chem. Theo. Comput. 2012, 8, 2844–2867. (55) Manz, T. A.; Limas, N. G. Introducing DDEC6 Atomic Population Analysis: Part 1. Charge Partitioning Theory and Methodology. RSC Adv. 2016, 6, 47771–47801. (56) Goldberg, D. E. Genetic Algorithms in Search, Optimization, and Machine Learning; Artificial Intelligence; Addison-Wesley Publishing Company, 1989. (57) Gagn, C. DEAP : Evolutionary Algorithms Made Easy. Journal of Machine Learning Research 2012, 13, 2171–2175. (58) Plimpton, S. Fast Parallel Algorithms for Short–Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (59) Bristow, J. K.; Skelton, J. M.; Svane, K. L.; Walsh, A.; Gale, J. D. A General Force Field for Accurate Phonon Properties of Metal-Organic Frameworks. Phys. Chem. Chem. Phys. 2016, 18, 29316–29329. (60) Nye., J. F. Physical Properties of Crystals. Cryst. Res. Tech. 1986, 21, 1508–1508. (61) Ryder, M. R.; Tan, J.-C. Explaining the Mechanical Mechanisms of Zeolitic Metal–Organic Frameworks: Revealing Auxeticity and Anomalous Elasticity. Dalton Trans. 2016, 45, 4154–4161. (62) 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, 1–6. (63) Chapman, K. W.; Halder, G. J.; Chupas, P. J. Pressure-Induced Amorphization and Porosity Modification in a Metal-Organic Framework. J. Am. Chem. Soc. 2009, 131, 17546–17547. (64) Ortiz, A. U.; Boutin, A.; Fuchs, A. H.; Coudert, F.-X. X. Investigating the PressureInduced Amorphization of Zeolitic Imidazolate Framework ZIF-8: Mechanical Instability Due to Shear Mode Softening. J. Phys. Chem. Lett. 2013, 4, 1861–1865. (65) Verploegh, R. J.; Nair, S.; Sholl, D. S. Temperature and Loading-Dependent Diffusion of Light Hydrocarbons in ZIF-8 as Predicted Through Fully Flexible Molecular Simulations. J. Am. Chem. Soc. 2015, 137, 15760–15771. (66) Dubbeldam, D.; Calero, S.; Ellis, D. E.; Snurr, R. Q. RASPA: Molecular Simulation Software for Adsorption and Diffusion in Flexible Nanoporous Materials. Mol. Simul. 2016, 42, 81–101. 37 ACS Paragon Plus Environment


(67) Bennett, T. D.; Cheetham, A. K. Amorphous Metal–Organic Frameworks. Acc. Chem. Res. 2014, 47, 1555–1562. (68) Bennett, T. D.; Simoncic, P.; Moggach, S. A.; Gozzo, F.; Macchi, P.; Keen, D. A.; Tan, J.-C.; Cheetham, A. K. Reversible Pressure-Induced Amorphization of a Zeolitic Imidazolate Framework (ZIF-4). Chem. Commun. 2011, 47, 7983. (69) Bennett, T. D.; Saines, P. J.; Keen, D. A.; Tan, J. C.; Cheetham, A. K. Ball-MillingInduced Amorphization of Zeolitic Imidazolate Frameworks (ZIFs) for the Irreversible Trapping of Iodine. Chem. Eur. J. 2013, 19, 7049–7055. (70) Horcajada, P.; Serre, C.; Vallet-Reg´ı, M.; Sebban, M.; Taulelle, F.; Férey, G. MetalOrganic Frameworks as Efficient Materials for Drug Delivery. Angew. Chem. Int. Ed. 2006, 45, 5974–5978. (71) Cao, S.; Bennett, T. D.; Keen, D. A.; Goodwin, A. L.; Cheetham, A. K. Amorphization of the Prototypical Zeolitic Imidazolate Framework ZIF-8 by Ball-Milling. Chem. Commun. 2012, 48, 7805–7807. (72) Mantha, S.; McDaniel, J. G.; Perroni, D. V.; Mahanthappa, M. K.; Yethiraj, A. Electrostatic Interactions Govern “Odd/Even” Effects in Water-Induced Gemini Surfactant Self-Assembly. J. Phys. Chem. B 2017, 121, 565–576. (73) Moh, P. Y.; Cubillas, P.; Anderson, M. W.; Attfield, M. P. Revelation of the Molecular Assembly of the Nanoporous Metal Organic Framework ZIF-8. J. Am. Chem. Soc. 2011, 133, 13304–13307. (74) Tian, F.; Cerro, A. M.; Mosier, A. M.; Wayment-Steele, H. K.; Shine, R. S.; Park, A.; Webster, E. R.; Johnson, L. E.; Johal, M. S.; Benz, L. Surface and Stability Characterization of A Nanoporous ZIF-8 Thin Film. J. Phys. Chem. C 2014, 118, 14449–14456. (75) Zhu, Y.; Ciston, J.; Zheng, B.; Miao, X.; Czarnik, C.; Pan, Y.; Sougrat, R.; Lai, Z.; Hsiung, C. E.; Yao, K. et al. Unravelling surface and interfacial structures of a metalorganic framework by transmission electron microscopy. Nat. Mat. 2017, 16, 532–536.


Page 38 of 40

Page 39 of 40 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


Graphical TOC Entry ∆"# =

"(&# ) "(&)*+ ) − &# &)*+

ZIF-FF

ZIF-8 (SOD)

Amorphous ZIF-8



J. R. Schmidt is a Professor of Chemistry at the University of Wisconsin-Madison and a member of the Theoretical Chemistry Institute. He earned his B.S. from Hope College in 2001, and a Ph.D. in Physical Chemistry from the University of Wisconsin (with Jim Skinner) in 2006, focusing on simulations of the dynamics and non-linear spectroscopy of aqueous solutions. After postdoctoral work at Yale (with John Tully), he returned to UW in 2008. Prof. Schmidt's research interests include developing accurate, first-principles force fields, nano-porous materials, and computational heterogeneous catalysis. Photo note: I have used this photo with ACS before, and they already have a copyright release on file.

ACS Paragon Plus Environment

Page 40 of 40

Flexible and Transferable ab Initio Force Field for Zeolitic Imidazolate

Recommend Documents