On flexible force fields for metal-organic frameworks: Recent developments and future prospects.

Classical force field simulations can be used to study structural, diffusion, and adsorption properties of metal-organic frameworks (MOFs). To account for the dynamic behavior of the material, parameterization schemes have been developed to derive force constants and the associated reference values by fitting on ab initio energies, vibrational frequencies, and elastic constants. Here, we review recent developments in flexible force field models for MOFs. Existing flexible force field models are generally able to reproduce the majority of experimentally observed structural and dynamic properties of MOFs. The lack of efficient sampling schemes for capturing stimuli-driven phase transitions, however, currently limits the full predictive potential of existing flexible force fields from being realized. This article is categorized under: Structure and Mechanism > Computational Materials ScienceMolecular and Statistical Mechanics > Molecular Mechanics.

INTRODUCTION

Metal–organic frameworks (MOFs) are an emerging class of crystalline nanoporous materials consisting of metallic nodes interconnected by organic linkers (Zhou, Long, & Yaghi, 2012). MOFs have been studied extensively over the past decade due to their promising applications in areas such as gas storage (Zhao et al., 2016), adsorptive separations (Li, Kuppler, & Zhou, 2009; Li, Sculley, & Zhou, 2012), sensor (Kreno et al., 2012) and electronic devices (Stassen et al., 2017), and catalysis (Huang, Liang, Wang, & Cao, 2017). Typically, MOFs have exceptionally large surface areas and porosities, with NU‐110 (Farha et al., 2012) as the surface area record holder with a BET surface area of 7,140 m2/g. All MOFs exhibit some form of flexibility (Chang, Yang, Xu, Hu, & Bu, 2015; Coudert, 2015; Schneemann et al., 2014), ranging from lattice vibrations around equilibrium positions to large‐scale structural transformations upon external stimuli. These flexible behaviors can have important implications for MOF application prospects (Burtch, Heinen, Bennett, Dubbeldam, & Allendorf, 2017) by enabling, for instance, exceptional gas storage or separation performances that is not possible in traditional rigid materials.

Selected large‐scale flexible modes demonstrated in MOFs are shown in Figure 1. An anomalous structural transformation due to temperature is negative thermal expansion (NTE) (Figure 1a), which is the contraction of the unit cell upon increasing temperature (Miller, Smith, Mackenzie, & Evans, 2009). Examples of well‐studied MOFs exhibiting NTE are the isorecticular MOF series (IRMOFs; Eddaoudi et al., 2002; Li, Eddaoudi, O'Keeffe, & Yaghi, 1999; Rowsell & Yaghi, 2004; Yaghi et al., 2003) that consist of ZnO4 tetrahedral building blocks connected with benzene 1,4‐dicarboxylate linkers. Flexible force field calculations predict that the thermal expansion coefficients of these materials increase with linker length (Dubbeldam, Walton, Ellis, & Snurr, 2007), and recently it was also predicted that the presence of adsorbates can control their thermal expansion behavior (Balestra et al., 2016). Another material that exhibits NTE is HKUST‐1 (Wu et al., 2008). HKUST‐1 has paddlewheel nodes interconnected by benzene‐1,3,5‐tricarboxylate ligands (Chui, Lo, Charmant, Orpen, & Williams, 1999). Rotation of the paddlewheel node and “trampoline‐like” motions of the organic linker at frequencies of 58 cm−1 (1.7 THz) and 81 cm−1 (2.4 THz), respectively, were proposed to explain the origin of its NTE (Ryder, Civalleri, Cinque, & Tan, 2016).

Figure 1

Illustration of different flexibility behaviors reported in metal–organic frameworks (MOFs). (a) Negative thermal expansion in isorecticular MOF series I (IRMOF‐1), (b) breathing in MIL‐53(Cr), (c) swelling in MIL‐88D and (d) negative gas adsorption in DUT‐49

Guest adsorption can induce phase transitions such as breathing (Figure 1b): a transition from a large pore (lp) to narrow pore structure (np), as observed in MIL‐53(Cr) (Serre et al., 2002) due to stimuli such as CO2 sorption. A swelling behavior has also been observed in MIL‐88 (Serre et al., 2007) due to pyridine uptake (Figure 1c). The breathing transition of MIL‐53(Cr) is responsible for the distinctive stepped curvature of its CO2 adsorption isotherm (Llewellyn et al., 2008; Serre et al., 2002) whereby, up until a loading of 2 mmol/g, the unit cell is in the lp phase and beyond that it is in the np phase. These structural changes are reversible and the material again transitions toward its lp phase at around 5 bar. Neimark, Coudert, Boutin, and Fuchs. (2010) argued that if a certain threshold of adsorption‐stress on the material has been reached then a structural transformation is induced, and the predicted magnitude of the adsorption stress needed to deform MIL‐53(Cr) was found to be in agreement with the experimentally applied external mercury pressure (Neimark et al., 2011). Furthermore, materials generally become stiffer upon guest adsorption (Canepa et al., 2015), as exhibited by the increased shear and bulk modulus for ZIF‐8 (Ortiz, Boutin, Fuchs, & Coudert, 2013a) and ZIF‐4 (Bennett et al., 2011), respectively. Interestingly, however, Fuchs and coworkers (2015) showed using classical force field models for two simple pore models that softening can occur at low loadings before stiffening then occurs at higher adsorbate loadings.

Another interesting phenomenon is the gate‐opening behavior observed in ZIF‐8, which can have important implications for molecular separations and sensing devices. The material’s N2 adsorption isotherm at 77 K shows a step at p/p0 around 0.2 (Fairen‐Jimenez et al., 2011), with the exact pressure depending on the ZIF‐8 particle size (Tanaka et al., 2015). This behavior has been attributed to the rotation of its methyl‐imidazolate linkers (Casco et al., 2016; Ryder et al., 2014). Recently, negative gas adsorption (Figure 1d), which is the spontaneous desorption of adsorbates at higher pressures, was observed in DUT‐49 (Krause et al., 2016). In this structure, methane adsorption activates the open to contracted phase transition which occurs due to deformation of the linker (Evans, Bocquet, & Coudert, 2016).

The origin of these phenomena, as well as flexibility in MOFs in general, remains a highly active field of research. Fundamental insight into these flexibility properties are usually obtained from a combination of experiments (e.g., tetrahertz vibration spectroscopy) or ab initio density functional theory (DFT) calculations (Ryder et al., 2014; Ryder, Civalleri, Cinque, & Tan, 2016; Tan et al., 2012). Tensorial analysis of the elastic constants has also shown to be a powerful approach for elucidating flexible behaviors in MOFs (Bennett, Cheetham, Fuchs, & Coudert, 2017; Ortiz, Boutin, Fuchs, & Coudert, 2012; Ortiz, Boutin, Fuchs, & Coudert, 2013b; Ryder, Civalleri, & Tan, 2016; Tan et al., 2012). The elastic constant tensor can be written, within the elastic limit, by Hooke’s law as (Nye, 1957)with σ being the stress and ɛ being the strain tensor using the Einstein summation and Voigt notation (Voigt, 1910). To the best of our knowledge, the only experimental data on the elastic tensor in the MOF‐community is available for single‐crystal ZIF‐8 obtained from Brillouin scattering (Tan et al., 2012), showing excellent agreement with DFT‐B3LYP calculations. The authors found an extremely low shear modulus of G min = 0.967 GPa, attributed to cooperative Zn‐imidazole‐Zn (bridging) and N‐Zn‐N (tetrahedral) bonding angles.

Softening of elastic constants can lead to negative eigenvalues of the elastic tensor and violation of the Born stability criteria (Born, 1940; Born & Huang, 1954) Crystal‐to‐amorphous transitions have been attributed to the reduction of shear moduli (Koike, 1993; Ortiz et al., 2013a; Rehn, Okamoto, Pearson, Bhadra, & Grimsditch, 1987; Tan et al., 2012). Note that during phase transitions chemical bonds do not necessarily break. This raises the question if chemical instability (e.g., weak chemical bonds) is caused by mechanical instability (lowering of elastic moduli) or vice versa. Missing linkers in UiO‐66 resulted in a decrease of the elastic moduli (Bennett et al., 2017). However, low elastic moduli imply that the material is more flexible and therefore can “open up” more for water to hydrolyze coordination bonds.

Classical force field calculations can also be used to understand MOF flexible properties. Flexible force field models have the advantage of making material predictions at significantly lower computational costs than ab initio methods, provided an accurate force field model is available for the structure of interest. This requirement motivates the discussion we provide in later sections regarding different approaches for the development of suitable flexible force fields in MOFs. Considerable advances have already been made in this field and, as an example, we show several force field‐calculated properties of IRMOF‐1 in Table 1 and their strong agreement with both experiments and ab initio calculations. For instance, Boyd et al. ( 2017 ) compared several generic force fields by computing bulk moduli and thermal expansion coefficients of various MOFs. All force fields reproduced the experimental or ab initio values within 5 GPa, and the expected NTE was also predicted, although the thermal expansion coefficients differ by up to 25 ppm/K.

Table 1

Lattice parameter a, volumetric thermal expansion coefficient α , bulk modulus K, and Young’s modulus E of isorecticular MOF series I (IRMOF‐1). Experimental values are reported at 300 K, density functional theory (DFT), and force field values are reported at 0 K, unless stated otherwise

	a (Å)	α V · 10−6 (K−1)	K (GPa)	E 100 (GPa)	E 111 (GPa)
Experimental	25.885 (Li et al., 1999)	−39 (Lock et al., 2010)	—	7.9 (Bahr et al., 2007)a	—
PW‐91 (Bahr et al., 2007)	26.04	—	16.33	21.95	10.06
PBE‐D3 (Banlusan & Strachan, 2017)	26.09	—	15.76	18.88 (17.7a)	2.91 (2.5a)
Greathouse and Allendorf (2006)and Greathouse and Allendorf (2008)	26.05	−36	20.0	35.5 (14.9a)	—
Dubbeldam et al. (2007)	25.965	−55	17.71	22.42b	2.90b
Han and Goddhard III (2007)c	25.291	−23.91	19.37	42.73b (31.14b , a)	5.29b (3.97b , a)
Tafipolsky, Amirjalayer, and Schmid (2007)	25.946	—	10.8a	—	—
Bristow, Tiana, and Walsh (2014)	25.901	−15.80	11.95	37.42 (Boyd, Moosavi, Witman, & Smit, 2017)b	1.19 (Boyd et al., 2017)b

300 K.

Young’s modulus calculated from the elastic tensor using the ELATE code (Gaillac, Pullumbi, & Coudert, 2016).

10 K.

The remainder of this review is organized as follows: we first summarize how unit cells of structural models are defined and modeled in terms of force field potentials, as well as the theoretical framework by which elasticity and adsorption properties are computed. The next section highlights various parameterization schemes for creating accurate flexible force field models, starting with the seminal work of Greathouse and Allendorf. The final last section then reviews the theoretical and practical challenges and prospects associated with force field development efforts into the future.

CLASSICAL FORCE FIELD SIMULATION METHODOLOGY

Intra‐ and intermolecular potentials

Classical or molecular mechanics force fields rely on user‐defined interaction potentials which have a chemically meaningful interpretation (Dubbeldam, Torres‐Knoop, & Walton, 2013; Frenkel & Smit, 2002; Leach, 2001). The total interaction energy is expressed as the sum of all these potentials, such aswhich is usually divided into bonding and nonbonding (Van der Waals and electrostatic) contributions. Equation (2) is generally referred to as a force field since most terms are written within the harmonic approximation and thus only require a force constant and a reference value. In general, one can add as many ad hoc potentials as one desires. The more terms that are present in the force field, however, the more computationally expensive the simulations become. There is always a trade‐off between accuracy and computational cost, although potentials that describe polarization effects (Becker, Heinen, Dubbeldam, Lin, & Vlugt, 2017) and donor–acceptor interactions (Campbell, Ferreiro‐Rangel, Fischer, Gomes, & Jorge, 2017; Dzubak et al., 2012; Heinen, Burtch, Walton, Fonseca‐Guerra, & Dubbeldam, 2016; Kulkarni & Sholl, 2016) are often crucial to describe the adsorption properties in strongly interacting systems such open‐metal site MOFs (Fang, Demir, Kamakoti, & Sholl, 2014; Sun & Deng, 2017).

To represent the system’s complex electrostatic potential classically, including multipole interactions, point charges must be assigned to the system. This can be done by minimizing an error function of the ab initio calculated electron density and atomic charges under the constraint that the total sum of the atomic charges equals the total charge of the system (Jensen, 2009). For periodic systems, this can be done using the repeating electrostatic potential extracted atomic (REPEAT) method (Campa, Mussard, & Woo, 2009). There is no unique method of constructing point charges as these are not quantum observables, and other commonly used periodic charge partitioning schemes include density derived electrostatic and chemical (DDEC; Manz & Sholl, 2010; Manz & Sholl, 2012) and Hirshfeld (1977). Most partition schemes assign partial charges to the atoms instead of formal charges. The partial atomic charges assigned to the framework can heavily affect adsorbate uptake and separation selectivities (Castillo, Vlugt, & Calero, 2008; Walton et al., 2008; Zheng, Liu, Yang, Zhong, & Mi, 2009), as well as thermal expansion coefficients (Hamad, Balestra, Bueno‐Perez, Calero, & Ruiz‐Salvador, 2015). For structural transformations of the unit cell, the use of fixed charges becomes questionable since the true electrostatic potential energy surface also changes as the atomic positions change. Charge equilibration schemes (CES) can provide an approximate solution by predicting atomic charges based on the current geometry and connectivity of the system (Rappe & Goddard, 1991). The efficiency of CES also make them popular for large‐scale screening studies of existing and hypothetical MOFS (Wilmer & Snurr, 2011), and, as such, these CES remain an active area of research in force field development (Bauer & Patel, 2012).

Developing flexible force field models for MOFs also involves the construction and parameterization of bonding, bending, and torsion potentials. The construction versus parameterization stages of determining these potentials should be distinguished. The functional form’s construction involves deciding which mathematical forms should be used for the bonds, bends, and torsions of the system to appropriately represent the material dynamics. For instance, this deals with the question: is a simple harmonic potential between the metal and the ligand atoms sufficient or is a more complex functional form such as the Morse potential needed. The parameterization stage deals with the separate question of: given a certain set a functional forms, what are the associated force constants and reference values that should be assigned?

For liquids and gases, force fields have been parameterized based on experimental phase equilibrium data, typically in the form of vapor–liquid equilibrium data (TraPPE; Martin & Siepmann, 1998; Martin & Siepmann, 1999), conformational energetics from ab initio calculations (OPLS‐AA; Jorgensen, Maxwell, & Tirado‐Rives, 1996; Jorgensen & Tirado‐Rives, 1988), heats of formation, and vibrational spectra (MM3; Allinger, Yuh, & Lii, 1989). For solids, the choice of experimental observables with which to parameterize the force field models is often less obvious. Zeolite force fields have been parameterized on experimental infrared spectra (Demontis, Suffritti, Quartieri, Fois, and Gamba (1988) and Nicholas, Hopfinger, Trouw, and Iton (1991)), and ab initio cluster calculations (Hill and Sauer (1995)). Even though the chemical environment of atoms in zeolites is less chemically diverse than in MOFs, it still remains challenging to produce accurate and transferable potentials in these systems (Bueno‐Pérez et al., 2012). The challenge is even greater in MOFs, including the challenge of more diverse chemical environments and organic–inorganic linkages.

Defining simulation cell, boundary conditions, and elasticity properties

The unit cell of a crystalline material is defined by its lattice parameters (cell lengths a, b, and c and cell angles α, β, and γ) and the fractional positions of the atoms (Jenkins & Snyder, 1996). A transformation matrix can then be defined to go from fractional coordinates to Cartesian coordinateswith basis vectors , , and . The transformation matrix can be separated into a symmetric part representing the strain and an antisymmetric part representing the rotation of the unit cell, respectively (Nos & Klein, 1983).

Since force field potentials are usually defined in Cartesian space, it is convenient to store positions in Cartesian space. Distance vectors are transformed to fractional space and periodic boundary conditions are applied to it, and then converted back to Cartesian spacewith the rint‐function returning the rounded integer value of its argument.

For an infinitesimal homogeneous strain, the strain tensor can be expressed in terms of the transformation matrix by (Dubbeldam, Krishna, & Snurr, 2009; Lutsko, 1989)with the unstrained transformation matrix, the inverse of the unstrained transformation matrix, the transpose of the strained transformation matrix, and the identity matrix. Expanding the internal energy in terms of the strain tensor results in the energy‐strain relations (Wallace, 1972)where are the adiabatic elastic constants. Heinen used this energy‐strain approach (Heinen et al. 2017) by straining IRMOF‐1 at 0 K along the a‐direction with an amount ε. These results are shown in Figure 2a. The asymmetry of the component reveals that contraction is more favorable than expansion. Figure 2b illustrates that too large strains (ε > 0.012) can occur outside the harmonic regime. Using the energy‐strain relation finite‐temperature elastic constants (T) can be computed. Finite‐temperature elastic constants can also be computed with two other methods using molecular dynamics (MD) or Monte Carlo (MC) (Lutsko, 1989; Ray, 1988; Workum & Pablo, 2003).

Figure 2

(a) Computing the elastic constants of isorecticular MOF series I (IRMOF‐1) from energy‐strain curves: Symmetric strain (in orange) and asymmetric strain (in green), (b) values as a function of polynomial fit range (the converged value is obtained for strains smaller than 1% here; the line denotes the value obtained from Equation 6). Inset shows the structure before and after applying the strain. (Reprinted with permission from Heinen, Burtch, Walton, and Dubbeldam. Copyright 2017 American Chemical Society.)

One method is via the strain‐fluctuation formula introduced by Rahman and Ray (Parrinello & Rahman, 1982; Ray, 1988)where 〈V〉 denotes the average unit cell volume and 〈ɛ〉 are the strain fluctuations in the (HtN) (Ray & Rahman, 1984) or (TtN) (Ray & Rahman, 1985) ensembles. The unit cell is allowed to change shape by the method of Parrinello and Rahman (1980). The other method is, at zero stress and constant volume, the stress‐fluctuation formula (Ray & Rahman, 1984; Ray & Rahman, 1985)where the first term is simply the Born term calculated at 0 K, the second term is the fluctuation term, and the last term is the temperature‐correction term. Generally, the stress‐fluctuation method converges faster than the strain‐fluctuation method since the Born term in the stress fluctuation method dominates and does not depend on fluctuations (Fay & Ray, 1992; Gao, Workum, Schall, & Harrison, 2006). All terms appearing in the stress‐fluctuation formula can be obtained from a single MD simulation run. Note, that there is still debate in literature on the correct and unique definition of elastic constants at finite pressures (Kamb, 1961; Marcus, Ma, & Qiu, 2002; Marcus & Qiu, 2004; Marcus & Qiu, 2009; Steinle‐Neumann & Cohen, 2004).

To simulate lattice dynamics at finite temperature and pressure, the isothermal‐isobaric (NpT)‐ensemble is often used. The first constant pressure–temperature simulation was conducted by Andersen (1980) using isotropic volume changes. For inhomogeneous systems (e.g., noncubic solids) isotropic volume changes might not be sufficient. For example, the large to narrow pore phase transition in MIL‐53 requires anisotropic unit cell shape changes. Methodology for unit cell shape changes was first pioneered by Parrinello and Rahman (1980, 1981) by extending the Lagrangian of the Andersen method. Note, that for Monte Carlo calculations, the isoenthalpic–isotension–isobar ensemble should be used (Fay & Ray, 1992). Recently, Rogge et al. (2015) compared three barostat coupling schemes for reproducing cell parameters and pressure–volume behavior. The Langevin (Feller, Zhang, Pastor, & Brooks, 1995; Quigley & Probert, 2004) and Martyna‐Tuckerman‐Tobias‐Klein (MTTK; Martyna, Tobias, & Klein, 1994; Martyna, Tuckerman, Tobias, & Klein, 1996) barostats gave similar results for reproducing lattice parameters. The bulk moduli of MIL‐53(Al) and IRMOF‐1 was an order of magnitude larger for the Berendsen barostat (Berendsen, Postma, van Gunsteren, DiNola, & Haak, 1984) as compared with the other barostats. Clavier et al. (2017) studied the effect of the relaxation time τ of the Nose‐Hoover barostat on the elastic constants computed from strain fluctuations. The authors observed a very strong dependence (changes of more than 50%) in the range 0.1–10.0 ps. At cryogenic temperature, acoustic modes are thermally excited which results in an T 3 dependency in the heat capacity and also in the volume, which is known as the Debye T 3 law (Kittel, 2005; Sayetat, Fertey, & Kessler, 1998; Takenaka, 2012). Typically, classical force fields do not account for this behavior.

If a comparison between a rigid and flexible force field is ought to be made, the reference structure should be obtained via a minimization procedure which utilizes the flexible force field model. The references bonds, bends, and torsion’s are imposed by the flexible force field mode and these do not have to be the same as those found by DFT minimized experimental determined structures (Leach, 2001).

Adsorption isotherms in the osmotic ensemble

A well‐known application for classical force field calculations is the computation of adsorption isotherms. Adsorption isotherms for rigid materials are computed in the grand‐canonical (μV T) ensemble (Dubbeldam et al., 2013; Nicholson & Parsonage, 1988). To account for framework flexibility, a different ensemble is needed. Graben and Ray (1991) and Ray and Graben (1990) proposed a new potential R(S, p, μ) in the early 1990s for a single‐component adiabatic open ensemble representing a system in contact with a pressure and chemical potential reservoir (Ray, 2005). Generalizing this potential for an isothermal two component, the probability distribution function is written as (Wolf, Lee, & Davis, 1993)with Λ being the thermal de Broglie wavelength defined by . This probability distribution function is equivalent to that of the osmotic ensemble Ω introduced by Mehta and Kofke in 1994 and has been used by many others (Banaszak, Faller, & de Pablo, 2004; Brennan & Madden, 2002; Coudert, Boutin, Jeffroy, Mellot‐Draznieks, & Fuchs, 2011; Coudert, Jeffroy, Fuchs, Boutin, & Mellot‐Draznieks, 2008). The osmotic ensemble is the Legendre transform of the two‐component Ray potential with respect to the entropy Ω = R − TS.

The acceptance probability for any of the four MC moves (insertion, deletion, translation, and volume change) is given by with being defined in Equation (9) as the trial state and the probability distribution for the current state.

It is known that the acceptance ratio of insertion moves involving high loadings and bulky adsorbates can drop significantly. Efficient schemes that tackle this problem are the continuous fractional component Monte Carlo (CFCMC) method by Shi and Maginn (2007) and the configuration bias continuous fractional component Monte Carlo (CBCFMC) by Torres‐Knoop, Balaji, Vlugt, and Dubbeldam (2014). For a mixture isotherm of xylenes in MTW, convergence is obtained using CBCFMC with about four times fewer production cycles that for configuration bias Monte Carlo (Torres‐Knoop et al., 2014). However, for flexible materials, the unit cell volume change MC‐move forms an additional problem. As the unit cell changes volume, all particles are displaced which is often energetic unfavorable (Coudert et al., 2011). The most common approach to allow for unit cell volume changes is by means of a hybrid MD/MC scheme (Chempath, Clark, & Snurr, 2003; Duane, Kennedy, Pendleton, & Roweth, 1987). Here, one of the Monte Carlo moves is a short molecular dynamics run that results in higher acceptance probability.

Dubbeldam et al. (2007) used the osmotic ensemble to assess the influence on the adsorption of CO2 in IRMOF‐1. Little to no influence was found, as the equilibrium positions of the atoms in IRMOF‐1 fluctuate around their equilibrium positions and large‐scale atomic rearrangements are absent in this structure. This was also found for CO2 adsorption in NH2‐MIL‐53(Al) (Garcia‐Perez, Serra‐Crespo, Hamad, Kapteijn, & Gascon, 2014) and HKUST‐1 (Zhao et al., 2011). Simulations in the osmotic ensemble show excellent agreement with experimental single‐component adsorption isotherm for IRMOF‐1 (Dubbeldam et al., 2007) and ZIF‐8 (Wu, Huang, Cai, & Jaroniec, 2014). Multicomponent adsorption of hydroxymethyl furfural‐water and furfural‐water mixtures in zeolite MFI in the osmotic ensemble showed better agreement of the experimental saturation capacity compared to simulations in the grand‐canonical ensemble (Santander, Tsapatsis, & Auerbach, 2013). Ghoufi et al. (2012) constructed a flexible force field that captured the lp to np and np to lp transition of carbon dioxide in MIL‐53(Cr). Also, the adsorption enthalpies in both phases were in good agreement with the experimental values.

PARAMETERIZATION SCHEMES

The first flexible MOF force field was reported in 2006 by Greathouse and Allendorf (2006) for the water decomposition of the prototypical IRMOF‐1. Structural collapse of the unit cell from 3.9% water content was observed and in agreement with experimental data (Huang et al., 2003). The authors argued that bonded Zn‐O interactions result in poor volume changes upon external stimuli due to bond and angle constraints. Therefore, Zn‐O interactions were modeled using nonbonding pair interactions. Organic linker parameters were taken from the generic CVFF force field (Dauber‐Osguthorpe et al., 1988); with slight modifications to better represent the experimental unit cell of IRMOF‐1. Inorganic parameters were obtained from DFT calculations on zincite (ZnO). The lattice parameter 25.61 Å was found to be in good agreement with those of reported x‐ray studies of 25.67–25.89 Å (Eddaoudi et al., 2002; Rowsell, Spencer, Eckert, Howard, & Yaghi, 2005). In 2008, the authors extensively validated the force field (Greathouse & Allendorf, 2008) and reported lattice parameters of 25.74 and 26.05 Å at 0 K by extrapolating finite‐temperature MD calculations and from energy minimizations. The overestimated lattice parameter from the energy minimization was solved by using the mode‐following (Baker, 1986) technique. This method eliminates negative eigenvalues to obtain a true minimum and resulted in a lattice parameter of a = 25.698 Å in agreement with the MD extrapolation. The space group was, however, reduced from Fm3 (255) to Fm (202) showing distortions of the metal‐linker node (Dubbeldam et al., 2009). Softening of the Young’s and bulk modulus occurred with increasing temperatures following general temperature‐dependent elasticity (Ledbetter, 2006). The 0 K force field values were 38.3% and 18.5% higher compared to the DFT values (Bahr et al., 2007).

Three other flexible force fields for IRMOF‐1 were introduced in 2007 by Schmid et al., Dubbeldam et al., and Han and Goddhard III. Schmid and coworkers (2007) used a building block approach for the parameterization of IRMOF‐1 based on the MM3(2000) force field. DFT‐B3LYP calculations on a tetranuclear zinc benzoate cluster (Zn4O[O2CC6H5]6) were used as a reference system. Selected bond distances and vibrational frequencies were compared with experimental data (Clegg et al., 1991; Ming‐cai, Chi‐wei, Chang‐chun, Liang‐jie, & Ju‐tang, 2004). Even though bond distances were slightly overestimated, the authors argued that the reproduction of the vibrational frequencies was more important for the lattice dynamics. By transforming the Hessian matrix, a set of force constants was obtained (excluding translation and rotational motion). Unlike, the Greathouse‐Allendorf force field, the Zn‐O bonds were considered partially bonded. To reproduce the characteristic asymmetric stretch of Zn4O at around 530 cm−1 (Hermes, Schroder, Amirjalayer, Schmid, & Fischer, 2006), the off‐diagonal terms (presenting coupling interactions) were considered for the bond‐bend Zn‐Ocent/Zn‐Ocent‐Zn interactions. An additional Zn4(O2CH)5‐BDC‐BDCZn4O(O2CH)5 cluster model was chosen for the internal torsion barrier of BDC linker.

Dubbeldam et al. (2007) reparameterized the Greathouse and Allendorf force field by reproducing the experimental lattice parameter and CO2 and methane adsorption isotherms. Here, the Zn‐O interactions are also considered nonbonding. A later refinement (Dubbeldam et al., 2009) mapped the 0 K elastic tensor on the ab initio calculated tensor (Samanta, Furuta, & Li, 2006). Dubbeldam et al. (2007) discovered NTE in MOFs and demonstrated NTE in MOFs as a general phenomenon associated with the struts‐and‐spring, large‐pore nature of MOFs. The classical model predicted large NTE which was later confirmed by experiments of Zhou, Wu, Yildirim, Simpson, and Walker (2008) using neutron powder diffraction and first‐principles calculations using minimization of the free energy. The underlying reason for NTE in nanoporous materials is different from condensed matter. In MOFs, the linker molecule undergoes transverse motions due to thermal vibrations such that at lower temperatures the linker becomes more rigid and stretches out. NTE is related to relatively rigid linkers connected to rigid metal clusters via flexible groups such as carboxylate groups and the porosity of the structure that allows adequate volume for motion. Therefore, it is expected that NTE could be a common phenomenon in MOFs.

Han and Goddard III (2007) used the DREIDING force field (Mayo, Olafson, & Goddard, 1990) to elucidate the mechanism of NTE. They argued that the amplitude of the rotations of the Zn‐O clusters is increased at higher temperatures as well as the rotations of the organic linker. Despite contraction of the unit cell upon heating, the finite‐temperature elastic constants, obtained from the energy‐strain relation, decreased. This is in agreement with findings of Greathouse and Allendorf and of recent DFT work (Banlusan & Strachan, 2017). The above mentioned force fields has been summarized in Table 1 and show good agreement with experiments and ab initio results.

Table 2 presents the method of parameterizing, charge scheme, and the metal‐linker interaction of various generic flexible force fields. The building block methodology is a popular approach utilizing ab initio calculations on nonperiodic clusters to fit parameters. QuickFF (Vanduyfhuys et al., 2015) is such an example. Figure 3 shows two clusters for the parameterization of MIL‐53(Al): one for the organic linker and one for inorganic node. Parameterization is divided per cluster into three steps: (a) determination of dihedral reference angles and multiplicities, (b) extraction of force field parameters for an internal coordinate, and (c) refinement of force constants and fitting of missing dihedral force constants. When functional forms occur in both clusters for one material, their parameters are averaged.

Table 2

Methods of parameterizing, atomic charge partition schemes, and metal‐linker interaction of generic parameterization schemes for flexible metal–organic frameworks

	Method	Charge scheme	Metal‐linker interaction	Reference
QuickFF	Clusters	Point and Gaussian chargesa	Bondingb	Vanduyfhuys et al. (2015)
MOF‐FF	Clusters	Gaussian charges	Morse and bendingc	Bureekaew et al. (2013)
UFF4MOF	Clusters	Universal force field (UFF; Rappe, Casewit, Colwell, Goddard, & Skiff, 1992)	Bonding	Addicoat, Vankova, Akter, and Heine (2014) and Coupry, Addicoat, and Heine (2016)
BTW‐FF	Periodic	Effective chargesd	Bonding and nonbonding	Bristow et al. (2014)
VMOF	Phonons	Charge equilibration and formale	Buckinghamf and Coulomb	Bristow, Skelton, Svane, Walsh, and Gale (2016)
HBWD	Elasticity	REPEAT	Bonding, bending, and torsion	Heinen et al. (2017)

point charges, else Gaussian charges.

Bonds, bends, and torsions.

Morse: , bend: .

Topological analysis of Bloch states.

Ligands charges: equilibration scheme of Gasteiger and Marsili (1978) and Gasteiger and Marsili (1980), metal nodes and inorganic oxygen: formal charges.

Modified MM3 Buckingham potential with A = 1.84·105, B = 12 and C = 2.25.

Figure 3

Cluster models representing MIL‐53(Al) used in QuickFF. (Reprinted with permission from Vanduyfhuys, Verstraelen, Vandichel, Waroquier, and Van Speybroeck. Copyright 2012 American Chemical Society.)

MOF‐FF, a third generation force field of Schmid and coworkers (2013), is also based on the building block methodology. The first generation (Tafipolsky et al., 2007) has been discussed on the preceding page. In the second generation (Tafipolsky & Schmid, 2009), a generic algorithm was used to derive force field potentials for the zinc benzoate and dilithium terephthalate. To keep the Zn atoms in the carboxylate plain, the Zn‐dependent torsion Zn‐Ocarb‐Ccarb‐Ocarb was included. In the third generation (Bureekaew et al., 2013), a new energy expression was proposed as well as a new parameterization scheme for Cu and Zn paddlewheel‐based structures and Zn4O and Zr6(OH)4‐based MOFs.

UFF4MOF (Addicoat et al., 2014; Coupry et al., 2016) is an extension of the original Universal force field (UFF; Rappe et al., 1992) that incorporates additional atom types found in the Computation‐Ready Experimental (CoRE) Database. (Chung et al., 2014) The only fitted parameter is the bond (covalent) radius of the atom types which is obtained by minimizing the residual error of a training set of secondary building units from gas phase ab initio calculations. Only 3 of the 19 minimized MOF structures had calculated lattice parameters of larger than 4% compared with the experimental values.

By extending the BTW‐FF model (Bristow et al., 2014), Gale and coworkers (2016) developed the vibrational metal–organic framework (VMOF) force field. Here, force field parameters are explicitly fitted on DFT‐PBEsol calculated phonon spectra of periodic binary oxides such as ZnO, ZrO2, TiO2, and Al2O3. A modified MM3 Buckingham potential for the metal‐linker interaction was needed to reproduce the ab initio and experimental structural and mechanical properties of the binary oxides more accurately. Inaccurate long‐range dispersion interactions resulted in considerable lower bulk moduli of various MOFs calculated with VMOF compared to DFT values. The lower frequency region could not be accurately reproduced due to discrepancies in the metal‐oxygen stretching modes and due to the use of formal charges as was argued by the authors (Bristow et al., 2016). It must be noted that it is challenging to obtain accurate THz data of the vibrational spectra using experimental techniques or ab initio calculations (Ryder et al., 2014).

Alternatively, one can fit functional forms of metal‐linker interface on the elastic tensor (Dubbeldam et al., 2007, 2009; Heinen et al., 2017). By exploring the parameter space of a defined force field, missing functional forms can be detected if component of the elastic tensor are of. By having a properly defined interaction for the metal‐linker interface that captures the elasticity of the entire unit cell, the force field extension toward functionalized MOFs (substituted linkers) should be straightforward under the assumption that the functionalization does not heavily affect the metal‐linker interface. The bonding parameters of organic functional groups are well parameterized in generic force fields, such as OPLS‐AA (Jorgensen et al., 1996; Jorgensen & Tirado‐Rives, 1988), AMBER (Weiner et al., 1984; Weiner & Kollman, 1981), and CVFF (Dauber‐Osguthorpe et al., 1988).

CHALLENGES AND PROSPECTS

Frequently, supercells are used in force field calculations in order to satisfy the minimum image convention. However, finite size effects regarding finding a global minimum and substitution patterns in functionalized linkers should in some cases also be considered (Heinen & Dubbeldam, 2016; Munn et al., 2016). Figure 4 shows classical optimized DMOF structure using the mode‐following minimization technique (Baker, 1986). It was found, after trial‐and‐error, that the lowest energy structure has two sequentially DABCO linkers in the staggered configuration along the b‐direction (Burtch, Jasuja, Dubbeldam, & Walton, 2013; Grosch & Paesani, 2012). The use of a single unit cell would artificially impose a strain on the unit cell.

Figure 4

Classically optimized DMOF with a = c = 15.483 Å and b = 19.283 Å (space group: I4/mcm). The DABCO linkers are along the b‐direction are in the staggered configuration

An important challenge is the development of efficient sampling schemes for large‐scale structural transformations (Demuynck et al., 2017). Currently, it takes on the order of hundreds of picoseconds before a lp to np phase transition is observed in MIL‐53(Al) using standard molecular dynamics (Rogge et al., 2015). This timescale is too long to sample the phase space under varying conditions, especially when adsorbates are present. For studying the lattice dynamics, umbrella sampling (Frenkel & Smit, 2002) is a useful approach to bias the free energy barrier. Here, an order parameter λ is defined that describes the transition from the lp to the np phase and vice versa, as is shown schematically in Figure 5. The free energy barrier is however loading dependent, meaning that for each loading a separate bias potential needs to be constructed.

Figure 5

Any movement from a large pore phase to a narrow pore phase can be described in terms of a parameter λ that ranges from 0 to 1 describing the progress of the transition. Using Umbrella sampling, the free energy barrier can be biased away

Another problem to be solved is the connection between atomistic simulations and continuum mechanical simulations. More specifically, which time and length scales are needed in order to recover results consistent with continuum mechanical theories? (Sengupta, Nielaba, Rao, & Binder, 2000) In continuum mechanics, the deformation gradient is a fundamental quantity that describes strain at large deformation but remains ill‐defined at the atomic level. The deformation gradient is defined aswith  = {X 1, X 2, X 3} being a point in the reference framework and  = () = (X 1, X 2, X 3) being a point on the deformation framework. Zimmerman, Bammann, and Gao (2009) proposed an atomistic‐scale mean‐value deformation gradient relation and showed that the established metric is consistent with continuum mechanics, namely that the curl of the deformation gradient is zero (known as the compatibility of the deformation field). Nonzero values of the curl of the deformation gradient were observed for materials containing defects but it remains unclear if this can be correlated to geometric information.

An interesting question which can be answered with force fields is if the elastic tensor predicts thermal expansion or vice versa? For axial solids, the principal coefficients of thermal expansion can be related to the elastic compliance ( = ) by (Munn, 1972).with C the heat capacity under constants stress, V the unit cell volume, and γ the directional Grüneisen functions. This suggests that a correct description of the elastic constants is crucial for capturing thermal expansion effects. However, the opposite cause–effect relation has also been suggested (Wang et al., 2010).

Last, we would like to point out the increasing use and development of machine learning techniques (McCall, 2005). Evolutionary approaches, such as genetic algorithms, are well suited for nonlinear high‐dimensional optimization issues that are involved in force field parameterization (Tafipolsky & Schmid, 2009). Other examples include the determination of polarizable force field parameters from ab initio data (Li et al., 2017), prediction of mechanical properties for zeolites (Evans & Coudert, 2017), identification of molecular characteristics that give rise of porous crystals (Evans et al., 2016), constructing atomic forces using Bayesian interference, or on‐the‐fly ab initio molecular dynamics simulations (Botu, Batra, Chapman, & Ramprasad, 2017; Huan et al., 2017; Li, Kermode, & De Vita, 2015).

CONFLICT OF INTEREST

The authors have declared no conflicts of interest for this article.

RELATED WIREs ARTICLES

https://doi.org/10.1002/wcms.1282