**2.2 Models for MOFs**

When the first MD simulations were performed to examine gas diffusion in MOFs at the beginning of 2004, there was no experimental data to validate the accuracy of MD studies. However, in general whenever experimental equilibrium properties such as adsorption isotherms have been reproduced by the molecular simulations, it has been observed that dynamic simulations based on the same interatomic potentials are also reliable. Therefore, many MD studies examining gas diffusion in MOFs first showed the good agreement between experiments and simulations for gas adsorption isotherms and then used the same potential models for gas diffusion simulations.(Skoulidas&Sholl, 2005) Here, it is useful to highlight that considering a wide gas loading range when comparing simulation results with experimental data is crucial.(Keskin et al., 2009b) It is unreasonable to compare outcome of simulations with the experimental measurements over a very narrow range of loading and assume that good (or poor) agreement with experiment will continue to high loadings.

The MD simulations have used general-purpose force fields such as the universal force field (UFF) (Rappe et al., 1992), DREIDING force field(Mayo et al., 1990) and optimized potential for liquid simulations all-atom (OPLS-AA)(Jorgensen et al., 1996) force field for representing the interactions between MOF atoms and adsorbates. A few studies have used quantum

used. In most of the MD simulations, spherical 12-6 Lennard-Jones (LJ) model (Buch, 1994) has been used for H2. The Buch potential is known to reproduce the experimental bulk equation of state accurately for H2. Two-site LJ models have also been used in the literature.(Yang&Zhong, 2005) The potential model of Darkrim and Levesque (Darkrim&Levesque, 1998) has been used to account for the quadrupole moment of H2 molecules. This potential consists of a LJ core placed at the center of mass of the molecule and point charges at the position of the two protons and the center of mass.(Liu et al., 2008b) Methane and argon diffusion simulations have been performed using single site-spherical LJ potentials. The CO2 potential consists of three LJ sites with charges located on each site to represent the quadrupole moment of CO2.(Potoff&Siepmann, 2001) In N2 diffusion simulations, N2 is generally represented as a three site model with two sites located at two N atoms and the third one located at its center of mass (COM) with partial point charges.(Potoff&Siepmann, 2001) For alkanes in MOFs, the TraPPE potential has been used.(Martin&Siepmann, 1997) As an example, the most widely used potential parameters

(ε: energy parameter, *kB*: Boltzmann constant, σ: size parameter) are listed in Table 1.

Atoms/Molecules ε/*kB* (K) σ (Å) References H2-H2 34.20 2.96 Buch, 1994 Ar-Ar 119.8 3.4 Clark, 1998 CH4-CH4 148.20 3.73 Martin&Siepmann, 1997 C-C (in CO2) 27.00 2.80 Potoff&Siepmann, 2001 O-O (in CO2) 79.00 3.05 Potoff&Siepmann, 2001 N (in N2) 36.4 3.32 Potoff&Siepmann, 2001 COM (N2) 0.00 3.32 Potoff&Siepmann, 2001

Table 1. Potential parameters used for adsorbate molecules in MD simulations.

When the first MD simulations were performed to examine gas diffusion in MOFs at the beginning of 2004, there was no experimental data to validate the accuracy of MD studies. However, in general whenever experimental equilibrium properties such as adsorption isotherms have been reproduced by the molecular simulations, it has been observed that dynamic simulations based on the same interatomic potentials are also reliable. Therefore, many MD studies examining gas diffusion in MOFs first showed the good agreement between experiments and simulations for gas adsorption isotherms and then used the same potential models for gas diffusion simulations.(Skoulidas&Sholl, 2005) Here, it is useful to highlight that considering a wide gas loading range when comparing simulation results with experimental data is crucial.(Keskin et al., 2009b) It is unreasonable to compare outcome of simulations with the experimental measurements over a very narrow range of loading and assume that good (or poor) agreement with experiment will continue to high

The MD simulations have used general-purpose force fields such as the universal force field (UFF) (Rappe et al., 1992), DREIDING force field(Mayo et al., 1990) and optimized potential for liquid simulations all-atom (OPLS-AA)(Jorgensen et al., 1996) force field for representing the interactions between MOF atoms and adsorbates. A few studies have used quantum

**2.2 Models for MOFs** 

loadings.

mechanical calculations to develop new potentials for specific MOF-adsorbate interactions.(Bordiga et al., 2005; Sagara et al., 2004) There are studies where the parameters of the force fields are refined to match the predictions of simulations with the experimental measurements (in most cases experimental adsorption isotherm data exist whereas experimental diffusion data do not exist) or using first principles calculations.(Sagara et al., 2004; Sagara et al., 2005a; Sagara et al., 2005b; Yang&Zhong, 2005; Yang&Zhong, 2006). Of course, one must be careful in refining force field parameters to match the results of simulations with the experimental data since the accuracy of the experiments are significantly affected by the defects of as-synthesized MOFs or trapped residual solvent molecules present in the samples.(Keskin et al., 2009b)

Most MD simulations performed to date have assumed rigid MOF structures which means the framework atoms are fixed at their crystallographic positions. Generally, the crystallographic data for MOFs are obtained from X-ray diffraction experiments. In rigid framework simulations, only the nonbonding parameters, describing the pair wise interactions between the adsorbate and the adsorbent atoms of the particular force field, were used. It can be anticipated that the assumption of a rigid framework brings a huge computational efficiency yet the inclusion of the lattice motion and deformation is crucial for an accurate description of diffusion of large gas molecules since they fit tightly in the MOF pores, forcing the MOF to deform in order to allow migration from pore to pore.

In order to include the lattice dynamics in MD simulations of gas diffusion in MOFs, Tafipolsky and coworkers extended the MM3 (Molecular Mechanics) force field (Lii et al., 1989), which is well known to accurately describe structures and conformational energies of organic molecules, with parameters for the Zn4O moiety based on the first principles of density functional theory (DFT) calculations of nonperiodic model systems.(Tafipolsky et al., 2007) After this force field accurately predicted the structure of MOF-5 (also known as isoreticular MOF-1, IRMOF-1), it was used in MD simulations to investigate the self diffusion of benzene in MOF-5.(Amirjalayer et al., 2007) The self diffusivity calculated from MD simulations, 2.49×10-9 m2/s, corresponds well to the experimental value determined by Stallmach et al., who found values between 1.8-2×10-9 m2/s.(Stallmach et al., 2006) Under identical conditions, MD simulations performed with a rigid MOF lattice gave substantially higher diffusion coefficient, 19.5×10-9 m2/s, which is almost one order of magnitude larger than the value obtained from MD simulations with flexible lattice.

Amirjalayer et al. later investigated the diffusion mechanism of benzene in MOF-5 using MD simulations based on fully flexible MM3 force field and computed the self diffusivity of benzene as a function of loading in MOF-5.(Amirjalayer&Schmid, 2009) The results were close to the experimentally determined self diffusivity for liquid benzene, therefore they concluded that MOF-5 is a very open structure with liquid like mobilities for benzene, ~2.49×10-9 m2/s. The third study including flexibility of MOFs in MD simulations used the modified MM3 force field to study hexane's self diffusion mechanism in IRMOF-1 and IRMOF-16 and concluded that the flexibility of IRMOF-16 is much larger than that of IRMOF-1 due to the nature of the organic linkers.(Xue&Zhong, 2009)

Greathouse and Allendorf developed a flexible hybrid force field for MD simulations of IRMOF-1(Greathouse&Allendorf, 2006) and the activation energy for benzene self diffusion calculated at low loadings using this force field was found to be in good agreement with previous MD simulations and nuclear magnetic resonance (NMR) results. (Greathouse

Recent Advances in Molecular Dynamics

and bulk phase fugacity, *f*

Simulations of Gas Diffusion in Metal Organic Frameworks 261

thermodynamic correction factor, a partial derivative relating the adsorbate concentration, *c*

ln () () . ln *t o <sup>f</sup> Dc Dc*

The thermodynamic correction factor is fully defined once the single component adsorption isotherm is known. Well developed approaches exist for calculating the corrected diffusion coefficient from MD simulations.(Kärger&Ruthven, 1992; Keil et al., 2000; Skoulidas&Sholl, 2003; Skoulidas&Sholl, 2005) For systems with a single adsorbed component, the corrected diffusivity is equivalent to the Maxwell-Stefan diffusion coefficient.(Kapteijn et al., 2000; Ruthven, 1984; Sholl, 2006) The corrected diffusivity includes information on the collective motion of multiple adsorbed molecules that is relevant to net mass transport and can be

> 1 <sup>1</sup> lim ( ) (0) <sup>6</sup> *Ni*

Here, *N* is the number of molecules, r*il*(*t*) is the three dimensional position vector of molecule *l* of species *i* at time *t* and the angular brackets denote that the ensemble average. Using MD simulations, one can record the trajectory of the gas molecules in the pores of MOFs and calculate the corrected diffusivity. A more microscopic measure of diffusion is the self diffusion coefficient which describes the motion of individual, tagged particles. In an isotropic three dimensional material, the self diffusivity is related to the mean squared

calculated using the following expression:(Kärger&Ruthven, 1992; Keil et al., 2000)

*o i il il <sup>t</sup> <sup>l</sup> <sup>D</sup> rt r Nt*

,

displacement of tagged particles by the Einstein relation:

,

*c*

1 1 1 lim ( ) (0) <sup>6</sup> *Ni*

*self i il il <sup>t</sup> <sup>t</sup> <sup>l</sup> <sup>D</sup> rt r t N*

simulations based on the method by Theodorou et al. (Theodorou et al., 1996):

This definition of self diffusivity is applicable to both single component and multicomponent systems.(Sanborn&Snurr, 2000) In general, all three diffusion coefficients described here, transport, corrected and self diffusivities are the functions of concentration and they are only equal in the limit of dilute concentrations.(Sholl, 2006) In some extreme cases, the self and corrected diffusivities vary by orders of magnitude.(Ackerman et al., 2003; Skoulidas et al., 2002) This observation sometimes underscores the value of characterizing these two diffusivities independently. Applications such as modeling of membranes, pressure swing adsorption require the accurate description of net mass transfer and in these processes generally the transport diffusivity is of greatest interest.(Sholl, 2006) Almost all applications of nanoporous materials in gas separations involve chemical mixtures; therefore it is important to describe the multi-component gas transport in nanopores. There are several mathematically equivalent formalisms such as Onsager, Fickian and Maxwell-Stefan to describe multi-component gas transport through nanoporous materials.(Krishna&van den Broeke, 1995; Wesselingh&Krishna, 2000) The Onsager formulation is based on irreversible thermodynamics and expresses the flux of each species in terms of chemical potentials. One can calculate the Onsager coefficient using MD

2

2

(3)

(2)

(1)

& Allendorf, 2008) In this force field, only nonbonded parameters were used to describe Zn−O interactions and the CVFF (consistent valence force field) was used with slight modifications to describe the benzene dicarboxylate linker. The magnitude of the diffusion constant was underestimated and this was attributed to the deficiencies in the CVFF portion of the force field.

The literature summary presented so far indicates that the number of MD simulation studies with flexible MOFs and flexible force fields is very limited. More research will sure be helpful to understand the importance of lattice dynamics on diffusivity of gas molecules in MOFs. Studies to date indicated that the lattice dynamics are specifically important in computing diffusivity of large gas molecules (such as benzene) in MOFs having relatively narrow pores. Studies on flexible force fields also suggested that a force field developed for a specific MOF can be adapted to similar MOF structures (as in the case of IRMOFs) with slight modifications for doing comparative studies to provide a comprehensive understanding of gas diffusion in flexible MOFs.

One major issue in carrying out MD simulations for MOFs is to assign partial charges to the MOF atoms that are required to calculate adsorbate-adsorbent interactions for some polar (including quadrupolar) adsorbates. Several MD studies computed the diffusivity of CO2 in MOFs' pores and to do this, partial charges must be assigned to MOF atoms. Recent studies showed that the effects of inclusion of framework charges are crucial at low loadings. If the charge-quadrupolar interactions are not taken into account in MD simulations then the diffusivities can be significantly overestimated.(Rankin et al., 2009) Force field-based classical MD simulations of MOFs typically treat electrostatic interactions between adsorbates and MOF atoms by assigning fixed point charges to each atom. In this context an important role for quantum mechanics (QM) calculations is to assign the point charges that can later be used in force field calculations. Unfortunately, multiple methods exist for partitioning the net electron density determined in a QM calculation(Keskin et al., 2009b) and none of these methods give an unambiguous definition of the resulting point charges. Keskin and coworkers reviewed the partial charges assigned to IRMOF-1 on the basis of QM calculations and showed that there is a significant variation in the charge values based on the method used.(Keskin et al., 2009b) This variation may have a significant impact on the outcome of classical force field calculations in examples where electrostatic interactions are important. Since QM calculations are time consuming and the charges obtained from these calculations are method sensitive, a strategy called connectivity based atom contribution method (CBAC) with which the partial charges of framework atoms can be estimated easily was proposed.(Xu&Zhong, 2010) A recent study on two different MOFs showed that CO2 adsorption isotherms and diffusivities computed using the charges from QM methods based on the ChelpG(Francl et al., 1996) DFT calculations are very similar to the ones computed using charges from CBAC method.(Keskin, 2011a)
