# Intermolecular Force Field Parameters Optimization for Computer Simulations of CH4 in ZIF-8

Authors: Phannika Kanthima; Pikul Puphasuk; Tawun Remsungnen
Journal: Journal of Nanotechnology
 (2016)

---

## Abstract

The differential evolution (DE) algorithm is applied for obtaining the optimized intermolecular interaction parameters between CH4 and 2-methylimidazolate ([C4N2H5]−) using quantum binding energies of CH4-[C4N2H5]− complexes. The initial parameters and their upper/lower bounds are obtained from the general AMBER force field. The DE optimized and the AMBER parameters are then used in the molecular dynamics (MD) simulations of CH4 molecules in the frameworks of ZIF-8. The results show that the DE parameters are better for representing the quantum interaction energies than the AMBER parameters. The dynamical and structural behaviors obtained from MD simulations with both sets of parameters are also of notable differences.


---

## Body

## 1. Introduction

Methane (CH4) is one of the greenhouse gases caused by farming, livestock, burning biofuel, and burning fuels like coal, oil, and natural gas. In order to reduce this gas in the atmosphere, efficient and cost saving methods are studied worldwide. The application of porous material for high capacity and cost-effective gas separation and storage is one intensively investigated way [1–7]. Zeolitic imidazolate frameworks (ZIFs), a subclass of metal organic frameworks (MOFs), which are composed of positive transition metal ions such as Zn, Co, and Cu are linked by ditopic imidazolate ligands. Because of their high porosity, thermal and chemical stability, and the exact structure, they are used in the greenhouse gases separation and storage [1–5, 8–11]. ZIF-8 [8–10, 12] is one of the ZIFs, which is composed of Zn2+ ion linked with ditopic 2-methylimidazolate ([C4N2H5]−) ligands with a chemical formula of Zn[C4N2H5]2. It is applied for the separation of CH4 from natural gas and the reduction of CH4 in the atmosphere. Structures of CH4 molecule, ZIF-8 framework, and [C4N2H5]− molecule are shown in Figure 1.Figure 1
Structures of CH4 molecule, ZIF-8 framework, and [C4N2H5]− molecule.The knowledge of chemical processes such as diffusion, adsorption, and reaction at the molecular level can help to design and improve the properties of these materials for efficient use in gas separation and storage [1–5, 8, 10, 11]. The key to understanding the behavior of the adsorption and diffusion of CH4 in ZIF-8 is the intermolecular interaction. Since the quantum method needs large resources and time consumption, it is not suitable for computer simulation which is composed of several atoms. Therefore, the development of intermolecular interactions functions is important. These functions are generated by the experimental data or the computational quantum data. Moreover, the optimization method is used to find the suitable forms and parameters of such functions. The famous form of the intermolecular interactions is the Lennard-Jones function [5, 8]. In this paper, the parameters of Lennard-Jones are optimized based onab initio using the heuristic method called the differential evolution (DE) in order to obtain the specific parameters for use in the molecular dynamic simulations of CH4 in the ZIF-8.


## 2. Models and Calculations

The geometrical structure of [C4N2H5]− is cut directly from the framework of single-crystal XRD data [9] and the geometric model of CH4 molecule is taken from the interaction site model with the bond length of 1.089 and bond angle of 109.5 [13]. A total of 1377 single point energies of CH4-[C4N2H5]− complexes are obtained at B2PLYPD/6-31G(d) level using Gaussian 09 package [14]. Some configurations of CH4-[C4N2H5]− complexes are illustrated in Figure 2.Figure 2
Some configurations of CH4-[C4N2H5]− complexes. In each complex, a methane molecule is started from the initial point and moved out along an arrow line.The interaction energy with the corrected basis sets superposition errorΔ
E
C
C is represented by the following equation [15]:(1)
Δ
E
C
C
=
E
A
B
A
B
-
E
A
B
A
-
E
A
B
B
,where E
A
B
A
B, E
A
B
A, and E
A
B
B are the total energy of CH4-[C4N2H5]− complex, the energy of CH4 molecule, and the energy of [C4N2H5]− molecule with the basis sets of both molecules, respectively [10]. The Lennard-Jones function is used to estimate the intermolecular interaction between CH4 with zero charges and [C4N2H5]− molecules:(2)
U
r
=
∑
i
=
1
m
∑
j
=
1
n
A
i
j
r
i
j
12
-
B
i
j
r
i
j
6
,where A
i
j
=
4
ε
i
j
σ
i
j
12, B
i
j
=
4
ε
i
j
σ
i
j
6, σ
i
j
=
σ
i
+
σ
j
/
2, and ε
i
j
=
ε
i
ε
j. σ and ε used for computing the interactions are shown in Table 1.Table 1
The AMBER force field for atoms of ZIF-8 framework and CH4 molecules [8].
Molecules
Atoms
σ (Å)
ε (kcal/mol)
ZIF-8
Zn
1.960
0.013
CR
3.400
0.086
N
3.250
0.170
CC
3.400
0.086
H4
2.421
0.015
CT
3.400
0.110
HT
2.650
0.016
CH4
H
2.650
0.016
C
3.400
0.110Although the general AMBER force field is famous for use in molecular dynamics, it is not suitable for some specific models. Thus, we use differential evolution algorithm [16] with the interaction energies from quantum mechanics method with zero charges to optimize the specific intermolecular interaction parameters A and B. The objective function is the following chi-square equation:(3)
χ
2
=
∑
k
=
1
N
w
k
E
C
C
k
-
U
k
2
d
f
,where k is the index of configuration, E
C
C
(
k
) is the counterpoise corrected interaction energy, U
(
k
) is the Lennard-Jones energy, w
(
k
)
=
1
/
E
C
C
(
k
)
-
w
E is the weighted value with respect to w
E which is set very close to the minimum energy of E
C
C, d
f
=
N
-
D is the degree of freedom, N is the total number of energies, and D is the total number of parameters.The optimized molecular interaction parameters and the original AMBER [8] parameters are used in computer simulations for studying the adsorption and diffusion of CH4 in porous material ZIF-8 by using DL_POLY program (version 2.20) [17, 18]. The flexible ZIF-8 framework is composed of eight (2
×
2
×
2) unit cells and loadings of CH4 molecules per unit cell are 1, 2, 4, 6, and 8, respectively. The NVT ensemble is set at 300 K and time step of 1.0 fs. The simulations are equilibrated for 1.0 ns (1,000,000 steps), and further 1.0 ns trajectories data are collected at every 200 time steps for studying the structural and dynamical properties.The site-site radial distribution functions (RDFs),g
α
β
(
r
), are a statistical analysis which is the probability of finding an atom β around a reference atom α:(4)
g
α
β
r
=
n
r
V
4
π
r
2
N
β
d
r
,where n
(
r
) is the number of β atoms in thickness of radius of circle d
r at distance r, V is the volume of box, and N
β is the number of all β atoms.A smaller particle will diffuse without direction. The self-diffusion coefficient,D
s (m2/s), is computed by mean squared displacement (MSD) method using the Einstein relation for a three-dimensional system [5, 19]:(5)
D
s
=
l
i
m
t
→
∞
1
6
t
1
N
∑
i
=
1
N
r
i
t
-
r
i
t
0
2
,where N is the total number of guest particles and r
i
(
t
) and r
i
(
t
0
) are the position vectors of the diffusing molecule in framework at time t and time origin t
0.


## 3. Results and Discussion

The intermolecular interaction parametersA
i
j and B
i
j between atom of CH4 and atom of [C4N2H5]− obtained from AMBER force field and DE method are presented in Table 2. It is clear that the optimized parameters provide the value of chi-square of 0.159 which is much smaller than 35.052 obtained with AMBER parameters. Thus, the energies with DE parameters are close to the QM energies more than those obtained with the AMBER parameters (see Figure 3).Table 2
ParametersA
i
j and B
i
j (kcal/mol) obtained with AMBER force field and DE method.
Atom
Force field
AMBER
DE
i
j
A
i
j
B
i
j
A
i
j
B
i
j
H
CR
86889
113
19016
104
N
90349
137
17399
257
CC
86889
113
19016
104
H4
4245
15
431
24
CT
98631
128
13590
26
HT
7675
22
768
40
C
CR
925931
599
1677902
762
N
1000645
740
1457350
1262
CC
925931
599
1677902
762
H4
60722
99
32285
37
CT
1050025
679
1439945
343
HT
98631
128
29852
36Figure 3
Comparison betweenΔ
E
C
C (red line), Δ
E
A
M
B
E
R (dash-purple), and Δ
E
D
E (dot-black) energies of some CH4-[C4N2H5]− complexes.The adsorption and diffusion of CH4 molecules in ZIF-8 framework are investigated by molecular dynamic simulations using the DE and AMBER parameters for comparison. The site-site RDFs between all atoms in framework to atoms H and C of CH4 molecules are presented in Figure 4. The results of RDFs obtained by all simulations indicate that atoms H4 and HT are favored sites for adsorption. The parameters of DE yield the average minimum distances from atoms H and C of CH4 molecules to atoms of the framework as 1.45 Å and 2.15 Å, which are shorter than those of 1.77 Å and 2.33 Å obtained with the original AMBER parameters. These obtained distances agree well with the interaction profiles (see Figure 3) where both minimum binding energies and minimum binding distances of QM and DE profiles are stronger and shorter than the AMBER profiles. The simulations with AMBER parameters show that CH4 molecules are distributed within the cavities and less diffuse between the cavities, while the CH4 molecules are distributed within the cavities more close to atoms H4 and HT and more diffuse between the cavities with the optimized DE parameters in the simulations (see Figure 5).Figure 4
RDFs of simulations between ZIF-8-H and ZIF-8-C using (a) AMBER parameters and (b) DE parameters.
(a)
(b)Figure 5
Diffusion of CH4 molecules in ZIF-8 using (a) AMBER parameters and (b) DE parameters.
(a)
(b)The self-diffusion coefficients for loadings 1, 2, 4, 6, and 8 of CH4 molecules per unit cell obtained by simulations are shown in Table 3. The D
s values obtained from AMBER parameters are in agreement with previous works using the same parameters set. For all loadings, the D
s values obtained by the DE parameters are approximately one magnitude higher than those obtained by AMBER parameters. The higher distributions and diffusion obtained from DE parameters can be explained by interaction profiles and the obtained RDFs. Since the CH4 molecules are more close to the imidazolate ligands which are flexible, the movements of CH4 molecules are additionally induced by the ligands, and locating around the cage windows increases the chance of cage-to-cage migration. By the way, these obtained D
s values agree well with some experimental values that are in a range in the order between 10−11 and 10−9 m2/s [2, 8, 20–22]. More comparative works between experiments and simulations are required for more discussions.Table 3
The self-diffusion coefficients (10−9 m2/s) of CH4-ZIF-8 obtained with AMBER and DE parameters.
Loadings of CH4 per unit cell
Self-diffusion coefficients
AMBER
DE
1
0.29
3.35
2
0.93
4.12
4
0.61
4.21
6
1.02
6.88
8
1.01
5.41


## 4. Conclusions

The molecular dynamics simulations are performed for investigating the structural and dynamical properties of CH4 molecules in the ZIF-8 framework. A total of 1377 configurations of CH4-[C4N2H5]− complexes are used for obtaining the binding energies with BSSE correction at a level of B2PLYPD/6-31G(d). Those energies are used as the data in parameters optimization. Since the DE method is based on random searching in the bounded space, the optimal parameters can be controlled to be positive values to archive the physical meaning of Lennard-Jones model. The energies obtained with DE parameters are more close to the QM data than those obtained with AMBER parameters which are not specifically optimized for these systems of interest.The site-site RDFs between atoms of the framework to atoms H and C of CH4 molecule indicate that atoms H4 and HT are favored sites of adsorption obtained by all simulations. The parameters of DE give shorter average minimum distances from atoms H and C of CH4 to all atoms of the framework than those obtained with AMBER parameters. These obtained results show correspondence between the structural properties and the energies profiles. The molecular distributions obtained by the simulations with AMBER parameters show that the CH4 molecules are distributed within the cavities and there is little diffusion between the cavities, while there is more diffusion between the cavities when using DE parameters. Moreover, the self-diffusion coefficients obtained by the optimized parameters are higher than those obtained by the AMBER parameters for all loadings. These indicate the effects of flexible imidazolate ligands on the movements of their bounded CH4 molecules.



---

*Source: 3926089-2016-05-24.xml*