Reactive Coarse Grained Force Field for Metal-Organic Frameworks applied to Modeling ZIF-8 Self-Assembly

TL;DR

This work tackles the challenge of modeling MOF self-assembly under experimentally relevant conditions by developing nb-CG-ZIF-FF, a reactive coarse-grained force field learned via MS-CG from atomistic ZIF-8 benchmarks. The model encodes Zn connectivity without explicit bonds and reproduces both crystalline structure and nucleation dynamics in DMSO, achieving roughly two orders of magnitude speed-up over all-atom simulations. It captures key features of the amorphous intermediate and ring formation during self-assembly, while implicitly learning Zn tetrahedral coordination. The method offers a general, data-driven route to study MOF formation, decomposition, and defect dynamics across other MOFs and synthesis conditions, enabling systematic exploration of concentration and composition effects at larger scales than accessible atomistically.

Abstract

Decoding the self-assembly mechanism of metal-organic frameworks is a crucial step in reducing trial-and-error tests in their synthesis protocols. Atomistic simulations have proven essential in revealing molecular-level features of MOF nucleation, but they still exhibit limitations in the simulation setups due to size constraints (inability of reaching realistic concentrations or exploring non-stoichiometric metal:ligand ratios). In this contribution, we develop a methodology to derive reactive coarse grained force fields based on multiscale coarse graining methods. We apply our novel methodology to the case of the archetypal zeolitic-imidazolate framework ZIF-8. Our coarse grained force field, which we call nb-CG-ZIF-FF, does not contain any explicit connectivity information, but learns the tetrahedral Zn-connectivity from many body correlations within an atomistic benchmark. nb-CG-ZIF-FF quantitatively reproduces the features of bulk, crystalline ZIF-8 as well as the structural evolution of pre-nucleation species in terms of Zn n-fold coordination populations from the atomistic benchmark. While the range of rings that are formed along the synthesis process are well captured by nb-CG-ZIF-FF, the model cannot exactly reproduce ring populations. Our reactive CG force field fitting approach can be applied to any MOF, opening new research avenues in modeling MOF formation, decomposition, defect dynamics and phase transition processes.

Figures (9)

• Figure 1: Mapping of the coarse grained reactive force field developed in this work: each Zn$^{2+}$ ion together with its four associated dummy atoms[small dark red spheres] is represented as a single bead [light red sphere]; each 2-methylimidazolate (mIm$^-$) ligand along with its two dummy atoms [small light blue spheres] is modeled as an individual bead [dark blue spheres]; and each DMSO molecule is treated as a single bead [gray sphere].
• Figure 2: Schematic representation of the CG methodology employed in this work. AA models of crystalline ZIF-8 and the Zn$^{2+}$ and mIm$^{-}$ ions solution in DMSO are converted to CG models and a reactive CG force field is obtained through MS-CG and parameter optimization. Then, a pressure matching stage is performed. The CG force field is validated by comparing observables computed from CG and AA simulations.
• Figure 3: Radial distribution functions g(r) for Zn-Zn (blue), Zn-ligand (green), and ligand-ligand (red) interactions in crystalline ZIF-8, obtained from the nb-CG-ZIF-FF model (dashed lines) and from coarsening the AA trajectory (solid lines).
• Figure 4: Snapshots of a molecular dynamics simulation of ZIF-8 self-assembly made with nb-CG-ZIF-FF. (a) Sequential snapshots (t = 0–200 ps) showing the progression from dissolved metal and ligand beads through the formation of linear chains, branched oligomers, and interconnected networks, up to an amorphous ZIF-8 structure. (b) View of a representative crystalline ZIF-8 structure at the CG mapping chosen in this work.
• Figure 5: Zn-ligand coordination number evolution: (a) NVT simulation (b) NVT simulation carried out with the fix deform command, where the system volume is allowed to gradually relax to match the final volume from the reference AA NPT simulation.