Recent Advances in Molecular Dynamics Simulations of Gas Diffusion in Metal Organic Frameworks

Over approximately the last decade, metal organic framework (MOF) materials have attracted a great deal of attention as a new addition to the classes of nanoporous materials. MOFs, also known as porous coordination polymers (PCPs) or porous coordination networks (PCNs), are hybrid materials composed of single metal ions or polynuclear metal clusters linked by organic ligands through strong coordination bonds. Due to these strong coordination bonds, MOFs are crystallographically well defined structures that can keep their permanent porosity and crystal structure after the removal of the guest species used during synthesis.(Eddaoudi et al., 2000; Li et al., 1999; Rowsell et al., 2005; Yaghi et al., 2003) MOFs typically have low densities (0.2-1 g/cm3), high surface areas (500-4500 m2/g), high porosities and reasonable thermal and mechanical stabilities. This combination of properties has made MOFs interesting materials for a wide range of potential applications, including gas storage, gas separation, catalysis and biomedical applications.(Eddaoudi et al., 2002; KeskinK MillwardY Mueller et al., 2006; Pan et al., 2004)


Introduction
Over approximately the last decade, metal organic framework (MOF) materials have attracted a great deal of attention as a new addition to the classes of nanoporous materials.MOFs, also known as porous coordination polymers (PCPs) or porous coordination networks (PCNs), are hybrid materials composed of single metal ions or polynuclear metal clusters linked by organic ligands through strong coordination bonds.Due to these strong coordination bonds, MOFs are crystallographically well defined structures that can keep their permanent porosity and crystal structure after the removal of the guest species used during synthesis.(Eddaoudi et al., 2000;Li et al., 1999;Rowsell et al., 2005;Yaghi et al., 2003) MOFs typically have low densities (0.2-1 g/cm 3 ), high surface areas (500-4500 m 2 /g), high porosities and reasonable thermal and mechanical stabilities.This combination of properties has made MOFs interesting materials for a wide range of potential applications, including gas storage, gas separation, catalysis and biomedical applications.(Eddaoudi et al., 2002;Keskin&Kizilel, 2011;Millward&Yaghi, 2005;Mueller et al., 2006;Pan et al., 2004) MOFs have become attractive alternatives to traditional nanoporous materials specifically in gas storage and gas separation since their synthesis can be readily adapted to control pore connectivity, structure and dimension by varying the linkers, ligands and metals in the material.(Düren et al., 2004;Eddaoudi et al., 2002;El-Kaderi et al., 2007) Hundreds of MOF materials with various physical and chemical characteristics have been synthesized to date.(James, 2003;Kitagawa et al., 2004;Uemura et al., 2005;Yaghi et al., 2003) Most of the studies in the literature have focused on a few specific MOF groups such as IRMOFs (isoreticular MOFs) (Eddaoudi et al., 2002), ZIFs (zeolite imidazolate frameworks) (Park et al., 2006), CPOs (coordination polymers of Oslo) (Dietzel et al., 2005;Dietzel et al., 2006), MILs (Materials of the Institute Lavoisier) (Loiseau et al., 2004), CuBTC (Copper 1,3,5benzenetricarboxylate) (Chui et al., 1999) and Zn(bdc)(ted) 0.5 (Zinc 1,4-benzenedicarboxylic acid-triethylenediamine) (Li et al., 1998).As an example Figure 1 shows the unit cell structure of one of the most widely studied MOFs, CuBTC (also known as HKUST-1).The figures from left to right represent an empty CuBTC structure, a CuBTC structure with CH 4 molecules in the pores at 100 bar, 298 K and a CuBTC structure with adsorbed CH 4 and H 2 molecules in the pores at 100 bar, 298 K for a bulk gas composition of CH 4 /H 2 :5/95.The enormous number of different possible MOFs indicates that purely experimental means for designing optimal MOFs for targeted applications is inefficient at best.Efforts to predict the performance of MOFs using molecular modeling play an important role in selecting materials for specific applications.In many applications that are envisioned for MOFs, diffusion behavior of gases is of paramount importance.Applications such as catalysis, membranes and sensors cannot be evaluated for MOFs without information on gas diffusion rates.Most of the information on gas diffusion in MOFs has been provided by molecular dynamics (MD) studies.Figure 2 indicates that the idea of using MD simulations to assess the diffusivity of gases in MOFs is a new area and there is a rapid growth in the number of publications featuring the terms 'MOFs', 'MD' and 'diffusion' over the past decade.used.In most of the MD simulations, spherical 12-6 Lennard-Jones (LJ) model (Buch, 1994) has been used for H 2 .The Buch potential is known to reproduce the experimental bulk equation of state accurately for H 2 .T w o -s i t e L J m o d e l s h a v e a l s o b e e n u s e d i n t h e literature.(Yang&Zhong, 2005) The potential model of Darkrim and Levesque (Darkrim&Levesque, 1998) has been used to account for the quadrupole moment of H 2 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 CO 2 potential consists of three LJ sites with charges located on each site to represent the quadrupole moment of CO 2 .(Potoff&Siepmann, 2001) In N 2 diffusion simulations, N 2 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)

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 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 Zn 4 O 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 m 2 /s, corresponds well to the experimental value determined by Stallmach et al., who found values between 1.8-2×10 -9 m 2 /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 m 2 /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 m 2 /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 & 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 CO 2 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 CO 2 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)

Predicting gas diffusivity using molecular dynamics simulations
In the literature, MD simulations have been used to predict three different types of gas diffusivities in MOFs.These are transport diffusivity, corrected diffusivity and self diffusivity.The transport diffusivity, which is also known as Fickian diffusivity or chemical diffusivity, can be defined without approximation in terms of corrected diffusivity, D o and a 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 calculated using the following expression: (Kärger&Ruthven, 1992;Keil et al., 2000)  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 displacement of tagged particles by the Einstein relation: 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&vanden 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 simulations based on the method by Theodorou et al. (Theodorou et al., 1996): In this formulation, V is the subsystem volume, k B is the Boltzmann constant, T is temperature, r il (t) is the three-dimensional position vector of molecule l of species i at time t and N i is the number of molecules of species i.The Onsager coefficients and the matrix of Fickian coefficients are mathematically equivalent and they are related to each other without approximation by expressions involving derivatives of the mixture adsorption isotherm for the adsorbed species.(Skoulidas et al., 2003) The Onsager coefficients from MD simulations can be converted to Fickian diffusion coefficients using the followings: In these equations, T is temperature, c i is the concentration of species i, f j is the fugacity of species j, k B is Boltzmann constant, L ij and L ik are the Onsager coefficients and D ii and D ij are the Fickian diffusivity coefficients.Using Onsager or Fickian diffusivities, one can calculate the flux (J) of a binary gas mixture through a membrane as follows: In these expressions, i c  and i   represent concentration gradient and chemical potential gradient of species i through the membrane, respectively.As Equation 9 suggests, gas fluxes in a MOF membrane can be calculated based on either of the formulations (Onsager or Fickian).

Single component diffusion
The transport rates of single component gas molecules inside the materials' pores are important in many potential applications of MOFs.For example, in equilibrium-based separations such as pressure swing adsorption, transport rates define limits on the cycle times that can be achieved.In these cases, molecular transport rates are mainly important if they are very slow.Since accurate characterization of molecular transport inside nanoporous materials using experiments is very challenging, most of the information that is currently available about single component gas diffusion in MOFs has obtained from MD simulations.
Skoulidas performed the first study of gas diffusion in a MOF material in the literature using equilibrium MD simulations and calculated the self, corrected and transport diffusivities of argon at 298 K as a function of pressure.(Skoulidas, 2004) Results showed that diffusion in CuBTC MOF is an activated process as in zeolites.The calculated diffusivities of Ar in CuBTC were similar to the diffusion in zeolites both in magnitude and concentration dependence.Sarkisov et al. used equilibrium MD simulations to calculate the self diffusivities of methane, n-pentane, n-hexane, n-heptane, cyclohexane and benzene in MOF-5 at 300 K at dilute loadings.(Sarkisov et al., 2004) They found that self diffusivities of n-alkanes in MOF-5 are comparable to those in the crystalline bipyridine system (0.1-3×10 -8 m 2 /s), but they show a stronger dependence on chain length because of the more open structure of MOF-5.Skoulidas and Sholl(Skoulidas&Sholl, 2005) then used equilibrium MD simulations to probe the self, corrected and transport diffusivity of a number of gas species, Ar, CH 4 , CO 2 , N 2 and H 2 in MOF-5 as a function of pore loading at room temperature.They also calculated self, corrected and transport diffusivities of Ar in MOF-2, MOF-3, MOF-5, CuBTC and MFI to make a comparison among different MOFs and a zeolite.They concluded that diffusion of gas molecules in MOFs is mostly dominated by motions where the adsorbed species remain in close contact with the surfaces defined by the pore structure throughout their diffusion.At the time of their study, there was no experimental data for gas diffusion in MOFs.Therefore, Skoulidas and Sholl could not directly comment on the accuracy of the MD simulations, however, they pointed out that using similar molecular simulation methods for gas diffusion in zeolites they got excellent agreement with the experiments.
Yang and Zhong performed constant temperature equilibrium MD simulations by a momentum scaling method to calculate the self diffusivity of H 2 in isoreticular MOFs, IRMOF-1, IRMOF-8 and IRMOF-18 as a function of pressure.(Yang&Zhong, 2005) Their results showed that the diffusivity of H 2 in IRMOFs is slightly larger than in zeolites due to the larger pore volume of IRMOFs.The self diffusivity of H 2 in IRMOFs at 77 K at low pressures was around 1-3×10 -8 m 2 /s whereas the self diffusivity of H 2 in various zeolites were experimentally measured to be around 0.1-1×10 -8 m 2 /s.The activation energy of H 2 in IRMOFs was between 2-3 kJ/mol which is close to the values measured in zeolite NaX (4 kJ/mol (Bär et al., 1999)) and single walled carbon nanotubes (1.12 kJ/mol (Narehood et al., 2003)).This MD study examined the effects of framework topology on the diffusivity of H 2 .For example, H 2 diffuses more rapidly in IRMOF-8 than that in IRMOF-1 because of the relatively larger pore sizes of the former.The diffusivity of H 2 in IRMOF-18 is much slower than diffusion in IRMOF-1 and IRMOF-8 due to the steric hindrance effects of the pendant CH 3 groups in IRMOF-18.
IRMOFs can be further categorized as catenated and non-catenated structures.In catenation, two or more identical frameworks are intergrown at the expense of pore volume.Early studies showed that catenated MOFs can give better adsorption properties compared to their counterparts.(Ryan et al., 2008) The first study about the effects of catenation on the gas diffusion used equilibrium MD simulations in the canonical ensemble to investigate H 2 diffusion.(Liu et al., 2008a) Nosé-Hoover chain thermostat as formulated by Martyna et al.(Martyna et al., 1996) was used to calculate room temperature self diffusivities of H 2 in catenated and non-catenated MOFs.The results showed that H 2 self diffusivity in the IRMOFs without catenation such as IRMOF-10, IRMOF-12, IRMOF-14, IRMOF-16 (30-90×10 - Lee and coworkers also investigated the diffusion of H 2 in catenated MOFs, IRMOF-9, IRMOF-11 and IRMOF-13 at 77 K. (Lee et al., 2009) Diffusivities reported by Liu et al. were larger than the ones reported by Lee et al. by one order of magnitude since the former group performed the MD simulations at room temperature.The results of two studies were consistent; the diffusion rate of H 2 is dramatically reduced by the catenation of IRMOFs due to the interpenetrated chains of the catenated structures and/or by tighter binding of the H 2 molecules in catenated structures.Equilibrium MD simulation studies showed that the effect of catenation on CH 4 diffusivity is much larger than that on H 2 diffusivity at room temperature.(Xue et al., 2009) Xue and coworkers discussed that the motion of both CH 4 and H 2 is restricted by the catenated structures of IRMOF-11 and IRMOF-13 while the stronger interactions between CH 4 and atoms of the catenated frameworks lead to stronger confinement effects than that of H 2 in these IRMOFs.
Liu and coworkers (Liu et al., 2008b) investigated the influence of quantum effects on H 2 diffusivity using MD simulations.They used both the classical and the Feynman-Hibbs (FH) (Feynman&Hibbs, 1965) effective Buch potentials with UFF in their MD simulations to calculate self, corrected and transport diffusivities of H 2 in a MOF called Zn(bdc)(ted) 0.5 at 77 K.The inclusion of quantum effects increased the self diffusivity of H 2 at zero loading which was explained by the decrease in the diffusion energy barrier due to a non-uniform smearing of solid-fluid potential within the FH formalism.At higher loadings, inclusion of quantum effects decreased H 2 diffusivity which was attributed to the steric hindrance in narrow pores due to the increase in the effective size parameter for the solid-fluid and fluidfluid interactions.In contrast to self diffusivity, transport diffusivity is not strongly influenced by the quantum effects at 77 K.
In order to compare the diffusivities of gases in MOFs with those in zeolites, MD simulations were performed to calculate self, corrected and transport diffusivities of CH 4 and CO 2 in silicalite, IRMOF-1 and C 168 schwarzite.(Babarao&Jiang, 2008) The simulations were carried out in a canonical ensemble with a Nosé-Hoover thermostat and the equations of motion were integrated using a sixth order Gear predictor-corrector algorithm.(Allen&Tildesley, 1987) Both self and corrected diffusivities of CH 4 and CO 2 were found to be larger in IRMOF-1 (D self-CH4 :4-5×10 -8 m 2 /s, D self-CO2 :2-3×10 -8 m 2 /s) compared to the diffusivities in MFI and C 168 .This was attributed to the large pore volume of the IRMOF-1.This work also showed that in the limit of infinite dilution the diffusivities at various temperatures exhibit a good Arrhenius relationship.In another MD study, NVT ensemble with Berendsen(Frenkel&Smit, 2002) thermostat was used to examine the self diffusivity of CH 4 in alkoxy functionalized IRMOF-1.(Jhon et al., 2007) As expected, CH 4 diffusion was hindered due to the constriction of the pores as the length of the alkoxy chains increases.Comparison of the results with the early MD studies of Sarkisov et al. (Sarkisov et al., 2004) revealed good agreement whereas there is an unexplained small discrepancy between the results of Jhon et al. and Skoulidas et al.(Skoulidas&Sholl, 2005) at higher loadings.
Zeolite imidazolate frameworks (ZIFs) are a subclass of MOFs with their tetrahedral networks that resemble those of zeolites with transition metals linked by imidazolate ligands.(Banerjee et al., 2008;Banerjee et al., 2009;Hayashi et al., 2007;Park et al., 2006) www.intechopen.com Recent Advances in Molecular Dynamics Simulations of Gas Diffusion in Metal Organic Frameworks 265 Zeolites are known with the Al(Si)O 2 unit formula, whereas ZIFs are recognized by M(Im) 2 where M is the transition metal (zinc, cobalt, copper, etc.) and Im is the imidazolate-type linker.Recent MD simulations focused on gas diffusion in ZIFs.For example, self and corrected diffusivities of CO 2 , CH 4 and H 2 were simulated using equilibrium MD in ZIF-68 and ZIF-70.(Rankin et al., 2009) That study underlined the importance of including chargequadrupole interactions on the diffusivity of CO 2 .Simulation results clearly revealed that addition of charge-quadrupole interaction terms results in almost one order of magnitude drop in the self and transport diffusivities of CO 2 at low loadings.At high loadings diffusivities calculated from MD simulations with or without charge-quadrupole interaction terms converge towards the same values.The diffusivities of CO 2 , CH 4 and H 2 in ZIF-68 were found to be lower than the ones in ZIF-70 since ZIF-68 has narrower pores hence provides stronger confinement of the adsorbate molecules in the pores.Self diffusivity of CO 2 in ZIF-68 and ZIF-69 was also computed by MD simulations.(Liu et al., 2009) The diffusion of CO 2 in ZIF-68 and ZIF-69 was found to be nearly an order of magnitude slower than that in IRMOF-10 and IRMOF-14.This was attributed to the smaller pores of ZIFs and their structural characteristic that causes larger steric hindrance.Pantatosaki and coworkers computed H 2 self diffusion in ZIF-8 using both LJ and FH potentials at 77 and 300 K. (Pantatosaki et al., 2010) The diffusivity predictions showed that quantum mechanical description of H 2 at ambient temperatures is unimportant whereas MD simulations showed a marked difference between the values obtained from the classical and quantum mechanical description at 77 K.A recent MD study computed self diffusivities of H 2 , CO 2 , CH 4 and N 2 in ZIF-2, ZIF-4, ZIF-5, ZIF-8 and ZIF-9.(Battisti et al., 2011) Results showed that gases except H 2 do not diffuse appreciably in ZIF-5 at least within the time interval of the MD calculations which makes ZIF-5 promising in H 2 separations as a molecular sieve.
Self diffusivities of H 2 , CH 4 and CO 2 in bioMOF-11 were computed from canonical ensemble MD simulations at 298 K. (Atci et al., 2011) BioMOFs are another subclass of MOFs that have been recently discovered.They incorporate simple biomolecules and biocompatible metal cations in their structures as linkers and metals.(An et al., 2009a;An et al., 2009b) Gas diffusion in bioMOFs was found to be similar to IRMOFs in terms of magnitude and loading dependence.As can be seen from the literature reviewed so far, most of the MD studies on MOFs computed self diffusivity of gases rather than corrected diffusivities since the calculation of the latter is computationally demanding.Keskin computed both single component self and corrected diffusivities of CH 4 and H 2 as a function of fugacity and pore loading in CPO-27-Ni.(Keskin, 2010a) The diffusivity of H 2 (4×10 -3 cm 2 /s) was faster than CH 4 (6×10 -4 cm 2 /s) as expected.Single component corrected diffusivities were found to be higher than the self diffusivities, since corrected diffusivity by definition includes information on the collective motion of multiple adsorbed molecules that is relevant to net mass transport.
Figure 3 represents the self diffusivity of CO 2 computed from MD simulations in the widely studied MOFs at room temperature.Gas diffusion in MOFs having large pores (IRMOF-1, CuBTC, Zn(bdc)(ted) 0.5 ) is higher than the one in MOFs having narrow pores (Cu(hfipbb)(H 2 hfipbb) 0.5 , MMIF).The CO 2 self diffusivity decreases with increased adsorbed loading in bioMOF-11, Cu(hfipbb)(H 2 hfipbb) 0.5 and MMIF since CO 2 reaches saturation in these MOFs due to their small pore volumes.The diffusivities in large pore MOFs do not change significantly with increased loadings since CO 2 is further away from the saturation loading in MOFs having large pore volumes.
As the literature cited so far indicated most MD studies focused on diffusion of small gas molecules such as Ar, CH 4 , CO 2 , N 2 , H 2 in MOFs.Limited number of studies investigated diffusion of larger gases.The self diffusivities of hexahydro-1,3,5-trinitro-1,3,5-triazine (RDX) were generated by MD in IRMOF-1, IRMOF-3, IRMOF-10.(Xiong et al., 2010) The trend for the self diffusivities of RDX in MOFs followed the pore sizes, highest in IRMOF-10 and lowest in IRMOF-3.The self diffusivity of ethane, n-butane, n-hexane and cyclohexane in a MOF with the organic linker tetrakis[4-(carboxyphenyl)oxamethyl]methane was studied using MD simulations.(Sun et al., 2011) For linear alkanes, the diffusivities decreased dramatically with increased chain length.The specific MOF studied in this work exhibited high selectivity towards n-hexane as a result of kinetics.(Atci et al., 2011;Erucar&Keskin, 2011;Keskin, 2011b;Keskin&Sholl, 2009b;Watanabe et al., 2009) In some cases, the adsorbate molecules cannot move in the MOF pores at a rate that can be measured by MD simulations.This case is generally observed when the kinetic diameter of the gas molecule is very similar in size to the pore diameter of the MOF.For example, initial MD simulations of adsorbed CH 4 in a rigid Cu(hfipbb)(H 2 hfipbb) 0.5 indicated that CH 4 can not move between adjacent cages on the nanosecond time scales accessible using MD due to the large energy barrier.(Watanabe et al., 2009)

267
where υ is the pre-exponential factor (10 12 -10 13 s -1 ), E trans is the transition energy barrier computed using DFT calculations for a flexible MOF, k B is the Boltzmann constant, T is temperature.These calculations resulted in a one-dimensional diffusivity, where a is the cage-to-cage distance along the pore.Keskin studied diffusion of CH 4 and CO 2 in a microporous metal-imidazolate framework and similarly observed that CH 4 diffusion in this MOF is not accessible by MD due to a very large energy barrier (95 kJ/mol) that exists for the adsorbate to move through the pore.(Keskin, 2011b) It is important to note that both Watanabe et al. and Keskin concluded that Cu(hfipbb)(H 2 hfipbb) 0.5 and MMIF are very promising materials for separation of CO 2 from CH 4 since CO 2 diffuses several orders of magnitude faster than CH 4 in these MOFs.

Mixture diffusion
In most practical applications, gases exist as mixtures rather than single components.For example, in membrane-based separations, at least two gas components exist.The relative transport rates of these components inside the material of interest are crucial in determining the overall performance of a material.Therefore, understanding mixture diffusion in MOFs is essential to design these materials as separation devices.In this section, MD simulations which predicted multi-component mixture diffusion in MOFs will be reviewed.These studies have mostly focused on self diffusivities of CO 2 /CH 4 , CO 2 /H 2 , CO 2 /N 2 , CO 2 /CO and CH 4 /H 2 mixtures.The diffusion selectivities of MOFs for these gas mixtures have been computed using MD to understand the potential of MOFs in kinetic-based separations.It is important to highlight the fact that number of mixture MD simulations is limited compared to the number of single component MD simulations since characterizing diffusivity of gas mixtures is harder than studying a single species.
Keskin and coworkers provided the first gas mixture diffusivity data in a MOF material using MD simulations.(Keskin et al., 2008) They computed self and Fickian diffusivities of CH 4 /H 2 mixtures at various compositions in CuBTC.Theoretical correlations that can estimate mixture self and Fickian diffusivities based on single component data were tested in that work and the predictions of the correlations were compared with the results of direct MD simulations.This will be discussed in detail in Section 5. Keskin and Sholl later showed that if MD (GCMC) simulations are used to compute the mixture diffusivities (adsorbed amounts of each species) in a MOF, one can simply estimate the selectivity of this MOF as a membrane using:(Keskin&Sholl, 2009b) where α perm,1/2 , α diff,1/2 and α sorp,1/2 represent permeation selectivity, diffusion selectivity and sorption selectivity of species 1 over species 2, respectively.In this approximate expression, the diffusion selectivity is defined as the ratio of self diffusivities in a binary mixture (D i,self (q i )) and the sorption selectivity is described as the ratio of adsorbed molar loadings, q i .This expression predicts a membrane's selectivity at a specified feed pressure and composition based on a single mixture GCMC simulation and an MD simulation performed at the loadings determined from this GCMC calculation.Keskin and Sholl computed self diffusivities for CH 4 /H 2 , CO 2 /CH 4 and CO 2 /H 2 mixtures in several MOFs, IRMOF-1, IRMOF-8, IRMOF-9, IRMOF-10, IRMOF-14, COF-102, Zn(bdc)(ted) 0.5 using MD simulations and based on these diffusivities they estimated the membrane selectivity of these MOFs.
Mixture self diffusivities of CH 4 /H 2 in CPO-27-Ni and CPO-27-Co were computed at 298 K for a wide pressure range and selectivity of these MOFs in CH 4 /H 2 separations were predicted.(Keskin, 2010a) MD calculations were carried out to evaluate diffusion selectivities and permeation selectivities of ZIF-3 and ZIF-10 for CH 4 /H 2 , CO 2 /CH 4 and CO 2 /H 2 mixtures.(Keskin, 2011a) Figure 4 shows the diffusion selectivities for these mixtures in the pores of ZIF-3 and ZIF-10 as a function of pressure at room temperature.Using the same approach, Krishna and   The self diffusivities of adsorbed CH 4 /H 2 mixtures were examined at different compositions in Zn(bdc)(ted) 0.5 .(Keskin, 2010b) The self diffusivities of CH 4 (H 2 ) in the CH 4 /H 2 mixture were larger (smaller) than pure component CH 4 (H 2 ) self diffusivity at the same loading.This observation is natural since the fast diffusing H 2 molecules in the mixture speeds up the slowly diffusing CH 4 molecules.Self diffusivities of CO 2 /CH 4 and CH 4 /H 2 mixtures were computed in ZIF-68 and ZIF-70 using NVT-MD simulations.(Liu et al., 2011) Results indicated that diffusion of CH 4 is increased with increasing concentration of H 2 in the CH 4 /H 2 mixture, while the diffusivity of H 2 decreases with increasing CH 4 concentration.In contrast, the diffusivity of CH 4 was essentially independent of the concentration of CO 2 in the CO 2 /CH 4 mixture, while CO 2 diffusivity decreases with increased CH 4 loading, even though the diffusivity of CH 4 is substantially larger than that of CO 2 .This unusual behavior was explained in terms of differences in adsorption site preferences due to chargequadrupole interactions.Another recent MD study examined the self diffusivities of equimolar CO 2 /ethane, CH 4 /ethane and CO 2 /methanol mixtures in Zn(tbip) including the flexibility effects.(Seehamart et al., 2011) Similar to previous observations, faster diffusing molecules accelerate the slower diffusing molecules whereas the slower ones slow down the faster ones through the channel of Zn(tbip).

Comparison of diffusivities from molecular dynamics with experiments
Measuring diffusivities of gas molecules in nanoporous materials is a challenging process, therefore experimentally measured diffusion data for gases in the pores of MOFs is still very limited.Stallmach and co-workers (Stallmach et al., 2006) carried out the first experimental study in the literature for diffusivity of hydrocarbons in MOF-5.They measured diffusion of methane, ethane, n-hexane, benzene by pulsed field gradient-nuclear magnetic resonance (PFG-NMR) which is a well-established technique for intra-crystalline diffusion studies in nanoporous materials.Diffusion of methane and ethane in MOF-5 was found to be faster than in NaX which was attributed to the larger pores of the former.This study supplied the first experimental data points for gas diffusion in MOFs for direct comparison between experiments and MD simulations.The measured diffusivity of n-hexane, 3.2-4.1×10 - m 2 /s, was found to be in a good agreement with the value of 2.2×10 -9 m 2 /s predicted by earlier MD simulations (Sarkisov et al., 2004) for a slightly higher loading.However, the self diffusivity of CH 4 measured by PFG-NMR was about one order of magnitude higher than the value of 3.1×10 -8 m 2 /s reported in MD simulations.(Sarkisov et al., 2004;Skoulidas&Sholl, 2005) Stallmach and coworkers attributed this discrepancy to the imperfections that may exist in the MOF structure and loadings used in MD simulations which were lower than the ones considered in the experiments.Zhao and coworkers (Zhao et al., 2009) measured diffusivity of CO 2 in MOF-5 and reported a value (8×10 -13 m 2 /s) which is several orders of magnitude smaller than the one obtained by the MD simulations (4×10 -9 m 2 /s) (Skoulidas&Sholl, 2005) and also significantly smaller than the diffusivity of larger adsorbates such as n-hexane, benzene measured by other groups.This large difference between experiments and simulations can be again attributed to the imperfections in the synthesized MOF structure.
The first experimental exploration of the H 2 self diffusivity in MOFs was performed by quasielastic neutron scattering (QENS) measurements.(Salles et al., 2008) The QENS technique has proved to be very powerful to extract the loading dependence of the diffusivities for a wide range of adsorbates including H 2 diffusivity in zeolites.(Jobic et al., 1999) Combining QENS technique with molecular simulations has been successful in the past to characterize the diffusion mechanism of various adsorbates in nanoporous materials.(Jobic&Theodorou, 2007) The self diffusivities of H 2 in MOFs, MIL-47(V) and MIL-53(Cr) were extracted from QENS measurements and compared with the ones predicted by MD simulations performed in the NVT ensemble using the Evans isokinetic www.intechopen.comthermostat.(Frenkel&Smit, 2002) Simulated data was in a good agreement with the experimentally measured data for both MILs.Experiments measured a diffusivity of 9×10 -8 m 2 /s (1.65×10 -7 m 2 /s) and simulations predicted 4.5×10 -8 m 2 /s (1.5×10 -7 m 2 /s) at a loading of 0.5 H 2 molecules per unit cell of MIL-53(Cr) (MIL-47(V)).In a similar study, QENS measurements were combined with MD simulations in NVT ensemble using either Berendsen or Evans thermostat to determine the self diffusivity of H 2 in the same MILs.(Salles et al., 2008) Two different force fields, spherical one site model (Frost et al., 2006) and explicit two atoms model (Yang&Zhong, 2005) were used in MD simulations of H 2 .Comparisons between QENS data and MD simulations clearly showed that the two force fields lead to very similar diffusivity values that produce the experimental value.This observation suggests that H 2 diffusion is not significantly affected by the potential model.
A combination of MD and QENS measurements were used to examine the diffusivity of water in MIL-53(Cr).(Salles et al., 2011) The breathing of this MOF upon water adsorption induces a structural transition between narrow pore (NP) and large pore (LP) forms.The self diffusivity of water was faster in LP form (8×10 -10 m 2 /s) compared to the one in NP form (2.5×10 -11 m 2 /s) since the confinement degree was much higher in NP structure.As an extension of this work, self, corrected and transport diffusivities of CO 2 in MIL-47(V) were determined using MD and QENS.(Salles et al., 2010) While self and corrected diffusivities exhibited a decreasing profile with increased loading as expected, transport diffusivity presented an unexpected trend with a decrease at low loadings.This behavior was attributed to the unusual evolution of thermodynamic correction factor.This work was a good example of probing the transport diffusivity of gases in MOFs by combining MD and QENS.
Two experiments studied diffusion of alkanes in MOFs: The diffusivity of n-butane, isobutane, 2-methylbutane and 2,2-dimethylpropane in CuBTC was investigated using infrared microscopy and MD simulations.(Chmelik et al., 2009) In another work, intracrystalline self diffusivities of propane, propene, n-butane, 1-butene, n-pentane and n-hexane in CuBTC were assessed using PFG-NMR and MD simulations.(Wehring et al., 2010) For the nalkanes, measured diffusivities within the experimental uncertainty agreed with the values from the MD simulations.The different trends observed in diffusivities of alkanes remained as an unsolved issue.

Comparison of diffusivities from molecular dynamics with theories
The prediction of transport properties in chemical mixtures from data taken from single component studies has been a long standing goal in describing mass transport in nanoporous materials.The validation of methods for this task can have great practical significance, but this type of validation can only be considered when high quality mixture diffusion data is available from MD.This section will present the validity of theoretical correlations by comparing their predictions for self diffusivities and Fickian diffusivities with the ones derived from MD.In the literature, Krishna-Paschek's (KP) correlation (Krishna&Paschek, 2002) and Skoulidas, Sholl and Krishna (SSK) correlation (Skoulidas et al., 2003) have been widely used to predict the self diffusivities and Fickian diffusivities of a binary gas mixture, respectively.
Keskin and Sholl studied MOF-5 as a membrane for separation of CO 2 /CH 4 mixtures.(Keskin&Sholl, 2007) As discussed in previous sections, in order to study transport of gas mixtures in membranes mixture diffusivity data is required.However, at the time of that study there was no binary diffusion data available for MOF-5.Keskin and Sholl applied the SSK approach to quantify mixture diffusion of CO 2 /CH 4 in MOF-5.This approach combines information from the loading dependence of the single component self diffusivities and corrected diffusivities (computed from MD simulations) with the binary adsorption isotherms (computed from GCMC simulations) to predict the loading and composition dependent matrix of binary diffusion coefficients.The SSK approach defines the mixture diffusivities for all loadings and compositions, an important feature of any description that will be used in examining a wide range of potential membrane operating conditions.Prior tests of this method by comparison with detailed atomic simulations of binary diffusion in silica zeolites and carbon nanotubes indicated that this approach is accurate for a wide variety of adsorbed mixtures.(Sholl, 2006) A year later, Keskin and coworkers presented the validity of SSK approach in a MOF.(Keskin et al., 2008) They examined both KP and SSK approaches by comparing predictions of these methods with the results of MD simulations for mixture transport of H 2 /CH 4 in CuBTC.In order to use SSK correlation, continuous functions describing the pure component self and corrected diffusivities were required.The self and corrected diffusivities of each species in H 2 /CH 4 mixture were calculated by MD simulations.Based on these single component diffusivities, the SSK approach predicted the Fickian diffusivities.Mixture MD simulations in a Nosé-Hoover thermostat in the NVT ensemble calculated Onsager coefficients (Equation 4) for H 2 /CH 4 mixture and these values were converted to Fickian diffusivities (Equations 5 and 6).The predictions of the SSK approach for the Fickian diffusivities were in good agreement with the direct MD simulations of binary diffusion, suggesting that this approach may be a powerful one for examining multicomponent diffusion in MOFs.Mixture self diffusivities were predicted using KP correlations based on single component self diffusivities, corrected diffusivities and fractional loadings.Comparison between KP predictions and mixture MD simulations were also found to be in a good agreement.The SSK approach was also used to obtain

Conclusion and outlook
Because of the large number of different MOFs that exist, efforts to predict the performance of MOFs using molecular modeling play an important role in selecting materials for specific applications.The high number of publications on MOFs and the dense interest of academy and industry on these new nanoporous materials hint that MOFs have numerous potential applications.Since almost all of these applications require the knowledge of molecular transport rates, MD studies become one of the most beneficial methods in studying MOFs.
As is evident from the volume of literature cited, this area is growing rapidly.The development of quantitative information about mixture diffusion in MOFs is just beginning (section 3.2) whereas a significant number of studies have already considered single component gas diffusion in MOFs (section 3.1).Detailed understanding of mixture diffusion in MOFs will be very beneficial for design of MOF membranes, adsorbents, catalysts and sensors.Current opportunities and challenges of using MD simulations in assessing transport rates of gases in MOFs will be addressed below.

Opportunities
The most significant opportunity of employing MD simulations for obtaining gas diffusivity in MOFs lies in areas where experiments for transport property of interest (transport diffusivities, energy barrier to diffusion) are challenging, not in reiterating properties that have already been addressed experimentally.Measuring diffusivity at a wide range of loadings in the pores at extreme conditions such as infinite dilute loading and/or saturation loading is experimentally difficult.MD simulations can provide information about gas diffusion in MOFs' pores under these conditions.Getting diffusivity data as a function of gas loading is crucial to design membranes, adsorbents, catalysts from MOFs that will work under a wide range of operating conditions.
As addressed in Section 5, the development of quantitative information about mixture diffusion in MOFs is just beginning.Since performing mixture MD simulations for MOFs with large frameworks and for gas mixtures at high adsorbed loadings are computationally demanding, theoretical correlations that predict mixture diffusion based on single component diffusion data are very useful.Recent research search showed that these models yield accurate results for at least simple chemical mixtures in MOFs.Testing and validation of theoretical correlations for predicting gas diffusivity in various subclasses of MOFs will be useful to widely utilize these correlations for different structures.
A great advantage of using MD simulations is to test hypothetical MOF structures for particular applications if the metric describing the performance of a material for the application can be directly calculated.For example, Düren and coworkers used GCMC simulations to design materials with large adsorption capacities for CH 4 .(Düren et al., 2004) In a similar way, MD simulations can be used to design materials with slow diffusivities for CH 4 and fast diffusivities for CO 2 to identify materials that will be promising in kinetic separation of CO 2 from CO 2 /CH 4 mixtures.

Challenges
The development of accurate classical interatomic potentials for describing gas diffusion in MOFs remains challenging.From the modeling perspective, it is important to use experimental diffusion data from a broad range of conditions to parameterize interatomic potentials whenever this is practical.However, as discussed in Section 4, the number of experimental data on gas diffusion in MOFs is very limited.Furthermore, developing potentials specific to a MOF structure is not the solution since hundreds of different MOF structures are available.Therefore, efforts to test and improve the transferability of potentials among related families of MOFs will have a great value.One of the major challenges in using MD simulations for MOFs was addressed in Section 2.2: absence of fully flexible force fields.Rigid framework assumption creates tremendous savings in computational effort.A handful of studies used flexible force fields to include the lattice dynamics effects on gas diffusivity in MOFs.These studies showed that there can be orders of magnitude difference between the diffusivity data from MD simulations using rigid framework and the one using flexible framework, specifically for large adsorbates.This issue indeed turns to be related with the challenge listed above, having accurate flexible interatomic potentials which can be applied to a family of MOF structures in a computationally meaningful time scale.
Another major challenge, especially in diffusivity simulations of CO 2 and N 2 , is the choice of method to assign partial charges to MOF atoms.The QM calculations were used to define partial point charges in the literature, however there is no unique way to accomplish this task and different charge decomposition methods can give rather different results.Studies have shown that charge effects are important especially for computing diffusivities at low loadings.Careful studies that establish reliable approaches in charge assignment will be very useful in employing MD simulations for diffusion of polar and quadrupolar molecules in MOFs.
To date MD simulations have been used to compute the transport rates of adsorbates in MOFs.One remaining challenge is to predict the long term stability of MOFs since this is a serious issue in practical applications of these materials.Although stability issue sounds to be most likely addressed by experimental studies, one recent MD study which investigated the mechanism of water induced decomposition of IRMOF-1(Greathouse&Allendorf, 2006) showed that molecular simulations can be also helpful in this area.

Fig. 1 .
Fig. 1.Unit cell representation of a widely studied MOF, CuBTC.From left to right: Empty CuBTC, CuBTC with CH 4 molecules (blue spheres) in the pores, CuBTC with CH 4 and H 2 (orange spheres) molecules in the pores.The atoms in the unit cell are copper (green), oxygen (red), carbon (gray) and hydrogen (white).
Fig. 2. Open bars represent the number of publications featuring the term 'metal organic framework', closed bars represent the number of publications featuring the terms 'metal organic framework' and 'molecular dynamics' and 'diffusion'.(Source: ISI Web of Science, retrieved August, 8 2011).

Table 1 .
As an example, the most widely used potential parameters (ε: energy parameter, k B : Boltzmann constant, σ: size parameter) are listed in Table1.Potential parameters used for adsorbate molecules in MD simulations.
(Atci et al., 2011)diffusion selectivities and permeation selectivities for equimolar CO 2 /CH 4 and CO 2 /H 2 mixtures as a function of total pore loading in CPO-27-Zn, CPO-27-Mg, IRMOF-1.(Krishna&vanBaten,2011) Atci and coworkers evaluated the mixture self diffusivities of CH 4 /H 2 , CO 2 /CH 4 and CO 2 /H 2 in bioMOF-11 at the adsorbed loadings calculated from GCMC simulations.(Atciet al., 2011) (Liu et al., 2011)calculated self diffusivities of CH 4 and CO 2 in IRMOF-1 as a function of total loading based on the adsorption of an equimolar mixture using MD simulations and compared their results with the predictions of KP correlation.(Babarao&Jiang,2008)Theorypredictionswere found to be in a fairly good agreement with MD simulations particularly for CH 4 diffusivity in IRMOF-1 whereas the CO 2 diffusivity was slightly overestimated by the theory.No certain reasoning was given for this overestimation.The predictions of KP correlations for mixture self diffusivities of CH 4 and H 2 were in reasonable agreement with the results of MD simulations for ZIF-68 and ZIF-70.(Liuet al., 2011)