Johannes P Dürholt

Guillaume Fraux

François-Xavier Coudert
email: fx.coudert@chimieparistech.psl.eu

Rochus Schmid
email: rochus.schmid@rub.de

Ab Initio Derived Force Fields for Zeolitic Imidazolate Frameworks: MOF-FF for ZIFs

come   L'archive ouverte pluridisciplinaire

Introduction

Zeolitic imidazolate frameworks (ZIFs) are a subclass of metal-organic frameworks, [START_REF] Horike | Soft Porous Crystals[END_REF][START_REF] Furukawa | The Chemistry and Applications of Metal-Organic Frameworks[END_REF][START_REF] Férey | Hybrid Porous Solids: Past, Present, Future[END_REF][START_REF] Farha | Rational Design, Synthesis, Purification, and Activation of Metal-Organic Framework Materials[END_REF] consisting of imidazolate linkers that bridge metal cations to form three-dimensional porous crystalline solids, which are isomorphous to zeolitic frameworks. [START_REF] Park | Exceptional Chemical and Thermal Stability of Zeolitic Imidazolate Frameworks[END_REF] Like other MOFs, ZIFs are neither static nor rigid, instead they exhibit different types of flexibility. Flexibility in MOFs means that, upon external stimuli such as temperature, mechanical pressure or guest molecule adsorption, the cell size and shape, and therefore the pore size and geometry can change drastically -yet reversibly. [START_REF] Schneemann | Flexible Metal-Organic Frameworks[END_REF][START_REF] Bennett | Flexibility and Disorder in Metal-Organic Frameworks[END_REF] The most prominent example for flexibility in ZIFs is the so-called swing effect in ZIF-8. ZIF-8 is a low density porous framework with the sodalite (sod) topology and chemical formula Zn(mim) 2 , where mim = 2-methylimidazolate. The sod topology features large spherical cages separated by 6-ring windows of small aperture, and 4-rings that connect different cages (see Fig. 5). By torsional motions of its imidazolate linkers, the frameworks can adsorb molecules with a kinetic diameter larger than its geometric window size. [START_REF] Wu | Force Field for ZIF-8 Flexible Frameworks: Atomistic Simulation of Adsorption, Diffusion of Pure Gases as CH4, H2, CO2 and N2[END_REF][START_REF] Verploegh | Temperature and Loading-Dependent Diffusion of Light Hydrocarbons in ZIF-8 as Predicted Through Fully Flexible Molecular Simulations[END_REF] It was recently shown by ab initio molecular dynamic simulations that functionalization of the organic linker can have a substantial influence on the swing effect. [START_REF] Coudert | Molecular Mechanism of Swing Effect in Zeolitic Imidazolate Framework ZIF-8: Continuous Deformation upon Adsorption[END_REF] Furthermore, Mortada and coworkers were able to synthesize ZIF-8(Cl) and ZIF-8(Br), both in the sod topology employing 2-bromo and 2-chloroimidazolate as organic linkers.

Interestingly they found that ZIF-8(Cl) exhibits a spring behavior with the highest amount of energy stored ever in high pressure intrusion-extrusion experiments. [START_REF] Mortada | Energetic Performances of ZIF-8 Derivatives: Impact of the Substitution (Me, Cl, or Br) on Imidazolate Linker[END_REF] These findings underline the potential of ZIFs for technical application. To exploit this potential, it is pivotal to gain an atomistic understanding of the underlying mechanisms, making molecular simulations a valuable tool, and the accurate description of the frameworks' flexibility a crucial goal.

However the use of periodic Density Functional Theory (DFT) calculations, and in particular DFT-based MD (also called ab initio MD) for such studies is limited to comparatively small length and time scales -because the computational cost increases very fast with both the size and number of atoms in the unit cell, and the necessary sampling time. As long as bond breaking is not involved, less accurate molecular mechanics methods can be used to investigate larger systems for longer time scales. The difficulty is to define an energy expression that describes the relevant part of the potential energy surface with good accuracy and to determine the corresponding parameters. For computing the conformational flexibility of the porous MOF or ZIF matrix (as well as for the host-guest interactions, dominated by physisorption) conventional non-reactive force fields (FFs) that employ a separation in bonded and non-bonded terms, are a sufficiently good approximation. [START_REF] Heinen | On Flexible Force Fields for Metal-Organic Frameworks: Recent Developments and Future Prospects[END_REF] However, the determination of parameters for the hybrid organic-inorganic part remains a challenging problem.

A frequently employed solution is to use so-called generic force fields, like UFF, where the parameters are generated by a rule based system from a much smaller set of atomic parameters. [START_REF] Rappe | a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations[END_REF] This allows the consistent treatment of a wide range of systems, including a number of MOFs, however, with a very limited and uncontrolled accuracy. This approach has recently been extended to MOFs, with the UFF4MOF extension of the UFF atomic parameter set. [START_REF] Addicoat | Extension of the Universal Force Field to Metal-Organic Frameworks[END_REF][START_REF] Coupry | Extension of the Universal Force Field for Metal-Organic Frameworks[END_REF][START_REF] Boyd | Force-Field Prediction of Materials Properties in Metal-Organic Frameworks[END_REF] For most of the computational studies dealing with ZIFs, the same approach was followed, by combining generic force field parameters from UFF and AMBER, modified in an ad hoc manner to obtain experimentally observed properties. The drawback of this approach is twofold: first, the force fields obtained may not describe physical properties that were not considered in the adjustment procedure. Second, the manual adjustment may lead to the correct macroscopic observables, at the price of an unphysical microscopic picture. A brief overview of the published force fields dealing with ZIFs is given in the following paragraph.

In 2012 Jiang et al. developed a flexible force field for ZIF-8. [START_REF] Hu | Development of a Force Field for Zeolitic Imidazolate Framework-8 with Structural Flexibility[END_REF] The equilibrium bond lengths and angles were set to the experimentally measured average values. Force constants for the organic linkers were adopted from the Amber force field. [START_REF] Cornell | A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules[END_REF] The parameters involv-ing the Zn atoms were derived by fitting to experimental lattice constants. In 2013 the same group published a refined version of their original FF in order to be able to model the sorption-induced structural transition by hybrid MC/MD simulations. [START_REF] Zhang | Sorption-Induced Structural Transition of Zeolitic Imidazolate Framework-8: A Hybrid Molecular Simulation Study[END_REF] The crucial torsions around the Zn atom were fitted to a single experimental N 2 isotherm. In a parallel effort, several FFs for ZIF-8 were developed by Demontis et al. . In 2011 they published partial charges for several ZIF systems derived from cluster calculation and periodic DFT calculations. [START_REF] Rana | Estimation of Partial Charges in Small Zeolite Imidazolate Frameworks from Density Functional Theory Calculations[END_REF] In 2012 the first FF based on these partial charges was published, [START_REF] Zheng | Force Field for Molecular Dynamics Computations in Flexible ZIF-8 Framework[END_REF] which is again based on Amber. Missing parameters for the organic linker were obtained via the parmcal software, [START_REF]Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations[END_REF] which calculates bond length and bond angle parameters based on empirical rules. The parameters necessary to describe the tetrahedral ZnN 4 were taken from cluster based quantum-chemical calculations [START_REF] Hoops | Force Field Design for Metalloproteins[END_REF][START_REF] Lin | Systematic Derivation of AMBER Force Field Parameters Applicable to Zinc-Containing Systems[END_REF][START_REF] Peters | Structural Survey of Zinc-Containing Proteins and Development of the Zinc AMBER Force Field (ZAFF)[END_REF] performed in order to be able to describe Zn containing biomolecules by Amber. In 2014 this group published a new FF based on a force-matching parameterization scheme. [START_REF] Gabrieli | Fast and Efficient Optimization of Molecular Dynamics Force Fields for Microporous Materials: Bonded Interactions via Force Matching[END_REF] Zheng et al. published an additional FF for ZIF-8 in 2013. [START_REF] Zheng | Molecular Dynamics Simulations on Gate Opening in ZIF-8: Identification of Factors for Ethane and Propane Separation[END_REF] They added an artificial long range bond between the carbon atoms in neighboring imidazolate linker in a 6-ring window in order to be able to describe the "swing effect" in the correct way. Wu et al. also developed a force field for ZIF-8, again based on Amber, UFF and experimental data. The parameters for the tetrahedral ZnN 4 were adopted from the already mentioned force field by Jiang et al. from 2013. [START_REF] Zhang | Sorption-Induced Structural Transition of Zeolitic Imidazolate Framework-8: A Hybrid Molecular Simulation Study[END_REF] In 2018 Verploegh et al.

published a study on the molecular diffusion in binary mixed linker ZIFs using a flexible FF called intraZIF-FF, whose parameterization methodology has not been published to date. [START_REF] Verploegh | Quantitative Predictions of Molecular Diffusion in Binary Mixed-Linker Zeolitic Imidazolate Frameworks Using Molecular Simulations[END_REF] It is thus clear that the majority of these force fields have been parameterized in a manually involved and ad hoc manner, often mixing different sources, ranging from different other force fields to experimental and theoretical reference data, which makes it difficult to transfer these parameterization methodologies to new systems for which less reference data is available (e.g in high-throughput screening of hypothetical systems).

On the contrary, we have developed over the years a consistent, transferable and automatic parameterization strategy with the intention to trade transferability for accuracy. It is based on a machine learning approach using evolutionary algorithms to derive all bonded parameters at once, relying only on a small set of first principles reference data, namely ab initio calculated structure and Hessian H. [START_REF] Tafipolsky | Systematic First Principles Parameterization of Force Fields for Metal-Organic Frameworks Using a Genetic Algorithm Approach[END_REF] This methodology is called MOF-FF [START_REF] Bureekaew | MOF-FF -A Flexible First-Principles Derived Force Field for Metal-Organic Frameworks[END_REF] and and was recently extended to parameterize also coarse-grained FFs for MOFs. [START_REF] Dürholt | Coarse Graining of Force Fields for Metal-Organic Frameworks[END_REF][START_REF] Dürholt | The Impact of Mesopores on the Mechanical Stability of HKUST-1: A Multiscale Investigation[END_REF][START_REF] Semino | Multiscale Modeling of the HKUST-1/Poly(Vinyl Alcohol) Interface: From an Atomistic to a Coarse Graining Approach[END_REF] In the meantime, other groups have published related approaches to derive FFs for MOFs in a consistent fashion, [START_REF] Bristow | Transferable Force Field for Metal-Organic Frameworks from First-Principles: BTW-FF[END_REF][START_REF] Bristow | A General Forcefield for Accurate Phonon Properties of Metal-Organic Frameworks[END_REF] including the Quick-FF methodology, which uses a different approach to derive the parameters from the same type of reference data as in MOF-FF. [START_REF] Louis | QuickFF: A Program for a Quick and Easy Derivation of Force Fields for Metal-organic Frameworks from Ab Initio Input[END_REF] MOF-FF and Quick-FF force fields are available for a variety of different MOF families.

Yet, force fields for the important family of ZIFs are still missing. A reason for this could be that until recently both methodologies were not capable of treating periodic reference data, which is almost necessary for deriving a ZIF force field from scratch since it is very difficult to construct a representative cluster model for ZIFs. Rana and coworkers showed, for example, that one needs a cluster with a size of 400 to 500 atoms to converge partial charges on the core atoms of the cluster, because a ZIF is made up by charged fragments and charge neutrality is only achieved in the periodic system. [START_REF] Rana | Estimation of Partial Charges in Small Zeolite Imidazolate Frameworks from Density Functional Theory Calculations[END_REF] Recently Vanduyfhuys et al. extended the Quick-FF methodology in a way that also periodic reference data can be treated. [START_REF] Vanduyfhuys | Extension of the QuickFF Force Field Protocol for an Improved Accuracy of Structural, Vibrational, Mechanical and Thermal Properties of Metal-Organic Frameworks[END_REF] In this study we improved the MOF-FF methodology in several aspects to be able to parameterize force fields in respect to periodic reference data. We demonstrate its capabilities by parameterizing force fields for ZIFs of varying topology and chemical composition (ZIF-8, ZIF-8(H), ZIF-8(Cl) and ZIF-8(Br)). We thus implemented the possibility of handling periodic reference data and introduced a more efficient optimizer. The obtained FFs were then validated in great detail against experimental and ab initio calculated data. We focused especially on the question how well the ZIF flexibility (swing effect) is reproduced by our FFs, as it is crucial to their description and rather difficult to capture. Furthermore, we questioned the transferability of our FFs by applying them to polymorphs for which they were not parameterized.

Methods

The basic idea behind the MOF-FF parameterization procedure is visualized in Fig. 1. The reference information i.e. DFT optimized structure and the curvature information represented in the the matrix of the second derivatives of the energy in respect to the coordinates -→ ∇ 2 E called Hessian H are calculated for a given reference system and serve as input for the FFgen code (written in Python) . The code is then used to find the best matching bonded parameter set P of the predefined force field energy expression in respect to the provided reference information for a given set of a priori defined van der Waals (vdW) potentials and charges. Whereas charges were calculated by a fit to the electrostatic potential of the optimized reference structure, vdW parameters were taken from the well known MM3 force field, [START_REF] Allinger | The MM3 Force Field for Hydrocarbons. 1[END_REF][START_REF] Allinger | Molecular Mechanics Parameters[END_REF] and used within the dispersion damped Buckingham potential, as implemented in MOF-FF. [START_REF] Bureekaew | MOF-FF -A Flexible First-Principles Derived Force Field for Metal-Organic Frameworks[END_REF] x 1 ,y 1 ,z 1 ,...

Hessian

FFgen FF DFT

MM3-vdW FF Setup Charges

Figure 1: General scheme of the MOF-FF parameterization methodology. First the reference information (optimized structure, Hessian H and atomic charges) is calculated by DFT for a given reference system, which could be either a representative cluster or a periodic system. This information is then used together with the corresponding van der Waals parameters -taken from the MM3 force field -as input for FFgen, which tries to identify the best matching bonded parameters for the actual FF setup with respect to the provided reference information.

During the parameterization process, a so called objective function Z is used to measure how well a set of parameters P reproduces the reference data. This objective function is then minimized by a suitable numerical optimizer. An optimization cycle consists of the following three steps: first, the atomic positions of the reference system are relaxed. Second, the Hessian H is calculated by a double sided finite difference approach. Third, the objective function Z is evaluated. For this purpose geometry and Hessian are projected from Cartesian coordinates to redundant internal coordinates (RICs). The RICs comprise more than the usual 3N -6 coordinates and are made up by all bonds (str), angles (ibe), dihedrals (tor) and improper dihedrals (obe) of the system, which are also used to calculate the bonded energy of the system in the FF. The reformulated objective function Z MOF-FF is displayed below.

Z MOF-FF (P ) = Z str (P ) + Z ibe (P ) + Z obe (P ) + Z tor (P ) + Z hes (P )

It is composed of four parts, measuring the difference for a specific RIC type between reference and force field, and a fifth contribution determined by the difference in the diagonal terms of the projected Hessians. Every term is formulated as weighted mean square deviation between FF and reference as shown for the example of Z str (P ) below:

Z str (P ) = w str M str Nstr i=1 ω i (r i (P ) -r ref i ) 2 , (2)
where the sum runs over all bonds r i in the system and the individual weights ω i per redundant internal coordinate q are assigned based on the atomtypes in the system. If a RIC i, defined by its atom types, occurs n i times in the system, it gets a weight of ω i = n -1 i . The condition, therefore, is such that the parameters belonging to the FF term describing the RIC of interest are in the set of variable parameters P . Furthermore, since M str = Nstr i=1 ω i , all single contributions to Z MOF-FF are weighted in the same manner. Of course this is somehow arbitrary, since different RIC types have different units. Bond lengths are compared in Å, angles in rad and force constants in mdyn/ Å (bonds) or mdyn Å/rad 2 (angles). With this choice all contributions to Z MOF-FF are roughly in the same order of magnitude. The weight w str can be used to increase the importance of Z str to the overall objective function Z M OF -F F . In addition one could include the off-diagonal elements of the Hessian in the objective function if also the corresponding cross terms are fitted.

The advantage of this reformulated objective function is a proper weighting, since in the original formulation the RICs were weighted all equally with the consequence that RICs get a higher weight in the overall objective function when they occur more often in the system.

Fitting to Periodic Reference Data

Besides the reformulated objective function, the first main innovation affects the input required by FFgen: To parameterize FFs for ZIFs and other materials, in which it is not trivial to derive a representative zero-dimenesional cluster we extended the parameterization strategy to use reference data that is periodic in one, two or three dimensions.

Technically, the support for periodic reference data was implemented by including periodic boundary conditions (PBCs) in the calculation of Wilson's B-matrix, by augmenting the objective function with a stress dependent term and by coupling FFgen with the well known LAMMPS molecular dynamics simulator. [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF] Wilson's B-matrix is needed to for the projection into redundant internal coordinates and is defined in Eq. 3, where the q i are the internal coordinates and the x j are the Cartesian ones.

The matrix measures how an internal coordinates q i changes when a Cartesian coordinate x j has changed. In the course of this study we have implemented that PBCs are considered when B is calculated.

B ij = ij ∂q i ∂x j , (3)
As already mentioned only the atomic positions are relaxed in every optimization cycle whereas the lattice is kept fixed. To ensure that the force field also describes the lattice dimensions in the correct way, we introduced an additional term Z lattice which depends on the stress tensor S and penalizes parameter sets P which cause a large stress tensor for the optimized geometry. Z lattice is defined as follows:

Z lattice = 1 9 (C -1 • (S • V )) 2 2 (4)
where V is the cell volume and C is the cell tensor. This additional term was then weighted by a factor of 0.1 and added to Z M OF -F F . Another option to include the cell shape in the objective function would be to run, in addition to the atomic relaxation, an additional lattice optimization; however we did not follow this approach due to its higher computational cost.

To handle periodic input data it is necessary to use a molecular mechanics back-end which is able to handle relatively small unit cells and is therefore not subject to the minimum image convention. For this reason the MOF-FF total energy expression was implemented in LAMMPS and coupled to FFgen. An additional advantage of this coupling is the ability to set up the total energy expression using any potential implemented in LAMMPS. This makes FFgen a versatile tool for the parameterization of classical FFs in general and not only in the MOF-FF formulation.

A new Optimizer: CMA-ES

The second main innovation affects the algorithm, which is used to optimize the objective function. Optimization of the objective function is not trivial because the search landscape is relatively bumpy with a lot of local minima. For this purpose we rely on a stochastic zeroorder optimizer which is able to escape from a local minimum. Such optimizers are also often referred to as black-box optimizers because they do not need any further information besides the actual value of the objective function at a given search point. In the original MOF-FF parameterization the PIKAIA optimizer was used, [START_REF] Charbonneau | Genetic Algorithms in Astronomy and Astrophysics[END_REF] which is a genetic algorithm specifically designed for use on continuous variables. Generally, genetic algorithms are employed for optimization on discrete variables, thus in order to represent the continuous variables in a genome of discrete numbers, one has to predefine for every parameter p i a range consisting of an upper p max i and a lower bound p min i so that it holds p min i ≤ p i ≤ p max i . One limitation of this approach is thus that it is necessary to define these ranges a priori. If the solution is outside of the range it must be readjusted manually.

Furthermore, the periodic reference systems used in this study are substantially larger than the cluster models employed before. The primitive cell of ZIF-8 consists of 138 atoms and is more than a factor of two larger than e.g. a benzoate paddle-wheel unit. We thus needed to change the optimizer to reduce the computational cost of the parameterization process. We implemented in our FFgen code an algorithm called Covariance Matrix Adaption Evolution Strategy (CMA-ES), which was developed by Hansen for optimization problems on continuous variables. [START_REF] Hansen | Towards a New Evolutionary Computation[END_REF][START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] CMA-ES has previously been applied for force field development, namely for deriving coarse grained FFs for MOFs [START_REF] Semino | Multiscale Modeling of the HKUST-1/Poly(Vinyl Alcohol) Interface: From an Atomistic to a Coarse Graining Approach[END_REF] and recently for the parameterization of ReaxFF. [START_REF] Shchygol | Systematic Comparison of Monte Carlo Annealing and Covariance Matrix Adaptation for the Optimization of ReaxFF Parameters[END_REF] In a nutshell, CMA-ES iteratively adopts a multivariate normal distribution (starting from an isotropic normal distribution with an initial width σ) in the parameter space to find a distribution whose random samples minimize the objective function. Thus, the individuals of a generation of size λ are sampled from a multivariate normal distribution where recombination is done by selecting a new mean for the distribution. During the optimization process, the covariance matrix is updated in order to optimize the shape of the multivariate normal distribution with respect to the search landscape defined by the objective function.

This amounts to learning a second order model of the objective function similar to the approximation of the inverse Hessian in quasi-Newton optimization methods, which are known for their quadratic convergence near the optimum. 45 A more elaborate description of the optimizer can be found in the supporting information. CMA-ES is implemented in FFgen using the ask-and-tell interface of the pycma library. [START_REF] Hansen | Python Implementation of CMA-ES[END_REF] We accelerated the parameterization process by distributing the evaluation of the fitness function over parallel processes.

Using CMA-ES instead of PIKAIA leads to a substantially faster convergence. To validate this claim we fitted a force field for ZIF-8(Br), using the same setup as described in Sec. 4.1.3, with both optimizers. The convergence behavior is visualized in a logarithmic scale in Fig. 2a.

One evaluation of the objective function takes around 2 s on a single core of a current desktop GNU/LINUX workstation. A further advantage of CMA-ES is that it is not necessary to define a priori any ranges for the parameters as needed in PIKAIA. But defining ranges can be very helpful for the following reasons: at the startup of the algorithm an isotropic normal distribution with a predefined stepsize σ is initialized. Since the parameters differ by their units, they have different orders of magnitude. For this reason we define every parameter in reduced parameter units by the help of an upper and a lower bound:

pi = p i p max i -p min i (5)
However, in contrast to PIKAIA the ranges are in general not hard constraints, so the parameter is allowed to escape them. Hard ranges are only applied to prevent unphysical parameters like negative bond lengths, bond angles or force constants.

The optimizer is started from an educated guess for the parameters and their appropriate ranges. For bond lengths r and bond angles θ it is started from the geometrically measured values of the reference system. Upper and lower bound were by default set to ±10 % from the geometrically measured average value. For bond and angle force constants initial ranges are usually defined as 0 mdyn Å-1 ≤ p i ≤ 8 mdyn Å-1 and 0 mdyn Å/rad 2 ≤ p i ≤ 2 mdyn Å/rad 2 .

For dihedral potentials 0 kcal mol -1 ≤ p i ≤ 20 kcal mol -1 was used as range for the barrier and for out-of-plane potentials 0 mdyn Å/rad 2 ≤ p i ≤ 1 mdyn Å/rad 2 was the default. The initial value was set to pi = 0.5.

3 Computational Details

Obtaining the Reference Information

The reference information needed for the FF optimization was obtained with periodic DFT calculations using the QUICKSTEP/CP2K package. [START_REF] Vandevondele | Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach[END_REF] QUICKSTEP/CP2K is based on a hybrid Gaussian plane-wave approach combining a Gaussian basis for the wave-functions with an auxiliary plane wave basis set for the representation of the density.

We found that it is crucial to use a high plane-wave cutoff E cut to obtain accurate Hessians, since in the QUICKSTEP module, the computation of Coulomb and exchange-correlation energies is performed on a real space grid. This representation breaks the translational invariance of the system, which can lead to spurious forces on the atoms (egg box effect), which can have a large influence on the Hessian calculated by a finite difference approach based on the atomic forces. [START_REF] Haigis | Challenges in First-Principles NPT Molecular Dynamics of Soft Porous Crystals: A Case Study on MIL-53(Ga)[END_REF] The gradient-corrected PBE functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] was used with an empirical correction for the dispersive interactions using the "D3" method by Grimme et al. . [START_REF] Grimme | A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu[END_REF] Double-ζ valence polarized Gaussian basis sets were employed for all atoms. For C, H, and N basis sets optimized for usage with the PBE functional were employed, whereas on Zn, Cl and Br basis sets, optimized for molecules (MOLOPT) were employed. The interaction between ions and valence electrons was represented by Goedecker-Teter-Hutter (GTH) type pseudo-potentials. [START_REF] Goedecker | Separable Dual-Space Gaussian Pseudopotentials[END_REF][START_REF] Hartwigsen | Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn[END_REF][START_REF] Krack | Pseudopotentials for H to Kr Optimized for Gradient-Corrected Exchange-Correlation Functionals[END_REF] Given the systems' size, the Brillouin zone was sampled at the Γ point only.

In order to obtain accurate reference data, strict convergence criteria had to be chosen both for the SCF and for the geometry optimizations. For the SCF a convergence criterion 10 -10 Hartree was applied, whereas for the optimizations the RMS force has to be lower than 10 -7 Hartree/bohr. Wherever possible, primitive cells were used for all DFT calculations.

The reference data generation can be divided into four steps:

1. Atomic coordinates and the lattice dimensions were optimized using a plane-wave cutoff for the density E cut of 600 Ry (systems without Cl or Br) and 700 Ry (systems with Cl or Br) together with a relative cutoff E rel cut of 40 Ry.

2. Afterwards only the atomic coordinates were re-optimized using a cutoff E cut of 2500 Ry together with a relative cutoff of E rel cut of 100 Ry. These cutoffs were used for all subsequent calculations.

3. The Hessian of the optimized structure was calculated by the help of a double sided finite difference scheme using a distortion of at least 0.001 bohr.

4. Charges were calculated by the REPEAT method 54 using its implementation in CP2K.

For this purpose two type of constraints were employed: The total charge of the system has to be zero and atoms with equal atom-types get same charges.

Ab initio Molecular Dynamics Simulations

We performed ab initio molecular dynamics simulations (AIMD) of the ZIF-8 isomorphs with Born-Oppenheimer dynamics using DFT for the calculation of the energy and the atomic forces using the QUICKSTEP/CP2K package. We used the same functional, dispersion correction and basis sets as for the static calculations described in Sec. 3.1. The cutoff for the density plane-wave basis set was set to 600 Ry.

We performed all AIMD simulations using periodic boundary conditions on a single crystallographic unit cell of the material (a ZIF-8 = 16.991 Å, a ZIF-8(Br) = 16.985 Å, a ZIF-8(Cl) = 16.998 Å). We used deuterated hydrogen atoms to allow for a larger time step (1 fs) in the integration of the equations of motion. Simulations were run in the canonical ensemble (N, V, T ) using a CSVR [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] thermostat with a time constant of 1 ps to control the temperature.

The total simulation time was between 30 ps and 40 ps depending on the isomorph.

For comparison we used also data based on BLYP [START_REF] Becke | Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density[END_REF] AIMD simulations published in a previous study. For further details on these simulations it is referred to the original publication. [START_REF] Coudert | Molecular Mechanism of Swing Effect in Zeolitic Imidazolate Framework ZIF-8: Continuous Deformation upon Adsorption[END_REF]

Classical Molecular Dynamics Simulations

Classical molecular dynamics simulations were performed using our in-house developed PYDLPOLY code on 2 × 2 × 2 crystallographic unit cells of the material, employing the same lattice constants as used for the AIMD simulations. Simulations were run in the canonical ensemble (N, V, T ) using a Berendsen thermostat for equilibration runs of 0.1 ns and a Nosé-Hoover thermostat for sampling runs of 1 ns. Thermostat relaxation times were set to 0.2 ps and 2.0 ps.

Calculation of Elastic Constants

Elastic constants C ij , which are the coefficients of the second-order elastic tensor C, were computed by using the numerical first derivative of the cell gradients corresponding to the each elastic coefficient C ij . For this purpose the optimized structure was deformed in each possible direction, applying both negative and positive strain (corresponding to compression and tension), and for each deformation the atomic coordinates were energy-minimized. For deformations along the normal coordinates strains of -1.0 % to 1.0 % in steps of 0.5 % were applied, whereas for shear deformations, strains from -4.0 % to 4.0 % in steps of 2.0 %. The PYMATGEN package was used to deform the structures and to perform the final analysis, [START_REF] Ong | Python Materials Genomics (Pymatgen): A Robust, Open-Source Python Library for Materials Analysis[END_REF] using the crystallographic unit cell as a reference configuration. Elastic constants at the ab initio level of theory were computed using the QUICKSTEP/CP2K package employing the same setup as used for the production of the reference data. Elastic constants at the FF level were computed by using PYDLPOLY and a 2 × 2 × 2 supercell.

Results and Discussion

4.1 FF Parameterization

Systems of Interest

We focused in this study on the development of four different FFs for ZIF-8, ZIF-8(H), ZIF-8(Br) and ZIF-8(Cl). All are composed of zinc and differently functionalized imidazolate linkers, and form 3D crystalline networks with the sod topology (chemical composition and atom types shown in Fig. 3).

Furthermore we calculated also the reference data for six other polymorphs with the ZIF-8(H) composition, which we used to probe into the transferability of the FFs between different topologies. Crystallographic and chemical information about the investigated systems are listed altogether in Table 1.

Partial Charges

The atomic partial charges obtained for ZIF-8, listed in Table 2, are in good agreement with those published by Rana et al. [START_REF] Rana | Estimation of Partial Charges in Small Zeolite Imidazolate Frameworks from Density Functional Theory Calculations[END_REF] These authors computed both single-point charges for a high-symmetry relaxed structure of ZIF-8 and charges as average over several snapshots obtained from ab initio MD simulations in the Born-Oppenheimer scheme and in the Car-Parrinello scheme. From the difference between them, they concluded that single-point charges are in agreement with the average value during the time evolution of the system and are thus sufficient. Since they did not calculate charges for ZIF-8(H), ZIF-8(Cl) and ZIF-8(Br), we demonstrate by this comparison the validity of the methodology to derive charges for our FF and used it also for ZIF-8(H), ZIF-8(Cl) and ZIF-8(Br).

For the seven hydrogen substituted ZIF polymorphs the charges are shown in Table S1.

Depending on the topology the charges on the Zn range from 0.4645 e to 0.5971 e with a standard deviation of 0.04 e, indicating that the impact of the topology on the charges is relatively minor, compared to the natural methodology-related uncertainty. We thus chose as charges used for our FF a unique set of partial charges for this class of materials, obtained by averaging over the seven structures.

In addition the charges for the halogen-substituted structures are listed in Table S2. The charges in the Zn atoms for the four chemically distinct species in the sod topology differ up to 50 % (ZIF-8: 0.7290 e, ZIF-8(H): 0.5118 e, ZIF-8(Cl): 0.4986 e, ZIF-8(Br): 0.4714 e). This demonstrates the considerable impact of the substituent on the whole framework and emphasizes the need for distinct parameter sets for every chemically distinct species, especially in the case of the non-polarizable MOF-FF where no Coulomb exclusions are applied.

Intramolecular Parameterization

The force field energy expression was set up for the four systems in the way that only diagonal terms were applied and no cross terms like stretch-bend potentials were used. All dihedral potentials were set up by imposing a multiplicity of two, besides the C2-N-Zn-N and C1-N-Zn-N where a multiplicity of three was used. Due to the disorder of the methyl groups, no dihedral potential was applied for N-C2-C3-H3 dihedral in ZIF-8(H) and the involved parameters were predefined for the actual parameterization. The comparison of the actual dihedral angle values was not included in the objective function (except in the case of ZIF-8, where w tor = 0.1 was employed), since the geometry is already imposed by the chosen multiplicity of the potential. Charges were chosen as described above and MM3 vdW parameters were employed. Equal weights of one were assigned to all different contributions (w str , w ibe , w wdiag ) of the objective function Z.

When using heuristic optimizers like CMA-ES on bumpy search landscapes, it is recommended to run the optimizer several times with increasing population sizes. [START_REF] Auger | A Restart CMA Evolution Strategy with Increasing Population Size[END_REF] For this reason we started always several runs increasing the populations size λ from the default size (λ = 4 + (3 ln N par ) = 16) up to 32 individuals per generation. To prove that λ = 32 is sufficient, we performed additional parameterization runs with population sizes from 8 to 128 individuals for ZIF-8(Br). Furthermore we conducted additional runs with 128 individuals and a doubled stepsize σ = 0.6. The results shown in Fig. 2b demonstrate that λ = 32 is an appropriate population size for our problem. In addition we perturbed during our restart strategy also the parameter ranges, if we felt that they were too loose or too strict.

Force Field Validation

Structural Properties

As a first validation step, we verified how well the parameterized FFs reproduce the structural properties of the materials, with respect to DFT-optimized structures. Fig. 4 shows a comparison between the final bond lengths, bending angles and dihedral angles obtained at the two levels of theory. In case of the FFs, we have also relaxed the lattice for this comparison. Note that since no Coulomb exclusion is used in MOF-FF and due to the complex periodic structures it is not clear beforehand that these structural parameters are always reproduced correctly. Overall agreement is excellent, with the biggest deviations observed for ZIF-8, due to the disorder introduced by the methyl groups' free rotation and the fact that no dihedral potential is included for the N-C2-C3-H3 torsion. We also investigated the lattice dimensions at zero Kelvin (full cell energy optimization) in comparison to the DFT results; see Table 3. The lattice dimensions of the structures which served as reference systems for the parameterization coincide almost perfectly. The biggest difference of 0.05 Å arises for ZIF-8, which is again due to the disorder introduced by the methyl groups. The MOF-FF parameterization methodology relies on a single-structure fit, which works best as long as structures are highly ordered and only one isomeric form exists.

a[ Å] b[ Å] c[ Å] α[ • ] β[ • ] γ[ • ] ZIF-

Deformations: Vibrations and Elasticity

Next we checked the accuracy of the FFs in representing deformations from the relaxed structures, comparing normal mode frequencies and elastic constants, which are directly linked to the second derivative of the energy in respect to the atomic coordinates, and unit cell parameters, respectively.

Fig. 4d shows the comparison of the FF normal modes against DFT data. The agreement is good, in particular for the low-frequency vibrations modes. Those modes are maximally delocalized and are mainly responsible for the lattice vibrations and flexibility of the framework, and consequently crucial to reproduce correctly. The modes between 500 cm -1 and 2000 cm -1 are more localized and involve especially distortions of the aromatic imidazolate rings. A possibility of achieving a better accuracy in this region in the future could be the incorporation of cross-terms into the FF.

Elastic constants calculated by different DFT methods -some published in the literature, and some computed as part of this work -are compared to those obtained by our FFs in Table 4. Zheng et al. predicted recently the elastic constants of differently functionalized ZIFs in the sod topology. They found that electron withdrawing groups improve the mechanical stability of the materials (ZIF-8(H) < ZIF-8(Cl) < ZIF-8(Br)). [START_REF] Zheng | Theoretical Prediction of the Mechanical Properties of Zeolitic Imidazolate Frameworks (ZIFs)[END_REF] Although the absolute numbers of our DFT calculations differ up to a few GPa this trend is also reproduced by our calculations. Furthermore this trend is even reproduced in our force field calculations. Only exception is C 44 , but already the differences from the reference calculations are here very subtle in comparison to those published in the literature. In general our FF systematically overestimates C 11 and C 12 whereas it underestimates C 44 . To improve here, one could either incorporate cross-terms or one could adjust the anharmonicities of the bonded potentials in order to tune the curvature of the effective potential.

Molecular Mechanism of the Swing Effect

The flexibility in ZIF-8 analogues is governed by intraframework dynamics and involves, as detailed in the introduction, the so-called swing effect of the imidazolate linkers that allows molecules larger than its window size to diffuse into the framework. Coudert investigated this effect in detail by ab initio MD simulations. [START_REF] Coudert | Molecular Mechanism of Swing Effect in Zeolitic Imidazolate Framework ZIF-8: Continuous Deformation upon Adsorption[END_REF] The swinging motion of the imidazolate linkers was characterized by the dihedral angle φ Zn 3 -Zn 2 -Zn 1 -CH 3 of the imidazolate around the Zn 1 -Zn 2 axis, where the "reference" of 0 • is the 6-ring of Zn (i.e. the window connection the cages) as shown in Fig. 5.

We used our force fields to run classical MD simulations in the (N, V, T ) ensemble using Histograms of the swinging angle φ for the four different systems are shown in Fig. 6. For ZIF-8, data from both PBE and BLYP is available (see Fig. 6a). From the difference between them one can estimate the errors made by different DFT exchange-correlation functionals which is an important measure for the accuracy of FFs. As one can see the differences are quite substantial, the histogram calculated with PBE has a larger spread, the "thermal" swing motion goes up to 20 • , whereas in the case of BLYP it is only 15 • . The histogram predicted by MOF-FF lies in between the two curves predicted by the two functionals having a cumulative overlap to the PBE functional of S = 0.84.

For ZIF-8(H) again data from both functionals is available (see Fig. 6c). Here the PBE and BLYP curves are more similar as in the case of ZIF-8. As already predicted by Coudert ZIF-8(H) exhibits a swing motion which much larger amplitude than ZIF-8. This effect is nicely resembled by our force field, resulting in a cumulative overlap S of 0.97. The thermal swing motion goes all the way up to 35 • (instead of 20 • for ZIF-8).

For the halogenated species (ZIF-8(Cl) and ZIF-8(Br)) the agreement between FF and DFT calculated data is less good than for ZIF-8 and ZIF-8(H), resulting for both systems in a cumulative overlap of S = 0.76. For ZIF-8(Br) the general shape of the curve is reproduced well but it is stretched along the x-axis towards higher swing angles by 5 • (see Fig. 6b).

In case of ZIF-8(Cl), the FF predicted curve is compressed to lower swing angles by 5 • degrees (see Fig. 6d). We emphasize, however, that these differences are of the same order of magnitude as the difference between the curves for ZIF-8 predicted at the DFT level of theory by the different exchange-correlation functionals.

Transferability and Overfitting

In general, the MOF-FF strategy of parameterizing a FF is to use non-periodic reference information if possible. This has the advantage of fitting in respect to sterically relaxed building blocks, which should increase the transferability of the FFs. For example, our

Copper paddle-wheel (PW) force field is not parameterized in respect to HKUST- For every chemically distint system the overlap between the FF and our reference DFT method is plotted togther with the histogram intersection S, giving the cumulative overlap. Uncertainties on the histograms were determined by dividing each trajectory into 10 evenly spaced and randomly arranged subtrajectories, and calculating for each the histogram. We plotted always ±2σ as error.

strained environments like MOF-14 and HKUST-1, i.e. in the tbo and pto topology. [START_REF] Amirjalayer | Exploring Network Topologies of Copper Paddle Wheel Based Metal-Organic Frameworks with a First-Principles Derived Force Field[END_REF] However, in the current work we chose a different strategy, since we fitted directly to the strained periodic structure. In the following, the transferability of the FF to other topologies and the amount of overfitting is assessed.

We thus optimized the structure of six ZIF-8(H) polymorphs and compared the energies to DFT results, obtained in this work as well as previously by Baburin et al. [START_REF] Baburin | Enumeration of Not-Yet-Synthesized Zeolitic Zinc Imidazolate MOF Networks: A Topological and DFT Approach[END_REF] (who did, however, not include any dispersion corrections). Is it possible to reproduce the relative stabilities of the isomorphs and predict the correct energetic ordering? The results are shown in Fig. 7a. DFT predicts the following energetic ordering: coi < zni < nog < ZIF-4 < ZIF-1 < ZIF-8(H) < ZIF-6, whereas the FF predicts coi < zni < ZIF-1 < nog < ZIF-4 < ZIF-8(H) < ZIF-6. Thus, the only difference is that the stability of ZIF-1 is overestimated.

The reason for this can be found in the lower symmetry of the other polymorphs compared to ZIF-8(H) which was used as reference system. Fig. 7b shows the distribution of the N -Zn -N angles for the other investigated polymorphs based on the DFT optimized structures. They show substantially larger deviations from the ideal tetrahedron angle as compared to ZIF-8(H). So they are further away from the reference structure which results in a worse performance. The same behavior can also be observed for the lattice constants of the other polymorphs as shown in Table 3. In comparison to the systems which served as reference information they show an inferior performance, however they are still well within the range of lattice constants predicted by different DFT methods. An exception is ZIF- qualitative agreement is achieved, but less quantitative agreement is observed compared to the accuracy on the actual reference system. We note, however, that the predictive power in respect to the other polymorphs could be increased by fitting a FF for each of them, but by this one would loose the possibility to compare the relative stabilities between different polymorphs. A solution would be to perform multi-structure fits using different polymorphs as reference system at once.

Conclusions

The aim of this study was to exploit the systematic and consistent MOF-FF force field parameterization methodology to arrive at an accurate and efficient potential for a range of ZIFs. By trading transferability for accuracy we parameterized explicitly for the chemically distinct systems, namely ZIF-8, ZIF-8(H), ZIF-8(Br) and ZIF-8(Cl), using the sod topology as a reference. For this purpose we needed to improve and extend the original methodology in several aspects. First, we added the possibility to employ also periodic reference systems, which enables us to use our approach also for those systems which can not be easily truncated to a cluster representations, as for example rod based MOFs like MIL-53 or MOF-74, or ZIFs in general. Second, we replaced the original genetic algorithm based global optimizer by another more efficient evolutionary strategy that is better suited for continuous variables.

The CMA-ES converges substantially faster and without constraining to parameter ranges, which is extremely beneficial since the numerical effort for the evaluation of the target function is numerically much more involved in case of a fit to periodic reference systems like in the case of ZIFs. Furthermore, it also paves the way towards a completely automated black box algorithm for FF parameterization.

As expected, the force fields are able to well reproduce structure and lattice parameters as well as dynamic properties like vibrational normal modes and elastic constants in comparison to available experimental results as well as the computed reference data at DFT level of theory.

A much more subtle property of ZIFs is their inherent flexibility, namely the so called swing effect, which allows molecules larger than the geometric window size to diffuse into the framework by a slight rotation of the imidazolate linkers. The ability of a force field to reproduce this behavior in an accurate way is crucial for its use in simulating guest molecule adsorption or for example heat conduction. We find that our force field is able to reproduce this flexibility within the same range of accuracy that is obtained in ab initio MD simulations using different exchange-correlation functionals.

In order to validate the transferability of the parameter set, also the energetic ranking of other ZIF topologies was tested. Interestingly, we find that the energetic ordering of the polymorphs ZIF-1, ZIF-4, ZIF-6, coi, nog and zni computed with the force field fitted to ZIF-8(H) is nicely reproduced as compared to periodic DFT calculations. A further improvement could be achieved by fitting to several reference structures at once to avoid overfitting. Future development efforts will likely aim at an efficient implementation of these multistructure fits.

All in all, this study demonstrates the potential of the here introduced extended MOF-FF parameterization methodology. Based on a single periodic reference structure (geometry and curvature) it is possible to derive an accurate force field in a consistent and systematic way which can be used to substantially extend length and time scales in MD simulations, within the constraint of no bond breaking, in an accuracy close to periodic DFT.

  , σ = 0. 3 λ = 16, σ = 0. 3 λ = 32, σ = 0. 3 λ = 64, σ = 0. 3 λ = 128, σ = 0. 3 λ = 128, σ = 0. 8 (b)

Figure 2 :

 2 Figure 2: Logarithmic convergence behavior of (a) the originally employed PIKAIA genetic algorithm vs. CMA-ES for a parameterization run of ZIF-8(Br) and (b) for several parameterization runs of ZIF-8(Br) using CMA-ES but different population sizes λ and initial step sizes σ.

Figure 3 :

 3 Figure 3: Chemical systems investigated in this study, showing the atom typing used throughout the text.

Figure 4 :

 4 Figure 4: Scatter plots visualizing the performance of our FFs with respect to the ab initio reference data. (a) Comparison of bond lengths (b) Comparison of bond angles (c) Comparison fo dihedral angles (d) Comparison of vibrational normal modes.

Figure 5 :

 5 Figure 5: Sodalite topology sod of ZIF-8 and view of the 6-ring window of ZIF-8, with the swing dihedral angle φ marked in light red. The reference angle of 0 • is the plane of the 6-ring.

Figure 6 :

 6 Figure 6: Histograms of the swing angle of the imidazolate linkers, for ZIF-8 (a), ZIF-8(Br) (b), ZIF-8(H) (c) and(ZIF-8(Cl) (d) computed at different levels of theory.For every chemically distint system the overlap between the FF and our reference DFT method is plotted togther with the histogram intersection S, giving the cumulative overlap. Uncertainties on the histograms were determined by dividing each trajectory into 10 evenly spaced and randomly arranged subtrajectories, and calculating for each the histogram. We plotted always ±2σ as error.

4 .

 4 The reason for this difference is not only the lower symmetry of the Zn coordination environment, but also the existence of several different phases of ZIF, which bedevils the situation for a FF fitted in respect to the sod topology. Finally we checked how well our FF describes the elastic tensor of a different polymorph, namely ZIF-6. DFT predicts the following matrix: c 11 = 8.33 GPa, c 12 = 7.58 GPa, c 13 = 10.99 GPa, c 33 = 16.65 GPa, c 44 = 0.67 GPa, c 66 = 0.30 GPa. The FF predicts: c 11 = 6.79 GPa, c 12 = 5.68 GPa, c 13 = 7.62 GPa, c 33 = 12.68 GPa, c 44 = 0.96 GPa, c 66 = 0.55 GPa. Also here the same trend is observed,

Figure 7 :

 7 Figure 7: (a) Scatter plot of the total energies of the 7 ZIF-8(H) polymorphs computed by the FF and two different DFT methods, namely by Baburin et al. 66 and us as described in Sec. 3.1. (b) Violin plot of the N-Zn-N angle distributions in the seven ZIF-8(H) polymorphs, based on the DFT optimized structures. Minimal and maximal values are shown together with the mean of the distributions.

Table 1 :

 1 Summary of the systems investigated in this study.

		Space group	Bravais lattice	Net topology	Functional group
	ZIF-1	P 2 1 /c	mono	crb	H
	ZIF-4	P bca	ortho	cag	H
	ZIF-6	I4 1 /amd	tetra	gls	H
	ZIF-8(H)	I43m	cubic	sod	H
	coi	I4 1	tetra	coi	H
	nog	P 2 1 c	mono	nog	H
	zni	I4 1 cd	tetra	zni	H
	ZIF-8(Br)	I43m	cubic	sod	Br
	ZIF-8(Cl)	I43m	cubic	sod	Cl
	ZIF-8	I43m	cubic	sod	Met

Table 2 :

 2 REPEAT charges (units in partial electron charges) for ZIF-8 computed by us in respect to the optimized structure in comparison to those published by Rana et al. based on snapshots extracted from BOMD and CPMD simulations.[START_REF] Rana | Estimation of Partial Charges in Small Zeolite Imidazolate Frameworks from Density Functional Theory Calculations[END_REF]

		Zn	N	C1	C2	C3	H1	H3
	SP	0.7290 -0.3417 -0.1923	0.4937 -0.5276	0.1627	0.1373
	BOMD 0.7362 -0.3008 -0.1924	0.4339 -0.6042	0.1585	0.1572
	CPMD	0.6894 -0.2800 -0.1910	0.4184 -0.5726	0.1536	0.1481

Table 3 :

 3 Lattice dimensions for the investigated ZIF systems computed computed by the FFs and different DFT methods. DFT results from the literature are listed only in extracts.

Table 4 :

 4 Elastic constants of ZIF-8, ZIF-8(H), ZIF-8(Br) and ZIF-8(Cl) computed by the FFs and different DFT methods. DFT results from the literature are listed only in extracts.

		C 11 [GPa]	C 12 [GPa]	C 44 [GPa]
	ZIF-8 FF	8.54	6.55	0.62
	ZIF-8 63	11.04	8.32	0.94
	ZIF-8(H) FF	6.65	4.95	1.12
	ZIF-8(H) cp2k	5.39	4.51	0.32
	ZIF-8(H) 62	8.95	7.59	2.36
	ZIF-8(Cl) FF	9.92	7.84	0.46
	ZIF-8(Cl) cp2k	9.23	7.35	0.86
	ZIF-8(Cl) 62	12.30	9.98	3.58
	ZIF-8(Br) FF	10.51	8.65	0.19
	ZIF-8(Br) cp2k	10.33	8.31	0.88
	ZIF-8(Br) 62	15.92	11.57	6.56
	the same lattice constants as in the corresponding AIMD for the different functionalizations
	to compare with the DFT results. Originally only ZIF-8 and ZIF-8(H) were investigated
	using the BLYP functional. 56,57 We performed further AIMD simulations, in the course of
	the present work and the work by Chaiplas et al. 64 , for ZIF-8, ZIF-8(H), ZIF-8(Br) and
	ZIF-8(Cl) using the PBE functional.

Acknowledgement

J.P.D. is grateful for the financial support by the Fonds der Chemischen Industrie (FCI) and by the Research School Plus from the Ruhr-Universität Bochum. Further financial support from the Deutsche Forschungsgemeinschaft (DFG) is acknowledged (grants SCHM1389/8-1 and SCHM1389/10-1). Access to HPC platforms was provided by a GENCI grant (A0050807069).

We thank Hendrik Heenen for implementing MOF-FF in LAMMPS.