**2. Molecular dynamics models**

256 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Fig. 1. Unit cell representation of a widely studied MOF, CuBTC. From left to right: Empty CuBTC, CuBTC with CH4 molecules (blue spheres) in the pores, CuBTC with CH4 and H2 (orange spheres) molecules in the pores. The atoms in the unit cell are copper (green),

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.

2002 2003 2004 2005 2006 2007 2008 2009 2010 2011

Years

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,

oxygen (red), carbon (gray) and hydrogen (white).

1

10

100

Number of publications

retrieved August, 8 2011).

1000

Gas diffusion is an observable consequence of the motion of atoms and molecules as a response to external force such as temperature, pressure or concentration change. Molecular dynamics (MD) is a natural method to simulate the motion and dynamics of atoms and molecules. The main concept in an MD simulation is to generate successive configurations of a system by integrating Newton's law of motion.(Frenkel&Smit, 2002) Using MD simulations, various diffusion coefficients can be measured from the trajectories showing how the positions and velocities of the particles vary with time in the system. Several different types of gas diffusion coefficients and the methods to measure them will be addressed in the next section in details.

In accessing the gas diffusion in nanoporous materials, equilibrium MD simulations which model the behavior of the system in equilibrium have been very widely utilized. In equilibrium MD simulations, first a short grand canonical Monte Carlo (GCMC) simulation is applied to generate the initial configurations of the atoms in the nanopores. Initial velocities are generally randomly assigned to each particle (atom) based on Maxwell-Boltzmann velocity distribution.(Allen&Tildesley, 1987; Frenkel&Smit, 2002) An initial NVT-MD (NVT: constant number of molecules, constant volume, constant temperature) simulation is performed to equilibrate the system. After the equilibration, Newton's equation is integrated and the positions of each particle in the system are recorded at a pre-specified rate. Nosé-Hoover thermostat is very widely applied to keep the desired temperature and the integration of the system dynamics is based on the explicit N-V-T chain integrator by Martyna et al.(Martyna et al., 1992; Martyna et al., 1996) By keeping temperature constant, Newton's equations are integrated in a canonical ensemble (NVT) instead of a microcanonical ensemble (NVE: constant number of molecules, constant volume and constant energy). To describe the dynamics of rigid-linear molecules such as carbon dioxide the MD algorithm of Ciccotti et al.(Ciccotti et al., 1982) is widely used. The so-called order *N* algorithm(Frenkel&Smit, 2002) is implemented to calculate the diffusivities from the saved trajectories.

In order to perform classical MD simulations to measure gas diffusion in MOFs' pores, force fields defining interactions between gas molecules-gas molecules and gas molecules-MOF's atoms are required. Once these force fields are specified, dynamical properties of the gases in the simulated material can be probed. These force fields will be studied in two parts: models for gas molecules (adsorbates) and models for MOFs (adsorbents).

#### **2.1 Models for gases**

Diffusion of hydrogen, methane, argon, carbon dioxide and nitrogen are very widely studied in MOFs. For H2, three different types of fluid-fluid potential models have been

Recent Advances in Molecular Dynamics

molecules present in the samples.(Keskin et al., 2009b)

than the value obtained from MD simulations with flexible lattice.

IRMOF-1 due to the nature of the organic linkers.(Xue&Zhong, 2009)

Simulations of Gas Diffusion in Metal Organic Frameworks 259

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

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

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

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

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.


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